"""This is where the chemistry scheme is implemented"""
import lir_achem.compute_ionisation as ci
import numpy as np
[docs]
def chemistry_eq(
t,
densities,
coefficients,
rad_here,
n_here,
compute_source,
EUV_file,
XR_file,
today,
):
"""Computes the derivative of the different species densities according to the chesmitry scheme
:param t: Time (in s) since the start of the computation
:param densities: Numpy array of species densities (cm-3). This is in the same (12, len(altitudes_D)). It is orders as
[Ne, O2m, Xm, NOp, O4p, Yp, Np, O2p, N2p, Op4S, Op2P, Op2D]
:param coefficients: Chemistry coefficients, as returned by the get_all_coefficients function
:param rad_here: Radiation class instance
:param n_here: Neutrals class instance
:param compute_source: Boolean. If True, the ionisation source will be recomputed at each time-step
:param EUV_file: Path to the EUV file of interest (from GOES EUV daily average)
:param XR_file: Path to the XR file of interest (from GOES 1s XR flux)
:param today: Datetime of the start of the computation
:returns: Derivatives of the densities, in the same format as densities"""
# =========================== Update the flux if needed =========================================
if compute_source: # Compute ionisation at that time
rad_here.update_flux(EUV_file, XR_file, t, today)
# Compute the absorbed fluxes
Phi_SXR, Phi_HXR, Phi_EUV = ci.compute_photon_flux(rad_here, HXR_bins=False)
# =========================== Compute the different ionisation sources =========================================
# O+ from O2
ionisation_Op = ci.ion_O2_to_Op(Phi_SXR, Phi_HXR, rad_here, n_here, HXR_bins=False)
# We differentiate the different energy levels of O+ using branching factors from Solomon & Qian (2005)
ionisation_Op4S = ionisation_Op * 0.390
ionisation_Op2D = ionisation_Op * 0.378
ionisation_Op2P = ionisation_Op * 0.224
# N2+ from N2
ionisation_N2p = ci.ion_N2_to_N2p(
Phi_SXR, Phi_HXR, rad_here, n_here, HXR_bins=False
)
# N+ from N2
ionisation_Np = ci.ion_N2_to_Np(Phi_SXR, Phi_HXR, rad_here, n_here, HXR_bins=False)
# NO+ from NO
ionisation_NOp = ci.ion_NO_to_NOp(Phi_EUV, rad_here, n_here)
# O2+ from O2(1Dg)
ionisation_O2p = ci.ionisation_O2_from_O21Dg(rad_here, n_here)
# Ne
ionisation_Ne = (
ionisation_N2p
+ ionisation_Np
+ ionisation_NOp
+ ionisation_Op4S
+ ionisation_Op2D
+ ionisation_Op2P
+ ionisation_O2p
)
# =========================== Densities =========================================
# Explicit the densities (so that's easier to debug)
length_one_density = int(np.size(densities) / 12)
Ne = densities[0:length_one_density]
O2m = densities[length_one_density : 2 * length_one_density]
Xm = densities[2 * length_one_density : 3 * length_one_density]
NOp = densities[3 * length_one_density : 4 * length_one_density]
O4p = densities[4 * length_one_density : 5 * length_one_density]
Yp = densities[5 * length_one_density : 6 * length_one_density]
Np = densities[6 * length_one_density : 7 * length_one_density]
O2p = densities[7 * length_one_density : 8 * length_one_density]
N2p = densities[8 * length_one_density : 9 * length_one_density]
Op4S = densities[9 * length_one_density : 10 * length_one_density]
Op2P = densities[10 * length_one_density : 11 * length_one_density]
Op2D = densities[11 * length_one_density :]
# =========================== Implement the chemistry scheme =========================================
# Electrons
dNe = (
ionisation_Ne
+ coefficients[6, :] * O2m
+ coefficients[7, :] * Xm
- (
coefficients[0, :] * N2p
+ coefficients[1, :] * O2p
+ coefficients[2, :] * NOp
+ coefficients[3, :] * O4p
+ coefficients[4, :] * Yp
)
* Ne
- coefficients[5, :] * Ne
)
# O2-
dO2m = (
coefficients[5, :] * Ne
- coefficients[6, :] * O2m
+ coefficients[29, :] * Xm
- coefficients[30, :] * O2m
- coefficients[37, :]
* (NOp + O4p + Yp + Np + O2p + N2p + Op4S + Op2P + Op2D)
* O2m
)
# X-
dXm = (
coefficients[30, :] * O2m
- coefficients[29, :] * Xm
- coefficients[7, :] * Xm
- coefficients[37, :]
* (NOp + O4p + Yp + Np + O2p + N2p + Op4S + Op2P + Op2D)
* Xm
)
# NO+
dNOp = (
-coefficients[2, :] * Ne * NOp
- coefficients[27, :] * NOp
+ coefficients[28, :] * Yp
+ coefficients[21, :] * O2p
+ coefficients[19, :] * N2p
+ coefficients[15, :] * Np
+ coefficients[8, :] * Op4S
+ ionisation_NOp
- coefficients[37, :] * (O2m + Xm) * NOp
)
# O4+
dO4p = (
-coefficients[3, :] * Ne * O4p
- coefficients[26, :] * O4p
- coefficients[23, :] * O4p
+ coefficients[22, :] * O2p
- coefficients[37, :] * (O2m + Xm) * O4p
)
# Y+
dYp = (
-coefficients[25, :] * Yp
+ coefficients[24, :] * O2p
+ coefficients[26, :] * O4p
- coefficients[28, :] * Yp
+ coefficients[27, :] * NOp
- coefficients[4, :] * Ne * Yp
- coefficients[37, :] * (O2m + Xm) * Yp
)
# N+
dNp = (
ionisation_Np
- (
coefficients[15, :]
+ coefficients[14, :]
+ coefficients[17, :]
+ coefficients[16, :]
+ coefficients[37, :] * (O2m + Xm)
)
* Np
)
# N2+
dN2p = (
ionisation_N2p
- coefficients[19, :] * N2p
- coefficients[0, :] * Ne * N2p
- coefficients[20, :] * N2p
+ coefficients[12, :] * Op2D
+ coefficients[17, :] * Np
+ coefficients[10, :] * Op2P
- coefficients[18, :] * N2p
- coefficients[37, :] * (O2m + Xm) * N2p
)
# O+(4S)
dOp4S = (
ionisation_Op4S
+ coefficients[16, :] * Np
- coefficients[9, :] * Op4S
+ (coefficients[31, :] + coefficients[32] * Ne) * Op2P
+ coefficients[20, :] * N2p
+ (coefficients[35, :] + coefficients[36] * Ne) * Op2D
- coefficients[8, :] * Op4S
- coefficients[37, :] * (O2m + Xm) * Op4S
)
# O+2P
dOp2P = (
ionisation_Op2P
- (
coefficients[11, :]
+ coefficients[10, :]
+ coefficients[33, :]
+ coefficients[34, :] * Ne
+ coefficients[31, :]
+ coefficients[32, :] * Ne
+ coefficients[37, :] * (O2m + Xm)
)
* Op2P
)
# O+2D
dOp2D = (
ionisation_Op2D
- (
coefficients[12, :]
+ coefficients[35, :]
+ coefficients[36, :] * Ne
+ coefficients[13, :]
+ coefficients[37, :] * (O2m + Xm)
)
* Op2D
+ (coefficients[33, :] + coefficients[34, :] * Ne) * Op2P
)
# O2+
# O2+ From chemistry (but ionisation doesn't countains UV)
dO2p = (
-coefficients[21, :] * O2p
- coefficients[22, :] * O2p
+ coefficients[23, :] * O4p
- coefficients[24, :] * O2p
+ coefficients[18, :] * N2p
+ coefficients[25, :] * Yp
+ coefficients[13, :] * Op2D
+ coefficients[11, :] * Op2P
+ coefficients[9, :] * Op4S
- coefficients[1, :] * Ne * O2p
+ coefficients[7, :] * Np
- coefficients[37, :] * (O2m + Xm) * O2p
+ ionisation_O2p
)
# O2+ From neutral condiction
# dO2p = (dNe + dO2m + dXm - dNOp - dO4p - dYp - dNp - dN2p - dOp4S - dOp2P - dOp2D)
return np.hstack(
(dNe, dO2m, dXm, dNOp, dO4p, dYp, dNp, dO2p, dN2p, dOp4S, dOp2P, dOp2D)
)