|
1
|
|
|
# -*- coding: utf-8 -*- |
|
2
|
|
|
|
|
3
|
|
|
from __future__ import division, print_function, absolute_import |
|
4
|
|
|
|
|
5
|
|
|
import datetime as dt |
|
6
|
|
|
import numpy as np |
|
7
|
|
|
import os |
|
8
|
|
|
import warnings |
|
9
|
|
|
|
|
10
|
|
|
from . import helpers |
|
11
|
|
|
|
|
12
|
|
|
# Below try..catch required for autodoc to work on readthedocs |
|
13
|
|
|
try: |
|
14
|
|
|
from . import fortranapex as fa |
|
15
|
|
|
except ImportError as err: |
|
16
|
|
|
warnings.warn("".join(["fortranapex module could not be imported, so ", |
|
17
|
|
|
"apexpy probably won't work. Make sure you have ", |
|
18
|
|
|
"a gfortran compiler. Wheels installation ", |
|
19
|
|
|
"assumes your compiler lives in /opt/local/bin"])) |
|
20
|
|
|
raise err |
|
21
|
|
|
|
|
22
|
|
|
# make sure invalid warnings are always shown |
|
23
|
|
|
warnings.filterwarnings('always', message='.*set to NaN where*', |
|
24
|
|
|
module='apexpy.apex') |
|
25
|
|
|
|
|
26
|
|
|
|
|
27
|
|
|
class ApexHeightError(ValueError): |
|
28
|
|
|
"""Specialized error type definition |
|
29
|
|
|
""" |
|
30
|
|
|
pass |
|
31
|
|
|
|
|
32
|
|
|
|
|
33
|
|
|
class Apex(object): |
|
34
|
|
|
"""Performs coordinate conversions, field-line mapping and base vector |
|
35
|
|
|
calculations. |
|
36
|
|
|
|
|
37
|
|
|
Parameters |
|
38
|
|
|
---------- |
|
39
|
|
|
date : float, :class:`dt.date`, or :class:`dt.datetime`, optional |
|
40
|
|
|
Determines which IGRF coefficients are used in conversions. Uses |
|
41
|
|
|
current date as default. If float, use decimal year. |
|
42
|
|
|
refh : float, optional |
|
43
|
|
|
Reference height in km for apex coordinates (the field lines are mapped |
|
44
|
|
|
to this height) |
|
45
|
|
|
datafile : str, optional |
|
46
|
|
|
Path to custom coefficient file |
|
47
|
|
|
|
|
48
|
|
|
Attributes |
|
49
|
|
|
---------- |
|
50
|
|
|
year : float |
|
51
|
|
|
Decimal year used for the IGRF model |
|
52
|
|
|
refh : float |
|
53
|
|
|
Reference height in km for apex coordinates |
|
54
|
|
|
datafile : str |
|
55
|
|
|
Path to coefficient file |
|
56
|
|
|
|
|
57
|
|
|
Notes |
|
58
|
|
|
----- |
|
59
|
|
|
The calculations use IGRF-13 with coefficients from 1900 to 2025 [1]_. |
|
60
|
|
|
|
|
61
|
|
|
The geodetic reference ellipsoid is WGS84. |
|
62
|
|
|
|
|
63
|
|
|
References |
|
64
|
|
|
---------- |
|
65
|
|
|
|
|
66
|
|
|
.. [1] Thébault, E. et al. (2015), International Geomagnetic Reference |
|
67
|
|
|
Field: the 12th generation, Earth, Planets and Space, 67(1), 79, |
|
68
|
|
|
:doi:`10.1186/s40623-015-0228-9`. |
|
69
|
|
|
|
|
70
|
|
|
""" |
|
71
|
|
|
|
|
72
|
|
|
def __init__(self, date=None, refh=0, datafile=None, fortranlib=None): |
|
73
|
|
|
|
|
74
|
|
|
if datafile is None: |
|
75
|
|
|
datafile = os.path.join(os.path.dirname(__file__), 'apexsh.dat') |
|
76
|
|
|
|
|
77
|
|
|
if fortranlib is None: |
|
78
|
|
|
fortranlib = fa.__file__ |
|
79
|
|
|
|
|
80
|
|
|
self.RE = 6371.009 # mean Earth radius |
|
81
|
|
|
self.set_refh(refh) # reference height |
|
82
|
|
|
|
|
83
|
|
|
if date is None: |
|
84
|
|
|
self.year = helpers.toYearFraction(dt.datetime.utcnow()) |
|
85
|
|
|
else: |
|
86
|
|
|
try: |
|
87
|
|
|
# Convert date/datetime object to decimal year |
|
88
|
|
|
self.year = helpers.toYearFraction(date) |
|
89
|
|
|
except AttributeError: |
|
90
|
|
|
# Failed while finding datetime attribute, so |
|
91
|
|
|
# date is probably an int or float; use directly |
|
92
|
|
|
self.year = date |
|
93
|
|
|
|
|
94
|
|
|
if not os.path.isfile(datafile): |
|
95
|
|
|
raise IOError('Data file does not exist: {}'.format(datafile)) |
|
96
|
|
|
|
|
97
|
|
|
if not os.path.isfile(fortranlib): |
|
98
|
|
|
raise IOError('Fortran library does not exist: {}'.format( |
|
99
|
|
|
fortranlib)) |
|
100
|
|
|
|
|
101
|
|
|
self.datafile = datafile |
|
102
|
|
|
self.fortranlib = fortranlib |
|
103
|
|
|
|
|
104
|
|
|
self.set_epoch(self.year) |
|
105
|
|
|
|
|
106
|
|
|
# vectorize fortran functions |
|
107
|
|
|
self._geo2qd = np.frompyfunc( |
|
108
|
|
|
lambda glat, glon, height: fa.apxg2q(glat, (glon + 180) % 360 - 180, |
|
109
|
|
|
height, 0)[:2], 3, 2) |
|
110
|
|
|
self._geo2apex = np.frompyfunc( |
|
111
|
|
|
lambda glat, glon, height: fa.apxg2all(glat, (glon + 180) % 360 |
|
112
|
|
|
- 180, height, self.refh, |
|
113
|
|
|
0)[2:4], 3, 2) |
|
114
|
|
|
self._geo2apexall = np.frompyfunc( |
|
115
|
|
|
lambda glat, glon, height: fa.apxg2all(glat, (glon + 180) % 360 |
|
116
|
|
|
- 180, height, self.refh, |
|
117
|
|
|
1), 3, 14) |
|
118
|
|
|
self._qd2geo = np.frompyfunc( |
|
119
|
|
|
lambda qlat, qlon, height, precision: fa.apxq2g(qlat, (qlon + 180) |
|
120
|
|
|
% 360 - 180, height, |
|
121
|
|
|
precision), 4, 3) |
|
122
|
|
|
self._basevec = np.frompyfunc( |
|
123
|
|
|
lambda glat, glon, height: fa.apxg2q(glat, (glon + 180) % 360 - 180, |
|
124
|
|
|
height, 1)[2:4], 3, 2) |
|
125
|
|
|
|
|
126
|
|
|
# vectorize other nonvectorized functions |
|
127
|
|
|
self._apex2qd = np.frompyfunc(self._apex2qd_nonvectorized, 3, 2) |
|
128
|
|
|
self._qd2apex = np.frompyfunc(self._qd2apex_nonvectorized, 3, 2) |
|
129
|
|
|
self._get_babs = np.frompyfunc(self._get_babs_nonvectorized, 3, 1) |
|
130
|
|
|
|
|
131
|
|
|
def convert(self, lat, lon, source, dest, height=0, datetime=None, |
|
132
|
|
|
precision=1e-10, ssheight=50 * 6371): |
|
133
|
|
|
"""Converts between geodetic, modified apex, quasi-dipole and MLT. |
|
134
|
|
|
|
|
135
|
|
|
Parameters |
|
136
|
|
|
---------- |
|
137
|
|
|
lat : array_like |
|
138
|
|
|
Latitude |
|
139
|
|
|
lon : array_like |
|
140
|
|
|
Longitude/MLT |
|
141
|
|
|
source : {'geo', 'apex', 'qd', 'mlt'} |
|
142
|
|
|
Input coordinate system |
|
143
|
|
|
dest : {'geo', 'apex', 'qd', 'mlt'} |
|
144
|
|
|
Output coordinate system |
|
145
|
|
|
height : array_like, optional |
|
146
|
|
|
Altitude in km |
|
147
|
|
|
datetime : :class:`datetime.datetime` |
|
148
|
|
|
Date and time for MLT conversions (required for MLT conversions) |
|
149
|
|
|
precision : float, optional |
|
150
|
|
|
Precision of output (degrees) when converting to geo. A negative |
|
151
|
|
|
value of this argument produces a low-precision calculation of |
|
152
|
|
|
geodetic lat/lon based only on their spherical harmonic |
|
153
|
|
|
representation. |
|
154
|
|
|
A positive value causes the underlying Fortran routine to iterate |
|
155
|
|
|
until feeding the output geo lat/lon into geo2qd (APXG2Q) reproduces |
|
156
|
|
|
the input QD lat/lon to within the specified precision (all |
|
157
|
|
|
coordinates being converted to geo are converted to QD first and |
|
158
|
|
|
passed through APXG2Q). |
|
159
|
|
|
ssheight : float, optional |
|
160
|
|
|
Altitude in km to use for converting the subsolar point from |
|
161
|
|
|
geographic to magnetic coordinates. A high altitude is used |
|
162
|
|
|
to ensure the subsolar point is mapped to high latitudes, which |
|
163
|
|
|
prevents the South-Atlantic Anomaly (SAA) from influencing the MLT. |
|
164
|
|
|
|
|
165
|
|
|
Returns |
|
166
|
|
|
------- |
|
167
|
|
|
lat : ndarray or float |
|
168
|
|
|
Converted latitude (if converting to MLT, output latitude is apex) |
|
169
|
|
|
lat : ndarray or float |
|
170
|
|
|
Converted longitude/MLT |
|
171
|
|
|
|
|
172
|
|
|
""" |
|
173
|
|
|
|
|
174
|
|
|
if datetime is None and ('mlt' in [source, dest]): |
|
175
|
|
|
raise ValueError('datetime must be given for MLT calculations') |
|
176
|
|
|
|
|
177
|
|
|
lat = helpers.checklat(lat) |
|
178
|
|
|
|
|
179
|
|
|
if source == dest: |
|
180
|
|
|
return lat, lon |
|
181
|
|
|
# from geo |
|
182
|
|
|
elif source == 'geo' and dest == 'apex': |
|
183
|
|
|
lat, lon = self.geo2apex(lat, lon, height) |
|
184
|
|
|
elif source == 'geo' and dest == 'qd': |
|
185
|
|
|
lat, lon = self.geo2qd(lat, lon, height) |
|
186
|
|
|
elif source == 'geo' and dest == 'mlt': |
|
187
|
|
|
lat, lon = self.geo2apex(lat, lon, height) |
|
188
|
|
|
lon = self.mlon2mlt(lon, datetime, ssheight=ssheight) |
|
189
|
|
|
# from apex |
|
190
|
|
|
elif source == 'apex' and dest == 'geo': |
|
191
|
|
|
lat, lon, _ = self.apex2geo(lat, lon, height, precision=precision) |
|
192
|
|
|
elif source == 'apex' and dest == 'qd': |
|
193
|
|
|
lat, lon = self.apex2qd(lat, lon, height=height) |
|
194
|
|
|
elif source == 'apex' and dest == 'mlt': |
|
195
|
|
|
lon = self.mlon2mlt(lon, datetime, ssheight=ssheight) |
|
196
|
|
|
# from qd |
|
197
|
|
|
elif source == 'qd' and dest == 'geo': |
|
198
|
|
|
lat, lon, _ = self.qd2geo(lat, lon, height, precision=precision) |
|
199
|
|
|
elif source == 'qd' and dest == 'apex': |
|
200
|
|
|
lat, lon = self.qd2apex(lat, lon, height=height) |
|
201
|
|
|
elif source == 'qd' and dest == 'mlt': |
|
202
|
|
|
lat, lon = self.qd2apex(lat, lon, height=height) |
|
203
|
|
|
lon = self.mlon2mlt(lon, datetime, ssheight=ssheight) |
|
204
|
|
|
# from mlt (input latitude assumed apex) |
|
205
|
|
|
elif source == 'mlt' and dest == 'geo': |
|
206
|
|
|
lon = self.mlt2mlon(lon, datetime, ssheight=ssheight) |
|
207
|
|
|
lat, lon, _ = self.apex2geo(lat, lon, height, precision=precision) |
|
208
|
|
|
elif source == 'mlt' and dest == 'apex': |
|
209
|
|
|
lon = self.mlt2mlon(lon, datetime, ssheight=ssheight) |
|
210
|
|
|
elif source == 'mlt' and dest == 'qd': |
|
211
|
|
|
lon = self.mlt2mlon(lon, datetime, ssheight=ssheight) |
|
212
|
|
|
lat, lon = self.apex2qd(lat, lon, height=height) |
|
213
|
|
|
# no other transformations are implemented |
|
214
|
|
|
else: |
|
215
|
|
|
estr = 'Unknown coordinate transformation: ' |
|
216
|
|
|
estr += '{} -> {}'.format(source, dest) |
|
217
|
|
|
raise NotImplementedError(estr) |
|
218
|
|
|
|
|
219
|
|
|
return lat, lon |
|
220
|
|
|
|
|
221
|
|
|
def geo2apex(self, glat, glon, height): |
|
222
|
|
|
"""Converts geodetic to modified apex coordinates. |
|
223
|
|
|
|
|
224
|
|
|
Parameters |
|
225
|
|
|
---------- |
|
226
|
|
|
glat : array_like |
|
227
|
|
|
Geodetic latitude |
|
228
|
|
|
glon : array_like |
|
229
|
|
|
Geodetic longitude |
|
230
|
|
|
height : array_like |
|
231
|
|
|
Altitude in km |
|
232
|
|
|
|
|
233
|
|
|
Returns |
|
234
|
|
|
------- |
|
235
|
|
|
alat : ndarray or float |
|
236
|
|
|
Modified apex latitude |
|
237
|
|
|
alon : ndarray or float |
|
238
|
|
|
Modified apex longitude |
|
239
|
|
|
|
|
240
|
|
|
""" |
|
241
|
|
|
|
|
242
|
|
|
glat = helpers.checklat(glat, name='glat') |
|
243
|
|
|
|
|
244
|
|
|
alat, alon = self._geo2apex(glat, glon, height) |
|
245
|
|
|
|
|
246
|
|
|
if np.any(alat == -9999): |
|
247
|
|
|
warnings.warn('Apex latitude set to NaN where undefined ' |
|
248
|
|
|
'(apex height may be < reference height)') |
|
249
|
|
|
if np.isscalar(alat): |
|
250
|
|
|
alat = np.nan |
|
251
|
|
|
else: |
|
252
|
|
|
alat[alat == -9999] = np.nan |
|
253
|
|
|
|
|
254
|
|
|
# if array is returned, dtype is object, so convert to float |
|
255
|
|
|
return np.float64(alat), np.float64(alon) |
|
256
|
|
|
|
|
257
|
|
|
def apex2geo(self, alat, alon, height, precision=1e-10): |
|
258
|
|
|
"""Converts modified apex to geodetic coordinates. |
|
259
|
|
|
|
|
260
|
|
|
Parameters |
|
261
|
|
|
---------- |
|
262
|
|
|
alat : array_like |
|
263
|
|
|
Modified apex latitude |
|
264
|
|
|
alon : array_like |
|
265
|
|
|
Modified apex longitude |
|
266
|
|
|
height : array_like |
|
267
|
|
|
Altitude in km |
|
268
|
|
|
precision : float, optional |
|
269
|
|
|
Precision of output (degrees). A negative value of this argument |
|
270
|
|
|
produces a low-precision calculation of geodetic lat/lon based only |
|
271
|
|
|
on their spherical harmonic representation. A positive value causes |
|
272
|
|
|
the underlying Fortran routine to iterate until feeding the output |
|
273
|
|
|
geo lat/lon into geo2qd (APXG2Q) reproduces the input QD lat/lon to |
|
274
|
|
|
within the specified precision. |
|
275
|
|
|
|
|
276
|
|
|
Returns |
|
277
|
|
|
------- |
|
278
|
|
|
glat : ndarray or float |
|
279
|
|
|
Geodetic latitude |
|
280
|
|
|
glon : ndarray or float |
|
281
|
|
|
Geodetic longitude |
|
282
|
|
|
error : ndarray or float |
|
283
|
|
|
The angular difference (degrees) between the input QD coordinates |
|
284
|
|
|
and the qlat/qlon produced by feeding the output glat and glon |
|
285
|
|
|
into geo2qd (APXG2Q) |
|
286
|
|
|
|
|
287
|
|
|
""" |
|
288
|
|
|
|
|
289
|
|
|
alat = helpers.checklat(alat, name='alat') |
|
290
|
|
|
|
|
291
|
|
|
qlat, qlon = self.apex2qd(alat, alon, height=height) |
|
292
|
|
|
glat, glon, error = self.qd2geo(qlat, qlon, height, precision=precision) |
|
293
|
|
|
|
|
294
|
|
|
return glat, glon, error |
|
295
|
|
|
|
|
296
|
|
|
def geo2qd(self, glat, glon, height): |
|
297
|
|
|
"""Converts geodetic to quasi-dipole coordinates. |
|
298
|
|
|
|
|
299
|
|
|
Parameters |
|
300
|
|
|
---------- |
|
301
|
|
|
glat : array_like |
|
302
|
|
|
Geodetic latitude |
|
303
|
|
|
glon : array_like |
|
304
|
|
|
Geodetic longitude |
|
305
|
|
|
height : array_like |
|
306
|
|
|
Altitude in km |
|
307
|
|
|
|
|
308
|
|
|
Returns |
|
309
|
|
|
------- |
|
310
|
|
|
qlat : ndarray or float |
|
311
|
|
|
Quasi-dipole latitude |
|
312
|
|
|
qlon : ndarray or float |
|
313
|
|
|
Quasi-dipole longitude |
|
314
|
|
|
|
|
315
|
|
|
""" |
|
316
|
|
|
|
|
317
|
|
|
glat = helpers.checklat(glat, name='glat') |
|
318
|
|
|
|
|
319
|
|
|
qlat, qlon = self._geo2qd(glat, glon, height) |
|
320
|
|
|
|
|
321
|
|
|
# if array is returned, dtype is object, so convert to float |
|
322
|
|
|
return np.float64(qlat), np.float64(qlon) |
|
323
|
|
|
|
|
324
|
|
|
def qd2geo(self, qlat, qlon, height, precision=1e-10): |
|
325
|
|
|
"""Converts quasi-dipole to geodetic coordinates. |
|
326
|
|
|
|
|
327
|
|
|
Parameters |
|
328
|
|
|
---------- |
|
329
|
|
|
qlat : array_like |
|
330
|
|
|
Quasi-dipole latitude |
|
331
|
|
|
qlon : array_like |
|
332
|
|
|
Quasi-dipole longitude |
|
333
|
|
|
height : array_like |
|
334
|
|
|
Altitude in km |
|
335
|
|
|
precision : float, optional |
|
336
|
|
|
Precision of output (degrees). A negative value of this argument |
|
337
|
|
|
produces a low-precision calculation of geodetic lat/lon based only |
|
338
|
|
|
on their spherical harmonic representation. A positive value causes |
|
339
|
|
|
the underlying Fortran routine to iterate until feeding the output |
|
340
|
|
|
geo lat/lon into geo2qd (APXG2Q) reproduces the input QD lat/lon to |
|
341
|
|
|
within the specified precision. |
|
342
|
|
|
|
|
343
|
|
|
Returns |
|
344
|
|
|
------- |
|
345
|
|
|
glat : ndarray or float |
|
346
|
|
|
Geodetic latitude |
|
347
|
|
|
glon : ndarray or float |
|
348
|
|
|
Geodetic longitude |
|
349
|
|
|
error : ndarray or float |
|
350
|
|
|
The angular difference (degrees) between the input QD coordinates |
|
351
|
|
|
and the qlat/qlon produced by feeding the output glat and glon |
|
352
|
|
|
into geo2qd (APXG2Q) |
|
353
|
|
|
|
|
354
|
|
|
""" |
|
355
|
|
|
|
|
356
|
|
|
qlat = helpers.checklat(qlat, name='qlat') |
|
357
|
|
|
|
|
358
|
|
|
glat, glon, error = self._qd2geo(qlat, qlon, height, precision) |
|
359
|
|
|
|
|
360
|
|
|
# if array is returned, dtype is object, so convert to float |
|
361
|
|
|
return np.float64(glat), np.float64(glon), np.float64(error) |
|
362
|
|
|
|
|
363
|
|
View Code Duplication |
def _apex2qd_nonvectorized(self, alat, alon, height): |
|
|
|
|
|
|
364
|
|
|
"""Convert from apex to quasi-dipole (not-vectorised) |
|
365
|
|
|
|
|
366
|
|
|
Parameters |
|
367
|
|
|
----------- |
|
368
|
|
|
alat : (float) |
|
369
|
|
|
Apex latitude in degrees |
|
370
|
|
|
alon : (float) |
|
371
|
|
|
Apex longitude in degrees |
|
372
|
|
|
height : (float) |
|
373
|
|
|
Height in km |
|
374
|
|
|
|
|
375
|
|
|
Returns |
|
376
|
|
|
--------- |
|
377
|
|
|
qlat : (float) |
|
378
|
|
|
Quasi-dipole latitude in degrees |
|
379
|
|
|
qlon : (float) |
|
380
|
|
|
Quasi-diplole longitude in degrees |
|
381
|
|
|
""" |
|
382
|
|
|
|
|
383
|
|
|
alat = helpers.checklat(alat, name='alat') |
|
384
|
|
|
|
|
385
|
|
|
# convert modified apex to quasi-dipole: |
|
386
|
|
|
qlon = alon |
|
387
|
|
|
|
|
388
|
|
|
# apex height |
|
389
|
|
|
hA = self.get_apex(alat) |
|
390
|
|
|
|
|
391
|
|
|
if hA < height: |
|
392
|
|
|
if np.isclose(hA, height, rtol=0, atol=1e-5): |
|
393
|
|
|
# allow for values that are close |
|
394
|
|
|
hA = height |
|
395
|
|
|
else: |
|
396
|
|
|
estr = 'height {:.3g} is > apex height'.format(np.max(height))\ |
|
397
|
|
|
+ ' {:.3g} for alat {:.3g}'.format(hA, alat) |
|
398
|
|
|
raise ApexHeightError(estr) |
|
399
|
|
|
|
|
400
|
|
|
salat = np.sign(alat) if alat != 0 else 1 |
|
401
|
|
|
qlat = salat * np.degrees(np.arccos(np.sqrt((self.RE + height) / |
|
402
|
|
|
(self.RE + hA)))) |
|
403
|
|
|
|
|
404
|
|
|
return qlat, qlon |
|
405
|
|
|
|
|
406
|
|
|
def apex2qd(self, alat, alon, height): |
|
407
|
|
|
"""Converts modified apex to quasi-dipole coordinates. |
|
408
|
|
|
|
|
409
|
|
|
Parameters |
|
410
|
|
|
---------- |
|
411
|
|
|
alat : array_like |
|
412
|
|
|
Modified apex latitude |
|
413
|
|
|
alon : array_like |
|
414
|
|
|
Modified apex longitude |
|
415
|
|
|
height : array_like |
|
416
|
|
|
Altitude in km |
|
417
|
|
|
|
|
418
|
|
|
Returns |
|
419
|
|
|
------- |
|
420
|
|
|
qlat : ndarray or float |
|
421
|
|
|
Quasi-dipole latitude |
|
422
|
|
|
qlon : ndarray or float |
|
423
|
|
|
Quasi-dipole longitude |
|
424
|
|
|
|
|
425
|
|
|
Raises |
|
426
|
|
|
------ |
|
427
|
|
|
ApexHeightError |
|
428
|
|
|
if `height` > apex height |
|
429
|
|
|
|
|
430
|
|
|
""" |
|
431
|
|
|
|
|
432
|
|
|
qlat, qlon = self._apex2qd(alat, alon, height) |
|
433
|
|
|
|
|
434
|
|
|
# if array is returned, the dtype is object, so convert to float |
|
435
|
|
|
return np.float64(qlat), np.float64(qlon) |
|
436
|
|
|
|
|
437
|
|
View Code Duplication |
def _qd2apex_nonvectorized(self, qlat, qlon, height): |
|
|
|
|
|
|
438
|
|
|
|
|
439
|
|
|
qlat = helpers.checklat(qlat, name='qlat') |
|
440
|
|
|
|
|
441
|
|
|
alon = qlon |
|
442
|
|
|
hA = self.get_apex(qlat, height) # apex height |
|
443
|
|
|
|
|
444
|
|
|
if hA < self.refh: |
|
445
|
|
|
if np.isclose(hA, self.refh, rtol=0, atol=1e-5): |
|
446
|
|
|
# allow for values that are close |
|
447
|
|
|
hA = self.refh |
|
448
|
|
|
else: |
|
449
|
|
|
estr = 'apex height ({:.3g}) is < reference height '.format(hA) |
|
450
|
|
|
estr += '({:.3g}) for qlat {:.3g}'.format(self.refh, qlat) |
|
451
|
|
|
raise ApexHeightError(estr) |
|
452
|
|
|
|
|
453
|
|
|
sqlat = np.sign(qlat) if qlat != 0 else 1 |
|
454
|
|
|
alat = sqlat * np.degrees(np.arccos(np.sqrt((self.RE + self.refh) / |
|
455
|
|
|
(self.RE + hA)))) |
|
456
|
|
|
|
|
457
|
|
|
return alat, alon |
|
458
|
|
|
|
|
459
|
|
|
def qd2apex(self, qlat, qlon, height): |
|
460
|
|
|
"""Converts quasi-dipole to modified apex coordinates. |
|
461
|
|
|
|
|
462
|
|
|
Parameters |
|
463
|
|
|
---------- |
|
464
|
|
|
qlat : array_like |
|
465
|
|
|
Quasi-dipole latitude |
|
466
|
|
|
qlon : array_like |
|
467
|
|
|
Quasi-dipole longitude |
|
468
|
|
|
height : array_like |
|
469
|
|
|
Altitude in km |
|
470
|
|
|
|
|
471
|
|
|
Returns |
|
472
|
|
|
------- |
|
473
|
|
|
alat : ndarray or float |
|
474
|
|
|
Modified apex latitude |
|
475
|
|
|
alon : ndarray or float |
|
476
|
|
|
Modified apex longitude |
|
477
|
|
|
|
|
478
|
|
|
Raises |
|
479
|
|
|
------ |
|
480
|
|
|
ApexHeightError |
|
481
|
|
|
if apex height < reference height |
|
482
|
|
|
|
|
483
|
|
|
""" |
|
484
|
|
|
|
|
485
|
|
|
alat, alon = self._qd2apex(qlat, qlon, height) |
|
486
|
|
|
|
|
487
|
|
|
# if array is returned, the dtype is object, so convert to float |
|
488
|
|
|
return np.float64(alat), np.float64(alon) |
|
489
|
|
|
|
|
490
|
|
|
def mlon2mlt(self, mlon, datetime, ssheight=50 * 6371): |
|
491
|
|
|
"""Computes the magnetic local time at the specified magnetic longitude |
|
492
|
|
|
and UT. |
|
493
|
|
|
|
|
494
|
|
|
Parameters |
|
495
|
|
|
---------- |
|
496
|
|
|
mlon : array_like |
|
497
|
|
|
Magnetic longitude (apex and quasi-dipole longitude are always |
|
498
|
|
|
equal) |
|
499
|
|
|
datetime : :class:`datetime.datetime` |
|
500
|
|
|
Date and time |
|
501
|
|
|
ssheight : float, optional |
|
502
|
|
|
Altitude in km to use for converting the subsolar point from |
|
503
|
|
|
geographic to magnetic coordinates. A high altitude is used |
|
504
|
|
|
to ensure the subsolar point is mapped to high latitudes, which |
|
505
|
|
|
prevents the South-Atlantic Anomaly (SAA) from influencing the MLT. |
|
506
|
|
|
|
|
507
|
|
|
Returns |
|
508
|
|
|
------- |
|
509
|
|
|
mlt : ndarray or float |
|
510
|
|
|
Magnetic local time [0, 24) |
|
511
|
|
|
|
|
512
|
|
|
Notes |
|
513
|
|
|
----- |
|
514
|
|
|
To compute the MLT, we find the apex longitude of the subsolar point at |
|
515
|
|
|
the given time. Then the MLT of the given point will be computed from |
|
516
|
|
|
the separation in magnetic longitude from this point (1 hour = 15 |
|
517
|
|
|
degrees). |
|
518
|
|
|
|
|
519
|
|
|
""" |
|
520
|
|
|
ssglat, ssglon = helpers.subsol(datetime) |
|
521
|
|
|
ssalat, ssalon = self.geo2apex(ssglat, ssglon, ssheight) |
|
522
|
|
|
|
|
523
|
|
|
# np.float64 will ensure lists are converted to arrays |
|
524
|
|
|
return (180 + np.float64(mlon) - ssalon) / 15 % 24 |
|
525
|
|
|
|
|
526
|
|
|
def mlt2mlon(self, mlt, datetime, ssheight=50 * 6371): |
|
527
|
|
|
"""Computes the magnetic longitude at the specified magnetic local time |
|
528
|
|
|
and UT. |
|
529
|
|
|
|
|
530
|
|
|
Parameters |
|
531
|
|
|
---------- |
|
532
|
|
|
mlt : array_like |
|
533
|
|
|
Magnetic local time |
|
534
|
|
|
datetime : :class:`datetime.datetime` |
|
535
|
|
|
Date and time |
|
536
|
|
|
ssheight : float, optional |
|
537
|
|
|
Altitude in km to use for converting the subsolar point from |
|
538
|
|
|
geographic to magnetic coordinates. A high altitude is used |
|
539
|
|
|
to ensure the subsolar point is mapped to high latitudes, which |
|
540
|
|
|
prevents the South-Atlantic Anomaly (SAA) from influencing the MLT. |
|
541
|
|
|
|
|
542
|
|
|
Returns |
|
543
|
|
|
------- |
|
544
|
|
|
mlon : ndarray or float |
|
545
|
|
|
Magnetic longitude [0, 360) (apex and quasi-dipole longitude are |
|
546
|
|
|
always equal) |
|
547
|
|
|
|
|
548
|
|
|
Notes |
|
549
|
|
|
----- |
|
550
|
|
|
To compute the magnetic longitude, we find the apex longitude of the |
|
551
|
|
|
subsolar point at the given time. Then the magnetic longitude of the |
|
552
|
|
|
given point will be computed from the separation in magnetic local time |
|
553
|
|
|
from this point (1 hour = 15 degrees). |
|
554
|
|
|
""" |
|
555
|
|
|
|
|
556
|
|
|
ssglat, ssglon = helpers.subsol(datetime) |
|
557
|
|
|
ssalat, ssalon = self.geo2apex(ssglat, ssglon, ssheight) |
|
558
|
|
|
|
|
559
|
|
|
# np.float64 will ensure lists are converted to arrays |
|
560
|
|
|
return (15 * np.float64(mlt) - 180 + ssalon + 360) % 360 |
|
561
|
|
|
|
|
562
|
|
|
def map_to_height(self, glat, glon, height, newheight, conjugate=False, |
|
563
|
|
|
precision=1e-10): |
|
564
|
|
|
"""Performs mapping of points along the magnetic field to the closest |
|
565
|
|
|
or conjugate hemisphere. |
|
566
|
|
|
|
|
567
|
|
|
Parameters |
|
568
|
|
|
---------- |
|
569
|
|
|
glat : array_like |
|
570
|
|
|
Geodetic latitude |
|
571
|
|
|
glon : array_like |
|
572
|
|
|
Geodetic longitude |
|
573
|
|
|
height : array_like |
|
574
|
|
|
Source altitude in km |
|
575
|
|
|
newheight : array_like |
|
576
|
|
|
Destination altitude in km |
|
577
|
|
|
conjugate : bool, optional |
|
578
|
|
|
Map to `newheight` in the conjugate hemisphere instead of the |
|
579
|
|
|
closest hemisphere |
|
580
|
|
|
precision : float, optional |
|
581
|
|
|
Precision of output (degrees). A negative value of this argument |
|
582
|
|
|
produces a low-precision calculation of geodetic lat/lon based only |
|
583
|
|
|
on their spherical harmonic representation. A positive value causes |
|
584
|
|
|
the underlying Fortran routine to iterate until feeding the output |
|
585
|
|
|
geo lat/lon into geo2qd (APXG2Q) reproduces the input QD lat/lon to |
|
586
|
|
|
within the specified precision. |
|
587
|
|
|
|
|
588
|
|
|
Returns |
|
589
|
|
|
------- |
|
590
|
|
|
newglat : ndarray or float |
|
591
|
|
|
Geodetic latitude of mapped point |
|
592
|
|
|
newglon : ndarray or float |
|
593
|
|
|
Geodetic longitude of mapped point |
|
594
|
|
|
error : ndarray or float |
|
595
|
|
|
The angular difference (degrees) between the input QD coordinates |
|
596
|
|
|
and the qlat/qlon produced by feeding the output glat and glon |
|
597
|
|
|
into geo2qd (APXG2Q) |
|
598
|
|
|
|
|
599
|
|
|
Notes |
|
600
|
|
|
----- |
|
601
|
|
|
The mapping is done by converting glat/glon/height to modified apex |
|
602
|
|
|
lat/lon, and converting back to geographic using newheight (if |
|
603
|
|
|
conjugate, use negative apex latitude when converting back) |
|
604
|
|
|
|
|
605
|
|
|
""" |
|
606
|
|
|
|
|
607
|
|
|
alat, alon = self.geo2apex(glat, glon, height) |
|
608
|
|
|
if conjugate: |
|
609
|
|
|
alat = -alat |
|
610
|
|
|
try: |
|
611
|
|
|
newglat, newglon, error = self.apex2geo(alat, alon, newheight, |
|
612
|
|
|
precision=precision) |
|
613
|
|
|
except ApexHeightError: |
|
614
|
|
|
raise ApexHeightError("newheight is > apex height") |
|
615
|
|
|
|
|
616
|
|
|
return newglat, newglon, error |
|
617
|
|
|
|
|
618
|
|
|
def _map_EV_to_height(self, alat, alon, height, newheight, X, EV): |
|
619
|
|
|
|
|
620
|
|
|
# make sure X is array of correct shape |
|
621
|
|
|
if (not (np.ndim(X) == 1 and np.size(X) == 3) and not ( |
|
622
|
|
|
np.ndim(X) == 2 and np.shape(X)[0] == 3)): |
|
623
|
|
|
# raise ValueError because if passing e.g. a (6,) ndarray the |
|
624
|
|
|
# reshape below will work even though the input is invalid |
|
625
|
|
|
raise ValueError(EV + ' must be (3, N) or (3,) ndarray') |
|
626
|
|
|
X = np.reshape(X, (3, np.size(X) // 3)) |
|
627
|
|
|
|
|
628
|
|
|
_, _, _, _, _, _, d1, d2, _, e1, e2, _ = self.basevectors_apex( |
|
629
|
|
|
alat, alon, height, coords='apex') |
|
630
|
|
|
|
|
631
|
|
|
if EV == 'E': |
|
632
|
|
|
v1 = e1 |
|
633
|
|
|
v2 = e2 |
|
634
|
|
|
else: |
|
635
|
|
|
v1 = d1 |
|
636
|
|
|
v2 = d2 |
|
637
|
|
|
|
|
638
|
|
|
# make sure v1 and v2 have shape (3, N) |
|
639
|
|
|
v1 = np.reshape(v1, (3, v1.size // 3)) |
|
640
|
|
|
v2 = np.reshape(v2, (3, v2.size // 3)) |
|
641
|
|
|
|
|
642
|
|
|
X1 = np.sum(X * v1, axis=0) # E dot e1 or V dot d1 |
|
643
|
|
|
X2 = np.sum(X * v2, axis=0) # E dot e2 or V dot d2 |
|
644
|
|
|
|
|
645
|
|
|
_, _, _, _, _, _, d1, d2, _, e1, e2, _ = self.basevectors_apex( |
|
646
|
|
|
alat, alon, newheight, coords='apex') |
|
647
|
|
|
|
|
648
|
|
|
if EV == 'E': |
|
649
|
|
|
v1 = d1 |
|
650
|
|
|
v2 = d2 |
|
651
|
|
|
else: |
|
652
|
|
|
v1 = e1 |
|
653
|
|
|
v2 = e2 |
|
654
|
|
|
|
|
655
|
|
|
# make sure v1 and v2 have shape (3, N) |
|
656
|
|
|
v1 = np.reshape(v1, (3, v1.size // 3)) |
|
657
|
|
|
v2 = np.reshape(v2, (3, v2.size // 3)) |
|
658
|
|
|
|
|
659
|
|
|
X_mapped = X1[np.newaxis, :] * v1 + X2[np.newaxis, :] * v2 |
|
660
|
|
|
|
|
661
|
|
|
return np.squeeze(X_mapped) |
|
662
|
|
|
|
|
663
|
|
|
def map_E_to_height(self, alat, alon, height, newheight, E): |
|
664
|
|
|
"""Performs mapping of electric field along the magnetic field. |
|
665
|
|
|
|
|
666
|
|
|
It is assumed that the electric field is perpendicular to B. |
|
667
|
|
|
|
|
668
|
|
|
Parameters |
|
669
|
|
|
---------- |
|
670
|
|
|
alat : (N,) array_like or float |
|
671
|
|
|
Modified apex latitude |
|
672
|
|
|
alon : (N,) array_like or float |
|
673
|
|
|
Modified apex longitude |
|
674
|
|
|
height : (N,) array_like or float |
|
675
|
|
|
Source altitude in km |
|
676
|
|
|
newheight : (N,) array_like or float |
|
677
|
|
|
Destination altitude in km |
|
678
|
|
|
E : (3,) or (3, N) array_like |
|
679
|
|
|
Electric field (at `alat`, `alon`, `height`) in geodetic east, |
|
680
|
|
|
north, and up components |
|
681
|
|
|
|
|
682
|
|
|
Returns |
|
683
|
|
|
------- |
|
684
|
|
|
E : (3, N) or (3,) ndarray |
|
685
|
|
|
The electric field at `newheight` (geodetic east, north, and up |
|
686
|
|
|
components) |
|
687
|
|
|
|
|
688
|
|
|
""" |
|
689
|
|
|
return self._map_EV_to_height(alat, alon, height, newheight, E, 'E') |
|
690
|
|
|
|
|
691
|
|
|
def map_V_to_height(self, alat, alon, height, newheight, V): |
|
692
|
|
|
"""Performs mapping of electric drift velocity along the magnetic field. |
|
693
|
|
|
|
|
694
|
|
|
It is assumed that the electric field is perpendicular to B. |
|
695
|
|
|
|
|
696
|
|
|
Parameters |
|
697
|
|
|
---------- |
|
698
|
|
|
alat : (N,) array_like or float |
|
699
|
|
|
Modified apex latitude |
|
700
|
|
|
alon : (N,) array_like or float |
|
701
|
|
|
Modified apex longitude |
|
702
|
|
|
height : (N,) array_like or float |
|
703
|
|
|
Source altitude in km |
|
704
|
|
|
newheight : (N,) array_like or float |
|
705
|
|
|
Destination altitude in km |
|
706
|
|
|
V : (3,) or (3, N) array_like |
|
707
|
|
|
Electric drift velocity (at `alat`, `alon`, `height`) in geodetic |
|
708
|
|
|
east, north, and up components |
|
709
|
|
|
|
|
710
|
|
|
Returns |
|
711
|
|
|
------- |
|
712
|
|
|
V : (3, N) or (3,) ndarray |
|
713
|
|
|
The electric drift velocity at `newheight` (geodetic east, north, |
|
714
|
|
|
and up components) |
|
715
|
|
|
|
|
716
|
|
|
""" |
|
717
|
|
|
|
|
718
|
|
|
return self._map_EV_to_height(alat, alon, height, newheight, V, 'V') |
|
719
|
|
|
|
|
720
|
|
|
def basevectors_qd(self, lat, lon, height, coords='geo', precision=1e-10): |
|
721
|
|
|
"""Returns quasi-dipole base vectors f1 and f2 at the specified |
|
722
|
|
|
coordinates. |
|
723
|
|
|
|
|
724
|
|
|
The vectors are described by Richmond [1995] [2]_ and |
|
725
|
|
|
Emmert et al. [2010] [3]_. The vector components are geodetic east and |
|
726
|
|
|
north. |
|
727
|
|
|
|
|
728
|
|
|
Parameters |
|
729
|
|
|
---------- |
|
730
|
|
|
lat : (N,) array_like or float |
|
731
|
|
|
Latitude |
|
732
|
|
|
lon : (N,) array_like or float |
|
733
|
|
|
Longitude |
|
734
|
|
|
height : (N,) array_like or float |
|
735
|
|
|
Altitude in km |
|
736
|
|
|
coords : {'geo', 'apex', 'qd'}, optional |
|
737
|
|
|
Input coordinate system |
|
738
|
|
|
precision : float, optional |
|
739
|
|
|
Precision of output (degrees) when converting to geo. A negative |
|
740
|
|
|
value of this argument produces a low-precision calculation of |
|
741
|
|
|
geodetic lat/lon based only on their spherical harmonic |
|
742
|
|
|
representation. |
|
743
|
|
|
A positive value causes the underlying Fortran routine to iterate |
|
744
|
|
|
until feeding the output geo lat/lon into geo2qd (APXG2Q) reproduces |
|
745
|
|
|
the input QD lat/lon to within the specified precision (all |
|
746
|
|
|
coordinates being converted to geo are converted to QD first and |
|
747
|
|
|
passed through APXG2Q). |
|
748
|
|
|
|
|
749
|
|
|
Returns |
|
750
|
|
|
------- |
|
751
|
|
|
f1 : (2, N) or (2,) ndarray |
|
752
|
|
|
f2 : (2, N) or (2,) ndarray |
|
753
|
|
|
|
|
754
|
|
|
References |
|
755
|
|
|
---------- |
|
756
|
|
|
.. [2] Richmond, A. D. (1995), Ionospheric Electrodynamics Using |
|
757
|
|
|
Magnetic Apex Coordinates, Journal of geomagnetism and |
|
758
|
|
|
geoelectricity, 47(2), 191–212, :doi:`10.5636/jgg.47.191`. |
|
759
|
|
|
|
|
760
|
|
|
.. [3] Emmert, J. T., A. D. Richmond, and D. P. Drob (2010), |
|
761
|
|
|
A computationally compact representation of Magnetic-Apex |
|
762
|
|
|
and Quasi-Dipole coordinates with smooth base vectors, |
|
763
|
|
|
J. Geophys. Res., 115(A8), A08322, :doi:`10.1029/2010JA015326`. |
|
764
|
|
|
|
|
765
|
|
|
""" |
|
766
|
|
|
|
|
767
|
|
|
glat, glon = self.convert(lat, lon, coords, 'geo', height=height, |
|
768
|
|
|
precision=precision) |
|
769
|
|
|
|
|
770
|
|
|
f1, f2 = self._basevec(glat, glon, height) |
|
771
|
|
|
|
|
772
|
|
|
# if inputs are not scalar, each vector is an array of arrays, |
|
773
|
|
|
# so reshape to a single array |
|
774
|
|
|
if f1.dtype == object: |
|
775
|
|
|
f1 = np.vstack(f1).T |
|
776
|
|
|
f2 = np.vstack(f2).T |
|
777
|
|
|
|
|
778
|
|
|
return f1, f2 |
|
779
|
|
|
|
|
780
|
|
|
def basevectors_apex(self, lat, lon, height, coords='geo', precision=1e-10): |
|
781
|
|
|
"""Returns base vectors in quasi-dipole and apex coordinates. |
|
782
|
|
|
|
|
783
|
|
|
The vectors are described by Richmond [1995] [4]_ and |
|
784
|
|
|
Emmert et al. [2010] [5]_. The vector components are geodetic east, |
|
785
|
|
|
north, and up (only east and north for `f1` and `f2`). |
|
786
|
|
|
|
|
787
|
|
|
Parameters |
|
788
|
|
|
---------- |
|
789
|
|
|
lat : (N,) array_like or float |
|
790
|
|
|
Latitude |
|
791
|
|
|
lon : (N,) array_like or float |
|
792
|
|
|
Longitude |
|
793
|
|
|
height : (N,) array_like or float |
|
794
|
|
|
Altitude in km |
|
795
|
|
|
coords : {'geo', 'apex', 'qd'}, optional |
|
796
|
|
|
Input coordinate system |
|
797
|
|
|
precision : float, optional |
|
798
|
|
|
Precision of output (degrees) when converting to geo. A negative |
|
799
|
|
|
value of this argument produces a low-precision calculation of |
|
800
|
|
|
geodetic lat/lon based only on their spherical harmonic |
|
801
|
|
|
representation. |
|
802
|
|
|
A positive value causes the underlying Fortran routine to iterate |
|
803
|
|
|
until feeding the output geo lat/lon into geo2qd (APXG2Q) reproduces |
|
804
|
|
|
the input QD lat/lon to within the specified precision (all |
|
805
|
|
|
coordinates being converted to geo are converted to QD first and |
|
806
|
|
|
passed through APXG2Q). |
|
807
|
|
|
|
|
808
|
|
|
Returns |
|
809
|
|
|
------- |
|
810
|
|
|
f3, g1, g2, g3, d1, d2, d3, e1, e2, e3 : (3, N) or (3,) ndarray |
|
811
|
|
|
|
|
812
|
|
|
Notes |
|
813
|
|
|
----- |
|
814
|
|
|
`f3`, `g1`, `g2`, and `g3` are not part of the Fortran code |
|
815
|
|
|
by Emmert et al. [2010] [5]_. They are calculated by this |
|
816
|
|
|
Python library according to the following equations in |
|
817
|
|
|
Richmond [1995] [4]_: |
|
818
|
|
|
|
|
819
|
|
|
* `g1`: Eqn. 6.3 |
|
820
|
|
|
* `g2`: Eqn. 6.4 |
|
821
|
|
|
* `g3`: Eqn. 6.5 |
|
822
|
|
|
* `f3`: Eqn. 6.8 |
|
823
|
|
|
|
|
824
|
|
|
References |
|
825
|
|
|
---------- |
|
826
|
|
|
|
|
827
|
|
|
.. [4] Richmond, A. D. (1995), Ionospheric Electrodynamics Using |
|
828
|
|
|
Magnetic Apex Coordinates, Journal of geomagnetism and |
|
829
|
|
|
geoelectricity, 47(2), 191–212, :doi:`10.5636/jgg.47.191`. |
|
830
|
|
|
|
|
831
|
|
|
.. [5] Emmert, J. T., A. D. Richmond, and D. P. Drob (2010), |
|
832
|
|
|
A computationally compact representation of Magnetic-Apex |
|
833
|
|
|
and Quasi-Dipole coordinates with smooth base vectors, |
|
834
|
|
|
J. Geophys. Res., 115(A8), A08322, :doi:`10.1029/2010JA015326`. |
|
835
|
|
|
|
|
836
|
|
|
""" |
|
837
|
|
|
|
|
838
|
|
|
glat, glon = self.convert(lat, lon, coords, 'geo', height=height, |
|
839
|
|
|
precision=precision) |
|
840
|
|
|
|
|
841
|
|
|
returnvals = self._geo2apexall(glat, glon, height) |
|
842
|
|
|
qlat = np.float64(returnvals[0]) |
|
843
|
|
|
alat = np.float64(returnvals[2]) |
|
844
|
|
|
f1, f2 = returnvals[4:6] |
|
845
|
|
|
d1, d2, d3 = returnvals[7:10] |
|
846
|
|
|
e1, e2, e3 = returnvals[11:14] |
|
847
|
|
|
|
|
848
|
|
|
# if inputs are not scalar, each vector is an array of arrays, |
|
849
|
|
|
# so reshape to a single array |
|
850
|
|
|
if f1.dtype == object: |
|
851
|
|
|
f1 = np.vstack(f1).T |
|
852
|
|
|
f2 = np.vstack(f2).T |
|
853
|
|
|
d1 = np.vstack(d1).T |
|
854
|
|
|
d2 = np.vstack(d2).T |
|
855
|
|
|
d3 = np.vstack(d3).T |
|
856
|
|
|
e1 = np.vstack(e1).T |
|
857
|
|
|
e2 = np.vstack(e2).T |
|
858
|
|
|
e3 = np.vstack(e3).T |
|
859
|
|
|
|
|
860
|
|
|
# make sure arrays are 2D |
|
861
|
|
|
f1 = f1.reshape((2, f1.size // 2)) |
|
862
|
|
|
f2 = f2.reshape((2, f2.size // 2)) |
|
863
|
|
|
d1 = d1.reshape((3, d1.size // 3)) |
|
864
|
|
|
d2 = d2.reshape((3, d2.size // 3)) |
|
865
|
|
|
d3 = d3.reshape((3, d3.size // 3)) |
|
866
|
|
|
e1 = e1.reshape((3, e1.size // 3)) |
|
867
|
|
|
e2 = e2.reshape((3, e2.size // 3)) |
|
868
|
|
|
e3 = e3.reshape((3, e3.size // 3)) |
|
869
|
|
|
|
|
870
|
|
|
# compute f3, g1, g2, g3 |
|
871
|
|
|
F1 = np.vstack((f1, np.zeros_like(f1[0]))) |
|
872
|
|
|
F2 = np.vstack((f2, np.zeros_like(f2[0]))) |
|
873
|
|
|
F = np.cross(F1.T, F2.T).T[-1] |
|
874
|
|
|
cosI = helpers.getcosIm(alat) |
|
875
|
|
|
k = np.array([0, 0, 1], dtype=np.float64).reshape((3, 1)) |
|
876
|
|
|
g1 = ((self.RE + np.float64(height)) |
|
877
|
|
|
/ (self.RE + self.refh)) ** (3 / 2) * d1 / F |
|
878
|
|
|
g2 = -1.0 / (2.0 * F * np.tan(np.radians(qlat))) * ( |
|
879
|
|
|
k + ((self.RE + np.float64(height)) |
|
880
|
|
|
/ (self.RE + self.refh)) * d2 / cosI) |
|
881
|
|
|
g3 = k * F |
|
882
|
|
|
f3 = np.cross(g1.T, g2.T).T |
|
883
|
|
|
|
|
884
|
|
|
if np.any(alat == -9999): |
|
885
|
|
|
warnings.warn(''.join(['Base vectors g, d, e, and f3 set to NaN ', |
|
886
|
|
|
'where apex latitude is undefined (apex ', |
|
887
|
|
|
'height may be < reference height)'])) |
|
888
|
|
|
mask = alat == -9999 |
|
889
|
|
|
f3 = np.where(mask, np.nan, f3) |
|
890
|
|
|
g1 = np.where(mask, np.nan, g1) |
|
891
|
|
|
g2 = np.where(mask, np.nan, g2) |
|
892
|
|
|
g3 = np.where(mask, np.nan, g3) |
|
893
|
|
|
d1 = np.where(mask, np.nan, d1) |
|
894
|
|
|
d2 = np.where(mask, np.nan, d2) |
|
895
|
|
|
d3 = np.where(mask, np.nan, d3) |
|
896
|
|
|
e1 = np.where(mask, np.nan, e1) |
|
897
|
|
|
e2 = np.where(mask, np.nan, e2) |
|
898
|
|
|
e3 = np.where(mask, np.nan, e3) |
|
899
|
|
|
|
|
900
|
|
|
return tuple(np.squeeze(x) for x in |
|
901
|
|
|
[f1, f2, f3, g1, g2, g3, d1, d2, d3, e1, e2, e3]) |
|
902
|
|
|
|
|
903
|
|
|
def get_apex(self, lat, height=None): |
|
904
|
|
|
""" Calculate apex height |
|
905
|
|
|
|
|
906
|
|
|
Parameters |
|
907
|
|
|
----------- |
|
908
|
|
|
lat : (float) |
|
909
|
|
|
Latitude in degrees |
|
910
|
|
|
height : (float or NoneType) |
|
911
|
|
|
Height above the surface of the earth in km or NoneType to use |
|
912
|
|
|
reference height (default=None) |
|
913
|
|
|
|
|
914
|
|
|
Returns |
|
915
|
|
|
---------- |
|
916
|
|
|
apex_height : (float) |
|
917
|
|
|
Height of the field line apex in km |
|
918
|
|
|
""" |
|
919
|
|
|
lat = helpers.checklat(lat, name='alat') |
|
920
|
|
|
if height is None: |
|
921
|
|
|
height = self.refh |
|
922
|
|
|
|
|
923
|
|
|
cos_lat_squared = np.cos(np.radians(lat)) ** 2 |
|
924
|
|
|
apex_height = (self.RE + height) / cos_lat_squared - self.RE |
|
925
|
|
|
|
|
926
|
|
|
return apex_height |
|
927
|
|
|
|
|
928
|
|
|
def set_epoch(self, year): |
|
929
|
|
|
"""Updates the epoch for all subsequent conversions. |
|
930
|
|
|
|
|
931
|
|
|
Parameters |
|
932
|
|
|
---------- |
|
933
|
|
|
year : float |
|
934
|
|
|
Decimal year |
|
935
|
|
|
|
|
936
|
|
|
""" |
|
937
|
|
|
# f2py |
|
938
|
|
|
self.year = np.float64(year) |
|
939
|
|
|
fa.loadapxsh(self.datafile, self.year) |
|
940
|
|
|
igrf_fn = os.path.join(os.path.dirname(__file__), 'igrf13coeffs.txt') |
|
941
|
|
|
if not os.path.exists(igrf_fn): |
|
942
|
|
|
raise OSError("File {} does not exist".format(igrf_fn)) |
|
943
|
|
|
fa.cofrm(self.year, igrf_fn) |
|
944
|
|
|
|
|
945
|
|
|
def set_refh(self, refh): |
|
946
|
|
|
"""Updates the apex reference height for all subsequent conversions. |
|
947
|
|
|
|
|
948
|
|
|
Parameters |
|
949
|
|
|
---------- |
|
950
|
|
|
refh : float |
|
951
|
|
|
Apex reference height in km |
|
952
|
|
|
|
|
953
|
|
|
Notes |
|
954
|
|
|
----- |
|
955
|
|
|
The reference height is the height to which field lines will be mapped, |
|
956
|
|
|
and is only relevant for conversions involving apex (not quasi-dipole). |
|
957
|
|
|
|
|
958
|
|
|
""" |
|
959
|
|
|
self.refh = refh |
|
960
|
|
|
|
|
961
|
|
|
def _get_babs_nonvectorized(self, glat, glon, height): |
|
962
|
|
|
bnorth, beast, bdown, babs = fa.feldg(1, glat, glon, height) |
|
963
|
|
|
# BABS is in guass, so convert to tesla |
|
964
|
|
|
return babs / 10000.0 |
|
965
|
|
|
|
|
966
|
|
|
def get_babs(self, glat, glon, height): |
|
967
|
|
|
"""Returns the magnitude of the IGRF magnetic field in tesla. |
|
968
|
|
|
|
|
969
|
|
|
Parameters |
|
970
|
|
|
---------- |
|
971
|
|
|
glat : array_like |
|
972
|
|
|
Geodetic latitude |
|
973
|
|
|
glon : array_like |
|
974
|
|
|
Geodetic longitude |
|
975
|
|
|
height : array_like |
|
976
|
|
|
Altitude in km |
|
977
|
|
|
|
|
978
|
|
|
Returns |
|
979
|
|
|
------- |
|
980
|
|
|
babs : ndarray or float |
|
981
|
|
|
Magnitude of the IGRF magnetic field |
|
982
|
|
|
|
|
983
|
|
|
""" |
|
984
|
|
|
|
|
985
|
|
|
babs = self._get_babs(glat, glon, height) |
|
986
|
|
|
|
|
987
|
|
|
# if array is returned, the dtype is object, so convert to float |
|
988
|
|
|
return np.float64(babs) |
|
989
|
|
|
|
|
990
|
|
|
def bvectors_apex(self, lat, lon, height, coords='geo', precision=1e-10): |
|
991
|
|
|
"""Returns the magnetic field vectors in apex coordinates. |
|
992
|
|
|
|
|
993
|
|
|
The apex magnetic field vectors described by Richmond [1995] [4]_ and |
|
994
|
|
|
Emmert et al. [2010] [5]_, specfically the Be3 and Bd3 components. The |
|
995
|
|
|
vector components are geodetic east, north, and up. |
|
996
|
|
|
|
|
997
|
|
|
Parameters |
|
998
|
|
|
---------- |
|
999
|
|
|
lat : (N,) array_like or float |
|
1000
|
|
|
Latitude |
|
1001
|
|
|
lon : (N,) array_like or float |
|
1002
|
|
|
Longitude |
|
1003
|
|
|
height : (N,) array_like or float |
|
1004
|
|
|
Altitude in km |
|
1005
|
|
|
coords : {'geo', 'apex', 'qd'}, optional |
|
1006
|
|
|
Input coordinate system |
|
1007
|
|
|
precision : float, optional |
|
1008
|
|
|
Precision of output (degrees) when converting to geo. A negative |
|
1009
|
|
|
value of this argument produces a low-precision calculation of |
|
1010
|
|
|
geodetic lat/lon based only on their spherical harmonic |
|
1011
|
|
|
representation. |
|
1012
|
|
|
A positive value causes the underlying Fortran routine to iterate |
|
1013
|
|
|
until feeding the output geo lat/lon into geo2qd (APXG2Q) reproduces |
|
1014
|
|
|
the input QD lat/lon to within the specified precision (all |
|
1015
|
|
|
coordinates being converted to geo are converted to QD first and |
|
1016
|
|
|
passed through APXG2Q). |
|
1017
|
|
|
|
|
1018
|
|
|
Returns |
|
1019
|
|
|
------- |
|
1020
|
|
|
Be3: (1, N) or (1,) ndarray |
|
1021
|
|
|
e3 : (3, N) or (3,) ndarray |
|
1022
|
|
|
Bd3: (1, N) or (1,) ndarray |
|
1023
|
|
|
d3 : (3, N) or (3,) ndarray |
|
1024
|
|
|
|
|
1025
|
|
|
Notes |
|
1026
|
|
|
----- |
|
1027
|
|
|
Be3 is not equivalent to the magnitude of the IGRF magnitude, but is |
|
1028
|
|
|
instead equal to the IGRF magnitude divided by a scaling factor, D. |
|
1029
|
|
|
Similarly, Bd3 is the IGRF magnitude multiplied by D. |
|
1030
|
|
|
|
|
1031
|
|
|
See Richmond, A. D. (1995) [4]_ equations 3.13 and 3.14 |
|
1032
|
|
|
|
|
1033
|
|
|
References |
|
1034
|
|
|
---------- |
|
1035
|
|
|
Richmond, A. D. (1995) [4]_ |
|
1036
|
|
|
Emmert, J. T. et al. (2010) [5]_ |
|
1037
|
|
|
|
|
1038
|
|
|
""" |
|
1039
|
|
|
glat, glon = self.convert(lat, lon, coords, 'geo', height=height, |
|
1040
|
|
|
precision=precision) |
|
1041
|
|
|
|
|
1042
|
|
|
babs = self.get_babs(glat, glon, height) |
|
1043
|
|
|
|
|
1044
|
|
|
_, _, _, _, _, _, d1, d2, d3, _, _, e3 = self.basevectors_apex( |
|
1045
|
|
|
glat, glon, height, coords='geo') |
|
1046
|
|
|
d1_cross_d2 = np.cross(d1.T, d2.T).T |
|
1047
|
|
|
D = np.sqrt(np.sum(d1_cross_d2 ** 2, axis=0)) |
|
1048
|
|
|
|
|
1049
|
|
|
Be3 = babs / D |
|
1050
|
|
|
Bd3 = babs * D |
|
1051
|
|
|
|
|
1052
|
|
|
return Be3, e3, Bd3, d3 |
|
1053
|
|
|
|
aburrell /
apexpy
Loading history...