"""This is where the Mitra-Rowe scheme is solve"""
import lir_achem.compute_ionisation as ci
import numpy as np
[docs]
def chemistry_mr_eq(
t,
densities,
coefficients,
rad_here,
n_here,
Phi_EUV,
compute_source,
today,
compute_sza=False,
HXR_bins=True,
):
"""Computes the derivative of the different species densities according to the Mitra-Rowe chemistry 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 shape (7, len(altitudes_D)). It is ordered as
[Ne, O2m, Xm, NOp, Yp, O2p, O4p, Phi_SXR, Phi_HXR]
:param coefficients: Chemistry coefficients, as returned by the get_mr_coefficients function
:param rad_here: Radiation class instance
:param n_here: Neutrals class instance
:param Phi_EUV: Absorbed EUV flux at D-region altitudes
:param compute_source: Boolean. If True, the ionisation source will be recomputed at each time-step
:param today: Datetime of the start of the computation (in UT)
:param compute_sza: Boolean. If True, the solar zenith angle will be recomputed at each time step and the absorption, Ch and H recomputed. Default:False
:param HXR_bins: Bollean, Default=True. If True, the ionisation will be computed using the bins in rad_here. If False, it will only use the average GOES HXR flux.
:returns: Derivatives of the densities, in the same format as densities"""
length_one_density = int(np.size(densities) / 9)
if (
compute_sza
): # Recompute the solar zenith angle, the absorption of solar flux, Ch and H
rad_here.update_chi(today, n_here)
if compute_source: # Compute ionisation at that time
rad_here.get_flux_now(t, today)
if (compute_source) or (compute_sza):
# Compute the absorbed fluxes
Phi_SXR, Phi_HXR, Phi_EUV = ci.compute_photon_flux(rad_here, HXR_bins)
if HXR_bins:
Phi_HXR = rad_here.Phi_HXR_bins
Phi_SXR = densities[7 * length_one_density : 8 * length_one_density]
if (not compute_source) and (not compute_sza) and (not HXR_bins):
Phi_SXR = densities[7 * length_one_density : 8 * length_one_density]
Phi_HXR = densities[8 * length_one_density : 9 * length_one_density]
# =========================== Densities =========================================
# Explicit the densities (so that's easier to debug)
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]
Yp = densities[4 * length_one_density : 5 * length_one_density]
O2p = densities[5 * length_one_density : 6 * length_one_density]
O4p = densities[6 * length_one_density : 7 * length_one_density]
# =========================== Compute the different ionisation sources =========================================
# Cosmic rays
ionisation_cr = ci.ion_from_cosmicrays(n_here, rad_here)
# NO+ from NO
ionisation_NOp = ci.ion_NO_to_NOp(Phi_EUV, rad_here, n_here)
# NO+ from conversion of O+ and N+
ionisation_NOp = ionisation_NOp + ci.ion_NOp_from_XR(
Phi_SXR, Phi_HXR, rad_here, n_here, HXR_bins=HXR_bins
)
# O2+ from O2(1Dg) and O+, N+ and N2+ conversion
ionisation_O2p = (
ci.ionisation_O2_from_O21Dg(rad_here, n_here)
+ ci.ionisation_O2_from_XR(
Phi_SXR, Phi_HXR, rad_here, n_here, HXR_bins=HXR_bins
)
+ ionisation_cr
)
# Ne
ionisation_Ne = ionisation_NOp + ionisation_O2p
# =========================== Implement the chemistry scheme =========================================
# Electrons
dNe = (
coefficients[2, :] * Xm
+ coefficients[3, :] * O2m
- coefficients[4, :] * Ne
- coefficients[9, :] * Ne * Yp
- coefficients[8, :] * Ne * O4p
- coefficients[6, :] * Ne * NOp
- coefficients[7, :] * Ne * O2p
+ ionisation_Ne
)
# O2m
dO2m = (
-coefficients[3, :] * O2m
- coefficients[5, :] * O2m
- coefficients[1, :] * (NOp + Yp + O2p + O4p) * O2m
+ coefficients[4, :] * Ne
)
# Xm
dXm = (
-coefficients[2, :] * Xm
- coefficients[13, :] * O4p * Xm
- coefficients[1, :] * (NOp + Yp + O2p) * Xm
+ coefficients[5, :] * O2m
)
# NO+
dNOp = (
ionisation_NOp
- coefficients[6, :] * Ne * NOp
- coefficients[1, :] * (O2m + Xm) * NOp
- coefficients[0, :] * NOp
+ coefficients[14, :] * O2p
)
# Yp
dYp = (
-coefficients[9, :] * Ne * Yp
- coefficients[1, :] * (O2m + Xm) * Yp
+ coefficients[12, :] * O4p
+ coefficients[0, :] * NOp
)
# O2+
dO2p = (
ionisation_O2p
- coefficients[14, :] * O2p
- coefficients[1, :] * (O2m + Xm) * O2p
- coefficients[10, :] * O2p
+ coefficients[11, :] * O4p
- coefficients[7, :] * Ne * O2p
)
# O4p
dO4p = (
-coefficients[8, :] * Ne * O4p
- coefficients[1, :] * O2m * O4p
- coefficients[12, :] * O4p
- coefficients[13, :] * O4p * Xm
+ coefficients[10, :] * O2p
- coefficients[11, :] * O4p
)
zero_array = np.zeros(np.shape(Phi_SXR))
return np.hstack((dNe, dO2m, dXm, dNOp, dYp, dO2p, dO4p, zero_array, zero_array))