"""
Functions for calculating quantities associated with electromagnetic
"""

__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.frequencies import plasma_frequency
from plasmapy.particles.decorators import particle_input
from plasmapy.particles.particle_class import ParticleLike
from plasmapy.utils.decorators import validate_quantities
from plasmapy.utils.exceptions import PhysicsError

[docs]
@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.Quantity[u.Hz],
n_e: u.Quantity[u.m**-3],
T_e: u.Quantity[u.K],
n_i: u.Quantity[u.m**-3] = None,
ion: ParticleLike = "p+",
kmax: u.Quantity[u.m] = None,
) -> u.Quantity[u.kg * u.m**-1 * u.s**-2]:
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 quasineutral condition
:math:n_i = n_e / Z\ .

ion : |particle-like|, default: "p+"
An instance of ~plasmapy.particles.particle_class.Particle, or
a string convertible to |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 :cite:t:bekefi:1966\ .

Examples
--------
>>> import astropy.units as u
>>> import numpy as np
>>> thermal_bremsstrahlung(10**15 * u.Hz, 1e10 * u.cm**-3, 2e7 * u.K)  # solar flare
<Quantity 8.17560238e-23 kg / (m s2)>
>>> thermal_bremsstrahlung(
...     10 ** np.arange(15, 16, 0.1) * u.Hz, 1e22 * u.cm**-3, 1e2 * u.eV
... )
<Quantity [ 79.59052452, 117.73282254, 127.85119908, 127.12505588,
121.01549498, 112.02367743, 101.45553309,  90.04503155,
78.23475796,  66.32227273] kg / (m s2)>
>>> thermal_bremsstrahlung(
...     1e17 * u.Hz, 1e16 * u.cm**-3, 1e4 * u.eV, ion="Fe-56 12+"
... )
<Quantity 2.16932808e-10 kg / (m s2)>
"""

if n_i is None:  # default is quasineutrality
n_i = n_e / ion.charge_number

# Default value of kmax is the electron 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)

ω_pe = plasma_frequency(n=n_e, particle="e-")

# Check that all ω < ω_pe (this formula is only valid in this limit)
if np.min(ω) < np.max(ω_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: "
f"ℏω/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.charge_number
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)