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 as in Aksnes et al., 2006 [0] |
34
|
|
|
|
35
|
|
|
Defined according to Aksnes et al., JGR 2006, Eq. (1), |
36
|
|
|
normalized to unit number flux, i.e. |
37
|
|
|
:math:`\int_0^\infty \phi(E) \text{d}E = 1`. |
38
|
|
|
|
39
|
|
|
Parameters |
40
|
|
|
---------- |
41
|
|
|
en: float or array_like (N,) |
42
|
|
|
Energy in [keV] |
43
|
|
|
en_0: float, optional |
44
|
|
|
Characteristic energy in [keV] of the distribution. |
45
|
|
|
Default: 10 keV |
46
|
|
|
|
47
|
|
|
Returns |
48
|
|
|
------- |
49
|
|
|
phi: float or array_like (N,) |
50
|
|
|
Normalized differential hemispherical number flux at `en` in [keV-1 cm-2 s-1] |
51
|
|
|
([keV] or scaled by 1 keV-2 cm-2 s-1, e.g. ). |
52
|
|
|
""" |
53
|
|
|
return 1. / en_0 * np.exp(-en / en_0) |
54
|
|
|
|
55
|
|
|
|
56
|
|
|
def gaussian_general(en, en_0=10., w=1.): |
57
|
|
|
r"""Gaussian number flux spectrum as in Fang2008 [1] |
58
|
|
|
|
59
|
|
|
Standard normal distribution with mu = en_0 and sigma = w / sqrt(2) |
60
|
|
|
for use in Fang et al., JGR 2008, Eq. (1). |
61
|
|
|
Almost normalized to unit number flux, i.e. |
62
|
|
|
:math:`\int_0^\infty \phi(E) \text{d}E = 1` |
63
|
|
|
(ignoring the negative tail). |
64
|
|
|
|
65
|
|
|
Parameters |
66
|
|
|
---------- |
67
|
|
|
en: float or array_like (N,) |
68
|
|
|
Energy in [keV] |
69
|
|
|
en_0: float, optional |
70
|
|
|
Characteristic energy in [keV], i.e. mode of the distribution. |
71
|
|
|
Default: 10 keV |
72
|
|
|
w: float, optional |
73
|
|
|
Width of the Gaussian distribution, in [keV]. |
74
|
|
|
|
75
|
|
|
Returns |
76
|
|
|
------- |
77
|
|
|
phi: float or array_like (N,) |
78
|
|
|
Normalized differential hemispherical number flux at `en` in [keV-1 cm-2 s-1] |
79
|
|
|
([keV] or scaled by 1 keV-2 cm-2 s-1, e.g.). |
80
|
|
|
""" |
81
|
|
|
return 1. / np.sqrt(np.pi * w**2) * np.exp(-(en - en_0)**2 / w**2) |
82
|
|
|
|
83
|
|
|
|
84
|
|
|
def maxwell_general(en, en_0=10.): |
85
|
|
|
r"""Maxwell number flux spectrum as in Fang2008 [1] |
86
|
|
|
|
87
|
|
|
Defined in Fang et al., JGR 2008, Eq. (1), |
88
|
|
|
normalized to :math:`\int_0^\infty \phi(E) \text{d}E = 1`. |
89
|
|
|
|
90
|
|
|
Parameters |
91
|
|
|
---------- |
92
|
|
|
en: float or array_like (N,) |
93
|
|
|
Energy in [keV] |
94
|
|
|
en_0: float, optional |
95
|
|
|
Characteristic energy in [keV], i.e. mode of the distribution. |
96
|
|
|
Default: 10 keV |
97
|
|
|
|
98
|
|
|
Returns |
99
|
|
|
------- |
100
|
|
|
phi: float or array_like (N,) |
101
|
|
|
Normalized differential hemispherical number flux at `en` in [keV-1 cm-2 s-1] |
102
|
|
|
([keV] or scaled by 1 keV-2 cm-2 s-1, e.g.). |
103
|
|
|
""" |
104
|
|
|
return en / en_0**2 * np.exp(-en / en_0) |
105
|
|
|
|
106
|
|
|
|
107
|
|
|
def pow_general(en, en_0=10., gamma=-3., het=True): |
108
|
|
|
r"""Power-law number flux spectrum as in Strickland1993 [3] |
109
|
|
|
|
110
|
|
|
Defined e.g. in Strickland et al., 1993, |
111
|
|
|
normalized to unit particle flux: |
112
|
|
|
:math:`\int_{E_0}^\infty \phi(E) \text{d}E = 1` |
113
|
|
|
for the high-energy tail version, and |
114
|
|
|
:math:`\int_0^{E_0} \phi(E) \text{d}E = 1` |
115
|
|
|
for the low-energy tail version. |
116
|
|
|
|
117
|
|
|
Parameters |
118
|
|
|
---------- |
119
|
|
|
en: float or array_like (N,) |
120
|
|
|
Energy in [keV] |
121
|
|
|
en_0: float, optional |
122
|
|
|
Characteristic energy in [keV], i.e. mode of the distribution. |
123
|
|
|
Default: 10 keV |
124
|
|
|
gamma: float, optional |
125
|
|
|
Exponent of the power-law distribution, in [keV]. |
126
|
|
|
het: bool, optional |
127
|
|
|
Return a high-energy tail (het, default: true) for en > en_0, |
128
|
|
|
or low-energy tail (false) for en < en_0. |
129
|
|
|
Adjusts the normalization accordingly. |
130
|
|
|
|
131
|
|
|
Returns |
132
|
|
|
------- |
133
|
|
|
phi: float or array_like (N,) |
134
|
|
|
Normalized differential hemispherical number flux at `en` in [keV-1 cm-2 s-1] |
135
|
|
|
([keV] or scaled by 1 keV-2 cm-2 s-1, e.g.). |
136
|
|
|
""" |
137
|
|
|
spec = (gamma + 1) / en_0 * (en / en_0)**gamma |
138
|
|
|
if het: |
139
|
|
|
spec[en < en_0] = 0. |
140
|
|
|
return -spec |
141
|
|
|
spec[en > en_0] = 0. |
142
|
|
|
return spec |
143
|
|
|
|
144
|
|
|
|
145
|
|
|
def pflux_exp(en, en_0=10.): |
146
|
|
|
r"""Exponential particle flux spectrum |
147
|
|
|
|
148
|
|
|
Normalized to unit energy flux: |
149
|
|
|
:math:`\int_0^\infty \phi(E) E \text{d}E = 1`. |
150
|
|
|
|
151
|
|
|
Parameters |
152
|
|
|
---------- |
153
|
|
|
en: float or array_like (N,) |
154
|
|
|
Energy in [keV] |
155
|
|
|
en_0: float, optional |
156
|
|
|
Characteristic energy in [keV], i.e. mode of the distribution. |
157
|
|
|
Default: 10 keV. |
158
|
|
|
|
159
|
|
|
Returns |
160
|
|
|
------- |
161
|
|
|
phi: float or array_like (N,) |
162
|
|
|
Hemispherical differential particle flux at `en` in [keV-1 cm-2 s-1] |
163
|
|
|
([kev-2] scaled by unit energy flux). |
164
|
|
|
""" |
165
|
|
|
return exp_general(en, en_0=en_0) / en_0 |
166
|
|
|
|
167
|
|
|
|
168
|
|
|
def pflux_gaussian(en, en_0=10., w=1): |
169
|
|
|
r"""Gaussian particle flux spectrum |
170
|
|
|
|
171
|
|
|
Defined in Fang et al., JGR 2008, Eq. (1). |
172
|
|
|
|
173
|
|
|
Normalized to :math:`\int_0^\infty \phi(E) E \text{d}E = 1`. |
174
|
|
|
(ignoring the negative tail). |
175
|
|
|
Parameters |
176
|
|
|
---------- |
177
|
|
|
en: float or array_like (N,) |
178
|
|
|
Energy in [keV] |
179
|
|
|
en_0: float, optional |
180
|
|
|
Characteristic energy in [keV], i.e. mode of the distribution. |
181
|
|
|
Default: 10 keV. |
182
|
|
|
|
183
|
|
|
Returns |
184
|
|
|
------- |
185
|
|
|
phi: float or array_like (N,) |
186
|
|
|
Hemispherical differential particle flux at `en` in [keV-1 cm-2 s-1] |
187
|
|
|
([kev-2] scaled by unit energy flux). |
188
|
|
|
""" |
189
|
|
|
return gaussian_general(en, en_0=en_0, w=w) / en_0 |
190
|
|
|
|
191
|
|
|
|
192
|
|
|
def pflux_maxwell(en, en_0=10.): |
193
|
|
|
r"""Maxwell particle flux spectrum as in Fang2008 [1] |
194
|
|
|
|
195
|
|
|
Defined in Fang et al., JGR 2008, Eq. (1). |
196
|
|
|
The total precipitating energy flux is fixed to 1 keV cm-2 s-1, |
197
|
|
|
multiply by Q_0 [keV cm-2 s-1] to scale the particle flux. |
198
|
|
|
|
199
|
|
|
Normalized to :math:`\int_0^\infty \phi(E) E \text{d}E = 1`. |
200
|
|
|
|
201
|
|
|
Parameters |
202
|
|
|
---------- |
203
|
|
|
en: float or array_like (N,) |
204
|
|
|
Energy in [keV] |
205
|
|
|
en_0: float, optional |
206
|
|
|
Characteristic energy in [keV], i.e. mode of the distribution. |
207
|
|
|
Default: 10 keV. |
208
|
|
|
|
209
|
|
|
Returns |
210
|
|
|
------- |
211
|
|
|
phi: float or array_like (N,) |
212
|
|
|
Hemispherical differential particle flux at `en` in [keV-1 cm-2 s-1] |
213
|
|
|
([kev-2] scaled by unit energy flux). |
214
|
|
|
""" |
215
|
|
|
return 0.5 / en_0 * maxwell_general(en, en_0) |
216
|
|
|
|
217
|
|
|
|
218
|
|
|
def pflux_pow(en, en_0=10., gamma=-3., het=True): |
219
|
|
|
r"""Power-law particle flux spectrum |
220
|
|
|
|
221
|
|
|
Defined e.g. in Strickland et al., 1993. |
222
|
|
|
Normalized to :math:`\int_{E_0}^\infty \phi(E) E \text{d}E = 1` |
223
|
|
|
for the high-energy tail version, and to |
224
|
|
|
:math:`\int_0^{E_0} \phi(E) E \text{d}E = 1` |
225
|
|
|
for the low-energy tail version. |
226
|
|
|
|
227
|
|
|
Parameters |
228
|
|
|
---------- |
229
|
|
|
en: float or array_like (N,) |
230
|
|
|
Energy in [keV] |
231
|
|
|
en_0: float, optional |
232
|
|
|
Characteristic energy in [keV], i.e. mode of the distribution. |
233
|
|
|
Default: 10 keV |
234
|
|
|
gamma: float, optional |
235
|
|
|
Exponent of the power-law distribution, in [keV]. |
236
|
|
|
het: bool, optional (default True) |
237
|
|
|
Return a high-energy tail (true) for en > en_0, |
238
|
|
|
or low-energy tail (false) for en < en_0. |
239
|
|
|
Adjusts the normalization accordingly. |
240
|
|
|
|
241
|
|
|
Returns |
242
|
|
|
------- |
243
|
|
|
phi: float or array_like (N,) |
244
|
|
|
Hemispherical differential particle flux at `en` in [keV-1 cm-2 s-1] |
245
|
|
|
([kev-2] scaled by unit energy flux). |
246
|
|
|
""" |
247
|
|
|
return (gamma + 2) / (gamma + 1) / en_0 * pow_general(en, en_0=en_0, gamma=gamma, het=het) |
248
|
|
|
|
249
|
|
|
|
250
|
|
|
def ediss_spec_int( |
251
|
|
|
ens, |
252
|
|
|
dfluxes, |
253
|
|
|
scale_height, |
254
|
|
|
rho, |
255
|
|
|
func, |
256
|
|
|
axis=-1, |
257
|
|
|
func_kws=None, |
258
|
|
|
): |
259
|
|
|
r"""Integrate over a given energy spectrum |
260
|
|
|
|
261
|
|
|
Integrates a mono-energetic parametrization `q`, e.g. from Fang et al., 2010 |
262
|
|
|
using the given differential particle spectrum `phi`: |
263
|
|
|
|
264
|
|
|
:math:`\int_\text{spec} \phi(E) q(E, Q) E \text{d}E` |
265
|
|
|
|
266
|
|
|
This function uses the differential spectrum evaluated at the given energy bins. |
267
|
|
|
|
268
|
|
|
Parameters |
269
|
|
|
---------- |
270
|
|
|
ens: array_like (M,...) |
271
|
|
|
Central (bin) energies of the spectrum |
272
|
|
|
dfluxes: array_like (M,...) |
273
|
|
|
Differential particle fluxes in the given bins |
274
|
|
|
scale_height: array_like (N,...) |
275
|
|
|
The atmospheric scale heights |
276
|
|
|
rho: array_like (N,...) |
277
|
|
|
The atmospheric densities, corresponding to the |
278
|
|
|
scale heights. |
279
|
|
|
func: callable |
280
|
|
|
Mono-energetic energy dissipation function to integrate. |
281
|
|
|
axis: int, optional |
282
|
|
|
The axis to use for integration, default: -1 (last axis). |
283
|
|
|
func_kws: dict-like, optional |
284
|
|
|
Optional keyword arguments to pass to the mono-energetic |
285
|
|
|
energy dissipation function. Default: `None` |
286
|
|
|
|
287
|
|
|
Returns |
288
|
|
|
------- |
289
|
|
|
en_diss: array_like (N) |
290
|
|
|
The dissipated energy profiles [keV]. |
291
|
|
|
|
292
|
|
|
See Also |
293
|
|
|
-------- |
294
|
|
|
ediss_specfun_int |
295
|
|
|
""" |
296
|
|
|
ens = np.atleast_1d(ens) |
297
|
|
|
dfluxes = np.atleast_1d(dfluxes) |
298
|
|
|
scale_height = np.atleast_1d(scale_height) |
299
|
|
|
rho = np.atleast_1d(rho) |
300
|
|
|
func_kws = func_kws or dict() |
301
|
|
|
ediss = func( |
302
|
|
|
ens[None, None, :], |
303
|
|
|
dfluxes, |
304
|
|
|
scale_height[..., None], |
305
|
|
|
rho[..., None], |
306
|
|
|
**func_kws, |
307
|
|
|
) |
308
|
|
|
return np.trapz(ediss * ens, ens, axis=axis) |
309
|
|
|
|
310
|
|
|
|
311
|
|
|
def ediss_specfun_int( |
312
|
|
|
energy, |
313
|
|
|
flux, |
314
|
|
|
scale_height, |
315
|
|
|
rho, |
316
|
|
|
ediss_func, |
317
|
|
|
ediss_kws=None, |
318
|
|
|
bounds=(0.1, 300.), |
319
|
|
|
nstep=128, |
320
|
|
|
spec_fun=pflux_maxwell, |
321
|
|
|
spec_kws=None, |
322
|
|
|
): |
323
|
|
|
"""Integrate mono-energetic parametrization over a spectrum |
324
|
|
|
|
325
|
|
|
Integrates the mono-energetic parametrization over a spectrum given by a |
326
|
|
|
functional dependence with characteristic energy `energy` and total energy |
327
|
|
|
flux `flux`. |
328
|
|
|
|
329
|
|
|
Parameters |
330
|
|
|
---------- |
331
|
|
|
energy: float or array_like (M,...) |
332
|
|
|
Characteristic energy E_0 [keV] of the spectral distribution. |
333
|
|
|
flux: float or array_like (M,...) |
334
|
|
|
Integrated energy flux Q_0 [keV / cm² / s¹] |
335
|
|
|
scale_height: float or array_like (N,...) |
336
|
|
|
The atmospheric scale heights [cm]. |
337
|
|
|
rho: float or array_like (N,...) |
338
|
|
|
The atmospheric mass density [g / cm³] |
339
|
|
|
ediss_func: callable |
340
|
|
|
Mono-energetic energy dissipation function to integrate. |
341
|
|
|
ediss_kws: dict-like, optional |
342
|
|
|
Optional keyword arguments to pass to the mono-energetic |
343
|
|
|
energy dissipation function. Default: `None` |
344
|
|
|
bounds: tuple, optional |
345
|
|
|
(min, max) [keV] of the integration range to integrate the Maxwellian. |
346
|
|
|
Make sure that this is appropriate to encompass the spectrum. |
347
|
|
|
Default: (0.1, 300.) |
348
|
|
|
nsteps: int, optional |
349
|
|
|
Number of integration steps, default: 128. |
350
|
|
|
spec_func: callable, optional, default :func:`pflux_maxwell` |
351
|
|
|
Spectral shape function, choices are: |
352
|
|
|
* :func:`pflux_exp` for a exponential spectrum |
353
|
|
|
* :func:`pflux_gaussian` for a Gaussian shaped spectrum |
354
|
|
|
* :func:`pflux_maxwell` for a Maxwellian shaped spectrum |
355
|
|
|
* :func:`pflux_pow` for a power-law |
356
|
|
|
spec_kws: dict-like, optional |
357
|
|
|
Optional keyword arguments to pass to the spectral function |
358
|
|
|
Default: `None` |
359
|
|
|
|
360
|
|
|
Returns |
361
|
|
|
------- |
362
|
|
|
en_diss: array_like (M,N) |
363
|
|
|
The dissipated energy profiles [keV]. |
364
|
|
|
|
365
|
|
|
See Also |
366
|
|
|
-------- |
367
|
|
|
ediss_spec_int |
368
|
|
|
""" |
369
|
|
|
energy = np.asarray(energy) |
370
|
|
|
flux = np.asarray(flux) |
371
|
|
|
bounds_l10 = np.log10(bounds) |
372
|
|
|
ens = np.logspace(*bounds_l10, num=nstep) |
373
|
|
|
ensd = np.reshape(ens, (-1,) + (1,) * energy.ndim) |
374
|
|
|
spec_kws = spec_kws or dict() |
375
|
|
|
# "overwrite" the characteristic energy |
376
|
|
|
spec_kws["en_0"] = energy.T |
377
|
|
|
dflux = flux.T * spec_fun(ensd, **spec_kws) |
378
|
|
|
return ediss_spec_int( |
379
|
|
|
ens, dflux.T, scale_height, rho, ediss_func, |
380
|
|
|
axis=-1, func_kws=ediss_kws, |
381
|
|
|
) |
382
|
|
|
|