eppaurora.brems.berger1974()

Complexity

Conditions

4

Size

Total Lines	101
Code Lines	33

Duplication

Lines	0
Ratio	0 %

Importance

Changes

0

# coding: utf-8
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