Source code for lir_achem.chemistry_coefficients

import numpy as np

# Here are computed the different coefficients for the chemical reactions

""" Each function computes the corresponding reaction rates, in units of cm^-3/s / charged_reactant

    We assume that the temperature stays below 300K in the D-region. References are given in each function
    The variable names match my own hand-written notes, but they are just used here and their names stay internal to each function
"""


[docs] def delta_1(N2, Tn): """O^+ (4S) -> NO+ This is the reaction rate of the reaction 1. in Table 6 of Pavlov, 2014 :param N2: N_2 neutral density (cm^-3) :param Tn: Neutral temperature (K)""" # Reaction rate A_1 = np.zeros(np.size(Tn)) Y = Tn / 300 A_1[Tn < 300] = (2.05 - 0.00308 * Tn[Tn < 300]) * 1e-12 A_1[Tn >= 300] = ( 1.71676e-12 - 7.19934e-13 * Y[Tn >= 300] + 1.33276e-13 * Y[Tn >= 300] ** 2 - 0.28213e-15 * Y[Tn >= 300] ** 3 + 6.39557e-16 * Y[Tn >= 300] ** 4 ) # Reaction rate * neutral reactant d1 = A_1 * N2 return d1
[docs] def delta_2(O2, Tn): """O^+ (4S) -> 0_2+ This is the reaction rate of the reaction 2. in Table 6 of Pavlov, 2014 :param O2: O_2 neutral density (cm^-3) :param Tn: Neutral temperature (K)""" # Reaction rate X = 300 / Tn A_2 = 1.6e-11 * X**0.52 + 5.5e-11 * np.exp(-6832 / Tn) # Reaction rate * Neutral reactant d2 = A_2 * O2 return d2
[docs] def epsilon_1(O): """O^+ (2P)) -> O^+ (4S) This is the reaction rate of the reactionS 3, 6 in Table 6 of Pavlov, 2014 :param O: O neutral density (cm^-3)""" # Reaction 3 B_1 = 5e-11 # Reaction 6 B_2 = 0.0775 # Reaction rate * reactant density eps_1 = B_1 * O + B_2 return eps_1
[docs] def epsilon_1p(Tn): """O^+ (2P)) -> O^+ (4S) This is the reaction rate of the reaction 8 in Table 6 of Pavlov, 2014 It is considered separately to the previous one because it doesn't involve neutral densities :param Tn: Neutral temperature (in K)""" X = 300 / Tn # Reaction 8 eps_1p = 2.5e-8 * X**0.5 return eps_1p
[docs] def delta_3(N2, Tn): """O+(2P) -> N2+ This is the reaction rate of the reaction 4 in Table 6 of Pavlov, 2014 :param N2: Density of N2 (cm^-3) :param Tn: Neutral temperature (in K)""" X = 300 / Tn # Reaction rate A_3 = 2e-10 * X**0.5 # Rate * density d3 = A_3 * N2 return d3
[docs] def delta_4(O2, Tn): """O+(2P) -> N2+ This is the reaction rate of the reaction 5 in Table 6 of Pavlov, 2014 :param O2: Density of O2 (cm^-3) :param Tn: Neutral temperature (in K)""" X = 300 / Tn # Reaction rate A_4 = 3.1e-10 * X**0.5 # Rate * density d4 = A_4 * O2 return d4
[docs] def epsilon_2(): """O^+ (2P)) -> O^+ (4S) This is the reaction rate of the reaction 7 in Table 6 of Pavlov, 2014 """ # Reaction 7 B_3 = 0.314 return B_3
[docs] def epsilon_2p(Tn): """O^+ (2P)) -> O^+ (4S) This is the reaction rate of the reaction 9 in Table 6 of Pavlov, 2014 :param Tn: Neutral temperature (in K)""" X = 300 / Tn # Reaction 9 B_5 = 7e-8 * X**0.5 # Rate * density eps_2p = B_5 return eps_2p
[docs] def epsilon_3(O): """O^+ (2D)) -> O^+ (4S) This is the reaction rate of the reaction 10 in Table 6 of Pavlov, 2014 :param O: Density of O (cm-3)""" # Reaction 10 B_6 = 5e-12 return B_6 * O
[docs] def epsilon_3p(Tn): """O^+ (2D)) -> O^+ (4S) This is the reaction rate of the reaction 13 in Table 6 of Pavlov, 2014 :param Tn: Neutral temperature (in K)""" X = 300 / Tn # Reaction 9 B_7 = 4e-8 * X**0.5 return B_7
[docs] def delta_5(Tn, N2): """O^+ (2D)) -> N2+ This is the reaction rate of the reaction 11 in Table 6 of Pavlov, 2014 :param Tn: Neutral temperature (in K) :param N2: N2 density (in cm-3)""" X = 300 / Tn # Reaction rate A_5 = 1.5e-10 * X**0.5 # Rate * density d5 = A_5 * N2 return d5
[docs] def delta_6(Tn, O2): """O^+ (2D)) -> O2+ This is the reaction rate of the reaction 12 in Table 6 of Pavlov, 2014 :param Tn: Neutral temperature (in K) :param O2: O2 density (in cm-3)""" X = 300 / Tn # Rate A_6 = 1e-10 * X**0.5 # Rate * density d6 = A_6 * O2 return d6
[docs] def delta_7(O2): """N+ -> O2+ This is the reaction rate of the reactions 14 and 15 in Table 6 of Pavlov, 2014 :param O2: O2 density (in cm-3)""" # Reaction 14 A_7 = 1.925e-10 # Reaction 15 A_8 = 8.25e-11 # Reaction rate * density d7 = (A_7 + A_8) * O2 return d7
[docs] def delta_8(O2, NO, Tn): """N+ -> NO+ This is the reaction rate of the reactions 16, 17 and 20 in Table 6 of Pavlov, 2014 :param O2: O2 density (in cm-3) :param NO: NO density (in cm-3) :param Tn: Neutral temperature (in K) """ # Reaction 16 A_9 = 4.95e-11 # Reaction 17 A_10 = 1.98e-10 # Reaction 20 A_13 = 5.72e-9 * Tn ** (-0.44) d8 = (A_9 + A_10) * O2 + A_13 * NO return d8
[docs] def delta_9(O, O2): """N+ -> O+(4S) This is the reaction rate of the reactions 18 and 19 in Table 6 of Pavlov, 2014 :param O2: O2 density (in cm-3) :param O: O density (in cm-3) """ # Reaction 18 A_11 = 2.75e-11 # Reaction 19 A_12 = 2.2e-12 d9 = A_11 * O2 + A_12 * O return d9
[docs] def delta_10(NO, Tn): """N+ -> N2+ This is the reaction rate of the reaction 21 in Table 6 of Pavlov, 2014 :param NO: NO density (in cm-3) :param Tn: Neutral temperature (in K) """ A_14 = 7.15e-10 * Tn ** (-0.44) return A_14 * NO
[docs] def alpha_1(Tn): """N2+ + Ne -> N2 This is the reaction rate of the reaction 22 in Table 6 of Pavlov, 2014 :param Tn: Neutral temperature (in K) """ X = 300 / Tn return 2.2e-7 * X ** (0.39)
[docs] def delta_11(O2, Tn): """N2+ -> O2+ This is the reaction rate of the reaction 23 in Table 6 of Pavlov, 2014 :param Tn: Neutral temperature (in K) :param O2: O2 density (in cm-3) """ Y = Tn / 300 A_15 = ( 1.285 - 1.428 * Y + 0.7656 * Y**2 - 0.2014 * Y**3 + 2.663e-2 * Y**4 - 1.389e-3 * Y**5 ) * 1e-10 return A_15 * O2
[docs] def delta_12(O, NO, Tn): """N2+ -> NO+ This is the reaction rate of the reactions 24 and 26 in Table 6 of Pavlov, 2014 :param Tn: Neutral temperature (in K) :param O: O density (in cm-3) :param NO: NO density (in cm-3) """ X = 300 / Tn # Reaction 24 A_16 = 1.4e-10 * X**0.44 - 7.0e-12 * X**0.23 # Reaction 26 A_18 = 7.5e-9 * Tn ** (-0.52) return A_16 * O + A_18 * NO
[docs] def delta_13(O, Tn): """N2+ -> O+(4S) This is the reaction rate of the reaction 25 in Table 6 of Pavlov, 2014 :param Tn: Neutral temperature (in K) :param O: O density (in cm-3) """ X = 300 / Tn A_17 = 7.0e-12 * X**0.23 return A_17 * O
[docs] def alpha_2(Tn): """O2+ + Ne -> O + O This is the reaction rate of the reaction 27 in Table 6 of Pavlov, 2014 :param Tn: Neutral temperature (in K) :param O: O density (in cm-3) """ X = 300 / Tn a2 = np.zeros(np.size(Tn)) a2[Tn < 1200] = 1.95e-7 * X[Tn < 1200] ** (0.7) a2[Tn >= 1200] = 1.93e-7 * X[Tn >= 1200] ** (0.61) return a2
[docs] def delta_14(NO, N): """O2+ -> NO+ This is the reaction rate of the reactions 28 and 29 in Table 6 of Pavlov, 2014 :param NO: NO density (in cm-3) :param N: N density (in cm-3) """ # Reaction 28 A_19 = 4.1e-10 # Reaction 19 A_20 = 1e-10 return A_19 * NO + A_20 * N
[docs] def alpha_3(Tn): """NO+ + Ne -> N + O This is the reaction rate of the reaction 30 in Table 6 of Pavlov, 2014 :param Tn: Neutral temperature (in K) """ X = 300 / Tn a3 = np.zeros(np.size(Tn)) a3[Tn < 1200] = 3.5e-7 * X[Tn < 1200] ** 0.69 a3[Tn >= 1200] = 3.02e-7 * X[Tn >= 1200] ** 0.56 return a3
[docs] def delta_15(O2, N2, Tn): """O2+ -> O4+ This is the reaction rate of the reaction 36 in Table 6 of Pavlov, 2014 :param O2: O2 density (in cm-3) :param N2: N2 density (in cm-3) :param Tn: Neutral temperature (in K) """ M = N2 + O2 X = 300 / Tn A_21 = 4e-30 * X**2.93 return A_21 * M * O2
[docs] def delta_16(O2, N2, O, Tn): """O4+ -> O2+ This is the reaction rate of the reactions 37 and 44 in Table 6 of Pavlov, 2014 :param O2: O2 density (in cm-3) :param N2: N2 density (in cm-3) :param Tn: Neutral temperature (in K) """ X = 300 / Tn M = O2 + N2 # Reaction 37 A_22 = 1.3e-6 * X ** (3.93) * np.exp(-4607 / Tn) # Reaction 44 A_29 = 3e-10 return A_22 * M + A_29 * O
[docs] def delta_17(N2, O2, H2O, Tn): """O2+ -> Y+ This is the reaction rate of the reactions 38, 39 and 41 in Table 6 of Pavlov, 2014 :param O2: O2 density (in cm-3) :param N2: N2 density (in cm-3) :param H2O: H2O density (in cm-3) :param Tn: Neutral temperature (in K) """ M = N2 + O2 X = 300 / Tn # Reaction 38 A_23 = 1e-30 * X**3.2 # Reaction 39 A_24 = 2.8e-28 # Reaction 41 A_26 = 2.3e-28 return A_23 * N2 * M + A_24 * H2O * N2 * A_26 * H2O * O2
[docs] def delta_18(N2, O2, Tn): """Y+ -> O2+ This is the reaction rate of the reactions 40 and 42 in Table 6 of Pavlov, 2014 :param O2: O2 density (in cm-3) :param N2: N2 density (in cm-3) :param Tn: Neutral temperature (in K) """ # Reaction 40 A_25 = 0.017 / Tn * np.exp(-8047 / Tn) # Reaction 42 A_27 = 0.014 / Tn * np.exp(-8047 / Tn) return A_25 * N2 + A_27 * O2
[docs] def delta_19(H2O, N2, O2): """O4+ -> Y+ This is the reaction rate of the reactions 43 and 46 in Table 6 of Pavlov, 2014 :param O2: O2 density (in cm-3) :param N2: N2 density (in cm-3) :param H2O: H2O density (in cm-3) """ M = O2 + N2 # Reaction 43 A_28 = 1.2e-9 # Reaction 46 A_31 = 4.1e-29 return A_28 * H2O + A_31 * N2 * M
[docs] def delta_20(O2, N2, H2O, Tn): """NO+ -> Y+ This is the reaction rate of the reactions 48, 49, 53, 55 in Table 6 of Pavlov, 2014 All three-bodies reaction involving CO2 are neglected, as they should negligible in front of those involving N2 and 02 From Rowe & Mitra, 1974 :param O2: O2 density (in cm-3) :param N2: N2 density (in cm-3) :param H2O: H2O density (in cm-3) :param Tn: Neutral temperature (in K) """ # X = 300 / Tn M = N2 + O2 # # Reaction 48 # A_48 = 3e-31 # # Reaction 49 # A_49 = 3e-31 * X**4.3 # # Reaction 53 # A_53 = 1.35e-28 * X ** (2.837) # # Reaction 55 # A_55 = 8.28e-29 * X**2.837 # return A_48 * O2 * M + A_49 * N2 * M + A_53 * H2O * N2 + A_55 * H2O * O2 # return 1e-31 * M**2 # Mitra, 1975 return 0.01 # Rowe @ Mitra, 1974
[docs] def delta_21(O2, N2, Tn): """Y+ -> NO+ This is the reaction rate of the reactions 50, 52, 54 and 56 in Table 6 of Pavlov, 2014 :param O2: O2 density (in cm-3) :param N2: N2 density (in cm-3) :param Tn: Neutral temperature (in K) """ X = 300 / Tn M = N2 + O2 # Reaction 50 A_50 = 3.8e-8 * X**5.3 * np.exp(-2321 / Tn) # Reaction 52 A_52 = 3.4e-7 * X**5 * np.exp(-3872 / Tn) # Reaction 54 A_54 = 3.5e-4 * X**3.837 * np.exp(-9316 / Tn) # Reaction 56 A_56 = 2.2e-4 * X**3.837 * np.exp(-9316 / Tn) return A_50 * M + A_52 * M + A_54 * N2 + A_56 * O2
[docs] def alpha_4(Tn): """O4+ + e -> 2 O2 This is the reaction rate of the reaction 110 in Table 6 of Pavlov, 2014 :param Tn: Neutral temperature (in K) """ X = 300 / Tn return 4.2e-6 * X**0.48
def alpha_5(): """Y+ + e -> products This value is taken from Mitra & Rowe, 1972""" return 1e-7
[docs] def beta_1(O2, N2, Tn): """e -> O2- This is the reaction rate of the reactions 1 and 2 of Table 10 of Pavlov, 2014 :param O2: O2 density (in cm-3) :param N2: N2 density (in cm-3) :param Tn: Neutral temperature (in K) """ X = 300 / Tn # Reaction 1 b_1 = 1.4e-29 * np.exp(-600 / Tn) * X # Reaction 2 b_2 = 1.07e-31 * X**2 * np.exp(-70 / Tn) return b_1 * O2**2 + b_2 * O2 * N2
[docs] def gamma_1(N2, O2, O, Tn): """O2- -> e This is the reaction rate of the reactions 26, 27, 31 of Table 10 of Pavlov, 2014 :param O2: O2 density (in cm-3) :param N2: N2 density (in cm-3) :param O: O density (in cm-3) :param Tn: Neutral temperature (in K) """ X = 300 / Tn # Reaction 26 g_1 = 1.9e-12 * X ** (-1.5) * np.exp(-4990 / Tn) # Reaction 27 g_2 = 2.7e-10 * X ** (-0.5) * np.exp(-5590 / Tn) # Reaction 31 g_3 = 2.1e-10 return g_1 * N2 + g_2 * O2 + g_3 * O
[docs] def delta_22(): """X- -> O2- This reaction is neglected for the moment (Mitra & Rowe, 1972)""" return 0
[docs] def delta_23(): """O2- -> X- From Mitra & Rowe, 1972""" return 0.5
[docs] def gamma_2(): """X- -> e From Mitra & Rowe, 1972""" return 5e-3
[docs] def alpha_5(): """X- + e -> products From Mitra & Rowe, 1972""" return 1e-7
[docs] def alpha_i(): """Recombination of each positive ion to the negative ones From Mitra & Rowe, 1972""" return 1e-7
[docs] def get_all_coefficients(n_here, e_here): """Get all the reactions coefficients in order :param n_here: Neutrals class instance :param e_here: Electrons class instance :returns: Array with all the coefficients. Ordered from alpha_1 to alpha_5, beta_1, gamma_1, gamma_2 and delta_1 to delta_23, epsilon_1, epsilon_1p, epsilon_2, epsilon_2p, epsilon_3, epsilon_3p """ result = np.zeros((38, len(e_here.altitudes))) start = int(e_here.altitudes[0]) stop = int(e_here.altitudes[-1]) + 1 Tn = n_here.Tn[start:stop] O2 = n_here.O2[start:stop] N2 = n_here.N2[start:stop] O = n_here.O[start:stop] NO = n_here.NO[start:stop] H2O = n_here.H2O[start:stop] N = n_here.N[start:stop] result[0, :] = alpha_1(Tn) result[1, :] = alpha_2(Tn) result[2, :] = alpha_3(Tn) result[3, :] = alpha_4(Tn) result[4, :] = alpha_5() result[5, :] = beta_1(O2, N2, Tn) result[6, :] = gamma_1(N2, O2, O, Tn) result[7, :] = gamma_2() result[8, :] = delta_1(N2, Tn) result[9, :] = delta_2(O2, Tn) result[10, :] = delta_3(N2, Tn) result[11, :] = delta_4(O2, Tn) result[12, :] = delta_5(Tn, N2) result[13, :] = delta_6(Tn, O2) result[14, :] = delta_7(O2) result[15, :] = delta_8(O2, NO, Tn) result[16, :] = delta_9(O, O2) result[17, :] = delta_10(NO, Tn) result[18, :] = delta_11(O2, Tn) result[19, :] = delta_12(O, NO, Tn) result[20, :] = delta_13(O, Tn) result[21, :] = delta_14(NO, N) result[22, :] = delta_15(O2, N2, Tn) result[23, :] = delta_16(O2, N2, O, Tn) result[24, :] = delta_17(N2, O2, H2O, Tn) result[25, :] = delta_18(N2, O2, Tn) result[26, :] = delta_19(H2O, N2, O2) result[27, :] = delta_20(O2, N2, H2O, Tn) result[28, :] = delta_21(O2, N2, Tn) result[29, :] = delta_22() result[30, :] = delta_23() result[31, :] = epsilon_1(O) result[32, :] = epsilon_1p(Tn) result[33, :] = epsilon_2() result[34, :] = epsilon_2p(Tn) result[35, :] = epsilon_3(O) result[36, :] = epsilon_3p(Tn) result[37, :] = alpha_i() return result