nrlmsise00.core._doc_param()   A
last analyzed

Complexity

Conditions 1

Size

Total Lines 5
Code Lines 5

Duplication

Lines 0
Ratio 0 %

Importance

Changes 0
Metric Value
eloc 5
dl 0
loc 5
rs 10
c 0
b 0
f 0
cc 1
nop 1
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