nrlmsise00.core.scale_height() - Code Metrics - Inspection of "dde1c2320e21e138d8d25a6428196c27ceb8b3fe" - st-bender/pynrlmsise00 - Measure and Improve Code Quality continuously with Scrutinizer

Passed

Branch — master (dde1c2)

by Stefan

created 2021-10-03 09:33 UTC

nrlmsise00.core.scale_height() A

↳ Parent: nrlmsise00.core

Complexity

Conditions

Size

Total Lines	34
Code Lines	6

Duplication

Lines	0
Ratio	0 %

Importance

Changes

Metric	Value
eloc	6
dl	0
loc	34
rs	10
c	0
b	0
f	0
cc	1
nop	4

# -*- coding: utf-8 -*-
# vim:fileencoding=utf-8
#
# Copyright (c) 2019 Stefan Bender
#
# This file is part of pynrlmsise00.
# pynrlmsise00 is free software: you can redistribute it or modify it
# under the terms of the GNU General Public License as published by
# the Free Software Foundation, version 2.
# See accompanying LICENSE file or http://www.gnu.org/licenses/gpl-2.0.html.
"""NRLMSISE-00 python wrapper for the C version [#]_

.. [#] https://www.brodo.de/space/nrlmsise
"""
from __future__ import absolute_import, division, print_function

from functools import wraps

import numpy as np

from ._nrlmsise00 import gtd7, gtd7d

__all__ = ["gtd7_flat", "gtd7d_flat", "msise_model", "msise_flat", "scale_height"]


def _doc_param(*sub):
	def dec(obj):
		obj.__doc__ = obj.__doc__.format(*sub)
		return obj
	return dec


def vectorize_function(pyfunc, **kwargs):
	"""Function decorator acting like :class:`numpy.vectorize`.

	All arguments are passed to :class:`numpy.vectorize`. Unlike
	:class:`numpy.vectorize`, this decorator does not convert the decorated
	function into an instance of the :class:`numpy.vectorize` class. As a
	result, the decorated function's docstring will be taken into account
	properly by doc generators (pydoc or Sphinx)
	"""
	wrapped = np.vectorize(pyfunc, **kwargs)
	@wraps(pyfunc)
	def run_wrapped(*args, **kwargs):
		return wrapped(*args, **kwargs)
	run_wrapped.__doc__ = wrapped.__doc__  # preserve np.vectorize's `doc` arg
	return run_wrapped


@_doc_param(gtd7.__doc__)
def _gtd7_flat(*args, **kwargs):
	"""Flattened variant of the MSIS `gtd7()` function

	Returns a single 11-element :class:`numpy.ndarray` instead of
	the two lists. All arguments except the keywords `flags` and
	`ap_a` can be :class:`numpy.ndarray` to facilitate calculations
	at many locations/times.

	{0}
	"""
	ds, ts = gtd7(*args, **kwargs)
	return np.asarray(ds + ts)


@_doc_param(gtd7d.__doc__)
def _gtd7d_flat(*args, **kwargs):
	"""Flattened variant of the MSIS `gtd7d()` function

	Returns a single 11-element :class:`numpy.ndarray` instead of
	the two lists. All arguments except the keywords `flags` and
	`ap_a` can be :class:`numpy.ndarray` to facilitate calculations
	at many locations/times.

	{0}
	"""
	ds, ts = gtd7d(*args, **kwargs)
	return np.asarray(ds + ts)


gtd7_flat = vectorize_function(_gtd7_flat,
		signature='(),(),(),(),(),(),(),(),(),()->(n)',
		excluded=["ap_a", "flags"])

gtd7d_flat = vectorize_function(_gtd7d_flat,
		signature='(),(),(),(),(),(),(),(),(),()->(n)',
		excluded=["ap_a", "flags"])


def msise_model(time, alt, lat, lon, f107a, f107, ap,
		lst=None, ap_a=None, flags=None, method="gtd7"):
	"""Interface to `gtd7()` [1]_ and `gtd7d()` [2]_

	Calls the C model function using a :class:`datetime.datetime`
	instance to calculate the day of year, seconds and so on.
	The `method` keyword decides which version to call, the difference
	being that `gtd7d()` includes anomalous oxygen in the total
	mass density (`d[5]`), `gtd7()` does not.

	.. [1] https://git.linta.de/?p=~brodo/nrlmsise-00.git;a=blob;f=nrlmsise-00.c#l916
	.. [2] https://git.linta.de/?p=~brodo/nrlmsise-00.git;a=blob;f=nrlmsise-00.c#l1044

	Parameters
	----------
	time: datetime.datetime
		Date and time as a `datetime.dateime`.
	alt: float
		Altitude in km.
	lat: float
		Latitude in degrees north.
	lon: float
		Longitude in degrees east.
	f107a: float
		The observed f107a (81-day running mean of f107) centred at date.
	f107: float
		The observed f107 value on the previous day.
	ap: float
		The ap value at date.
	lst: float, optional
		The local solar time, can be different from the calculated one.
	ap_a: list, optional
		List of length 7 containing ap values to be used when flags[9] is set
		to -1, otherwise no effect.
	flags: list, optional
		List of length 24 setting the NRLMSIS switches explicitly.
	method: string, optional
		Set to "gtd7d" to use `gtd7d()` (which includes anomalous oxygen
		in the total mass density) instead of the "standard" `gtd7()` function
		without it.

	Returns
	-------
	densities: list of floats
		0. He number density [cm^-3]
		1. O number density [cm^-3]
		2. N2 number density [cm^-3]
		3. O2 number density [cm^-3]
		4. AR number density [cm^-3]
		5. total mass density [g cm^-3] (includes d[8] in gtd7d)
		6. H number density [cm^-3]
		7. N number density [cm^-3]
		8. Anomalous oxygen number density [cm^-3]
	temperatures: list of floats
		0. Exospheric temperature [K]
		1. Temperature at `alt` [K]

	Note
	----
	No date and time conversion will be attempted, the input needs to be
	converted before, `dateutil` or `pandas` provide convenient functions.

	The local solar time is calculated from time and longitude, except when
	`lst` is set, then that value is used.

	The solar and geomagnetic indices have to be provided, so far the values
	are not included in the module.
	"""
	year = time.year
	doy = int(time.strftime("%j"))
	sec = (time.hour * 3600.
			+ time.minute * 60.
			+ time.second
			+ time.microsecond * 1e-6)
	if lst is None:
		lst = sec / 3600. + lon / 15.0

	kwargs = {}
	if ap_a is not None:
		kwargs.update({"ap_a": ap_a})
	if flags is not None:
		kwargs.update({"flags": flags})

	if method == "gtd7d":
		return gtd7d(year, doy, sec, alt, lat, lon, lst, f107a, f107, ap, **kwargs)
	return gtd7(year, doy, sec, alt, lat, lon, lst, f107a, f107, ap, **kwargs)


def _msise_flat(*args, **kwargs):
	ds, ts = msise_model(*args, **kwargs)
	return np.asarray(ds + ts)


_msise_flatv = np.vectorize(_msise_flat,
		signature='(),(),(),(),(),(),(),(),(),()->(n)',
		excluded=["method"])


@_doc_param(msise_model.__doc__.replace("Interface", "interface"))
def msise_flat(*args, **kwargs):
	"""Flattened {0}
	Additional note for the `flat` version
	--------------------------------------
	This flattened version returns a single 11-element array instead of two
	separate lists."""
	# Set keyword arguments to None if not given to make `np.vectorize` happy
	lst = kwargs.pop("lst", None)
	ap_a = kwargs.pop("ap_a", None)
	flags = kwargs.pop("flags", None)

	return _msise_flatv(*args, lst=lst, ap_a=ap_a, flags=flags, **kwargs)


def scale_height(alt, lat, molw, temp):
	"""Atmospheric scale height

	Extracted from the C-code for easy access, with
	constants updated to the standard SI values.
	It is reasonably fast and :mod:`numpy`-broadcasting compatible.

	Parameters
	----------
	alt: float or array_like
		Altitude in [km].
	lat: float or array_like
		Geodetic latitude in [degrees N].
	molw: float or array_like
		Molecular mass at alt in [kg / mol].
		Can be derived by dividing the total mass density by
		the sum of number densities, e.g. from :func:`gtd7()`,
		multiplied by Avogadro's constant (6.0221407e23 / mol).
	temp: float or array_like
		Temperature at alt in [K].

	Returns
	-------
	scale_height: float or array_like
		Scale height in [m].
	"""
	# Rgas = 8.31446261815324  # J / K / mol
	# dgtr = 1.7453292519943295e-2  # rad / deg
	c2 = np.cos(2. * 1.7453292519943295e-2 * lat)
	gsurf = 9.80665 * (1.0 - 0.0026373 * c2)  # m / s2
	# re in [km]
	re = 2. * gsurf / (3.085462e-6 + 2.27e-9 * c2) * 1.0e-3
	g = gsurf / (1 + alt / re)**2
	return 8.31446261815324 * temp / (g * molw)


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
151			The local solar time is calculated from time and longitude, except when
152			`lst` is set, then that value is used.
153
154			The solar and geomagnetic indices have to be provided, so far the values
155			are not included in the module.
156			"""
157			year = time.year
158			doy = int(time.strftime("%j"))
159			sec = (time.hour * 3600.
160			+ time.minute * 60.
161			+ time.second
162			+ time.microsecond * 1e-6)
163			if lst is None:
164			lst = sec / 3600. + lon / 15.0
165
166			kwargs = {}
167			if ap_a is not None:
168			kwargs.update({"ap_a": ap_a})
169			if flags is not None:
170			kwargs.update({"flags": flags})
171
172			if method == "gtd7d":
173			return gtd7d(year, doy, sec, alt, lat, lon, lst, f107a, f107, ap, **kwargs)
174			return gtd7(year, doy, sec, alt, lat, lon, lst, f107a, f107, ap, **kwargs)
175
176
177			def _msise_flat(args, *kwargs):
178			ds, ts = msise_model(args, *kwargs)
179			return np.asarray(ds + ts)
180
181
182			_msise_flatv = np.vectorize(_msise_flat,
183			signature='(),(),(),(),(),(),(),(),(),()->(n)',
184			excluded=["method"])
185
186
187			@_doc_param(msise_model.__doc__.replace("Interface", "interface"))
188			def msise_flat(args, *kwargs):
189			"""Flattened {0}
190			Additional note for the `flat` version
191			--------------------------------------
192			This flattened version returns a single 11-element array instead of two
193			separate lists."""
194			# Set keyword arguments to None if not given to make `np.vectorize` happy
195			lst = kwargs.pop("lst", None)
196			ap_a = kwargs.pop("ap_a", None)
197			flags = kwargs.pop("flags", None)
198
199			return _msise_flatv(args, lst=lst, ap_a=ap_a, flags=flags, *kwargs)
200
201
202			def scale_height(alt, lat, molw, temp):
203			"""Atmospheric scale height
204
205			Extracted from the C-code for easy access, with
206			constants updated to the standard SI values.
207			It is reasonably fast and :mod:`numpy`-broadcasting compatible.
208
209			Parameters
210			----------
211			alt: float or array_like
212			Altitude in [km].
213			lat: float or array_like
214			Geodetic latitude in [degrees N].
215			molw: float or array_like
216			Molecular mass at alt in [kg / mol].
217			Can be derived by dividing the total mass density by
218			the sum of number densities, e.g. from :func:`gtd7()`,
219			multiplied by Avogadro's constant (6.0221407e23 / mol).
220			temp: float or array_like
221			Temperature at alt in [K].
222
223			Returns
224			-------
225			scale_height: float or array_like
226			Scale height in [m].
227			"""
228			# Rgas = 8.31446261815324 # J / K / mol
229			# dgtr = 1.7453292519943295e-2 # rad / deg
230			c2 = np.cos(2. * 1.7453292519943295e-2 * lat)
231			gsurf = 9.80665 * (1.0 - 0.0026373 * c2) # m / s2
232			# re in [km]
233			re = 2. * gsurf / (3.085462e-6 + 2.27e-9 * c2) * 1.0e-3
234			g = gsurf / (1 + alt / re)**2
235			return 8.31446261815324 * temp / (g * molw)
236

st-bender / pynrlmsise00

Branch — master (dde1c2)

nrlmsise00.core.scale_height() A

Complexity

Size

Duplication

Importance

Duplication Side-by-Side

Filter issues like