-
Notifications
You must be signed in to change notification settings - Fork 306
/
radiation.py
157 lines (125 loc) · 4.84 KB
/
radiation.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
"""
Functions for calculating quantities associated with electromagnetic
radiation.
"""
__all__ = [
"thermal_bremsstrahlung",
]
import astropy.constants as const
import astropy.units as u
import numpy as np
from scipy.special import exp1
from plasmapy.formulary.parameters import plasma_frequency
from plasmapy.particles import Particle, particle_input
from plasmapy.utils.decorators import validate_quantities
from plasmapy.utils.exceptions import PhysicsError
@validate_quantities(
frequencies={"can_be_negative": False},
n_e={"can_be_negative": False},
n_i={"can_be_negative": False},
T_e={"can_be_negative": False, "equivalencies": u.temperature_energy()},
)
@particle_input
def thermal_bremsstrahlung(
frequencies: u.Hz,
n_e: u.m ** -3,
T_e: u.K,
n_i: u.m ** -3 = None,
ion_species: Particle = "H+",
kmax: u.m = None,
) -> np.ndarray:
r"""
Calculate the bremsstrahlung emission spectrum for a Maxwellian plasma
in the Rayleigh-Jeans limit :math:`ℏ ω ≪ k_B T_e`
.. math::
\frac{dP}{dω} = \frac{8 \sqrt{2}}{3\sqrt{π}}
\bigg ( \frac{e^2}{4 π ε_0} \bigg )^3
\bigg ( m_e c^2 \bigg )^{-\frac{3}{2}}
\bigg ( 1 - \frac{ω_{pe}^2}{ω^2} \bigg )^\frac{1}{2}
\frac{Z_i^2 n_i n_e}{\sqrt(k_B T_e)}
E_1(y)
where :math:`E_1` is the exponential integral
.. math::
E_1 (y) = - \int_{-y}^∞ \frac{e^{-t}}{t}dt
and :math:`y` is the dimensionless argument
.. math::
y = \frac{1}{2} \frac{ω^2 m_e}{k_{max}^2 k_B T_e}
where :math:`k_{max}` is a maximum wavenumber approximated here as
:math:`k_{max} = 1/λ_B` where :math:`λ_B` is the electron
de Broglie wavelength.
Parameters
----------
frequencies : `~astropy.units.Quantity`
Array of frequencies over which the bremsstrahlung spectrum will be
calculated (convertible to Hz).
n_e : `~astropy.units.Quantity`
Electron number density in the plasma (convertible to m\ :sup:`-3`\ ).
T_e : `~astropy.units.Quantity`
Temperature of the electrons (in K or convertible to eV).
n_i : `~astropy.units.Quantity`, optional
Ion number density in the plasma (convertible to m\ :sup:`-3`\ ). Defaults
to the quasi-neutral condition :math:`n_i = n_e / Z`\ .
ion : `str` or `~plasmapy.particles.Particle`, optional
An instance of `~plasmapy.particles.Particle`, or a string
convertible to `~plasmapy.particles.Particle`.
kmax : `~astropy.units.Quantity`
Cutoff wavenumber (convertible to radians per meter). Defaults
to the inverse of the electron de Broglie wavelength.
Returns
-------
spectrum : `~astropy.units.Quantity`
Computed bremsstrahlung spectrum over the frequencies provided.
Notes
-----
For details, see "Radiation Processes in Plasmas" by
Bekefi. `ISBN 978\\-0471063506`_.
.. _`ISBN 978\\-0471063506`: https://ui.adsabs.harvard.edu/abs/1966rpp..book.....B/abstract
"""
# Default n_i is n_e/Z:
if n_i is None:
n_i = n_e / ion_species.integer_charge
# Default value of kmax is the electrom thermal de Broglie wavelength
if kmax is None:
kmax = (np.sqrt(const.m_e.si * const.k_B.si * T_e) / const.hbar.si).to(1 / u.m)
# Convert frequencies to angular frequencies
ω = (frequencies * 2 * np.pi * u.rad).to(u.rad / u.s)
# Calculate the electron plasma frequency
ω_pe = plasma_frequency(n=n_e, particle="e-")
# Check that all ω < wpe (this formula is only valid in this limit)
if np.min(ω) < ω_pe:
raise PhysicsError(
"Lowest frequency must be larger than the electron "
f"plasma frequency {ω_pe:.1e}, but min(ω) = {np.min(ω):.1e}"
)
# Check that the parameters given fall within the Rayleigh-Jeans limit
# hω << kT_e
rj_const = (
np.max(ω) * const.hbar.si / (2 * np.pi * u.rad * const.k_B.si * T_e)
).to(u.dimensionless_unscaled)
if rj_const.value > 0.1:
raise PhysicsError(
"Rayleigh-Jeans limit not satisfied: "
"hbar*ω/kT_e = {rj_const.value:.2e} > 0.1. "
"Try lower ω or higher T_e."
)
# Calculate the bremsstrahlung power spectral density in several steps
c1 = (
(8 / 3)
* np.sqrt(2 / np.pi)
* (const.e.si ** 2 / (4 * np.pi * const.eps0.si)) ** 3
* 1
/ (const.m_e.si * const.c.si ** 2) ** 1.5
)
Zi = ion_species.integer_charge
c2 = (
np.sqrt(1 - ω_pe ** 2 / ω ** 2)
* Zi ** 2
* n_i
* n_e
/ np.sqrt(const.k_B.si * T_e)
)
# Dimensionless argument for exponential integral
arg = 0.5 * ω ** 2 * const.m_e.si / (kmax ** 2 * const.k_B.si * T_e) / u.rad ** 2
# Remove units, get ndarray of values
arg = (arg.to(u.dimensionless_unscaled)).value
return c1 * c2 * exp1(arg)