1
|
|
|
# -*- coding: utf-8 -*- |
2
|
|
|
# vim:fileencoding=utf-8 |
3
|
|
|
# |
4
|
|
|
# Copyright (c) 2019 Stefan Bender |
5
|
|
|
# |
6
|
|
|
# This file is part of pynrlmsise00. |
7
|
|
|
# pynrlmsise00 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
|
|
|
"""NRLMSISE-00 python wrapper for the C version [#]_ |
12
|
|
|
|
13
|
|
|
.. [#] https://www.brodo.de/space/nrlmsise |
14
|
|
|
""" |
15
|
|
|
from __future__ import absolute_import, division, print_function |
16
|
|
|
|
17
|
|
|
from functools import wraps |
18
|
|
|
|
19
|
|
|
import numpy as np |
20
|
|
|
|
21
|
|
|
from ._nrlmsise00 import gtd7, gtd7d |
22
|
|
|
|
23
|
|
|
__all__ = ["gtd7_flat", "gtd7d_flat", "msise_model", "msise_flat", "scale_height"] |
24
|
|
|
|
25
|
|
|
|
26
|
|
|
def _doc_param(*sub): |
27
|
|
|
def dec(obj): |
28
|
|
|
obj.__doc__ = obj.__doc__.format(*sub) |
29
|
|
|
return obj |
30
|
|
|
return dec |
31
|
|
|
|
32
|
|
|
|
33
|
|
|
def vectorize_function(pyfunc, **kwargs): |
34
|
|
|
"""Function decorator acting like :class:`numpy.vectorize`. |
35
|
|
|
|
36
|
|
|
All arguments are passed to :class:`numpy.vectorize`. Unlike |
37
|
|
|
:class:`numpy.vectorize`, this decorator does not convert the decorated |
38
|
|
|
function into an instance of the :class:`numpy.vectorize` class. As a |
39
|
|
|
result, the decorated function's docstring will be taken into account |
40
|
|
|
properly by doc generators (pydoc or Sphinx) |
41
|
|
|
""" |
42
|
|
|
wrapped = np.vectorize(pyfunc, **kwargs) |
43
|
|
|
@wraps(pyfunc) |
44
|
|
|
def run_wrapped(*args, **kwargs): |
45
|
|
|
return wrapped(*args, **kwargs) |
46
|
|
|
run_wrapped.__doc__ = wrapped.__doc__ # preserve np.vectorize's `doc` arg |
47
|
|
|
return run_wrapped |
48
|
|
|
|
49
|
|
|
|
50
|
|
|
@_doc_param(gtd7.__doc__) |
51
|
|
|
def _gtd7_flat(*args, **kwargs): |
52
|
|
|
"""Flattened variant of the MSIS `gtd7()` function |
53
|
|
|
|
54
|
|
|
Returns a single 11-element :class:`numpy.ndarray` instead of |
55
|
|
|
the two lists. All arguments except the keywords `flags` and |
56
|
|
|
`ap_a` can be :class:`numpy.ndarray` to facilitate calculations |
57
|
|
|
at many locations/times. |
58
|
|
|
|
59
|
|
|
{0} |
60
|
|
|
""" |
61
|
|
|
ds, ts = gtd7(*args, **kwargs) |
62
|
|
|
return np.asarray(ds + ts) |
63
|
|
|
|
64
|
|
|
|
65
|
|
|
@_doc_param(gtd7d.__doc__) |
66
|
|
|
def _gtd7d_flat(*args, **kwargs): |
67
|
|
|
"""Flattened variant of the MSIS `gtd7d()` function |
68
|
|
|
|
69
|
|
|
Returns a single 11-element :class:`numpy.ndarray` instead of |
70
|
|
|
the two lists. All arguments except the keywords `flags` and |
71
|
|
|
`ap_a` can be :class:`numpy.ndarray` to facilitate calculations |
72
|
|
|
at many locations/times. |
73
|
|
|
|
74
|
|
|
{0} |
75
|
|
|
""" |
76
|
|
|
ds, ts = gtd7d(*args, **kwargs) |
77
|
|
|
return np.asarray(ds + ts) |
78
|
|
|
|
79
|
|
|
|
80
|
|
|
gtd7_flat = vectorize_function(_gtd7_flat, |
81
|
|
|
signature='(),(),(),(),(),(),(),(),(),()->(n)', |
82
|
|
|
excluded=["ap_a", "flags"]) |
83
|
|
|
|
84
|
|
|
gtd7d_flat = vectorize_function(_gtd7d_flat, |
85
|
|
|
signature='(),(),(),(),(),(),(),(),(),()->(n)', |
86
|
|
|
excluded=["ap_a", "flags"]) |
87
|
|
|
|
88
|
|
|
|
89
|
|
|
def msise_model(time, alt, lat, lon, f107a, f107, ap, |
90
|
|
|
lst=None, ap_a=None, flags=None, method="gtd7"): |
91
|
|
|
"""Interface to `gtd7()` [1]_ and `gtd7d()` [2]_ |
92
|
|
|
|
93
|
|
|
Calls the C model function using a :class:`datetime.datetime` |
94
|
|
|
instance to calculate the day of year, seconds and so on. |
95
|
|
|
The `method` keyword decides which version to call, the difference |
96
|
|
|
being that `gtd7d()` includes anomalous oxygen in the total |
97
|
|
|
mass density (`d[5]`), `gtd7()` does not. |
98
|
|
|
|
99
|
|
|
.. [1] https://git.linta.de/?p=~brodo/nrlmsise-00.git;a=blob;f=nrlmsise-00.c#l916 |
100
|
|
|
.. [2] https://git.linta.de/?p=~brodo/nrlmsise-00.git;a=blob;f=nrlmsise-00.c#l1044 |
101
|
|
|
|
102
|
|
|
Parameters |
103
|
|
|
---------- |
104
|
|
|
time: datetime.datetime |
105
|
|
|
Date and time as a `datetime.dateime`. |
106
|
|
|
alt: float |
107
|
|
|
Altitude in km. |
108
|
|
|
lat: float |
109
|
|
|
Latitude in degrees north. |
110
|
|
|
lon: float |
111
|
|
|
Longitude in degrees east. |
112
|
|
|
f107a: float |
113
|
|
|
The observed f107a (81-day running mean of f107) centred at date. |
114
|
|
|
f107: float |
115
|
|
|
The observed f107 value on the previous day. |
116
|
|
|
ap: float |
117
|
|
|
The ap value at date. |
118
|
|
|
lst: float, optional |
119
|
|
|
The local solar time, can be different from the calculated one. |
120
|
|
|
ap_a: list, optional |
121
|
|
|
List of length 7 containing ap values to be used when flags[9] is set |
122
|
|
|
to -1, otherwise no effect. |
123
|
|
|
flags: list, optional |
124
|
|
|
List of length 24 setting the NRLMSIS switches explicitly. |
125
|
|
|
method: string, optional |
126
|
|
|
Set to "gtd7d" to use `gtd7d()` (which includes anomalous oxygen |
127
|
|
|
in the total mass density) instead of the "standard" `gtd7()` function |
128
|
|
|
without it. |
129
|
|
|
|
130
|
|
|
Returns |
131
|
|
|
------- |
132
|
|
|
densities: list of floats |
133
|
|
|
0. He number density [cm^-3] |
134
|
|
|
1. O number density [cm^-3] |
135
|
|
|
2. N2 number density [cm^-3] |
136
|
|
|
3. O2 number density [cm^-3] |
137
|
|
|
4. AR number density [cm^-3] |
138
|
|
|
5. total mass density [g cm^-3] (includes d[8] in gtd7d) |
139
|
|
|
6. H number density [cm^-3] |
140
|
|
|
7. N number density [cm^-3] |
141
|
|
|
8. Anomalous oxygen number density [cm^-3] |
142
|
|
|
temperatures: list of floats |
143
|
|
|
0. Exospheric temperature [K] |
144
|
|
|
1. Temperature at `alt` [K] |
145
|
|
|
|
146
|
|
|
Note |
147
|
|
|
---- |
148
|
|
|
No date and time conversion will be attempted, the input needs to be |
149
|
|
|
converted before, `dateutil` or `pandas` provide convenient functions. |
150
|
|
|
For example to convert :class:`pandas.DatetimeIndex` |
151
|
|
|
use ``<index>.to_pydatetime()``, |
152
|
|
|
for :class:`astropy.time.Time` use ``<Time>.to_datetime()``. |
153
|
|
|
|
154
|
|
|
The local solar time is calculated from time and longitude, except when |
155
|
|
|
`lst` is set, then that value is used. |
156
|
|
|
|
157
|
|
|
The solar and geomagnetic indices have to be provided, so far the values |
158
|
|
|
are not included in the module. |
159
|
|
|
""" |
160
|
|
|
year = time.year |
161
|
|
|
doy = int(time.strftime("%j")) |
162
|
|
|
sec = (time.hour * 3600. |
163
|
|
|
+ time.minute * 60. |
164
|
|
|
+ time.second |
165
|
|
|
+ time.microsecond * 1e-6) |
166
|
|
|
if lst is None: |
167
|
|
|
lst = sec / 3600. + lon / 15.0 |
168
|
|
|
|
169
|
|
|
kwargs = {} |
170
|
|
|
if ap_a is not None: |
171
|
|
|
kwargs.update({"ap_a": ap_a}) |
172
|
|
|
if flags is not None: |
173
|
|
|
kwargs.update({"flags": flags}) |
174
|
|
|
|
175
|
|
|
if method == "gtd7d": |
176
|
|
|
return gtd7d(year, doy, sec, alt, lat, lon, lst, f107a, f107, ap, **kwargs) |
177
|
|
|
return gtd7(year, doy, sec, alt, lat, lon, lst, f107a, f107, ap, **kwargs) |
178
|
|
|
|
179
|
|
|
|
180
|
|
|
def _msise_flat(*args, **kwargs): |
181
|
|
|
ds, ts = msise_model(*args, **kwargs) |
182
|
|
|
return np.asarray(ds + ts) |
183
|
|
|
|
184
|
|
|
|
185
|
|
|
_msise_flatv = np.vectorize(_msise_flat, |
186
|
|
|
signature='(),(),(),(),(),(),(),(),(),()->(n)', |
187
|
|
|
excluded=["method"]) |
188
|
|
|
|
189
|
|
|
|
190
|
|
|
@_doc_param(msise_model.__doc__.replace("Interface", "interface")) |
191
|
|
|
def msise_flat(*args, **kwargs): |
192
|
|
|
"""Flattened {0} |
193
|
|
|
Attention |
194
|
|
|
--------- |
195
|
|
|
This flattened version returns a single 11-element array instead of two |
196
|
|
|
separate lists. Additionally, it can take :class:`numpy.ndarray` as input. |
197
|
|
|
However, the `time` input needs to contain entries of :class:`datetime.datetime`, |
198
|
|
|
e.g. using :meth:`pandas.DatetimeIndex.to_pydatetime()` |
199
|
|
|
or :meth:`astropy.time.Time.to_datetime()`. |
200
|
|
|
""" |
201
|
|
|
# Set keyword arguments to None if not given to make `np.vectorize` happy |
202
|
|
|
lst = kwargs.pop("lst", None) |
203
|
|
|
ap_a = kwargs.pop("ap_a", None) |
204
|
|
|
flags = kwargs.pop("flags", None) |
205
|
|
|
|
206
|
|
|
return _msise_flatv(*args, lst=lst, ap_a=ap_a, flags=flags, **kwargs) |
207
|
|
|
|
208
|
|
|
|
209
|
|
|
def scale_height(alt, lat, molw, temp): |
210
|
|
|
"""Atmospheric scale height |
211
|
|
|
|
212
|
|
|
Extracted from the C-code for easy access, with |
213
|
|
|
constants updated to the standard SI values. |
214
|
|
|
It is reasonably fast and :mod:`numpy`-broadcasting compatible. |
215
|
|
|
|
216
|
|
|
Parameters |
217
|
|
|
---------- |
218
|
|
|
alt: float or array_like |
219
|
|
|
Altitude in [km]. |
220
|
|
|
lat: float or array_like |
221
|
|
|
Geodetic latitude in [degrees N]. |
222
|
|
|
molw: float or array_like |
223
|
|
|
Molecular mass at alt in [kg / mol]. |
224
|
|
|
Can be derived by dividing the total mass density by |
225
|
|
|
the sum of number densities, e.g. from :func:`gtd7()`, |
226
|
|
|
multiplied by Avogadro's constant (6.0221407e23 / mol). |
227
|
|
|
temp: float or array_like |
228
|
|
|
Temperature at alt in [K]. |
229
|
|
|
|
230
|
|
|
Returns |
231
|
|
|
------- |
232
|
|
|
scale_height: float or array_like |
233
|
|
|
Scale height in [m]. |
234
|
|
|
""" |
235
|
|
|
# Rgas = 8.31446261815324 # J / K / mol |
236
|
|
|
# dgtr = 1.7453292519943295e-2 # rad / deg |
237
|
|
|
c2 = np.cos(2. * 1.7453292519943295e-2 * lat) |
238
|
|
|
gsurf = 9.80665 * (1.0 - 0.0026373 * c2) # m / s2 |
239
|
|
|
# re in [km] |
240
|
|
|
re = 2. * gsurf / (3.085462e-6 + 2.27e-9 * c2) * 1.0e-3 |
241
|
|
|
g = gsurf / (1 + alt / re)**2 |
242
|
|
|
return 8.31446261815324 * temp / (g * molw) |
243
|
|
|
|