eppaurora.brems

Complexity

Total Complexity

4

Size/Duplication

Total Lines	149
Duplicated Lines	0 %

Importance

Changes

0

1 Function

Rating	Name	Duplication	Size	Complexity
A	berger1974()	0	101	4

# 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.
"""Atmospheric bremsstrahlung ionization parametrization [1]_

.. [1] Berger, M.J., Seltzer, S.M., Maeda, K.,
	Some new results on electron transport in the atmosphere,
	Journal of Atmospheric and Terrestrial Physics, v36, i4, pp. 591--617,
	April 1974,
	doi: 10.1016/0021-9169(74)90085-3
"""

import numpy as np
from scipy import interpolate

__all__ = ["berger1974"]

E_BR = [2., 5., 10., 20., 50., 100., 200., 500., 1000., 2000.]

Z_BR = [
	2e-6, 4e-6, 8e-6,
	2e-5, 4e-5, 8e-5,
	2e-4, 4e-4, 8e-4,
	2e-3, 4e-3, 8e-3,
	2e-2, 4e-2, 8e-2,
	2e-1, 4e-1,
]

A_BR = [
	[np.nan] * 9 + [8.6e-6, 5.1e-7] + [np.nan] * 6,
	[np.nan] * 8 + [7.2e-4, 1.4e-4, 3.3e-5, 5.2e-6, 1.1e-7] + [np.nan] * 4,
	[np.nan] * 7 + [2.0e-3, 8.8e-4, 2.7e-4, 9.4e-5, 2.8e-5, 3.0e-6, 3.8e-7, 3.9e-8] + [np.nan] * 2,
	[np.nan] * 6 + [3.4e-3, 1.7e-3, 8.0e-4, 3.0e-4, 1.3e-4, 5.7e-5, 1.5e-5, 3.3e-6, 5.7e-7, 2.9e-8] + [np.nan],
	[np.nan] * 5 + [7.3e-3, 2.2e-3, 1.1e-3, 5.7e-4, 2.3e-4, 1.1e-4, 5.4e-5, 2.0e-5, 8.3e-6, 2.6e-6, 3.1e-7, 2.8e-8],
	[np.nan] * 4 + [1.1e-2, 7.9e-3, 1.7e-3, 8.4e-4, 4.4e-4, 1.9e-4, 1.0e-4, 5.4e-5, 2.3e-5, 1.1e-5, 3.9e-6, 5.4e-7, 2.4e-8],
	[np.nan] * 3 + [5.0e-3, 5.3e-3, 5.0e-3, 1.8e-3, 6.5e-4, 3.3e-4, 1.5e-4, 8.8e-5, 5.2e-5, 2.4e-5, 1.1e-5, 3.8e-6, 1.7e-7, 5.9e-9],
	[np.nan] * 2 + [2.2e-3, 2.4e-3, 2.5e-3, 2.5e-3, 1.6e-3, 5.2e-4, 2.8e-4, 1.5e-4, 9.9e-5, 6.5e-5, 2.8e-5, 9.5e-6, 1.5e-6, 6.0e-9] + [np.nan],
	[np.nan] + [1.0e-3, 1.1e-3, 1.2e-3, 1.3e-3, 1.4e-3, 1.2e-3, 5.5e-4, 3.2e-4, 1.9e-4, 1.3e-4, 8.3e-5, 2.8e-5, 5.5e-6, 2.4e-7] + [np.nan] * 2,
	[6.8e-4, 7.6e-4, 8.4e-4, 9.3e-4, 9.9e-4, 1.1e-3, 1.1e-3, 6.8e-4, 4.5e-4, 2.9e-4, 1.9e-4, 9.2e-5, 1.7e-5, 1.2e-6] + [np.nan] * 3,
]


def berger1974(
	energy, flux,
	scale_height, rho,
	ens=None, zm_p_en=None, coeffs=None,
	fillna=None, log3=True,
	rbf="multiquadric",
):
	"""Bremsstrahlung ionization by secondary electrons

	Formulae and parameters as described in [#]_.

	By default, the `log(coefficients)` are interpolated wrt. `log(energy)`
	and `log(zm)` using :class:`scipy.interpolated.Rbf`.
	The default "multiquadric" should work fine, if not consider
	using "thin-plate" splines.

	Parameters
	----------
	energy: array_like (M, ...)
		Energy E_0 of the mono-energetic electron beam [keV].
		A scalar (0-D) value is promoted to 1-D with one element.
	flux: array_like (M,...)
		Energy flux Q_0 of the mono-energetic electron beam [keV / cm² / s¹].
	scale_height: array_like (N, ...)
		The atmospheric scale heights [cm].
	rho: array_like (N, ...)
		The atmospheric mass density [g / cm³].
	ens: array_like (I,), optional
		The energies (one axis) of the coefficient array,
		used to interpolate the coefficients to `energy`.
		Defaults to the Berger et al., 1974 coefficients.
	zm_p_en: array_like (J,), optional
		The atmospheric depth (the other axis) of the coefficient array,
		used to interpolate the coefficients to `z` = `scale_height` * `rhos`.
		Defaults to the Berger et al., 1974 coefficients.
	coeffs: array_like, (J, I), optional
		The bremsstrahlung energy dissipation coefficients.
		Defaults to the Berger et al., 1974 coefficients.
	fillna: float or None, optional (default `None`)
		Value to use for `nan` values in `coeffs`, `None` skips them.
	log3: bool, optional (default `True`)
		Interpolate the coefficients as log(ens)-log(zm)-log(coeff)
		instead of a linear variant.
	rbf: str or callable, optional (default "multiquadric")
		Radial basis functions to use for :class:`scipy.interpolate.Rbf`.

	Returns
	-------
	a_br: array_like (M, N)
		A scalar (0-D) `energy` is promoted to 1-D, and the result will
		have shape (1, N), *not* (N,).
		Energy dissipation rate, units: [keV cm⁻³ s⁻¹]

	References
	----------
	.. [#] Berger, M.J., Seltzer, S.M., Maeda, K.,
		Some new results on electron transport in the atmosphere,
		Journal of Atmospheric and Terrestrial Physics, v36, i4, pp. 591--617,
		April 1974,
		doi: 10.1016/0021-9169(74)90085-3

	See also
	--------
	scipy.interpolate.Rbf
	"""
	energy = np.atleast_1d(np.asarray(energy, dtype=float))
	ens = np.asarray(ens) or np.asarray(E_BR)
	zm_p_en = np.asarray(zm_p_en) or np.asarray(Z_BR)

	coeffs = (
		np.array(coeffs, copy=fillna is not None)
		or np.array(A_BR, copy=fillna is not None)
	)
	nans = np.isnan(coeffs)
	if fillna is not None:
		coeffs[nans] = fillna
		# update nan positions (should be all False)
		nans = np.isnan(coeffs)

	z = scale_height * rho / energy
	# reshape by `numpy`'s automatic broacdasting to the same shape as z
	enp = np.ones_like(z) * energy

	pts = [
		(_e, _z, coeffs[_i, _j])
		for _i, _e in enumerate(ens)
		for _j, _z in enumerate(zm_p_en)
		if not nans[_i, _j]
	]
	pts = np.array(pts, copy=False)

	if log3:
		enp = np.log(enp)
		z = np.log(z)
		pts = np.log(pts)
	intp = interpolate.Rbf(*(pts.T), function=rbf)
	abr_zm = intp(enp, z)

	if log3:
		abr_zm = np.exp(abr_zm)
	return abr_zm * rho * flux