eppaurora.models.zhangpaxton2008.zp2008() - Code Metrics - st-bender/pyeppaurora - Measure and Improve Code Quality continuously with Scrutinizer

eppaurora.models.zhangpaxton2008.zp2008() A
last analyzed 2024-04-05 18:05 UTC

↳ Parent: eppaurora.models.zhangpaxton2008

Complexity

Conditions

Size

Total Lines	54
Code Lines	19

Duplication

Lines	0
Ratio	0 %

Importance

Changes

Metric	Value
cc	3
eloc	19
nop	5
dl	0
loc	54
rs	9.45
c	0
b	0
f	0

How to fix Long Method

# coding: utf-8
# Copyright (c) 2023 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.
"""Empirical model for electron energy and flux

Implements the empirical Kp-driven model for auroral electrons,
providing mean energy and energy flux as described in [ZP08]_.

.. [ZP08] Zhang and Paxton, JASTP, 70, 1231--1242, 2008,
	https://doi.org/10.1016/j.jastp.2008.03.008
"""
from os import path
from pkg_resources import resource_filename

import numpy as np

__all__ = [
	"epstein_coeffs",
	"epstein_eval",
	"hemispheric_power",
	"read_zp2008_coeffs",
	"zp2008",
]

COEFF_FILE = "Zhang2008.txt"
COEFF_PATH = resource_filename(__name__, path.join("data", COEFF_FILE))
COEFF_NAMES = ["A", "B", "C", "D"]

# Kp_model bin edges and centres as in Zhang and Paxton, 2008
KP_BINE = [
	(0.0, 1.5), (1.5, 3.0), (3.0, 4.5), (4.5, 6.0), (6.0, 8.0), (8.0, 10.0)
]
KP_BINC = np.asarray(KP_BINE).mean(axis=1)


def hemispheric_power(Kp):
	"""Hemispheric Power in GW from Kp

	Zhang and Paxton, 2008, Eqs. (1) and (2)

	Parameters
	----------
	Kp: float or array_like
		Geomagnetic Kp index value(s).

	Returns
	-------
	HP: float or array_like
		Hemispheric power in [GW], same shape as `Kp`.
	"""
	Kp = np.asarray(Kp)
	return np.where(
		Kp <= 5.0,
		38.66 * np.exp(0.1967 * Kp) - 33.99,
		4.592 * np.exp(0.4731 * Kp) + 20.47,
	)


def read_zp2008_coeffs(file=None, nf=6, nKp=len(KP_BINC)):
	"""Read Epstein coefficient tables from file

	Parameters
	----------
	file: str, optional
		The text file containing the coefficient tables, with the same layout
		as in [ZP08]_. Defaults to the packaged coefficient file.
	nf: int, optional
		The number of harmonic (Fourier) terms, defaults to 6 as in [ZP08]_
	nKp: int, optional
		The number of Kp bins, defaults to 6 as in [ZP08]_

	Returns
	-------
	Q0_table, Em_table: tuple of ``numpy.recarray``
		The tables for the total energy flux Q0 and the electron mean energy Em,
		the row names indicate the frequency and the columns the Epstein
		parameters A, B, C, D.

	References
	----------
	.. [ZP08] Zhang and Paxton, JASTP, 70, 1231--1242, 2008,
		https://doi.org/10.1016/j.jastp.2008.03.008
	"""
	fdata = np.genfromtxt(
		file or COEFF_PATH,
		delimiter=" ",
		dtype=None,
		encoding="utf-8",
		names=["name"] + COEFF_NAMES,
	)
	# number of coefficients per Kp bin
	nc = 2 * nf + 1
	# Energy fluxes are in Table 1.
	Q0tab = [fdata[n * nc:(n + 1) * nc] for n in range(nKp)]
	# Mean energies are in Table 2.
	Emtab = [fdata[n * nc:(n + 1) * nc] for n in range(nKp, 2 * nKp)]
	return Q0tab, Emtab


def find_Kp_idx(Kp):
	# maximum avoids negative indices
	return np.maximum(0, np.searchsorted(KP_BINC, Kp, side="left") - 1)


def epstein_coeffs(angle, table):
	r"""Epstein coefficients from table

	Returns the Epstein coefficients as read from the table,
	ready for evaluation with :func:`epstein_eval()`.

	Parameters
	----------
	angle: float or array_like
		The magnetic local time hour angle: :math:`angle = \text{MLT} * 2\pi / 24`
	table: array_like
		Table of N harmonic (Fourier) coefficients.
		The 4 constant offsets are in the first row, then Nx4 cosine amplitudes
		followed by Nx4 sine amplitudes, each for the coefficients A, B, C, and D.

	Returns
	-------
	coeffs: array_like (4,)
		The Epstein coefficients for the MLT angle.

	See Also
	--------
	epstein_eval
	"""
	coeffs = np.array(table[COEFF_NAMES].tolist())
	nf = (len(coeffs) - 1) // 2
	if (2 * nf + 1) != len(coeffs):
		raise ValueError("Number of coefficients is inconsistent.")
	fs = np.arange(1, nf + 1)
	cos = np.cos(fs * angle).dot(coeffs[1:nf + 1])
	sin = np.sin(fs * angle).dot(coeffs[nf + 1:2 * nf + 1])
	return coeffs[0] + cos + sin


def epstein_eval(x, coeffs):
	r"""Epstein function evaluated at x

	The so-called Epstein function is defined by

	.. math:: E(x) = \frac{A\exp\{(x - B) / C\}}{(1 + \exp\{(x - B) / D\})^2}

	Parameters
	----------
	x: float or array_like
		Argument of the Epstein function, e.g. the magnetic co-latitude
		:math:`x = 90 - |Mlat|`.
	coeffs: array_like
		Epstein coefficients, e.g. from :func:`epstein_coeffs()`.

	Returns
	-------
	y: float or array_like
		The Epstein function with coefficients as given evaluated at x.
	"""
	a, b, c, d = coeffs
	loc = x - b
	return a * np.exp(loc / c) / (1 + np.exp(loc / d))**2


def zp2008(mlat, mlt, Kp, Q0table=None, Emtable=None):
	u"""Electron total energy flux and mean energy model

	Implements the model algorithm as given in Appendix A of [ZP08]_.
	Defaults to using the packaged coefficients taken from that reference,
	but custom tables for the Q0 and Em Fourier coefficients can be provided.

	Parameters
	----------
	mlat: float
		(Geo)Magnetic latitude in [degrees].
	mlt: float
		Magnetic local time in [hours].
	Kp: float
		Geomagnetic Kp index value(s).
	Q0table: np.recarray, optional
		Fourier coefficient table for the Epstein coefficients
		for the energy flux. E.g. as returned by `read_zp2008_coeffs()`.
	Emtable: np.recarray, optional
		Fourier coefficient table for the Epstein coefficients
		for the mean energy. E.g. as returned by `read_zp2008_coeffs()`.

	Returns
	-------
	(Q0, Em): tuple
		Electron energy flux Q0 in [mW m⁻²] (= [erg s⁻¹ cm⁻²]),
		and electron mean energy in [keV].

	References
	----------
	.. [ZP08] Zhang and Paxton, JASTP, 70, 1231--1242, 2008,
		https://doi.org/10.1016/j.jastp.2008.03.008
	"""
	if (Q0table is None) or (Emtable is None):
		Q0t, Emt = read_zp2008_coeffs()
		Q0table = Q0table or Q0t
		Emtable = Emtable or Emt

	angle = mlt * np.pi / 12.0
	x = 90.0 - np.abs(mlat)
	ix = find_Kp_idx(Kp)
	ixs = [ix, ix + 1]
	Kps = KP_BINC[ixs]
	Q0_epst = [
		epstein_eval(x, epstein_coeffs(angle, Q0table[ix]))
		for ix in ixs
	]
	Em_epst = [
		epstein_eval(x, epstein_coeffs(angle, Emtable[ix]))
		for ix in ixs
	]
	Q0 = np.interp(hemispheric_power(Kp), hemispheric_power(Kps), Q0_epst)
	Em = np.interp(Kp, Kps, Em_epst)
	return Q0, Em


1			# coding: utf-8
2			# Copyright (c) 2023 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			"""Empirical model for electron energy and flux
10
11			Implements the empirical Kp-driven model for auroral electrons,
12			providing mean energy and energy flux as described in [ZP08]_.
13
14			.. [ZP08] Zhang and Paxton, JASTP, 70, 1231--1242, 2008,
15			https://doi.org/10.1016/j.jastp.2008.03.008
16			"""
17			from os import path
18			from pkg_resources import resource_filename
19
20			import numpy as np
21
22			__all__ = [
23			"epstein_coeffs",
24			"epstein_eval",
25			"hemispheric_power",
26			"read_zp2008_coeffs",
27			"zp2008",
28			]
29
30			COEFF_FILE = "Zhang2008.txt"
31			COEFF_PATH = resource_filename(__name__, path.join("data", COEFF_FILE))
32			COEFF_NAMES = ["A", "B", "C", "D"]
33
34			# Kp_model bin edges and centres as in Zhang and Paxton, 2008
35			KP_BINE = [
36			(0.0, 1.5), (1.5, 3.0), (3.0, 4.5), (4.5, 6.0), (6.0, 8.0), (8.0, 10.0)
37			]
38			KP_BINC = np.asarray(KP_BINE).mean(axis=1)
39
40
41			def hemispheric_power(Kp):
42			"""Hemispheric Power in GW from Kp
43
44			Zhang and Paxton, 2008, Eqs. (1) and (2)
45
46			Parameters
47			----------
48			Kp: float or array_like
49			Geomagnetic Kp index value(s).
50
51			Returns
52			-------
53			HP: float or array_like
54			Hemispheric power in [GW], same shape as `Kp`.
55			"""
56			Kp = np.asarray(Kp)
57			return np.where(
58			Kp <= 5.0,
59			38.66 * np.exp(0.1967 * Kp) - 33.99,
60			4.592 * np.exp(0.4731 * Kp) + 20.47,
61			)
62
63
64			def read_zp2008_coeffs(file=None, nf=6, nKp=len(KP_BINC)):
65			"""Read Epstein coefficient tables from file
66
67			Parameters
68			----------
69			file: str, optional
70			The text file containing the coefficient tables, with the same layout
71			as in [ZP08]_. Defaults to the packaged coefficient file.
72			nf: int, optional
73			The number of harmonic (Fourier) terms, defaults to 6 as in [ZP08]_
74			nKp: int, optional
75			The number of Kp bins, defaults to 6 as in [ZP08]_
76
77			Returns
78			-------
79			Q0_table, Em_table: tuple of ``numpy.recarray``
80			The tables for the total energy flux Q0 and the electron mean energy Em,
81			the row names indicate the frequency and the columns the Epstein
82			parameters A, B, C, D.
83
84			References
85			----------
86			.. [ZP08] Zhang and Paxton, JASTP, 70, 1231--1242, 2008,
87			https://doi.org/10.1016/j.jastp.2008.03.008
88			"""
89			fdata = np.genfromtxt(
90			file or COEFF_PATH,
91			delimiter=" ",
92			dtype=None,
93			encoding="utf-8",
94			names=["name"] + COEFF_NAMES,
95			)
96			# number of coefficients per Kp bin
97			nc = 2 * nf + 1
98			# Energy fluxes are in Table 1.
99			Q0tab = [fdata[n * nc:(n + 1) * nc] for n in range(nKp)]
100			# Mean energies are in Table 2.
101			Emtab = [fdata[n * nc:(n + 1) * nc] for n in range(nKp, 2 * nKp)]
102			return Q0tab, Emtab
103
104
105			def find_Kp_idx(Kp):
106			# maximum avoids negative indices
107			return np.maximum(0, np.searchsorted(KP_BINC, Kp, side="left") - 1)
108
109
110			def epstein_coeffs(angle, table):
111			r"""Epstein coefficients from table
112
113			Returns the Epstein coefficients as read from the table,
114			ready for evaluation with :func:`epstein_eval()`.
115
116			Parameters
117			----------
118			angle: float or array_like
119			The magnetic local time hour angle: :math:`angle = \text{MLT} * 2\pi / 24`
120			table: array_like
121			Table of N harmonic (Fourier) coefficients.
122			The 4 constant offsets are in the first row, then Nx4 cosine amplitudes
123			followed by Nx4 sine amplitudes, each for the coefficients A, B, C, and D.
124
125			Returns
126			-------
127			coeffs: array_like (4,)
128			The Epstein coefficients for the MLT angle.
129
130			See Also
131			--------
132			epstein_eval
133			"""
134			coeffs = np.array(table[COEFF_NAMES].tolist())
135			nf = (len(coeffs) - 1) // 2
136			if (2 * nf + 1) != len(coeffs):
137			raise ValueError("Number of coefficients is inconsistent.")
138			fs = np.arange(1, nf + 1)
139			cos = np.cos(fs * angle).dot(coeffs[1:nf + 1])
140			sin = np.sin(fs * angle).dot(coeffs[nf + 1:2 * nf + 1])
141			return coeffs[0] + cos + sin
142
143
144			def epstein_eval(x, coeffs):
145			r"""Epstein function evaluated at x
146
147			The so-called Epstein function is defined by
148
149			.. math:: E(x) = \frac{A\exp\{(x - B) / C\}}{(1 + \exp\{(x - B) / D\})^2}
150
151			Parameters
152			----------
153			x: float or array_like
154			Argument of the Epstein function, e.g. the magnetic co-latitude
155			:math:`x = 90 - \|Mlat\|`.
156			coeffs: array_like
157			Epstein coefficients, e.g. from :func:`epstein_coeffs()`.
158
159			Returns
160			-------
161			y: float or array_like
162			The Epstein function with coefficients as given evaluated at x.
163			"""
164			a, b, c, d = coeffs
165			loc = x - b
166			return a * np.exp(loc / c) / (1 + np.exp(loc / d))**2
167
168
169			def zp2008(mlat, mlt, Kp, Q0table=None, Emtable=None):
170			u"""Electron total energy flux and mean energy model
171
172			Implements the model algorithm as given in Appendix A of [ZP08]_.
173			Defaults to using the packaged coefficients taken from that reference,
174			but custom tables for the Q0 and Em Fourier coefficients can be provided.
175
176			Parameters
177			----------
178			mlat: float
179			(Geo)Magnetic latitude in [degrees].
180			mlt: float
181			Magnetic local time in [hours].
182			Kp: float
183			Geomagnetic Kp index value(s).
184			Q0table: np.recarray, optional
185			Fourier coefficient table for the Epstein coefficients
186			for the energy flux. E.g. as returned by `read_zp2008_coeffs()`.
187			Emtable: np.recarray, optional
188			Fourier coefficient table for the Epstein coefficients
189			for the mean energy. E.g. as returned by `read_zp2008_coeffs()`.
190
191			Returns
192			-------
193			(Q0, Em): tuple
194			Electron energy flux Q0 in [mW m⁻²] (= [erg s⁻¹ cm⁻²]),
195			and electron mean energy in [keV].
196
197			References
198			----------
199			.. [ZP08] Zhang and Paxton, JASTP, 70, 1231--1242, 2008,
200			https://doi.org/10.1016/j.jastp.2008.03.008
201			"""
202			if (Q0table is None) or (Emtable is None):
203			Q0t, Emt = read_zp2008_coeffs()
204			Q0table = Q0table or Q0t
205			Emtable = Emtable or Emt
206
207			angle = mlt * np.pi / 12.0
208			x = 90.0 - np.abs(mlat)
209			ix = find_Kp_idx(Kp)
210			ixs = [ix, ix + 1]
211			Kps = KP_BINC[ixs]
212			Q0_epst = [
213			epstein_eval(x, epstein_coeffs(angle, Q0table[ix]))
214			for ix in ixs
215			]
216			Em_epst = [
217			epstein_eval(x, epstein_coeffs(angle, Emtable[ix]))
218			for ix in ixs
219			]
220			Q0 = np.interp(hemispheric_power(Kp), hemispheric_power(Kps), Q0_epst)
221			Em = np.interp(Kp, Kps, Em_epst)
222			return Q0, Em
223

st-bender / pyeppaurora

eppaurora.models.zhangpaxton2008.zp2008() A last analyzed 2024-04-05 18:05 UTC

Complexity

Size

Duplication

Importance

How to fix Long Method

Long Method

Duplication Side-by-Side

Filter issues like

eppaurora.models.zhangpaxton2008.zp2008() A
last analyzed 2024-04-05 18:05 UTC