|
1
|
|
|
# -*- coding: utf-8 -*- |
|
2
|
|
|
# vim:fileencoding=utf-8 |
|
3
|
|
|
# |
|
4
|
|
|
# Copyright (c) 2017-2018 Stefan Bender |
|
5
|
|
|
# |
|
6
|
|
|
# This file is part of sciapy. |
|
7
|
|
|
# sciapy is free software: you can redistribute it or modify it |
|
8
|
|
|
# under the terms of the GNU General Public License as published by |
|
9
|
|
|
# the Free Software Foundation, version 2. |
|
10
|
|
|
# See accompanying LICENSE file or http://www.gnu.org/licenses/gpl-2.0.html. |
|
11
|
1 |
|
"""IGRF geomagnetic coordinates |
|
12
|
|
|
|
|
13
|
|
|
This is a python (mix) version of GMPOLE and GMCOORD from |
|
14
|
|
|
http://www.ngdc.noaa.gov/geomag/geom_util/utilities_home.shtml |
|
15
|
|
|
to transform geodetic to geomagnetic coordinates. |
|
16
|
|
|
It uses the IGRF 2012 model and coefficients [#]_. |
|
17
|
|
|
|
|
18
|
|
|
.. [#] Thébault et al. 2015, |
|
19
|
|
|
International Geomagnetic Reference Field: the 12th generation. |
|
20
|
|
|
Earth, Planets and Space, 67 (79) |
|
21
|
|
|
http://nora.nerc.ac.uk/id/eprint/511258 |
|
22
|
|
|
https://doi.org/10.1186/s40623-015-0228-9 |
|
23
|
|
|
""" |
|
24
|
1 |
|
from __future__ import absolute_import, division, print_function |
|
25
|
|
|
|
|
26
|
1 |
|
import logging |
|
27
|
1 |
|
from collections import namedtuple |
|
28
|
1 |
|
from pkg_resources import resource_filename |
|
29
|
|
|
|
|
30
|
1 |
|
import numpy as np |
|
31
|
1 |
|
from scipy.interpolate import interp1d |
|
32
|
1 |
|
from scipy.special import lpmn |
|
33
|
|
|
|
|
34
|
1 |
|
__all__ = ['gmpole', 'gmag_igrf'] |
|
35
|
|
|
|
|
36
|
|
|
# The WGS84 reference ellipsoid |
|
37
|
1 |
|
Earth_ellipsoid = { |
|
38
|
|
|
"a": 6378.137, # semi-major axis of the ellipsoid in km |
|
39
|
|
|
"b": 6356.7523142, # semi-minor axis of the ellipsoid in km |
|
40
|
|
|
"fla": 1. / 298.257223563, # flattening |
|
41
|
|
|
"re": 6371.2 # Earth's radius in km |
|
42
|
|
|
} |
|
43
|
|
|
|
|
44
|
1 |
View Code Duplication |
def _ellipsoid(ellipsoid_data=Earth_ellipsoid): |
|
|
|
|
|
|
45
|
|
|
# extends the dictionary with the eccentricities |
|
46
|
1 |
|
ell = namedtuple('ellip', ["a", "b", "fla", "eps", "epssq", "re"]) |
|
47
|
1 |
|
ell.a = ellipsoid_data["a"] |
|
48
|
1 |
|
ell.b = ellipsoid_data["b"] |
|
49
|
1 |
|
ell.fla = ellipsoid_data["fla"] |
|
50
|
1 |
|
ell.re = ellipsoid_data["re"] |
|
51
|
|
|
# first eccentricity squared |
|
52
|
1 |
|
ell.epssq = 1. - ell.b**2 / ell.a**2 |
|
53
|
|
|
# first eccentricity |
|
54
|
1 |
|
ell.eps = np.sqrt(ell.epssq) |
|
55
|
1 |
|
return ell |
|
56
|
|
|
|
|
57
|
1 |
|
def _date_to_frac_year(year, month, day): |
|
58
|
|
|
# fractional year by dividing the day of year by the overall |
|
59
|
|
|
# number of days in that year |
|
60
|
1 |
|
extraday = 0 |
|
61
|
1 |
|
if ((year % 4 == 0) and (year % 100 != 0)) or (year % 400 == 0): |
|
62
|
|
|
extraday = 1 |
|
63
|
1 |
|
month_days = [0, 31, 28 + extraday, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31] |
|
64
|
1 |
|
doy = np.sum(month_days[:month]) + day |
|
65
|
1 |
|
return year + (doy - 1) / (365.0 + extraday) |
|
66
|
|
|
|
|
67
|
1 |
View Code Duplication |
def _load_igrf_file(filename="IGRF.tab"): |
|
|
|
|
|
|
68
|
|
|
"""Load IGRF coefficients |
|
69
|
|
|
|
|
70
|
|
|
Parameters |
|
71
|
|
|
---------- |
|
72
|
|
|
filename: str, optional |
|
73
|
|
|
The file with the IGRF coefficients. |
|
74
|
|
|
|
|
75
|
|
|
Returns |
|
76
|
|
|
------- |
|
77
|
|
|
interpol: `scipy.interpolate.interp1d` instance |
|
78
|
|
|
Interpolator instance, called with the fractional |
|
79
|
|
|
year to obtain the IGRF coefficients for the epoch. |
|
80
|
|
|
""" |
|
81
|
1 |
|
igrf_tab = np.genfromtxt(filename, skip_header=3, dtype=None) |
|
82
|
|
|
|
|
83
|
1 |
|
sv = igrf_tab[igrf_tab.dtype.names[-1]][1:].astype(np.float) |
|
84
|
|
|
|
|
85
|
1 |
|
years = np.asarray(igrf_tab[0].tolist()[3:-1]) |
|
86
|
1 |
|
years = np.append(years, [years[-1] + 5]) |
|
87
|
|
|
|
|
88
|
1 |
|
coeffs = [] |
|
89
|
1 |
|
for i in range(1, len(igrf_tab)): |
|
90
|
1 |
|
coeff = np.asarray(igrf_tab[i].tolist()[3:-1]) |
|
91
|
1 |
|
coeff = np.append(coeff, np.array([5]) * sv[0] + coeff[-1]) |
|
92
|
1 |
|
coeffs.append(coeff) |
|
93
|
|
|
|
|
94
|
1 |
|
return interp1d(years, coeffs) |
|
95
|
|
|
|
|
96
|
1 |
View Code Duplication |
def _geod_to_spher(phi, lon, Ellip, HeightAboveEllipsoid=0.): |
|
|
|
|
|
|
97
|
|
|
"""Convert geodetic to spherical coordinates |
|
98
|
|
|
|
|
99
|
|
|
Converts geodetic coordinates on the WGS-84 reference ellipsoid |
|
100
|
|
|
to Earth-centered Earth-fixed Cartesian coordinates, |
|
101
|
|
|
and then to spherical coordinates. |
|
102
|
|
|
""" |
|
103
|
1 |
|
CosLat = np.cos(np.radians(phi)) |
|
104
|
1 |
|
SinLat = np.sin(np.radians(phi)) |
|
105
|
|
|
|
|
106
|
|
|
# compute the local radius of curvature on the WGS-84 reference ellipsoid |
|
107
|
1 |
|
rc = Ellip.a / np.sqrt(1.0 - Ellip.epssq * SinLat**2) |
|
108
|
|
|
|
|
109
|
|
|
# compute ECEF Cartesian coordinates of specified point (for longitude=0) |
|
110
|
1 |
|
xp = (rc + HeightAboveEllipsoid) * CosLat |
|
111
|
1 |
|
zp = (rc * (1.0 - Ellip.epssq) + HeightAboveEllipsoid) * SinLat |
|
112
|
|
|
|
|
113
|
|
|
# compute spherical radius and angle phi of specified point |
|
114
|
1 |
|
rg = np.sqrt(xp**2 + zp**2) |
|
115
|
|
|
# geocentric latitude |
|
116
|
1 |
|
phig = np.degrees(np.arcsin(zp / rg)) |
|
117
|
1 |
|
return phig, lon, rg |
|
118
|
|
|
|
|
119
|
1 |
View Code Duplication |
def _igrf_model(coeffs, Lmax, r, theta, phi, R_E=Earth_ellipsoid["re"]): |
|
|
|
|
|
|
120
|
|
|
"""Evaluates the IGRF model function at the given location |
|
121
|
|
|
""" |
|
122
|
|
|
rho = R_E / r |
|
123
|
|
|
sin_theta = np.sin(theta) |
|
124
|
|
|
cos_theta = np.cos(theta) |
|
125
|
|
|
Plm, dPlm = lpmn(Lmax, Lmax, cos_theta) |
|
126
|
|
|
logging.debug("R_E: %s, r: %s, rho: %s", R_E, r, rho) |
|
127
|
|
|
logging.debug("rho: %s, theta: %s, sin_theta: %s, cos_theta: %s", |
|
128
|
|
|
rho, theta, sin_theta, cos_theta) |
|
129
|
|
|
Bx, By, Bz = 0., 0., 0. # Btheta, Bphi, Br |
|
130
|
|
|
idx = 0 |
|
131
|
|
|
rho_l = rho |
|
132
|
|
|
K_l1 = 1. |
|
133
|
|
|
for l in range(1, Lmax + 1): |
|
134
|
|
|
rho_l *= rho # rho^(l+1) |
|
135
|
|
|
# m = 0 values, h_l^m = 0 |
|
136
|
|
|
Bxl = rho_l * coeffs[idx] * dPlm[0, l] * (-sin_theta) |
|
137
|
|
|
Byl = 0. |
|
138
|
|
|
Bzl = -(l + 1) * rho_l * coeffs[idx] * Plm[0, l] |
|
139
|
|
|
if l > 1: |
|
140
|
|
|
K_l1 *= np.sqrt((l - 1) / (l + 1)) |
|
141
|
|
|
idx += 1 |
|
142
|
|
|
K_lm = -K_l1 |
|
143
|
|
|
for m in range(1, l + 1): |
|
144
|
|
|
cfi = K_lm * rho_l * (coeffs[idx] * np.cos(m * phi) |
|
145
|
|
|
+ coeffs[idx + 1] * np.sin(m * phi)) |
|
146
|
|
|
Bxl += cfi * dPlm[m, l] * (-sin_theta) |
|
147
|
|
|
Bzl += -(l + 1) * cfi * Plm[m, l] |
|
148
|
|
|
if sin_theta != 0: |
|
149
|
|
|
Byl += K_lm * rho_l * m * Plm[m, l] * ( |
|
150
|
|
|
- coeffs[idx] * np.sin(m * phi) + |
|
151
|
|
|
coeffs[idx + 1] * np.cos(m * phi)) |
|
152
|
|
|
else: |
|
153
|
|
|
Byl += 0. |
|
154
|
|
|
if m < l: |
|
155
|
|
|
# K_lm for the next m |
|
156
|
|
|
K_lm /= -np.sqrt((l + m + 1) * (l - m)) |
|
157
|
|
|
idx += 2 |
|
158
|
|
|
Bx += Bxl |
|
159
|
|
|
By += Byl |
|
160
|
|
|
Bz += Bzl |
|
161
|
|
|
Bx = rho * Bx |
|
162
|
|
|
By = -rho * By / sin_theta |
|
163
|
|
|
Bz = rho * Bz |
|
164
|
|
|
return Bx, By, Bz |
|
165
|
|
|
|
|
166
|
1 |
View Code Duplication |
def igrf_mag(date, lat, lon, alt, filename="IGRF.tab"): |
|
|
|
|
|
|
167
|
|
|
"""Evaluate the local magnetic field using the IGRF model |
|
168
|
|
|
|
|
169
|
|
|
Evaluates the IGRF coefficients to calculate the Earth's |
|
170
|
|
|
magnetic field at the given location. |
|
171
|
|
|
The results agree to within a few decimal places with |
|
172
|
|
|
https://ngdc.noaa.gov/geomag-web/ |
|
173
|
|
|
|
|
174
|
|
|
Parameters |
|
175
|
|
|
---------- |
|
176
|
|
|
date: `datetime.date` or `datetime.datetime` instance |
|
177
|
|
|
The date for the evaluation. |
|
178
|
|
|
lat: float |
|
179
|
|
|
Geographic latitude in degrees north |
|
180
|
|
|
lon: float |
|
181
|
|
|
Geographic longitude in degrees east |
|
182
|
|
|
alt: float |
|
183
|
|
|
Altitude above ground in km. |
|
184
|
|
|
filename: str, optional |
|
185
|
|
|
File containing the IGRF coefficients. |
|
186
|
|
|
|
|
187
|
|
|
Returns |
|
188
|
|
|
------- |
|
189
|
|
|
Bx: float |
|
190
|
|
|
Northward component of the magnetic field, B_N. |
|
191
|
|
|
By: float |
|
192
|
|
|
Eastward component of the magnetic field, B_E. |
|
193
|
|
|
Bz: float |
|
194
|
|
|
Downward component of the magnetic field, B_D. |
|
195
|
|
|
""" |
|
196
|
|
|
ellip = _ellipsoid() |
|
197
|
|
|
# date should be datetime.datetime or datetime.date instance, |
|
198
|
|
|
# or something else that provides .year, .month, and .day attributes |
|
199
|
|
|
frac_year = _date_to_frac_year(date.year, date.month, date.day) |
|
200
|
|
|
glat, glon, grad = _geod_to_spher(lat, lon, ellip, alt) |
|
201
|
|
|
sin_theta = np.sin(np.radians(90. - glat)) |
|
202
|
|
|
cos_theta = np.cos(np.radians(90. - glat)) |
|
203
|
|
|
|
|
204
|
|
|
rho = np.sqrt((ellip.a * sin_theta)**2 + (ellip.b * cos_theta)**2) |
|
205
|
|
|
r = np.sqrt(alt**2 + 2 * alt * rho + |
|
206
|
|
|
(ellip.a**4 * sin_theta**2 + ellip.b**4 * cos_theta**2) / rho**2) |
|
207
|
|
|
cd = (alt + rho) / r |
|
208
|
|
|
sd = (ellip.a**2 - ellip.b**2) / rho * cos_theta * sin_theta / r |
|
209
|
|
|
logging.debug("rho: %s, r: %s, (alt + rho) / r: %s, R_E / (R_E + h): %s, R_E / r: %s", |
|
210
|
|
|
rho, r, cd, ellip.re / (ellip.re + alt), ellip.re / r) |
|
211
|
|
|
|
|
212
|
|
|
cos_theta, sin_theta = cos_theta*cd - sin_theta*sd, sin_theta*cd + cos_theta*sd |
|
213
|
|
|
logging.debug("r: %s, spherical coordinates (radius, rho, theta, lat): %s, %s, %s, %s", |
|
214
|
|
|
r, grad, rho, np.degrees(np.arccos(cos_theta)), 90. - glat) |
|
215
|
|
|
|
|
216
|
|
|
# evaluate the IGRF model in spherical coordinates |
|
217
|
|
|
igrf_file = resource_filename(__name__, filename) |
|
218
|
|
|
igrf_coeffs = _load_igrf_file(igrf_file)(frac_year) |
|
219
|
|
|
Bx, By, Bz = _igrf_model(igrf_coeffs, 13, r, np.radians(90. - glat), np.radians(glon)) |
|
220
|
|
|
logging.debug("spherical geomagnetic field (Bx, By, Bz): %s, %s, %s", Bx, By, Bz) |
|
221
|
|
|
logging.debug("spherical dip coordinates: lat %s, lon %s", |
|
222
|
|
|
np.degrees(np.arctan(0.5 * Bz / np.sqrt(Bx**2 + By**2))), |
|
223
|
|
|
np.degrees(np.arctan(-By / Bz))) |
|
224
|
|
|
# back to geodetic coordinates |
|
225
|
|
|
Bx, Bz = cd * Bx + sd * Bz, cd * Bz - sd * Bx |
|
226
|
|
|
logging.debug("geodetic geomagnetic field (Bx, By, Bz): %s, %s, %s", Bx, By, Bz) |
|
227
|
|
|
logging.debug("geodetic dip coordinates: lat %s, lon %s", |
|
228
|
|
|
np.degrees(np.arctan(0.5 * Bz / np.sqrt(Bx**2 + By**2))), |
|
229
|
|
|
np.degrees(np.arctan(-By / Bz))) |
|
230
|
|
|
return Bx, By, Bz |
|
231
|
|
|
|
|
232
|
1 |
View Code Duplication |
def gmpole(date, r_e=Earth_ellipsoid["re"], filename="IGRF.tab"): |
|
|
|
|
|
|
233
|
|
|
"""Centered dipole geomagnetic pole coordinates |
|
234
|
|
|
|
|
235
|
|
|
Parameters |
|
236
|
|
|
---------- |
|
237
|
|
|
date: datetime.datetime or datetime.date |
|
238
|
|
|
r_e: float, optional |
|
239
|
|
|
Earth radius to evaluate the dipole's off-centre shift. |
|
240
|
|
|
filename: str, optional |
|
241
|
|
|
File containing the IGRF parameters. |
|
242
|
|
|
|
|
243
|
|
|
Returns |
|
244
|
|
|
------- |
|
245
|
|
|
(lat_n, phi_n): tuple of floats |
|
246
|
|
|
Geographic latitude and longitude of the centered dipole |
|
247
|
|
|
magnetic north pole. |
|
248
|
|
|
(lat_s, phi_s): tuple of floats |
|
249
|
|
|
Geographic latitude and longitude of the centered dipole |
|
250
|
|
|
magnetic south pole. |
|
251
|
|
|
(dX, dY, dZ): tuple of floats |
|
252
|
|
|
Magnetic variations in Earth-centered Cartesian coordinates |
|
253
|
|
|
for shifting the dipole off-center. |
|
254
|
|
|
B_0: float |
|
255
|
|
|
The magnitude of the magnetic field. |
|
256
|
|
|
""" |
|
257
|
1 |
|
igrf_file = resource_filename(__name__, filename) |
|
258
|
1 |
|
gh_func = _load_igrf_file(igrf_file) |
|
259
|
|
|
|
|
260
|
1 |
|
frac_year = _date_to_frac_year(date.year, date.month, date.day) |
|
261
|
1 |
|
logging.debug("fractional year: %s", frac_year) |
|
262
|
1 |
|
g10, g11, h11, g20, g21, h21, g22, h22 = gh_func(frac_year)[:8] |
|
263
|
|
|
|
|
264
|
|
|
# This function finds the location of the north magnetic pole in spherical coordinates. |
|
265
|
|
|
# The equations are from Wallace H. Campbell's "Introduction to Geomagnetic Fields" |
|
266
|
1 |
|
B_0_sq = g10**2 + g11**2 + h11**2 |
|
267
|
1 |
|
theta_n = np.arccos(-g10 / np.sqrt(B_0_sq)) |
|
268
|
1 |
|
phi_n = np.arctan2(h11, g11) |
|
269
|
1 |
|
lat_n = 0.5 * np.pi - theta_n |
|
270
|
1 |
|
logging.debug("centered dipole pole coordinates " |
|
271
|
|
|
"(lat, theta, phi): %s, %s, %s", |
|
272
|
|
|
np.degrees(lat_n), np.degrees(theta_n), np.degrees(phi_n)) |
|
273
|
|
|
|
|
274
|
|
|
# calculate dipole offset according to Fraser and Smith, 1987 |
|
275
|
1 |
|
L_0 = 2 * g10 * g20 + np.sqrt(3) * (g11 * g21 + h11 * h21) |
|
276
|
1 |
|
L_1 = -g11 * g20 + np.sqrt(3) * (g10 * g21 + g11 * g22 + h11 * h22) |
|
277
|
1 |
|
L_2 = -h11 * g20 + np.sqrt(3) * (g10 * h21 - h11 * g22 + g11 * h22) |
|
278
|
1 |
|
E = (L_0 * g10 + L_1 * g11 + L_2 * h11) / (4 * B_0_sq) |
|
279
|
|
|
|
|
280
|
1 |
|
xi = (L_0 - g10 * E) / (3 * B_0_sq) |
|
281
|
1 |
|
eta = (L_1 - g11 * E) / (3 * B_0_sq) |
|
282
|
1 |
|
zeta = (L_2 - h11 * E) / (3 * B_0_sq) |
|
283
|
1 |
|
dx = eta * r_e |
|
284
|
1 |
|
dy = zeta * r_e |
|
285
|
1 |
|
dz = xi * r_e |
|
286
|
|
|
|
|
287
|
1 |
|
delta = np.sqrt(dx**2 + dy**2 + dz**2) |
|
288
|
1 |
|
theta_d = np.arccos(dz / delta) |
|
289
|
1 |
|
lambda_d = 0.5 * np.pi - theta_d |
|
290
|
1 |
|
phi_d = np.pi + np.arctan(dy / dx) |
|
291
|
|
|
|
|
292
|
1 |
|
sin_lat_ed = (np.sin(lambda_d) * np.sin(lat_n) |
|
293
|
|
|
+ np.cos(lambda_d) * np.cos(lat_n) * np.cos(phi_d - phi_n)) |
|
294
|
1 |
|
lat_ed = np.arcsin(sin_lat_ed) |
|
295
|
1 |
|
theta_ed = 0.5 * np.pi - lat_ed |
|
296
|
|
|
|
|
297
|
1 |
|
sin_lon_ed = np.sin(theta_d) * np.sin(phi_d - phi_n) / np.sin(theta_ed) |
|
298
|
1 |
|
lon_ed = np.pi - np.arcsin(sin_lon_ed) |
|
299
|
1 |
|
logging.debug("eccentric dipole pole coordinates " |
|
300
|
|
|
"(lat, theta, lon): %s, %s, %s", |
|
301
|
|
|
np.degrees(theta_ed), np.degrees(lat_ed), np.degrees(lon_ed)) |
|
302
|
|
|
|
|
303
|
1 |
|
dX = delta * np.sin(theta_ed) * np.cos(lon_ed) |
|
304
|
1 |
|
dY = delta * np.sin(theta_ed) * np.sin(lon_ed) |
|
305
|
1 |
|
dZ = delta * np.cos(theta_ed) |
|
306
|
1 |
|
logging.debug("magnetic variations (dX, dY, dZ): %s, %s, %s", dX, dY, dZ) |
|
307
|
|
|
|
|
308
|
|
|
# North pole, south pole coordinates |
|
309
|
1 |
|
return ((np.degrees(lat_n), np.degrees(phi_n)), |
|
310
|
|
|
(-np.degrees(lat_n), np.degrees(phi_n + np.pi)), |
|
311
|
|
|
(dX, dY, dZ), |
|
312
|
|
|
np.sqrt(B_0_sq)) |
|
313
|
|
|
|
|
314
|
1 |
View Code Duplication |
def gmag_igrf(date, lat, lon, alt=0., |
|
|
|
|
|
|
315
|
|
|
centered_dipole=False, |
|
316
|
|
|
igrf_name="IGRF.tab"): |
|
317
|
|
|
"""Centered or eccentric dipole geomagnetic coordinates |
|
318
|
|
|
|
|
319
|
|
|
Parameters |
|
320
|
|
|
---------- |
|
321
|
|
|
date: datetime.datetime |
|
322
|
|
|
lat: float |
|
323
|
|
|
Geographic latitude in degrees north |
|
324
|
|
|
lon: float |
|
325
|
|
|
Geographic longitude in degrees east |
|
326
|
|
|
alt: float, optional |
|
327
|
|
|
Altitude in km. Default: 0. |
|
328
|
|
|
centered_dipole: bool, optional |
|
329
|
|
|
Returns the centered dipole geomagnetic coordinates |
|
330
|
|
|
if set to True, returns the eccentric dipole |
|
331
|
|
|
geomagnetic coordinates if set to False. |
|
332
|
|
|
Default: False |
|
333
|
|
|
igrf_name: str, optional |
|
334
|
|
|
Default: "IGRF.tab" |
|
335
|
|
|
|
|
336
|
|
|
Returns |
|
337
|
|
|
------- |
|
338
|
|
|
geomag_latitude: numpy.ndarray or float |
|
339
|
|
|
Geomagnetic latitude in eccentric dipole coordinates, |
|
340
|
|
|
centered dipole coordinates if `centered_dipole` is True. |
|
341
|
|
|
geomag_longitude: numpy.ndarray or float |
|
342
|
|
|
Geomagnetic longitude in eccentric dipole coordinates, |
|
343
|
|
|
centered dipole coordinates if `centered_dipole` is True. |
|
344
|
|
|
""" |
|
345
|
1 |
|
ellip = _ellipsoid() |
|
346
|
1 |
|
glat, glon, grad = _geod_to_spher(lat, lon, ellip, alt) |
|
347
|
1 |
|
(lat_GMP, lon_GMP), _, (dX, dY, dZ), B_0 = gmpole(date, ellip.re, igrf_name) |
|
348
|
1 |
|
latr, lonr = np.radians(glat), np.radians(glon) |
|
349
|
1 |
|
lat_GMPr, lon_GMPr = np.radians(lat_GMP), np.radians(lon_GMP) |
|
350
|
1 |
|
sin_lat_gmag = (np.sin(latr) * np.sin(lat_GMPr) |
|
351
|
|
|
+ np.cos(latr) * np.cos(lat_GMPr) * np.cos(lonr - lon_GMPr)) |
|
352
|
1 |
|
lon_gmag_y = np.cos(latr) * np.sin(lonr - lon_GMPr) |
|
353
|
1 |
|
lon_gmag_x = (np.cos(latr) * np.sin(lat_GMPr) * np.cos(lonr - lon_GMPr) |
|
354
|
|
|
- np.sin(latr) * np.cos(lat_GMPr)) |
|
355
|
1 |
|
lat_gmag = np.arcsin(sin_lat_gmag) |
|
356
|
1 |
|
lat_gmag_geod = np.arctan(np.tan(lat_gmag) / (1. - ellip.epssq)) |
|
357
|
|
|
|
|
358
|
1 |
|
B_r = -2. * B_0 * sin_lat_gmag / (grad / ellip.re)**3 |
|
359
|
1 |
|
B_th = -B_0 * np.cos(lat_gmag) / (grad / ellip.re)**3 |
|
360
|
1 |
|
logging.debug("B_r: %s, B_th: %s, dip lat: %s", |
|
361
|
|
|
-B_r, B_th, np.degrees(np.arctan(0.5 * B_r / B_th))) |
|
362
|
|
|
|
|
363
|
1 |
|
lon_gmag = np.arctan2(lon_gmag_y, lon_gmag_x) |
|
364
|
1 |
|
logging.debug("centered dipole coordinates: " |
|
365
|
|
|
"lat_gmag: %s, lon_gmag: %s", |
|
366
|
|
|
np.degrees(lat_gmag), np.degrees(lon_gmag)) |
|
367
|
1 |
|
logging.debug("lat_gmag_geod: %s, lat_GMPr: %s", |
|
368
|
|
|
np.degrees(lat_gmag_geod), np.degrees(lat_GMPr)) |
|
369
|
1 |
|
if centered_dipole: |
|
370
|
|
|
return (np.degrees(lat_gmag), np.degrees(lon_gmag)) |
|
371
|
|
|
|
|
372
|
|
|
# eccentric dipole coordinates (shifted dipole) |
|
373
|
1 |
|
phi_ed = np.arctan2((grad * np.cos(lat_gmag) * np.sin(lon_gmag) - dY), |
|
374
|
|
|
(grad * np.cos(lat_gmag) * np.cos(lon_gmag) - dX)) |
|
375
|
1 |
|
theta_ed = np.arctan2((grad * np.cos(lat_gmag) * np.cos(lon_gmag) - dX), |
|
376
|
|
|
((grad * np.sin(lat_gmag) - dZ) * np.cos(phi_ed))) |
|
377
|
1 |
|
lat_ed = 0.5 * np.pi - theta_ed |
|
378
|
1 |
|
lat_ed_geod = np.arctan(np.tan(lat_ed) / (1. - ellip.epssq)) |
|
379
|
1 |
|
logging.debug("lats ed: %s", np.degrees(lat_ed)) |
|
380
|
1 |
|
logging.debug("lats ed geod: %s", np.degrees(lat_ed_geod)) |
|
381
|
1 |
|
logging.debug("phis ed: %s", np.degrees(phi_ed)) |
|
382
|
|
|
return (np.degrees(lat_ed_geod), np.degrees(phi_ed)) |
|
383
|
|
|
|
st-bender /
sciapy
Loading history...