Source code for lir_achem.compute_ionisation

"""Here is computed the ionisation rate from the absorbed solar flux"""

import numpy as np


[docs] def ion_from_cosmicrays(n_here, rad_here): """Computes the ionisation of O2 by cosmic rays. We use the approximation described in Lehtinen et al., 2007. We will suppose that CR produce both O2+ and e- :param n_here: Neutrals class instance :param rad_here: Radiation class instance""" Qp = 10 # cm^-3 hp = 15 # km N = n_here.He + n_here.O + n_here.N2 + n_here.O2 + n_here.Ar + n_here.N + n_here.AO # Get Np, value of N at hp id = np.argmin(np.abs(n_here.altitudes - hp)) Np = N[id] N_over_Np = N / Np Qr = Qp * N_over_Np * np.exp(1 - N_over_Np) return Qr[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)]
[docs] def ion_N2_to_N2p( Phi_SXR, Phi_HXR, rad_here, n_here, sigma_SXR_N2p=5.019e-22, sigma_HXR_N2p=5.115e-23, HXR_bins=True, ): """Computes the ionisation of N2 (in cm-3/s-1) to give N2+. The ionisation cross section come from Table 3 of Pavlov, 2014 and its erratum If HXR_bin, the cross-sections for HXR come from Siskind et al., 2022 and are multiplied by the branching factors in Solomon & Qian, 2005. pe/pi is taken from Siskind et al., 2022 Following Nicolet & Aikin, 1960 and the results from Bourdeau et al., 1966, each ionisation by X-ray is multiplied by the number of created ion pair, which is 45 for SXR and 165 for HXR :param Phi_SXR, Phi_HXR: Absorbed SXR and HXR fluxes :param rad_here: Radiation class instance :param n_here: Neutrals class instance :param sigma_SXR_N2p: Ionisation cross-section of N2 to give N2+ in cm^2 due to SXR (Default: 5.019e-22 cm^2) :param sigma_HXR_N2p: Ionisation cross-section of N2 to give N2+ in cm^2 due to HXR (Default: 5.115e-23 cm^2) :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: Ionisation rate (in part/cm^3 /s) in the D-region altitudes (given by the altitudes_D in rad_here) """ # Ionisation by SXR ionisation_N2p = ( (sigma_SXR_N2p * Phi_SXR * 45) * n_here.N2[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 ) # in m-3 s-1 if HXR_bins: # We add the correct terms bin by bin ionisation_N2p += ( Phi_HXR[0, :] * 3.38e-23 * 0.040 * 555.2 * n_here.N2[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 ) ionisation_N2p += ( Phi_HXR[1, :] * 1.49e-22 * 0.040 * 313.1 * n_here.N2[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 ) ionisation_N2p += ( Phi_HXR[2, :] * 4.55e-22 * 0.040 * 220.0 * n_here.N2[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 ) ionisation_N2p += ( Phi_HXR[3, :] * 1.03e-21 * 0.040 * 169.3 * n_here.N2[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 ) ionisation_N2p += ( Phi_HXR[4, :] * 2.04e-21 * 0.040 * 137.9 * n_here.N2[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 ) ionisation_N2p += ( Phi_HXR[5, :] * 3.91e-21 * 0.040 * 109.7 * n_here.N2[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 ) ionisation_N2p += ( Phi_HXR[6, :] * 3.91e-21 * 0.040 * 109.7 * n_here.N2[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 ) else: ionisation_N2p = ( ionisation_N2p + sigma_HXR_N2p * Phi_HXR * 165 * n_here.N2[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 ) # Conversion ionisation_N2p = ionisation_N2p * 1e-6 # in cm-3 s-1 return ionisation_N2p
[docs] def ion_N2_to_Np( Phi_SXR, Phi_HXR, rad_here, n_here, sigma_SXR_Np=2.320e-20, sigma_HXR_Np=2.351e-21, HXR_bins=True, ): """Computes the ionisation of N2 (in cm-3/s-1) to give N+. The ionisation cross section come from Table 3 of Pavlov, 2014 and its erratum Following Nicolet & Aikin, 1960 and the results from Bourdeau et al., 1966, each ionisation by X-ray is multiplied by the number of created ion pair, which is 45 for SXR and 165 for HXR :param Phi_SXR, Phi_HXR: Absorbed SXR and HXR fluxes :param rad_here: Radiation class instance :param n_here: Neutrals class instance :param sigma_SXR_Np: Ionisation cross-section of N2 to give N+ in cm^2 due to SXR (Default: 2.320e-20 cm^2) :param sigma_HXR_Np: Ionisation cross-section of N2 to give N+ in cm^2 due to HXR (Default: 2.351e-21 cm^2) :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. If it is True, the ionisation of N2 to N+ is put to 0 :returns: Ionisation rate (in part/cm^3 /s) in the D-region altitudes (given by the altitudes_D in rad_here) """ ionisation_Np = ( (sigma_SXR_Np * Phi_SXR * 45) * n_here.N2[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 ) # in m-3 s-1 if HXR_bins: # We add the correct terms bin by bin ionisation_Np += ( Phi_HXR[0, :] * 3.38e-23 * 0.960 * 555.2 * n_here.N2[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 ) ionisation_Np += ( Phi_HXR[1, :] * 1.49e-22 * 0.960 * 313.1 * n_here.N2[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 ) ionisation_Np += ( Phi_HXR[2, :] * 4.55e-22 * 0.960 * 220.0 * n_here.N2[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 ) ionisation_Np += ( Phi_HXR[3, :] * 1.03e-21 * 0.960 * 169.3 * n_here.N2[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 ) ionisation_Np += ( Phi_HXR[4, :] * 2.04e-21 * 0.960 * 137.9 * n_here.N2[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 ) ionisation_Np += ( Phi_HXR[5, :] * 3.91e-21 * 0.960 * 109.7 * n_here.N2[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 ) ionisation_Np += ( Phi_HXR[6, :] * 3.91e-21 * 0.960 * 109.7 * n_here.N2[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 ) else: ionisation_Np += ( (sigma_HXR_Np * Phi_HXR * 165) * n_here.N2[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 ) # in m-3 s-1 ionisation_Np = ionisation_Np * 1e-6 # in cm-3 s-1 return ionisation_Np
[docs] def ion_O2_to_Op( Phi_SXR, Phi_HXR, rad_here, n_here, sigma_HXR_Op=2.436e-21, sigma_SXR_Op=2.890e-20, HXR_bins=True, ): """Computes the ionisation of O2 or O (in cm-3/s-1) to give O+. The ionisation cross section come from Table 3 of Pavlov, 2014 Following Nicolet & Aikin, 1960 and the results from Bourdeau et al., 1966, each ionisation by X-ray is multiplied by the number of created ion pair, which is 45 for SXR and 165 for HXR NOTE: for the first HXR bin, Siskind et al., do not indicated pe/pi. We take the average value in Solomon & Qian, 2005 instead. We also do not discriminate between the different excited states of O+ If HXR_bin, the cross-sections for HXR come from Siskind et al., 2022 and are multiplied by the branching factors in Solomon & Qian, 2005. pe/pi is taken from Siskind et al., 2022 NOTE: The Branching ratio for the ionisation of O2+ is 0 in Solomon & Qian, 2005, but it was underestimated (see Pavlov, 2014). Instead, we recompute it from Table 3 of Pavlov, 2014 and Table A2 of Solomon & Qian, 2005. This gives beta=0.011 for O2 -> 02+, SO beta=0.989 for O2 -> 2O+ :param Phi_SXR, Phi_HXR: Absorbed SXR and HXR fluxes :param rad_here: Radiation class instance :param n_here: Neutrals class instance :param sigma_SXR_Op: Ionisation cross-section of O2 to give O+ in cm^2 due to SXR (Default: 2.890e-20 cm^2) :param sigma_HXR_Op: Ionisation cross-section of O2 to give O+ in cm^2 due to HXR (Default: 2.436e-21 cm^2) :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: Ionisation rate (in part/cm^3 /s) in the D-region altitudes (given by the altitudes_D in rad_here) """ # Ionisation from SXR ionisation_Op = ( (sigma_SXR_Op * Phi_SXR * 45) * n_here.O[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 ) # in m-3 s-1 if HXR_bins: # here we compute the ionisation from O+ and O2+ together ionisation_Op += Phi_HXR[0, :] * ( 3.1e-23 * 218.12 * n_here.O[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 + 1 * 6.1e-23 * 333.9 * n_here.O2[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 ) # in m-3 s-1 ionisation_Op += Phi_HXR[1, :] * ( 1.4e-22 * 213.7 * n_here.O[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 + 1 * 2.73e-22 * 187.2 * n_here.O2[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 ) # in m-3 s-1 ionisation_Op += Phi_HXR[2, :] * ( 4.2e-22 * 150.3 * n_here.O[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 + 1 * 8.32e-22 * 131.6 * n_here.O2[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 ) ionisation_Op += Phi_HXR[3, :] * ( 9.3e-22 * 116.6 * n_here.O[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 + 1 * 1.87e-21 * 102.1 * n_here.O2[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 ) ionisation_Op += Phi_HXR[4, :] * ( 1.8e-21 * 95.5 * n_here.O[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 + 1 * 3.66e-21 * 83.6 * n_here.O2[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 ) ionisation_Op += Phi_HXR[5, :] * ( 3.5e-21 * 76.6 * n_here.O[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 + 1 * 6.94e-21 * 67.1 * n_here.O2[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 ) ionisation_Op += Phi_HXR[6, :] * ( 3.5e-21 * 76.6 * n_here.O[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 + 1 * 6.94e-21 * 67.1 * n_here.O2[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 ) else: ionisation_Op += ( sigma_HXR_Op * Phi_HXR * 165 * n_here.O2[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 ) # in m-3 s-1 # Conversion ionisation_Op = ionisation_Op * 1e-6 # in cm-3 s-1 return ionisation_Op
[docs] def ion_O2_to_O2p( Phi_SXR, Phi_HXR, rad_here, n_here, sigma_HXR_O2p=2.625e-23, sigma_SXR_O2p=3.174e-22, HXR_bins=True, ): """Computes the ionisation of O2 (in cm-3/s-1) to give O2+. The ionisation cross section come from Table 3 of Pavlov, 2014 Following Nicolet & Aikin, 1960 and the results from Bourdeau et al., 1966, each ionisation by X-ray is multiplied by the number of created ion pair, which is 45 for SXR and 165 for HXR If HXR_bin, the cross-sections for HXR come from Siskind et al., 2022 and are multiplied by the branching factors in Solomon & Qian, 2005. pe/pi is taken from Siskind et al., 2022 NOTE: The Branching ratio for the ionisation of O2+ is 0 in Solomon & Qian, 2005, but it was underestimated (see Pavlov, 2014). Instead, we recompute it from Table 3 of Pavlov, 2014 and Table A2 of Solomon & Qian, 2005. This gives beta=0.011 :param Phi_SXR, Phi_HXR: Absorbed SXR and HXR fluxes :param rad_here: Radiation class instance :param n_here: Neutrals class instance :param sigma_SXR_O2p: Ionisation cross-section of O2 to give O+ in cm^2 due to SXR (Default: 2.625e-23 cm^2) :param sigma_HXR_O2p: Ionisation cross-section of O2 to give O+ in cm^2 due to HXR (Default: 3.174e-22 cm^2) :param HXR_bins: Boolean, 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: Ionisation rate (in part/cm^3 /s) in the D-region altitudes (given by the altitudes_D in rad_here) """ # Ionisation from SXR ionisation_Op = ( (sigma_SXR_O2p * Phi_SXR * 45) * n_here.O2[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 ) # in m-3 s-1 if HXR_bins: ionisation_Op += ( Phi_HXR[0, :] * 6.21e-23 * 0 * 333.9 * n_here.O2[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 ) ionisation_Op += ( Phi_HXR[1, :] * 2.73e-22 * 0 * 187.2 * n_here.O2[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 ) ionisation_Op += ( Phi_HXR[2, :] * 8.32e-22 * 0 * 131.6 * n_here.O2[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 ) ionisation_Op += ( Phi_HXR[3, :] * 1.87e-21 * 0 * 102.1 * n_here.O2[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 ) ionisation_Op += ( Phi_HXR[4, :] * 3.66e-21 * 0 * 83.6 * n_here.O2[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 ) ionisation_Op += ( Phi_HXR[5, :] * 6.94e-21 * 0 * 67.1 * n_here.O2[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 ) ionisation_Op += ( Phi_HXR[6, :] * 6.94e-21 * 0 * 67.1 * n_here.O2[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 ) else: ionisation_Op += ( sigma_HXR_O2p * Phi_HXR * 165 * n_here.O2[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e2 ) # Conversion ionisation_Op = ionisation_Op * 1e-6 # in cm-3 s-1 return ionisation_Op
[docs] def ionisation_O2_from_O21Dg(rad_here, n_here): """Computes the ionisation of O2(1Dg) from UV radiation As an approximation, following Crutzen et al., 1970 (Figure 1), we will take [O_2 (1Dg)] = 1e10 cm-3 in the atmopshere The ionisation of O2(1Dg) is computed following Paulsen, 1972 and Swider, 1979, and Krivolutsky, 2015 The ionisation dependence on solar activity is then computed from Pavlov, 2014 Equation 17 and 18. Instead however of using A10.7 in Equation 18, we simply use F10.7 :param rad_here: Radiation class instance :param n_here: Neutral class instance""" # Put O2 in the right shape O2 = n_here.O2[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] ionisation_O2 = 1e10 * ( 5.49e-10 * np.exp(-2.406e-20 * O2 * rad_here.Ch * rad_here.H * 1e5) + 2.614e-9 * np.exp(-8.508e-20 * O2 * rad_here.Ch * rad_here.H * 1e5) ) # Equation 18 ionisation_O2 = ionisation_O2 * (0.764 + 3.693e-3 * (n_here.f107 - 80)) return ionisation_O2
[docs] def ionisation_O2_from_XR(Phi_SXR, Phi_HXR, rad_here, n_here, HXR_bins=True): """Computes the ionisation of O2 with resutls from the ionisation of O+, N+ and N2+ then converted into O2+ At first approximation, all N2+ are converted to O2+ (Rowe, 1974, Mitra, 1974) The XR also create O+ and N+, which are converted into NO+ and O2+ :param Phi_SXR, Phi_HXR: Absorbed SXR and HXR fluxes :param rad_here: Radiation class instance :param n_here: Neutrals class instance :param HXR_bins: Boolean, default:True. If True, will compute the ionisation from HXR using the bins in rad_here. """ ionisation_N2p = ion_N2_to_N2p( Phi_SXR, Phi_HXR, rad_here, n_here, HXR_bins=HXR_bins ) ionisation_O2p = ion_O2_to_O2p( Phi_SXR, Phi_HXR, rad_here, n_here, HXR_bins=HXR_bins ) # O+ and N+ b_Op, b_Np = conversion_OpNp(rad_here, n_here) ionisation_Op = ion_O2_to_Op(Phi_SXR, Phi_HXR, rad_here, n_here, HXR_bins=HXR_bins) ionisation_Np = ion_N2_to_Np(Phi_SXR, Phi_HXR, rad_here, n_here, HXR_bins=HXR_bins) conversion_ions = ionisation_Op * (1 - b_Op) + ionisation_Np * (1 - b_Np) return ionisation_N2p + ionisation_O2p + conversion_ions
[docs] def ion_NOp_from_XR(Phi_SXR, Phi_HXR, rad_here, n_here, HXR_bins=True): """Computes the creation of NO+ with resutls from the ionisation of O+, N+ converted into NO+ :param Phi_SXR, Phi_HXR: Absorbed SXR and HXR fluxes :param rad_here: Radiation class instance :param n_here: Neutrals class instance :param HXR_bins: Boolean, if True the ionisation from HXR is computed from the bins in rad_here. (default:True) """ # O+ and N+ b_Op, b_Np = conversion_OpNp(rad_here, n_here) ionisation_Op = ion_O2_to_Op(Phi_SXR, Phi_HXR, rad_here, n_here, HXR_bins=HXR_bins) ionisation_Np = ion_N2_to_Np(Phi_SXR, Phi_HXR, rad_here, n_here, HXR_bins=HXR_bins) conversion_ions = ionisation_Op * b_Op + ionisation_Np * b_Np return conversion_ions
[docs] def conversion_OpNp(rad_here, n_here): """Computes the fraction of O+ and N+ converted into NO+ We suppose that O+ is created into its ground state (which is a simplification) and that O+ and N+ do not interact with other species EXCEPT to be converted into either O2+ or NO+. All reactions and reaction rates are from Pavlov, 2014, Table 6 :param rad_here: Radiation class instance :param n_here: Neutrals class instance :returns: - b_Op Fraction of O+ ions converted into NO+ - b_Np: Fraction of N+ ions converted into O2+""" start = int(rad_here.altitudes_D[0]) stop = int(rad_here.altitudes_D[-1]) + 1 O2 = n_here.O2[start:stop] N2 = n_here.N2[start:stop] Tn = n_here.Tn[start:stop] X = 300 / Tn # -------------------------- O+ ions # O+ + N2 -> N + NO+ (Reaction 1 of Table 6) rate_1 = 1e-12 * (2.05 - 0.00308 * Tn) * N2 # O+ + O2 -> O + O2+ (Reaction 2 of Table 6) rate_2 = 1.6e-11 * X**0.52 + 5.5e-11 * np.exp(-6832 / Tn) * O2 b_Op = rate_1 / (rate_2 + rate_1) # -------------------------- N+ ions # N+ + O2 -> O2+ + N (Reactions 14 and 15 of Table 6) rate_2 = (1.925e-10 + 8.25e-11) * O2 # N+ + O2 -> NO+ + O (Reactions 16 and 17 of Table 6) rate_1 = (4.95e-11 + 1.98e-10) * O2 b_Np = rate_1 / (rate_2 + rate_1) return b_Op, b_Np
[docs] def ion_NO_to_NOp(Phi_EUV, rad_here, n_here, sigma_EUV_NOp=2e-18): """Computes the ionisation of NO (in cm-3/s-1) to give NO+. The ionisation cross section come from Watanabe et al., 1967 :param Phi_EUV: Absorbed EUV flux :param rad_here: Radiation class instance :param n_here: Neutrals class instance :param sigma_EUV_NOp: Ionisation cross-section of NO to give NO+ in cm^2 due to SXR (Default: 2e-18 cm^2) :returns: Ionisation rate (in part/cm^3 /s) in the D-region altitudes (given by the altitudes_D in rad_here) """ ionisation_NOp = ( sigma_EUV_NOp * 1e-4 * Phi_EUV * n_here.NO[rad_here.altitudes_D[0] : (rad_here.altitudes_D[-1] + 1)] * 1e6 ) # in m-3 s-1 ionisation_NOp = ionisation_NOp * 1e-6 # in cm-3 s-1 return ionisation_NOp
[docs] def compute_photon_flux(rad_here, HXR_bins=True): """Here we just convert the solar flux (in W/m^2) to the photon flux (particle/m^2/s) :param rad_here: Radiation class instance :param HXR_bins: Boolean (Default: True). If True, Phi_HXR is an array (size 7) with the fluxes for each bin in rad_here """ # Constants (in units of m and s) clight = 3e8 h = 6.62e-34 # Flux of photons at average lambda Phi_SXR = rad_here.SXR * np.exp(-rad_here.tau_SXR) * 0.6e-9 / (h * clight) Phi_EUV = rad_here.EUV * np.exp(-rad_here.tau_EUV) * 121.6e-9 / (h * clight) if HXR_bins: Phi_HXR = np.zeros((7, np.size(rad_here.tau_HXR))) Phi_HXR[0, :] = ( rad_here.bin_0 * np.exp(-rad_here.tau_HXR_0) * 0.075e-9 / (h * clight) ) Phi_HXR[1, :] = ( rad_here.bin_1 * np.exp(-rad_here.tau_HXR_1) * 0.125e-9 / (h * clight) ) Phi_HXR[2, :] = ( rad_here.bin_2 * np.exp(-rad_here.tau_HXR_2) * 0.175e-9 / (h * clight) ) Phi_HXR[3, :] = ( rad_here.bin_3 * np.exp(-rad_here.tau_HXR_3) * 0.225e-9 / (h * clight) ) Phi_HXR[4, :] = ( rad_here.bin_4 * np.exp(-rad_here.tau_HXR_4) * 0.275e-9 / (h * clight) ) Phi_HXR[5, :] = ( rad_here.bin_5 * np.exp(-rad_here.tau_HXR_5) * 0.325e-9 / (h * clight) ) Phi_HXR[6, :] = ( rad_here.bin_6 * np.exp(-rad_here.tau_HXR_6) * 0.375e-9 / (h * clight) ) # Store this in radiation class rad_here.Phi_HXR_bins = Phi_HXR else: # Average over 0.05-0.4 nm Phi_HXR = rad_here.HXR * np.exp(-rad_here.tau_HXR) * 0.2250e-9 / (h * clight) return Phi_SXR, Phi_HXR, Phi_EUV