|
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 [1]_. |
|
13
|
|
|
|
|
14
|
|
|
.. [1] Zhang and Paxton, JASTP, 70, 1231--1242, 2008 |
|
15
|
|
|
""" |
|
16
|
|
|
from os import path |
|
17
|
|
|
from pkg_resources import resource_filename |
|
18
|
|
|
|
|
19
|
|
|
import numpy as np |
|
20
|
|
|
|
|
21
|
|
|
__all__ = [ |
|
22
|
|
|
"epstein_coeffs", |
|
23
|
|
|
"epstein_eval", |
|
24
|
|
|
"hemispheric_power", |
|
25
|
|
|
"read_zp2008_coeffs", |
|
26
|
|
|
"zp2008", |
|
27
|
|
|
] |
|
28
|
|
|
|
|
29
|
|
|
COEFF_FILE = "Zhang2008.txt" |
|
30
|
|
|
COEFF_PATH = resource_filename(__name__, path.join("data", COEFF_FILE)) |
|
31
|
|
|
COEFF_NAMES = ["A", "B", "C", "D"] |
|
32
|
|
|
|
|
33
|
|
|
# Kp_model bin edges and centres as in Zhang and Paxton, 2008 |
|
34
|
|
|
KP_BINE = [ |
|
35
|
|
|
(0.0, 1.5), (1.5, 3.0), (3.0, 4.5), (4.5, 6.0), (6.0, 8.0), (8.0, 10.0) |
|
36
|
|
|
] |
|
37
|
|
|
KP_BINC = np.asarray(KP_BINE).mean(axis=1) |
|
38
|
|
|
|
|
39
|
|
|
|
|
40
|
|
|
def hemispheric_power(Kp): |
|
41
|
|
|
"""Hemispheric Power in GW from Kp |
|
42
|
|
|
|
|
43
|
|
|
Zhang and Paxton, 2008, Eqs. (1) and (2) |
|
44
|
|
|
|
|
45
|
|
|
Parameters |
|
46
|
|
|
----------- |
|
47
|
|
|
Kp: float or array_like |
|
48
|
|
|
Geomagnetic Kp index value(s). |
|
49
|
|
|
|
|
50
|
|
|
Returns |
|
51
|
|
|
------- |
|
52
|
|
|
HP: float or array_like |
|
53
|
|
|
Hemispheric power in [GW], same shape as `Kp`. |
|
54
|
|
|
""" |
|
55
|
|
|
Kp = np.asarray(Kp) |
|
56
|
|
|
return np.where( |
|
57
|
|
|
Kp <= 5.0, |
|
58
|
|
|
38.66 * np.exp(0.1967 * Kp) - 33.99, |
|
59
|
|
|
4.592 * np.exp(0.4731 * Kp) + 20.47, |
|
60
|
|
|
) |
|
61
|
|
|
|
|
62
|
|
|
|
|
63
|
|
|
def read_zp2008_coeffs(file=None, nf=6, nKp=len(KP_BINC)): |
|
64
|
|
|
fdata = np.genfromtxt( |
|
65
|
|
|
file or COEFF_PATH, |
|
66
|
|
|
delimiter=" ", |
|
67
|
|
|
dtype=None, |
|
68
|
|
|
names=["name"] + COEFF_NAMES, |
|
69
|
|
|
) |
|
70
|
|
|
# number of coefficients per Kp bin |
|
71
|
|
|
nc = 2 * nf + 1 |
|
72
|
|
|
# Energy fluxes are in Table 1. |
|
73
|
|
|
Q0tab = [fdata[n * nc:(n + 1) * nc] for n in range(nKp)] |
|
74
|
|
|
# Mean energies are in Table 2. |
|
75
|
|
|
Emtab = [fdata[n * nc:(n + 1) * nc] for n in range(nKp, 2 * nKp)] |
|
76
|
|
|
return Q0tab, Emtab |
|
77
|
|
|
|
|
78
|
|
|
|
|
79
|
|
|
def find_Kp_idx(Kp): |
|
80
|
|
|
# maximum avoids negative indices |
|
81
|
|
|
return np.maximum(0, np.searchsorted(KP_BINC, Kp, side="left") - 1) |
|
82
|
|
|
|
|
83
|
|
|
|
|
84
|
|
|
def epstein_coeffs(angle, table): |
|
85
|
|
|
r""" |
|
86
|
|
|
:math:`angle = MLT * 2\pi / 24` |
|
87
|
|
|
""" |
|
88
|
|
|
coeffs = np.array(table[COEFF_NAMES].tolist()) |
|
89
|
|
|
nf = (len(coeffs) - 1) // 2 |
|
90
|
|
|
if (2 * nf + 1) != len(coeffs): |
|
91
|
|
|
raise ValueError("Number of coefficients is inconsistent.") |
|
92
|
|
|
fs = np.arange(1, nf + 1) |
|
93
|
|
|
cos = np.cos(fs * angle) @ coeffs[1:nf + 1] |
|
94
|
|
|
sin = np.sin(fs * angle) @ coeffs[nf + 1:2 * nf + 1] |
|
95
|
|
|
return coeffs[0] + cos + sin |
|
96
|
|
|
|
|
97
|
|
|
|
|
98
|
|
|
def epstein_eval(x, coeffs): |
|
99
|
|
|
"""Epstein function evaluated at x |
|
100
|
|
|
|
|
101
|
|
|
x = 90 - |Mlat| |
|
102
|
|
|
coeffs = Epstein coefficients, e.g. from `epstein_coeffs()` |
|
103
|
|
|
""" |
|
104
|
|
|
a, b, c, d = coeffs |
|
105
|
|
|
loc = x - b |
|
106
|
|
|
return a * np.exp(loc / c) / (1 + np.exp(loc / d))**2 |
|
107
|
|
|
|
|
108
|
|
|
|
|
109
|
|
|
def zp2008(mlat, mlt, Kp, Q0table=None, Emtable=None): |
|
110
|
|
|
u""" |
|
111
|
|
|
Parameters |
|
112
|
|
|
---------- |
|
113
|
|
|
mlat: float |
|
114
|
|
|
(Geo)Magnetic latitude in [degrees]. |
|
115
|
|
|
mlt: float |
|
116
|
|
|
Magnetic local time in [hours]. |
|
117
|
|
|
Kp: float |
|
118
|
|
|
Geomagnetic Kp index value(s). |
|
119
|
|
|
Q0table: np.recarray, optional |
|
120
|
|
|
Fourier coefficient table for the Epstein coefficients |
|
121
|
|
|
for the energy flux. E.g. as returned by `read_zp2008_coeffs()`. |
|
122
|
|
|
Emtable: np.recarray, optional |
|
123
|
|
|
Fourier coefficient table for the Epstein coefficients |
|
124
|
|
|
for the mean energy. E.g. as returned by `read_zp2008_coeffs()`. |
|
125
|
|
|
|
|
126
|
|
|
Returns |
|
127
|
|
|
------- |
|
128
|
|
|
(Q0, Em): tuple |
|
129
|
|
|
Electron energy flux Q0 in [mW m⁻²] (or [erg s⁻¹ cm⁻²]), |
|
130
|
|
|
and electron mean energy in [keV]. |
|
131
|
|
|
""" |
|
132
|
|
|
if (Q0table is None) or (Emtable is None): |
|
133
|
|
|
Q0t, Emt = read_zp2008_coeffs() |
|
134
|
|
|
Q0table = Q0table or Q0t |
|
135
|
|
|
Emtable = Emtable or Emt |
|
136
|
|
|
|
|
137
|
|
|
angle = mlt * np.pi / 12.0 |
|
138
|
|
|
x = 90.0 - np.abs(mlat) |
|
139
|
|
|
ix = find_Kp_idx(Kp) |
|
140
|
|
|
ixs = [ix, ix + 1] |
|
141
|
|
|
Kps = KP_BINC[ixs] |
|
142
|
|
|
Q0_epst = [ |
|
143
|
|
|
epstein_eval(x, epstein_coeffs(angle, Q0table[ix])) |
|
144
|
|
|
for ix in ixs |
|
145
|
|
|
] |
|
146
|
|
|
Em_epst = [ |
|
147
|
|
|
epstein_eval(x, epstein_coeffs(angle, Emtable[ix])) |
|
148
|
|
|
for ix in ixs |
|
149
|
|
|
] |
|
150
|
|
|
Q0 = np.interp(hemispheric_power(Kp), hemispheric_power(Kps), Q0_epst) |
|
151
|
|
|
Em = np.interp(Kp, Kps, Em_epst) |
|
152
|
|
|
return Q0, Em |
|
153
|
|
|
|
st-bender /
pyeppaurora
