Passed
Push — master ( 01234e...e0a97d )
by Stefan
01:10
created

eppaurora.spectra.ediss_specfun_int()   A

Complexity

Conditions 1

Size

Total Lines 70
Code Lines 22

Duplication

Lines 0
Ratio 0 %

Importance

Changes 0
Metric Value
cc 1
eloc 22
nop 10
dl 0
loc 70
rs 9.352
c 0
b 0
f 0

How to fix   Long Method    Many Parameters   

Long Method

Small methods make your code easier to understand, in particular if combined with a good name. Besides, if your method is small, finding a good name is usually much easier.

For example, if you find yourself adding comments to a method's body, this is usually a good sign to extract the commented part to a new method, and use the comment as a starting point when coming up with a good name for this new method.

Commonly applied refactorings include:

Many Parameters

Methods with many parameters are not only hard to understand, but their parameters also often become inconsistent when you need more, or different data.

There are several approaches to avoid long parameter lists:

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