eppaurora.spectra.pflux_gaussian() - Code Metrics - st-bender/pyeppaurora - Measure and Improve Code Quality continuously with Scrutinizer

eppaurora.spectra.pflux_gaussian() A
last analyzed 2024-04-05 18:05 UTC

↳ Parent: eppaurora.spectra

Complexity

Conditions

Size

Total Lines	39
Code Lines	2

Duplication

Lines	0
Ratio	0 %

Importance

Changes

Metric	Value
cc	1
eloc	2
nop	3
dl	0
loc	39
rs	10
c	0
b	0
f	0

# coding: utf-8
# Copyright (c) 2020 Stefan Bender
#
# This file is part of pyeppaurora.
# pyeppaurora 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.
"""Particle precipitation spectra

Includes variants describing a normalized particle flux,
as well as variants describing a normalized energy flux.
"""

import numpy as np

__all__ = [
	"exp_general",
	"gaussian_general",
	"maxwell_general",
	"pow_general",
	"pflux_exp",
	"pflux_gaussian",
	"pflux_maxwell",
	"pflux_pow",
	"ediss_spec_int",
	"ediss_specfun_int",
]


# General normalized spectra, standard distributions
def exp_general(en, en_0=10.):
	r"""Exponential number flux spectrum

	.. math::
		\phi(E, E_0) = 1 / E_0 \cdot \exp\{-E / E_0\}

	Standard exponential distribution with
	:math:`\lambda` = 1 / ``en_0`` or :math:`\beta` = ``en_0``.
	normalized to unit number flux, i.e.
	:math:`\int_0^\infty \phi(E) \text{d}E = 1`.

	Parameters
	----------
	en: float or array_like (N,)
		Energy in [keV]
	en_0: float, optional
		Characteristic energy in [keV] of the distribution.
		Default: 10 keV

	Returns
	-------
	phi: float or array_like (N,)
		Normalized differential hemispherical number flux at `en` in [keV-1 cm-2 s-1]
		([keV] or scaled by 1 keV-2 cm-2 s-1, e.g.).
	"""
	return 1. / en_0 * np.exp(-en / en_0)


def gaussian_general(en, en_0=10., w=1.):
	r"""Gaussian number flux spectrum

	Standard normal distribution with
	:math:`\mu` = ``en_0`` and :math:`\sigma` = ``w`` / sqrt(2):

	.. math::
		\phi(E, E_0, W) = 1 / \sqrt{\pi}W \cdot \exp\{-(E - E_0)^2 / W^2\}

	Almost normalized to unit number flux
	:math:`\int_0^\infty \phi(E) \text{d}E = 1`
	(ignoring the negative tail for large ``en_0`` / ``w`` ratios).

	Parameters
	----------
	en: float or array_like (N,)
		Energy in [keV]
	en_0: float, optional
		Characteristic energy in [keV], i.e. mode of the distribution.
		Default: 10 keV
	w: float, optional
		Width of the Gaussian distribution, in [keV].
		Default: 1 keV

	Returns
	-------
	phi: float or array_like (N,)
		Normalized differential hemispherical number flux at `en` in [keV-1 cm-2 s-1]
		([keV] or scaled by 1 keV-2 cm-2 s-1, e.g.).
	"""
	return 1. / np.sqrt(np.pi * w**2) * np.exp(-(en - en_0)**2 / w**2)


def maxwell_general(en, en_0=10.):
	r"""Maxwell number flux spectrum

	.. math::
		\phi(E, E_0) = E / E_0^2 \cdot \exp\{-E / E_0\}

	Equal to a standard Gamma distribution with
	:math:`\alpha` = 2 and :math:`\beta` = 1 / ``en_0``,
	or
	:math:`k` = 2 and :math:`\theta` = ``en_0``.
	Normalized to :math:`\int_0^\infty \phi(E) \text{d}E = 1`.

	Parameters
	----------
	en: float or array_like (N,)
		Energy in [keV]
	en_0: float, optional
		Characteristic energy in [keV], i.e. mode of the distribution.
		Default: 10 keV

	Returns
	-------
	phi: float or array_like (N,)
		Normalized differential hemispherical number flux at `en` in [keV-1 cm-2 s-1]
		([keV] or scaled by 1 keV-2 cm-2 s-1, e.g.).
	"""
	return en / en_0**2 * np.exp(-en / en_0)


def pow_general(en, en_0=10., gamma=-3., het=True):
	r"""Power-law number flux spectrum

	.. math::
		\phi(E, E_0, \gamma) = \mp (\gamma + 1) / E_0 \cdot (E / E_0)^\gamma

	The minus-sign (-) and is used for the high-energy tail variant,
	and the plus-sign (+) for the low-energy tail variant.
	The exponent ``gamma`` needs to be set appropriately,
	< -1 for `het`, and > 1 for `let`.

	The "high-energy tail" version (`het` = `True`)
	resembles a Pareto distribution with
	scale parameter :math:`x_m` = ``en_0``
	and shape parameter :math:`\alpha` = -(``gamma`` + 1).

	Adapted from Strickland et al., 1993 [#]_ and
	normalized to unit particle flux:
	:math:`\int_{E_0}^\infty \phi(E) \text{d}E = 1`
	for the high-energy tail version, and
	:math:`\int_0^{E_0} \phi(E) \text{d}E = 1`
	for the low-energy tail version.

	Parameters
	----------
	en: float or array_like (N,)
		Energy in [keV]
	en_0: float, optional
		Characteristic energy in [keV], i.e. mode of the distribution.
		Default: 10 keV
	gamma: float, optional
		Exponent of the power-law distribution, in [keV].
	het: bool, optional
		Return a high-energy tail (het, default: true) for en > en_0,
		or low-energy tail (false) for en < en_0.
		Adjusts the normalization accordingly.

	Returns
	-------
	phi: float or array_like (N,)
		Normalized differential hemispherical number flux at `en` in [keV-1 cm-2 s-1]
		([keV] or scaled by 1 keV-2 cm-2 s-1, e.g.).

	References
	----------
	.. [#] D. J. Strickland, R. E. Daniell Jr, J. R. Jasperse, B. Basu
		J. Geophys. Res., 98(A12), pp. 21533--21548, 1993
		doi: `10.1029/93JA01645 <https://doi.org/10.1029/93JA01645>`_
	"""
	isscalar = (np.ndim(en) == 0)
	en = np.atleast_1d(en)
	spec = (gamma + 1) / en_0 * (en / en_0)**gamma
	if het:
		spec[en < en_0] = 0.
		return -spec[0] if isscalar else -spec
	spec[en > en_0] = 0.
	return spec[0] if isscalar else spec


def pflux_exp(en, en_0=10.):
	r"""Exponential particle flux spectrum

	.. math::
		\phi(E, E_0) = 1 / E_0^2 \cdot \exp\{-E / E_0\}

	Normalized to unit energy flux:
	:math:`\int_0^\infty \phi(E) E \text{d}E = 1`.

	Scales to arbitrary energy flux :math:`Q` via multiplication:
	:math:`\tilde\phi = Q \cdot \phi`.

	Parameters
	----------
	en: float or array_like (N,)
		Energy in [keV]
	en_0: float, optional
		Characteristic energy in [keV], i.e. mode of the distribution.
		Default: 10 keV.

	Returns
	-------
	phi: float or array_like (N,)
		Hemispherical differential particle flux at `en` in [keV-1 cm-2 s-1]
		([keV-2] scaled by unit energy flux).

	See Also
	--------
	exp_general
	"""
	return exp_general(en, en_0=en_0) / en_0


def pflux_gaussian(en, en_0=10., w=1):
	r"""Gaussian particle flux spectrum

	As used in, e.g., Strickland et al., 1993 [#]_

	.. math::
		\phi(E, E_0, W) = 1 / \sqrt{\pi}E_0W \cdot \exp\{-(E - E_0)^2 / W^2\}

	Normalized to :math:`\int_0^\infty \phi(E) E \text{d}E = 1`
	(ignoring the negative tail).

	Scales to arbitrary energy flux :math:`Q` via multiplication:
	:math:`\tilde\phi = Q \cdot \phi`.

	Parameters
	----------
	en: float or array_like (N,)
		Energy in [keV]
	en_0: float, optional
		Characteristic energy in [keV], i.e. mode of the distribution.
		Default: 10 keV.

	Returns
	-------
	phi: float or array_like (N,)
		Hemispherical differential particle flux at `en` in [keV-1 cm-2 s-1]
		([kev-2] scaled by unit energy flux).

	References
	----------
	.. [#] D. J. Strickland, R. E. Daniell, J. R. Jasperse, B. Basu
		J. Geophys. Res., 98(A12), pp. 21533--21548, 1993
		doi: 10.1029/93JA01645

	See Also
	--------
	gaussian_general
	"""
	return gaussian_general(en, en_0=en_0, w=w) / en_0


def pflux_maxwell(en, en_0=10.):
	r"""Maxwell particle flux spectrum

	As used in, e.g., Strickland et al., 1993 [#]_

	.. math::
		\phi(E, E_0) = E / 2E_0^3 \cdot \exp\{-E / E_0\}

	Equal to a standard Gamma distribution with
	:math:`\alpha` = 3 and :math:`\beta` = 1 / ``en_0``,
	or
	:math:`k` = 3 and :math:`\theta` = ``en_0``.
	The total precipitating energy flux is fixed to 1 keV cm-2 s-1,
	multiply by Q_0 [keV cm-2 s-1] to scale the particle flux.

	Normalized to :math:`\int_0^\infty \phi(E) E \text{d}E = 1`.

	Scales to arbitrary energy flux :math:`Q` via multiplication:
	:math:`\tilde\phi = Q \cdot \phi`.

	Parameters
	----------
	en: float or array_like (N,)
		Energy in [keV]
	en_0: float, optional
		Characteristic energy in [keV], i.e. mode of the distribution.
		Default: 10 keV.

	Returns
	-------
	phi: float or array_like (N,)
		Hemispherical differential particle flux at `en` in [keV-1 cm-2 s-1]
		([kev-2] scaled by unit energy flux).

	References
	----------
	.. [#] D. J. Strickland, R. E. Daniell, J. R. Jasperse, B. Basu
		J. Geophys. Res., 98(A12), pp. 21533--21548, 1993
		doi: 10.1029/93JA01645

	See Also
	--------
	maxwell_general
	"""
	return 0.5 / en_0 * maxwell_general(en, en_0)


def pflux_pow(en, en_0=10., gamma=-3., het=True):
	r"""Power-law particle flux spectrum

	As used in, e.g., Strickland et al., 1993 [#]_

	.. math::
		\phi(E, E_0, \gamma) = \mp (\gamma + 2) / E_0^2 \cdot (E / E_0)^\gamma

	The minus-sign (-) and is used for the high-energy tail variant,
	and the plus-sign (+) for the low-energy tail variant.
	The exponent ``gamma`` needs to be set appropriately,
	< -1 for `het`, and > 1 for `let`.

	Normalized to :math:`\int_{E_0}^\infty \phi(E) E \text{d}E = 1`
	for the high-energy tail version, and to
	:math:`\int_0^{E_0} \phi(E) E \text{d}E = 1`
	for the low-energy tail version.

	Scales to arbitrary energy flux :math:`Q` via multiplication:
	:math:`\tilde\phi = Q \cdot \phi`.

	Parameters
	----------
	en: float or array_like (N,)
		Energy in [keV]
	en_0: float, optional
		Characteristic energy in [keV], i.e. mode of the distribution.
		Default: 10 keV
	gamma: float, optional
		Exponent of the power-law distribution, in [keV].
	het: bool, optional (default True)
		Return a high-energy tail (true) for en > en_0,
		or low-energy tail (false) for en < en_0.
		Adjusts the normalization accordingly.

	Returns
	-------
	phi: float or array_like (N,)
		Hemispherical differential particle flux at `en` in [keV-1 cm-2 s-1]
		([keV-2] scaled by unit energy flux).

	References
	----------
	.. [#] D. J. Strickland, R. E. Daniell Jr, J. R. Jasperse, B. Basu
		J. Geophys. Res., 98(A12), pp. 21533--21548, 1993
		doi: `10.1029/93JA01645 <https://doi.org/10.1029/93JA01645>`_

	See Also
	--------
	pow_general
	"""
	return (gamma + 2) / (gamma + 1) / en_0 * pow_general(en, en_0=en_0, gamma=gamma, het=het)


def ediss_spec_int(
	ens,
	dfluxes,
	scale_height,
	rho,
	func,
	axis=-1,
	func_kws=None,
):
	r"""Integrate over a given energy spectrum

	Integrates a mono-energetic parametrization `q`, e.g. from Fang et al., 2010
	using the given differential particle spectrum `phi`:

	:math:`\int_\text{spec} \phi(E) q(E, Q) E \text{d}E`

	This function uses the differential spectrum evaluated at the given energy bins.

	Parameters
	----------
	ens: array_like (M,...)
		Central (bin) energies of the spectrum
	dfluxes: array_like (M,...)
		Differential particle fluxes in the given bins
	scale_height: array_like (N,...)
		The atmospheric scale heights
	rho: array_like (N,...)
		The atmospheric densities, corresponding to the
		scale heights.
	func: callable
		Mono-energetic energy dissipation function to integrate.
	axis: int, optional
		The axis to use for integration, default: -1 (last axis).
	func_kws: dict-like, optional
		Optional keyword arguments to pass to the mono-energetic
		energy dissipation function. Default: `None`

	Returns
	-------
	en_diss: array_like (N)
		The dissipated energy profiles [keV].

	See Also
	--------
	ediss_specfun_int
	"""
	ens = np.atleast_1d(ens)
	dfluxes = np.atleast_1d(dfluxes)
	scale_height = np.atleast_1d(scale_height)
	rho = np.atleast_1d(rho)
	func_kws = func_kws or dict()
	ediss = func(
		ens[None, None, :],
		dfluxes,
		scale_height[..., None],
		rho[..., None],
		**func_kws
	)
	return np.trapz(ediss * ens, ens, axis=axis)


def ediss_specfun_int(
	energy,
	flux,
	scale_height,
	rho,
	ediss_func,
	ediss_kws=None,
	bounds=(0.1, 300.),
	nstep=128,
	spec_fun=pflux_maxwell,
	spec_kws=None,
):
	"""Integrate mono-energetic parametrization over a spectrum

	Integrates the mono-energetic parametrization over a spectrum given by a
	functional dependence with characteristic energy `energy` and total energy
	flux `flux`.

	Parameters
	----------
	energy: float or array_like (M,...)
		Characteristic energy E_0 [keV] of the spectral distribution.
	flux: float or array_like (M,...)
		Integrated energy flux Q_0 [keV / cm² / s¹]
	scale_height: float or array_like (N,...)
		The atmospheric scale heights [cm].
	rho: float or array_like (N,...)
		The atmospheric mass density [g / cm³]
	ediss_func: callable
		Mono-energetic energy dissipation function to integrate.
	ediss_kws: dict-like, optional
		Optional keyword arguments to pass to the mono-energetic
		energy dissipation function. Default: `None`
	bounds: tuple, optional
		(min, max) [keV] of the integration range to integrate the Maxwellian.
		Make sure that this is appropriate to encompass the spectrum.
		Default: (0.1, 300.)
	nsteps: int, optional
		Number of integration steps, default: 128.
	spec_func: callable, optional, default :func:`pflux_maxwell`
		Spectral shape function, choices are:

		* :func:`pflux_exp` for a exponential spectrum
		* :func:`pflux_gaussian` for a Gaussian shaped spectrum
		* :func:`pflux_maxwell` for a Maxwellian shaped spectrum
		* :func:`pflux_pow` for a power-law
	spec_kws: dict-like, optional
		Optional keyword arguments to pass to the spectral function
		Default: `None`

	Returns
	-------
	en_diss: array_like (M,N)
		The dissipated energy profiles [keV].

	See Also
	--------
	ediss_spec_int
	"""
	energy = np.asarray(energy)
	flux = np.asarray(flux)
	bounds_l10 = np.log10(bounds)
	ens = np.logspace(*bounds_l10, num=nstep)
	ensd = np.reshape(ens, (-1,) + (1,) * energy.ndim)
	spec_kws = spec_kws or dict()
	# "overwrite" the characteristic energy
	spec_kws["en_0"] = energy.T
	dflux = flux.T * spec_fun(ensd, **spec_kws)
	return ediss_spec_int(
		ens, dflux.T, scale_height, rho, ediss_func,
		axis=-1, func_kws=ediss_kws,
	)


1			# coding: utf-8
2			# Copyright (c) 2020 Stefan Bender
3			#
4			# This file is part of pyeppaurora.
5			# pyeppaurora is free software: you can redistribute it or modify
6			# it under the terms of the GNU General Public License as published
7			# by the Free Software Foundation, version 2.
8			# See accompanying LICENSE file or http://www.gnu.org/licenses/gpl-2.0.html.
9			"""Particle precipitation spectra
10
11			Includes variants describing a normalized particle flux,
12			as well as variants describing a normalized energy flux.
13			"""
14
15			import numpy as np
16
17			__all__ = [
18			"exp_general",
19			"gaussian_general",
20			"maxwell_general",
21			"pow_general",
22			"pflux_exp",
23			"pflux_gaussian",
24			"pflux_maxwell",
25			"pflux_pow",
26			"ediss_spec_int",
27			"ediss_specfun_int",
28			]
29
30
31			# General normalized spectra, standard distributions
32			def exp_general(en, en_0=10.):
33			r"""Exponential number flux spectrum
34
35			.. math::
36			\phi(E, E_0) = 1 / E_0 \cdot \exp\{-E / E_0\}
37
38			Standard exponential distribution with
39			:math:`\lambda` = 1 / ``en_0`` or :math:`\beta` = ``en_0``.
40			normalized to unit number flux, i.e.
41			:math:`\int_0^\infty \phi(E) \text{d}E = 1`.
42
43			Parameters
44			----------
45			en: float or array_like (N,)
46			Energy in [keV]
47			en_0: float, optional
48			Characteristic energy in [keV] of the distribution.
49			Default: 10 keV
50
51			Returns
52			-------
53			phi: float or array_like (N,)
54			Normalized differential hemispherical number flux at `en` in [keV-1 cm-2 s-1]
55			([keV] or scaled by 1 keV-2 cm-2 s-1, e.g.).
56			"""
57			return 1. / en_0 * np.exp(-en / en_0)
58
59
60			def gaussian_general(en, en_0=10., w=1.):
61			r"""Gaussian number flux spectrum
62
63			Standard normal distribution with
64			:math:`\mu` = ``en_0`` and :math:`\sigma` = ``w`` / sqrt(2):
65
66			.. math::
67			\phi(E, E_0, W) = 1 / \sqrt{\pi}W \cdot \exp\{-(E - E_0)^2 / W^2\}
68
69			Almost normalized to unit number flux
70			:math:`\int_0^\infty \phi(E) \text{d}E = 1`
71			(ignoring the negative tail for large ``en_0`` / ``w`` ratios).
72
73			Parameters
74			----------
75			en: float or array_like (N,)
76			Energy in [keV]
77			en_0: float, optional
78			Characteristic energy in [keV], i.e. mode of the distribution.
79			Default: 10 keV
80			w: float, optional
81			Width of the Gaussian distribution, in [keV].
82			Default: 1 keV
83
84			Returns
85			-------
86			phi: float or array_like (N,)
87			Normalized differential hemispherical number flux at `en` in [keV-1 cm-2 s-1]
88			([keV] or scaled by 1 keV-2 cm-2 s-1, e.g.).
89			"""
90			return 1. / np.sqrt(np.pi * w*2) np.exp(-(en - en_0)2 / w2)
91
92
93			def maxwell_general(en, en_0=10.):
94			r"""Maxwell number flux spectrum
95
96			.. math::
97			\phi(E, E_0) = E / E_0^2 \cdot \exp\{-E / E_0\}
98
99			Equal to a standard Gamma distribution with
100			:math:`\alpha` = 2 and :math:`\beta` = 1 / ``en_0``,
101			or
102			:math:`k` = 2 and :math:`\theta` = ``en_0``.
103			Normalized to :math:`\int_0^\infty \phi(E) \text{d}E = 1`.
104
105			Parameters
106			----------
107			en: float or array_like (N,)
108			Energy in [keV]
109			en_0: float, optional
110			Characteristic energy in [keV], i.e. mode of the distribution.
111			Default: 10 keV
112
113			Returns
114			-------
115			phi: float or array_like (N,)
116			Normalized differential hemispherical number flux at `en` in [keV-1 cm-2 s-1]
117			([keV] or scaled by 1 keV-2 cm-2 s-1, e.g.).
118			"""
119			return en / en_0*2 np.exp(-en / en_0)
120
121
122			def pow_general(en, en_0=10., gamma=-3., het=True):
123			r"""Power-law number flux spectrum
124
125			.. math::
126			\phi(E, E_0, \gamma) = \mp (\gamma + 1) / E_0 \cdot (E / E_0)^\gamma
127
128			The minus-sign (-) and is used for the high-energy tail variant,
129			and the plus-sign (+) for the low-energy tail variant.
130			The exponent ``gamma`` needs to be set appropriately,
131			< -1 for `het`, and > 1 for `let`.
132
133			The "high-energy tail" version (`het` = `True`)
134			resembles a Pareto distribution with
135			scale parameter :math:`x_m` = ``en_0``
136			and shape parameter :math:`\alpha` = -(``gamma`` + 1).
137
138			Adapted from Strickland et al., 1993 [#]_ and
139			normalized to unit particle flux:
140			:math:`\int_{E_0}^\infty \phi(E) \text{d}E = 1`
141			for the high-energy tail version, and
142			:math:`\int_0^{E_0} \phi(E) \text{d}E = 1`
143			for the low-energy tail version.
144
145			Parameters
146			----------
147			en: float or array_like (N,)
148			Energy in [keV]
149			en_0: float, optional
150			Characteristic energy in [keV], i.e. mode of the distribution.
151			Default: 10 keV
152			gamma: float, optional
153			Exponent of the power-law distribution, in [keV].
154			het: bool, optional
155			Return a high-energy tail (het, default: true) for en > en_0,
156			or low-energy tail (false) for en < en_0.
157			Adjusts the normalization accordingly.
158
159			Returns
160			-------
161			phi: float or array_like (N,)
162			Normalized differential hemispherical number flux at `en` in [keV-1 cm-2 s-1]
163			([keV] or scaled by 1 keV-2 cm-2 s-1, e.g.).
164
165			References
166			----------
167			.. [#] D. J. Strickland, R. E. Daniell Jr, J. R. Jasperse, B. Basu
168			J. Geophys. Res., 98(A12), pp. 21533--21548, 1993
169			doi: `10.1029/93JA01645 <https://doi.org/10.1029/93JA01645>`_
170			"""
171			isscalar = (np.ndim(en) == 0)
172			en = np.atleast_1d(en)
173			spec = (gamma + 1) / en_0 * (en / en_0)**gamma
174			if het:
175			spec[en < en_0] = 0.
176			return -spec[0] if isscalar else -spec
177			spec[en > en_0] = 0.
178			return spec[0] if isscalar else spec
179
180
181			def pflux_exp(en, en_0=10.):
182			r"""Exponential particle flux spectrum
183
184			.. math::
185			\phi(E, E_0) = 1 / E_0^2 \cdot \exp\{-E / E_0\}
186
187			Normalized to unit energy flux:
188			:math:`\int_0^\infty \phi(E) E \text{d}E = 1`.
189
190			Scales to arbitrary energy flux :math:`Q` via multiplication:
191			:math:`\tilde\phi = Q \cdot \phi`.
192
193			Parameters
194			----------
195			en: float or array_like (N,)
196			Energy in [keV]
197			en_0: float, optional
198			Characteristic energy in [keV], i.e. mode of the distribution.
199			Default: 10 keV.
200
201			Returns
202			-------
203			phi: float or array_like (N,)
204			Hemispherical differential particle flux at `en` in [keV-1 cm-2 s-1]
205			([keV-2] scaled by unit energy flux).
206
207			See Also
208			--------
209			exp_general
210			"""
211			return exp_general(en, en_0=en_0) / en_0
212
213
214			def pflux_gaussian(en, en_0=10., w=1):
215			r"""Gaussian particle flux spectrum
216
217			As used in, e.g., Strickland et al., 1993 [#]_
218
219			.. math::
220			\phi(E, E_0, W) = 1 / \sqrt{\pi}E_0W \cdot \exp\{-(E - E_0)^2 / W^2\}
221
222			Normalized to :math:`\int_0^\infty \phi(E) E \text{d}E = 1`
223			(ignoring the negative tail).
224
225			Scales to arbitrary energy flux :math:`Q` via multiplication:
226			:math:`\tilde\phi = Q \cdot \phi`.
227
228			Parameters
229			----------
230			en: float or array_like (N,)
231			Energy in [keV]
232			en_0: float, optional
233			Characteristic energy in [keV], i.e. mode of the distribution.
234			Default: 10 keV.
235
236			Returns
237			-------
238			phi: float or array_like (N,)
239			Hemispherical differential particle flux at `en` in [keV-1 cm-2 s-1]
240			([kev-2] scaled by unit energy flux).
241
242			References
243			----------
244			.. [#] D. J. Strickland, R. E. Daniell, J. R. Jasperse, B. Basu
245			J. Geophys. Res., 98(A12), pp. 21533--21548, 1993
246			doi: 10.1029/93JA01645
247
248			See Also
249			--------
250			gaussian_general
251			"""
252			return gaussian_general(en, en_0=en_0, w=w) / en_0
253
254
255			def pflux_maxwell(en, en_0=10.):
256			r"""Maxwell particle flux spectrum
257
258			As used in, e.g., Strickland et al., 1993 [#]_
259
260			.. math::
261			\phi(E, E_0) = E / 2E_0^3 \cdot \exp\{-E / E_0\}
262
263			Equal to a standard Gamma distribution with
264			:math:`\alpha` = 3 and :math:`\beta` = 1 / ``en_0``,
265			or
266			:math:`k` = 3 and :math:`\theta` = ``en_0``.
267			The total precipitating energy flux is fixed to 1 keV cm-2 s-1,
268			multiply by Q_0 [keV cm-2 s-1] to scale the particle flux.
269
270			Normalized to :math:`\int_0^\infty \phi(E) E \text{d}E = 1`.
271
272			Scales to arbitrary energy flux :math:`Q` via multiplication:
273			:math:`\tilde\phi = Q \cdot \phi`.
274
275			Parameters
276			----------
277			en: float or array_like (N,)
278			Energy in [keV]
279			en_0: float, optional
280			Characteristic energy in [keV], i.e. mode of the distribution.
281			Default: 10 keV.
282
283			Returns
284			-------
285			phi: float or array_like (N,)
286			Hemispherical differential particle flux at `en` in [keV-1 cm-2 s-1]
287			([kev-2] scaled by unit energy flux).
288
289			References
290			----------
291			.. [#] D. J. Strickland, R. E. Daniell, J. R. Jasperse, B. Basu
292			J. Geophys. Res., 98(A12), pp. 21533--21548, 1993
293			doi: 10.1029/93JA01645
294
295			See Also
296			--------
297			maxwell_general
298			"""
299			return 0.5 / en_0 * maxwell_general(en, en_0)
300
301
302			def pflux_pow(en, en_0=10., gamma=-3., het=True):
303			r"""Power-law particle flux spectrum
304
305			As used in, e.g., Strickland et al., 1993 [#]_
306
307			.. math::
308			\phi(E, E_0, \gamma) = \mp (\gamma + 2) / E_0^2 \cdot (E / E_0)^\gamma
309
310			The minus-sign (-) and is used for the high-energy tail variant,
311			and the plus-sign (+) for the low-energy tail variant.
312			The exponent ``gamma`` needs to be set appropriately,
313			< -1 for `het`, and > 1 for `let`.
314
315			Normalized to :math:`\int_{E_0}^\infty \phi(E) E \text{d}E = 1`
316			for the high-energy tail version, and to
317			:math:`\int_0^{E_0} \phi(E) E \text{d}E = 1`
318			for the low-energy tail version.
319
320			Scales to arbitrary energy flux :math:`Q` via multiplication:
321			:math:`\tilde\phi = Q \cdot \phi`.
322
323			Parameters
324			----------
325			en: float or array_like (N,)
326			Energy in [keV]
327			en_0: float, optional
328			Characteristic energy in [keV], i.e. mode of the distribution.
329			Default: 10 keV
330			gamma: float, optional
331			Exponent of the power-law distribution, in [keV].
332			het: bool, optional (default True)
333			Return a high-energy tail (true) for en > en_0,
334			or low-energy tail (false) for en < en_0.
335			Adjusts the normalization accordingly.
336
337			Returns
338			-------
339			phi: float or array_like (N,)
340			Hemispherical differential particle flux at `en` in [keV-1 cm-2 s-1]
341			([keV-2] scaled by unit energy flux).
342
343			References
344			----------
345			.. [#] D. J. Strickland, R. E. Daniell Jr, J. R. Jasperse, B. Basu
346			J. Geophys. Res., 98(A12), pp. 21533--21548, 1993
347			doi: `10.1029/93JA01645 <https://doi.org/10.1029/93JA01645>`_
348
349			See Also
350			--------
351			pow_general
352			"""
353			return (gamma + 2) / (gamma + 1) / en_0 * pow_general(en, en_0=en_0, gamma=gamma, het=het)
354
355
356			def ediss_spec_int(
357			ens,
358			dfluxes,
359			scale_height,
360			rho,
361			func,
362			axis=-1,
363			func_kws=None,
364			):
365			r"""Integrate over a given energy spectrum
366
367			Integrates a mono-energetic parametrization `q`, e.g. from Fang et al., 2010
368			using the given differential particle spectrum `phi`:
369
370			:math:`\int_\text{spec} \phi(E) q(E, Q) E \text{d}E`
371
372			This function uses the differential spectrum evaluated at the given energy bins.
373
374			Parameters
375			----------
376			ens: array_like (M,...)
377			Central (bin) energies of the spectrum
378			dfluxes: array_like (M,...)
379			Differential particle fluxes in the given bins
380			scale_height: array_like (N,...)
381			The atmospheric scale heights
382			rho: array_like (N,...)
383			The atmospheric densities, corresponding to the
384			scale heights.
385			func: callable
386			Mono-energetic energy dissipation function to integrate.
387			axis: int, optional
388			The axis to use for integration, default: -1 (last axis).
389			func_kws: dict-like, optional
390			Optional keyword arguments to pass to the mono-energetic
391			energy dissipation function. Default: `None`
392
393			Returns
394			-------
395			en_diss: array_like (N)
396			The dissipated energy profiles [keV].
397
398			See Also
399			--------
400			ediss_specfun_int
401			"""
402			ens = np.atleast_1d(ens)
403			dfluxes = np.atleast_1d(dfluxes)
404			scale_height = np.atleast_1d(scale_height)
405			rho = np.atleast_1d(rho)
406			func_kws = func_kws or dict()
407			ediss = func(
408			ens[None, None, :],
409			dfluxes,
410			scale_height[..., None],
411			rho[..., None],
412			**func_kws
413			)
414			return np.trapz(ediss * ens, ens, axis=axis)
415
416
417			def ediss_specfun_int(
418			energy,
419			flux,
420			scale_height,
421			rho,
422			ediss_func,
423			ediss_kws=None,
424			bounds=(0.1, 300.),
425			nstep=128,
426			spec_fun=pflux_maxwell,
427			spec_kws=None,
428			):
429			"""Integrate mono-energetic parametrization over a spectrum
430
431			Integrates the mono-energetic parametrization over a spectrum given by a
432			functional dependence with characteristic energy `energy` and total energy
433			flux `flux`.
434
435			Parameters
436			----------
437			energy: float or array_like (M,...)
438			Characteristic energy E_0 [keV] of the spectral distribution.
439			flux: float or array_like (M,...)
440			Integrated energy flux Q_0 [keV / cm² / s¹]
441			scale_height: float or array_like (N,...)
442			The atmospheric scale heights [cm].
443			rho: float or array_like (N,...)
444			The atmospheric mass density [g / cm³]
445			ediss_func: callable
446			Mono-energetic energy dissipation function to integrate.
447			ediss_kws: dict-like, optional
448			Optional keyword arguments to pass to the mono-energetic
449			energy dissipation function. Default: `None`
450			bounds: tuple, optional
451			(min, max) [keV] of the integration range to integrate the Maxwellian.
452			Make sure that this is appropriate to encompass the spectrum.
453			Default: (0.1, 300.)
454			nsteps: int, optional
455			Number of integration steps, default: 128.
456			spec_func: callable, optional, default :func:`pflux_maxwell`
457			Spectral shape function, choices are:
458
459			* :func:`pflux_exp` for a exponential spectrum
460			* :func:`pflux_gaussian` for a Gaussian shaped spectrum
461			* :func:`pflux_maxwell` for a Maxwellian shaped spectrum
462			* :func:`pflux_pow` for a power-law
463			spec_kws: dict-like, optional
464			Optional keyword arguments to pass to the spectral function
465			Default: `None`
466
467			Returns
468			-------
469			en_diss: array_like (M,N)
470			The dissipated energy profiles [keV].
471
472			See Also
473			--------
474			ediss_spec_int
475			"""
476			energy = np.asarray(energy)
477			flux = np.asarray(flux)
478			bounds_l10 = np.log10(bounds)
479			ens = np.logspace(*bounds_l10, num=nstep)
480			ensd = np.reshape(ens, (-1,) + (1,) * energy.ndim)
481			spec_kws = spec_kws or dict()
482			# "overwrite" the characteristic energy
483			spec_kws["en_0"] = energy.T
484			dflux = flux.T * spec_fun(ensd, **spec_kws)
485			return ediss_spec_int(
486			ens, dflux.T, scale_height, rho, ediss_func,
487			axis=-1, func_kws=ediss_kws,
488			)
489

st-bender / pyeppaurora

eppaurora.spectra.pflux_gaussian() A last analyzed 2024-04-05 18:05 UTC

Complexity

Size

Duplication

Importance

Duplication Side-by-Side

Filter issues like

eppaurora.spectra.pflux_gaussian() A
last analyzed 2024-04-05 18:05 UTC