"""This is the module computing the source of the ionisation in the D-region"""
import numpy as np
import scipy as sp
[docs]
def compute_tau(rad_here, neutral_profiles, chi):
"""Computes the absorbed flux at D-regions altitudes
:param rad_here: radiation class instance
:param neutral_profiles: Class instance of neutrals
:param altitudes_D: Altitude of interest (km - Only one altitude here)
:param chi: Solar zenith angle (in °)
:returns: Values of tau for EUV, SXR and HXR wavelengths. tau_HXR_bins is an array with taus computed for bins of 0.05 nm length (e.g. 0.05-0.1 nm).
tau is returned because if the neutrals densities do not vary on solar flare timescales, then tau should stay the same. In this case, it is only necessary to compute it once, and not at each iteration
"""
# Definition of the cross-section (N2 and O2 are the only ones responsible for the absorption)
# Absorption of SXR
O2_SXR_cs = 3.4e-20 # Table 1, Pavlov, 2014; cm^2
O_SXR_cs = 1.7e-20
N2_SXR_cs = 2.01e-20
# Absorption of HXR
O2_HXR_cs = 4.5e-21 # Table 1, Pavlov, 2014
O_HXR_cs = 2.3e-21
N2_HXR_cs = 2.5e-21
tau_EUV = np.zeros(np.shape(rad_here.altitudes_D))
tau_SXR = np.zeros(np.shape(rad_here.altitudes_D))
tau_HXR_bins = np.zeros((7, np.size(rad_here.altitudes_D)))
tau_HXR = np.zeros(np.shape(rad_here.altitudes_D)) # Averaged HXR flux
# Compute Ch and H
Ch, H = chapman(chi, neutral_profiles)
# Absorption of Lyman alpha
# Only O2 absorbes the flux at that wavelength, with a cross section given by Huffman, 1969
O2_EUV_cs = compute_Lya_absorption(neutral_profiles, chi, Ch, H)
for i in range(len(rad_here.altitudes_D)):
altitude_here = rad_here.altitudes_D[i]
tau_EUV[i] = get_tau(
chi, neutral_profiles, 0, O2_EUV_cs, 0, altitude_here, Ch, H
) # Equation (4) of Pavlov, 2014
tau_SXR[i] = get_tau(
chi, neutral_profiles, O_SXR_cs, O2_SXR_cs, N2_SXR_cs, altitude_here, Ch, H
)
tau_0 = get_tau(
chi, neutral_profiles, 3.1e-23, 6.21e-23, 3.38e-23, altitude_here, Ch, H
)
tau_1 = get_tau(
chi, neutral_profiles, 1.4e-22, 2.73e-22, 1.49e-22, altitude_here, Ch, H
)
tau_2 = get_tau(
chi, neutral_profiles, 4.2e-22, 8.32e-22, 4.55e-22, altitude_here, Ch, H
)
tau_3 = get_tau(
chi, neutral_profiles, 9.3e-22, 1.87e-21, 1.03e-21, altitude_here, Ch, H
)
tau_4 = get_tau(
chi, neutral_profiles, 1.8e-21, 3.66e-21, 2.04e-21, altitude_here, Ch, H
)
tau_5 = get_tau(
chi, neutral_profiles, 3.5e-21, 6.94e-21, 3.91e-21, altitude_here, Ch, H
)
tau_6 = tau_5
# Put this into tau_HXR
tau_HXR_bins[0, i] = tau_0
tau_HXR_bins[1, i] = tau_1
tau_HXR_bins[2, i] = tau_2
tau_HXR_bins[3, i] = tau_3
tau_HXR_bins[4, i] = tau_4
tau_HXR_bins[5, i] = tau_5
tau_HXR_bins[6, i] = tau_6
tau_HXR[i] = get_tau(
chi,
neutral_profiles,
O_HXR_cs,
O2_HXR_cs,
N2_HXR_cs,
altitude_here,
Ch,
H,
)
return tau_EUV, tau_SXR, tau_HXR, tau_HXR_bins, Ch, H
[docs]
def compute_Lya_absorption(n_here, chi, Ch=None, H=None):
"""Compute the absorption of Lyman-alpha radiation from Reddmann and Uhl, 2003
:param n_here: Neutrals class instance
:param chi: Solar zenith angle (in °)
:param Ch, H: Chapman function and atmosphere scale height (optional). If they are not given, they are recomputed
:returns: Absorption cross-section (in cm2) of O2 for Lyman-alpha"""
# Compute Ch and H if not given
if Ch is None or H is None:
Ch, H = chapman(chi, n_here)
# Compute N(O2) = [O2] H Ch
NO2 = n_here.O2 * Ch * H
# Get F(O2)
F = F_O2(chi, NO2)
# Compute cross-section
sigma = (
7.65e-21
* (1 + 0.35 * np.exp(-5e-21 * NO2))
* F
* (1.16 - 0.0021 * n_here.Tn + 6e-6 * n_here.Tn**2)
)
return sigma
[docs]
def F_O2(chi, NO2):
"""Function necessary to compute the absorption cross-section of O2 for Lyman-a
Taken from Reddmann & Uhl, 2003 and Pavlov, 2014
:param chi: Solar zenith angle (in °)
:param NO2: N(O2)=[O2]HCh, which is computed in the compute_Lya_absorption
:return: F, value of the function for different solar zenith angles"""
cx = np.cos(chi * np.pi / 180)
if cx < 0.8:
F = 1.1 + 0.1 * np.tanh(3e-21 * NO2 * (0.8 - cx) - 2.4)
else:
F = 1.1 + 0.1 * np.tanh(-2.4)
return F
[docs]
def get_tau(chi, neutrals_here, cs_O, cs_O2, cs_N2, altitude_here, Ch=None, H=None):
"""Computes tau in Equation (4) of Pavlov, 2014
:param chi: Solar zenith angle
:param neutrals_here: Neutrals class instance
:param cs_O: Cross section of O at that wavelength in cm2
:param cs_O2: Cross section of O2 at that wavelength
:param cs_N2: Cross section of N2 at that wavelength
:param altitude_here: Altitude (in km) of the computation
:param Ch: Chapman function (if computed before). Default: None
:param H: Atmosphere scale-height (in km, if computed previously). Default: None
:returns: - tau"""
id = np.argmin(np.abs(neutrals_here.altitudes - altitude_here))
# For Lyman-a absorption, the parametrisation of the cross-section depends on height
if np.size(cs_O2) > 1:
cs_O2 = cs_O2[id]
# Get the correction from the Chapman function (curvature of the atmosphere)
if Ch is None:
Ch, H = chapman(chi, neutrals_here)
# Sum on the neutrals and their cross-sections
tau = (
neutrals_here.O[id] * cs_O
+ neutrals_here.O2[id] * cs_O2
+ neutrals_here.N2[id] * cs_N2
) # in cm-1
tau = tau * H[id] * 1e5 # 1e5 because H was in km (1 km = 1e5 cm)
# Get the values for the right altitude
Ch_here = Ch[id]
return tau * Ch_here
[docs]
def chapman(chi, neutrals_here):
"""Computes the integral of the Chapman function. It is based on Brasseur & Solomon (2005)
if the solar zenith angle is less than 90°, and on Smith & Smith (1972) if it is above
:param chi: Solar zenith angle (in °)
:param neutrals_here: neutrals class instance
:returns: Ch, value of the Chapman function and H, height scale (in km). Both are arrays, with the same size as the altitudes in neutrals_here
"""
# Scale-height
H = compute_H(neutrals_here)
# Constant Earth-radius
Re = 6367.088132098377173 # km
X = (Re + neutrals_here.altitudes) / H
Y = np.abs(np.cos(chi * np.pi / 180)) * np.sqrt(0.5 * X)
if chi <= 90:
Ch = np.sqrt(0.5 * np.pi * X) * f_chapman(Y)
else:
Ch = np.sqrt(2 * np.pi * X) * (
np.sqrt(
np.sin(chi * np.pi / 180) * np.exp(X * (1 - np.sin(chi * np.pi / 180)))
- 0.5 * f_chapman(Y)
)
)
return Ch, H
[docs]
def compute_H(neutrals_here):
"""Computes the height scale in the atmosphere
:param neutrals_here: 'neutrals' class instance
:returns: H height scale (in km). Array with the same size as the densities in neutrals_here
"""
# Constants
k = 1.3806505e-23
g = 9.81
# Compute total density (in cm-3)
tot_density = (
+ neutrals_here.N2
+ neutrals_here.O2
)
# Mean mass (in kg)
M = neutrals_here.Total / tot_density
# Scale height (in km)
H = k * neutrals_here.Tn / (M * g) * 1e-3 # in km
return H
[docs]
def f_chapman(Y):
"""Function 'f' necessary for the computation of the Chapman function according to Smith & Smith (1972)"""
f = np.exp(Y**2) * sp.special.erfc(Y)
return f