"""This modules computes the HF absorption"""
import igrf
import numpy as np
import scipy.integrate as sci
from scipy.integrate import solve_ivp
import pytz
from pysolar.solar import get_altitude
import matplotlib.pyplot as plt
import pickle
import lir_achem.compute_collision_frequency as ccf
import lir_achem.class_definition as cd
import lir_achem.compute_coefficients_mitra as ccm
import lir_achem.mitra_rowe_scheme as mrs
import lir_achem.compute_ionisation as ci
[docs]
def compute_total_absorption(
f, today, glat, glon, e_here, n_here, ions, pm, abs_coeff="AH", ignore_B=False
):
"""Compute the total HF absorption for a wave launched vertically at frequency f
as described in Davies, 1990 (around p215)
:param f: Frequency of the wave (Hz)
:param today: Datetime, time of the computation
:param glat/glon: Latitude/Longitude for the computation
:param e_here: Electron class instance
:param n_here: Neutrals class instance
:param ions: List of ions class instance (ions participating in the absorption)
:param pm: Defines whether we use '+' or '-' in the denominator of the Appleton-Hartree equation. It correspond to choosing a polarisation for the propagating wave. If pm = 1, we will choose '+', if not, we will chose '-'
:param abs_coeff: String, chooses which form of the absorbtion coefficent we use. 'HA'-> We solve the Appleton Hartree equation (default). 'Eccles'-< we use the form proposed in Eccles et al., 2005
:param ignore_B: Ignore B field. Default: False
"""
# Theta and B for different altitudes for phi = 0 and phi = 180
# Theta and B for different altitudes
theta_0 = np.zeros(np.shape(e_here.altitudes))
B_0 = np.zeros(np.shape(e_here.altitudes))
theta_180 = np.zeros(np.shape(e_here.altitudes))
B_180 = np.zeros(np.shape(e_here.altitudes))
for i in range(np.size(e_here.altitudes)):
theta_h, B_h = magnetic_theta_and_B(
today, glat, glon, e_here.altitudes[i], phi=0
)
theta_0[i] = theta_h
B_0[i] = B_h
theta_h, B_h = magnetic_theta_and_B(
today, glat, glon, e_here.altitudes[i], phi=180
)
theta_180[i] = theta_h
B_180[i] = B_h
# Absorption coefficients
if abs_coeff == "AH":
abs_coeff_up, fN, mu = absorption_coefficient(
today,
glat,
glon,
e_here.altitudes,
n_here,
ions,
e_here,
f,
pm,
0,
ignore_B,
theta=theta_0,
B=B_0,
)
abs_coeff_down, fN, mu = absorption_coefficient(
today,
glat,
glon,
e_here.altitudes,
n_here,
ions,
e_here,
f,
pm,
180,
ignore_B,
theta=theta_180,
B=B_180,
)
elif abs_coeff == "Eccles":
abs_coeff_up = absorption_coefficient_eccles(n_here, e_here, ions, f)
abs_coeff_down = abs_coeff_up
fN, dummy = compute_frequencies(e_here, B_0, e_here.altitudes)
else:
print("Unrecognised keyword for the absorption coefficient")
return np.zeros(np.shape(e_here.altitudes))
# Bounds for the integration
alt_reflection = highest_altitude_wave(fN, f, e_here.altitudes)
min_altitude = np.min(e_here.altitudes)
# Constrain the arrays between the bounds
id_min = np.argmin(np.abs(min_altitude - e_here.altitudes))
id_max = np.argmin(np.abs(alt_reflection - e_here.altitudes))
abs_coeff = (
abs_coeff_up[id_min : id_max + 1] + abs_coeff_down[id_min : id_max + 1]
) # dB/km
altitudes = e_here.altitudes[id_min : id_max + 1] # km
# Integration
L = sci.simpson(abs_coeff, altitudes) * (-8.68)
return L, abs_coeff, altitudes
[docs]
def highest_altitude_wave(fN, f, altitudes):
"""Computes the highest altitude reached by the wave. If its frequency is below the plasma frequency, it is reflected at this altitude. If not, the function returns the highest altitude of the D-region
:param fN: Array, plasma frequency (in Hz). Its size must match the one of e_here
:param f: Frequency of the wave (in Hz)
:param altitudes: Altitudes array (in km)
:returns: altitude_reflected: Altitude at which the wave is reflected"""
id = np.argmin(np.abs(f - fN))
altitude_ref = altitudes[id]
if f > np.max(fN):
altitude_ref = np.max(altitudes)
return altitude_ref
[docs]
def absorption_coefficient_eccles(n_here, e_here, ions, f):
"""Compute the absorption coefficients at the given altitude (from Davies, 1990)
The notations of the code are those of Eccles et al., 2005 (Equation 3)
:param n_here: Neutrals class instance
:param e_here: Electron class instance
:param ions: List of ions class instance, ions participating in the absorption
:param f: Frequency of the VLF wave (in Hz)
:returns: kappa, absorption coefficent (same shape than altitudes, in dB/km)
"""
vei = np.zeros(np.shape(e_here.altitudes))
for i in ions:
vei += ccf.ei_collisionfreq(e_here, i)
ve = ccf.en_collisionfreq(e_here, n_here) + vei
kappa = 4.6e-2 * (e_here.densities * 1e6 * ve / (ve**2 + (2 * np.pi * f) ** 2))
return kappa
[docs]
def absorption_coefficient(
today,
glat,
glon,
altitudes,
n_here,
ions,
e_here,
f,
pm,
phi,
ignore_B=False,
theta=None,
B=None,
):
"""Compute the absorption coefficients at the given altitude (from Davies, 1990)
The notations of the code are those of Zawdie et al., 2017 (Equation 2)
:param today: Datetime, date at the computation
:param glat/glon: Latitude/longitude of the point
:param altitudes: 1D Array, altitudes for the computation (in km)
:param n_here: Neutrals class instance
:param ions: List of ions class instance
:param e_here: Electron class instance
:param f: Frequency of the VLF wave (in Hz)
:param pm: Defines whether we use '+' or '-' in the denominator of the Appleton-Hartree equation. It correspond to choosing a polarisation for the propagating wave. If pm = 1, we will choose '+', if not, we will chose '-'
:param phi: Angle of the wave with the vertical (in °)
:param ignore_B: Ignore magnetic field (Boolean). Default: False
:param theta: Angle between the wave and the magnetic field for each altitude (default: None)
:param B: Norm of the magnetic field at each altitude (Default: None)
:returns: kappa, absorption coefficent (same shape than altitudes, in dB/km)
"""
if (theta is None) and (B is None):
# Theta and B for different altitudes
theta = np.zeros(np.shape(altitudes))
B = np.zeros(np.shape(altitudes))
for i in range(np.size(altitudes)):
theta_h, B_h = magnetic_theta_and_B(today, glat, glon, altitudes[i], phi)
theta[i] = theta_h
B[i] = B_h
# n and mu
n_2, vah, fecf, fN = appleton_hartree(
altitudes, n_here, ions, e_here, f, pm, theta, B
) # n^2
mu = np.real(np.sqrt(n_2 + 0j)) # Real part of n
# The two polarisations
sign = -1
if pm == 1:
sign = 1
if ignore_B:
sign = 0
kappa = np.zeros(np.shape(mu))
# kappa
kappa[mu > 0] = (
4.6e-2
/ mu[mu > 0]
* e_here.densities[mu > 0]
* 1e6
* vah[mu > 0]
/ (
vah[mu > 0] ** 2
+ (
2 * np.pi * f
+ sign
* np.abs(2 * np.pi * fecf[mu > 0] * np.cos(theta[mu > 0] * np.pi / 180))
)
** 2
)
)
# The absolute value in the formula above is present in the origin article by Appleton & Pigott (1954)
return kappa, fN, mu
[docs]
def appleton_hartree(altitudes, n_here, ions, e_here, f, pm, theta, B, ignore_B=True):
"""Compute the square of the index of refraction at the given altitude
The notations of the code are those of Zawdie et al., 2017 (Equations 3 to 6)
:param n_here: Neutrals class instance
:param ions: List of ions class instance
:param e_here: Electron class instance
:param f: Frequency of the VLF wave (in Hz)
:param pm: Defines whether we use '+' or '-' in the denominator of the Appleton-Hartree equation. It correspond to choosing a polarisation for the propagating wave. If pm = 1, we will choose '+', if not, we will chose '-'
:param theta: Angle of the wave with the B field (Array, same size as altitudes)
:param B: Magnetic field (in T). Array, same size as altitudes
:returns: n_2, square of the refractive index. n_2 has the same shape as altitudes. vah, total collision frequency. theta, angle (in °) of the wave compared to the B field (to avoid computing it again). fecf, electron cyclotron frequency (in s-1) to avoid recomputing
"""
# Initialisation of n_2 and fecf
n_2 = np.zeros(np.shape(altitudes))
fecf = np.zeros(np.shape(altitudes))
fN = np.zeros(np.shape(altitudes))
# Collision frequencies
# Compute electron-ions collision frequency
vei = np.zeros(np.shape(altitudes))
for ion in ions:
vei += ccf.ei_collisionfreq(e_here, ion)
# Compute electron-neutrals collision frequency
ven = ccf.en_collisionfreq(e_here, n_here)
# Total collision frequency
vah = vei + ven
# Denominator
sign = 1
if pm != 1:
sign = -1
for i in range(np.size(altitudes)):
id = np.argmin(np.abs(altitudes[i] - e_here.altitudes))
# Compute plasma and cyclotron frequencies
fN_1, fecf_1 = compute_frequencies(e_here, B[i], altitudes[i])
fecf[i] = fecf_1
fN[i] = fN_1
# X, Y, Z
X = (fN_1 / f) ** 2
Y = fecf_1 / f
Zah = vah[id] / (2 * np.pi * f)
if ignore_B:
Y = 0
# Numerator of Equation (3)
num = 1 - 1j * Zah - X
den = (
2 * (1 - 1j * Zah) * (1 - 1j * Zah - X)
- Y**2 * np.sin(theta[i] * np.pi / 180) ** 2
+ sign
* np.sqrt(
Y**4 * np.sin(theta[i] * np.pi / 180) ** 4
+ 4
* Y**2
* np.cos(theta[i] * np.pi / 180) ** 2
* (1 - 1j * Zah - X) ** 2
)
)
n_2[i] = 1 - 2 * X * num / den
return n_2, vah, fecf, fN
[docs]
def compute_frequencies(e_here, B, altitude=75):
"""Compute the plasma frequency and the cyclotron frequency at altitude 'altitude'
:param e_here: Electrons class instance
:param B: Norm of the magnetic field (may be obtained through magnetic_theta_and_B)
:param altitude: Altitude at the computation (in km). Default: 75 km
:return: fN, plasma frequency (s-1) and fecf, electron collision frequency (s-1)
"""
# Constants
e = 1.602e-19 # Charge electron (C)
me = 9.109e-31 # Mass electron (kg)
eO = 8.854e-12 # Vaccum permittivity
# Get the correct altitude
id = np.argmin(np.abs(altitude - e_here.altitudes))
# Plasma frequency
fN = np.sqrt(e_here.densities[id] * 1e6 * e**2 / (eO * me))
# Cyclotron frequency
fecf = e * B / me
return fN, fecf
[docs]
def magnetic_theta_and_B(today, glat, glon, altitude=75, phi=0):
"""Computes the angle of the field to the vertical and the total magnetic field
:param today: Datetime, time for the computation
:param glat/glon: Latitude and longitude for the computation
:param altitude: Altitude at the computation (in km). Default: 75 km
:param phi: Angle of the wave with the vertical (in °, 0° is upwards). Default: 0
:returns: theta, angle of the magnetic field and the vertical, B, norm of the field
"""
mag = igrf.igrf(today, glat, glon, altitude)
# Get the right values from mag
incl = mag.incl.values[
0
] # Inclination of the magnetic field (by convention, it is positive when the B field points downwards)
B = mag.total.values[0] * 1e-9 # Total norm, T
theta = incl + 90 - phi # Angle between the B field and a wave launched vertically
return theta, B
[docs]
def get_HF_absorption_from_scratch(
today,
glon,
glat,
EUV_file,
XR_file,
time_end,
f,
chi=None,
AH_O=True,
AH_X=True,
Ec=True,
ignore_B=False,
time_initialisation=4e3,
time_HF=60,
alpha=None,
compute_sza=False,
HXR_bins=True,
max_altitude=90,
file_manual_ini_neutrals=None,
file_manual_ini_ions=None,
time_first_ini = 40
):
"""Runs the entire pipeline to get an estimate of the HF absorption at that location
:param today: Datetime, must be timezone-blind BUT the time must correspond to the UT time
:param glon/glat: Geographic longitude and latitude (in °)
:param EUV_file: Filename with the daily EUV average (from GOES)
:param XR_file: Filename with the 1s XR flux (from GOES)
:param time_end: Maximum number of second after today when we run the simulation. today MUST be in quiet time since it is used for the initialisation
:param f: HF frequency of the waves (in Hz). Maybe an array
:param chi: If it is specified, will run everything as normal BUT the sza will be set to the specified value (Default: None)
:param AH_O: Boolean, if True computes the absorption of a HF vertical wave (O-mode) using Appleton-Hartree equation (Default:True)
:param AH_X: Boolean, if True computes the absorption of a HF vertical wave (X-mode) using Appleton-Hartree equation (Default:True)
:param Ec: Boolean, if True computes the absorption of a HF vertical wave using the equation in Eccles, 2005 (Default:True)
:param ignore_B: Boolean. If True, we ignore the contribution of the magnetic field. Default:False
:param time_initialisation: Number of seconds needed to initialise the ions the second time (Default:4000)
:param time_HF: Time resolution to compute the HF absorption (Default:60, meaning one point out of 60 in the output)
:param alpha: If specified, corrects [NO] with alpha, thus only initialising ions once (Default: None)
: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: Boolean, default:True. If True, the ionisation from HXR will be discretized in bins. If False, it will only be caused by the average HXR wavelength of GOES data.
:param max_altitude: Maximum altitude of the D-region (Default, 90). Physically, the only constraint is the ionisation of O2(1Dg), which is only valid below 90 km (Paulsen, 1972). However, we can set this to None, in which case it will be limited by FIRI, or to another value above 60 km
:param file_manual_ini_neutrals: Default: None. If a path to a pickle file is specified with precomputed neutral profiles, the neutrals will be initialised from the file. Note: This file must contain exactly the attributes of the neutrals class instance (see documentation)
:param file_manual_ini_ions: Default: None. If a path to a pickle file is specified with precomputed ions profiles, the ions will be initialised from the file. Note: This file must contain all ions and electron profiles (see documentation for the required format)
:param time_first_ini: Time (in s) for the first initialisation (without adjusted NO). Default: 40
:returns: A dictionnary, countaining the different ion and electron densities through time, the time and altitudes array abd the estimation of the HF absorption using the Appleton Hartree equation or the one from Eccles et al., (2005)
"""
# -------------------------- INITIALISATION OF NEUTRALS, ELECTRONS and FLUXES ---------------------------------------
print("Initialisation of the neutrals .... ")
# Altitudes arrays
altitudes_neutrals = np.arange(0, 2000, 1)
# Initialisation of profiles
if file_manual_ini_neutrals is None:
n_here = cd.neutrals(today, altitudes_neutrals, glat, glon) # Neutrals
else:
print("[Neutral initialisation from file]")
n_here = cd.neutrals(
today,
altitudes_neutrals,
glat,
glon,
manual_ini=True,
ini_file=file_manual_ini_neutrals,
)
print("Done !")
print("Initialisation of the electrons.....")
e_here = cd.electrons(glat, glon, today, n_here.f107, n_here) # Electrons
if max_altitude is not None:
# Restrict the electrons to 60 - 90 km (optionnal)
# Note: We should not go above 90 km because the computation of the production of O2p from Lyman-alpha is only valid below
id_min = np.argmin(np.abs(e_here.altitudes - 60))
id_max = np.argmin(np.abs(e_here.altitudes - int(max_altitude)))
e_here.altitudes = e_here.altitudes[id_min : id_max + 1]
e_here.densities = e_here.densities[id_min : id_max + 1]
e_here.temperature = e_here.temperature[id_min : id_max + 1]
print("Done !")
print("Initialisation of the radiation fluxes....")
# Get date in good format
time_here_aware = pytz.utc.localize(today)
# IMPORTANT: They must be integers, and match the e_here altitude array if we want to run the chemistry scheme
altitudes_D = e_here.altitudes.astype(np.int32)
if chi is None:
# Compute solar zenith angle
chi = float(90) - get_altitude(glat, glon, time_here_aware)
print("Solar zenith angle :", chi)
# Flux initialisation
rad_here = cd.radiation(today, EUV_file, XR_file, altitudes_D, n_here, chi)
# Compute photon fluxes, using the taus already computed
Phi_SXR, Phi_HXR, Phi_EUV = ci.compute_photon_flux(rad_here, HXR_bins)
print("Done !")
# ------------------------------ FIRST INITIALISATION OF IONS ------------------------------------------------------------
# Get average Phi_HXR even if HXR_bins is True (for size in hstack below, but it is recomputed anyways in chemistry_mr_eq)
dummy1, Phi_HXR_1dim, dummy2 = ci.compute_photon_flux(rad_here, HXR_bins=False)
print("Initialisation of the ions....")
if file_manual_ini_ions is None:
if alpha is None:
# Load the chemistry coefficients
coefficients = ccm.get_mr_coeffs(n_here, e_here)
# Initial conditions
zero_density = np.zeros(np.shape(e_here.densities))
densities = np.hstack(
(
zero_density,
zero_density,
zero_density,
zero_density,
zero_density,
zero_density,
zero_density,
Phi_SXR,
Phi_HXR_1dim,
)
)
# Run model
sol_0 = solve_ivp(
mrs.chemistry_mr_eq,
[0, time_first_ini],
densities,
args=(
coefficients,
rad_here,
n_here,
Phi_EUV,
False,
today,
False,
HXR_bins,
),
)
# Get everything in the right format
length_one_density = int(np.size(e_here.densities))
Ne = sol_0.y[0:length_one_density, :]
O2m = sol_0.y[length_one_density : 2 * length_one_density, :]
Xm = sol_0.y[2 * length_one_density : 3 * length_one_density, :]
NOp = sol_0.y[3 * length_one_density : 4 * length_one_density, :]
Yp = sol_0.y[4 * length_one_density : 5 * length_one_density, :]
O2p = sol_0.y[5 * length_one_density : 6 * length_one_density, :]
O4p = sol_0.y[6 * length_one_density : 7 * length_one_density, :]
# Check the error in electron density at 75 km
id = np.argmin(np.abs(e_here.altitudes - 75))
alpha = e_here.densities[id] / Ne[id, -1]
n_here.NO = n_here.NO * alpha
print("Alpha : ", alpha)
print("FIRI: ", e_here.densities[id])
print("Modeled: ", Ne[id, -1])
print("Done !")
else:
print("Alpha specified :", alpha)
# ------------------------------- SECOND INITIALISATION OF IONS WITH AJUSTED NO ---------------
print("Initialisation of ions with ajusted [NO]....")
# Load the chemistry coefficients
coefficients = ccm.get_mr_coeffs(n_here, e_here)
# Initial conditions
zero_density = np.zeros(np.shape(e_here.densities))
densities = np.hstack(
(
zero_density,
zero_density,
zero_density,
zero_density,
zero_density,
zero_density,
zero_density,
Phi_SXR,
Phi_HXR_1dim,
)
)
sol_1 = solve_ivp(
mrs.chemistry_mr_eq,
[0, time_initialisation],
densities,
args=(
coefficients,
rad_here,
n_here,
Phi_EUV,
False,
today,
False,
HXR_bins,
),
)
# Get everything in the right format
length_one_density = int(np.size(e_here.densities))
Ne = sol_1.y[0:length_one_density, :]
O2m = sol_1.y[length_one_density : 2 * length_one_density, :]
Xm = sol_1.y[2 * length_one_density : 3 * length_one_density, :]
NOp = sol_1.y[3 * length_one_density : 4 * length_one_density, :]
Yp = sol_1.y[4 * length_one_density : 5 * length_one_density, :]
O2p = sol_1.y[5 * length_one_density : 6 * length_one_density, :]
O4p = sol_1.y[6 * length_one_density : 7 * length_one_density, :]
Phi_SXR = sol_1.y[7 * length_one_density : 8 * length_one_density, :]
Phi_HXR = sol_1.y[8 * length_one_density : 9 * length_one_density, :]
# We take the densities given at the end of the initialisation
densities = np.hstack(
(
Ne[:, -1],
O2m[:, -1],
Xm[:, -1],
NOp[:, -1],
Yp[:, -1],
O2p[:, -1],
O4p[:, -1],
Phi_SXR[:, -1],
Phi_HXR[:, -1],
)
)
print("Done !")
else: # An input file has been specified for the initialisation of ions
print("[Ion initialisation from file]")
# Open file:
with open(file_manual_ini_ions, "rb") as file_here:
input = pickle.load(file_here)
densities = np.hstack(
(
input["Ne"],
input["O2m"],
input["Xm"],
input["NOp"],
input["Yp"],
input["O2p"],
input["O4p"],
Phi_SXR,
Phi_HXR_1dim,
)
)
# Put correct density in e_here
e_here.densities = input["Ne"]
length_one_density = int(np.size(e_here.densities))
# Load the chemistry coefficients
coefficients = ccm.get_mr_coeffs(n_here, e_here)
print("Done !")
# ------------------------------- DYNAMICS DURING THE SOLAR FLARES ----------------------------
print("Variations of densities during flares....")
# Implementation of the scheme (aroung 2.5 min for each 1/2 hour in the modelling)
sol = solve_ivp(
mrs.chemistry_mr_eq,
[0, time_end],
densities,
args=(
coefficients,
rad_here,
n_here,
Phi_EUV,
True,
today,
compute_sza,
HXR_bins,
),
method="BDF",
)
the_time = sol.t
print("Done !")
# -------------------------------- HF ABSORPTION --------------------------------------
# Initialisation of the ion classes
# O2m
O2m_c = cd.ion_species(e_here.altitudes, -1)
# Xm
Xm_c = cd.ion_species(e_here.altitudes, -1)
# NOp
NOp_c = cd.ion_species(e_here.altitudes, 1)
# O4p
O4p_c = cd.ion_species(e_here.altitudes, 1)
# Yp
Yp_c = cd.ion_species(e_here.altitudes, 1)
# O2p
O2p_c = cd.ion_species(e_here.altitudes, 1)
print("Computation of HF absorption....")
HF_time_array = the_time[0::time_HF]
# HF_absorption array
HF_abs_O = np.zeros((np.size(HF_time_array), np.size(f)))
HF_abs_X = np.zeros((np.size(HF_time_array), np.size(f)))
HF_abs_Ec = np.zeros((np.size(HF_time_array), np.size(f)))
# Get HF absorption for each time in the_time
for t in range(np.size(HF_time_array)):
i = t * time_HF
# Get everything in the right format
e_here.densities = sol.y[0:length_one_density, i]
O2m_c.density = sol.y[length_one_density : 2 * length_one_density, i]
Xm_c.density = sol.y[2 * length_one_density : 3 * length_one_density, i]
NOp_c.density = sol.y[3 * length_one_density : 4 * length_one_density, i]
Yp_c.density = sol.y[4 * length_one_density : 5 * length_one_density, i]
O2p_c.density = sol.y[5 * length_one_density : 6 * length_one_density, i]
O4p_c.density = sol.y[6 * length_one_density : 7 * length_one_density, i]
# Positive ions contributing to the collision frequencies
ions = [NOp_c, O4p_c, Yp_c, O2p_c]
for k in range(np.size(f)):
# O-mode
if AH_O:
HF, dummy1, dummy2 = compute_total_absorption(
f[k], today, glat, glon, e_here, n_here, ions, 1, "AH", ignore_B
)
HF_abs_O[t, k] = HF
# X-mode
if AH_X:
HF, dummy1, dummy2 = compute_total_absorption(
f[k], today, glat, glon, e_here, n_here, ions, 0, "AH", ignore_B
)
HF_abs_X[t, k] = HF
# Eccles
if Ec:
HF, dummy1, dummy2 = compute_total_absorption(
f[k], today, glat, glon, e_here, n_here, ions, 1, "Eccles", ignore_B
)
HF_abs_Ec[t, k] = HF
print("Done !")
# -------------------------------- OUPUT --------------------------------------
print("Preparation output....")
# Put everything in a dictionary to return
the_time = today.hour + today.minute / 60 + today.second / 3600 + the_time / 3600
la_date = today.strftime("%Y%m%d")
# Return the fluxes
time_XR = rad_here.XR_times
id_min = np.argmin(np.abs(the_time[0] - time_XR))
id_max = np.argmin(np.abs(the_time[-1] - time_XR))
SXR = rad_here.SXR_array[id_min : id_max + 1]
HXR = rad_here.HXR_array[id_min : id_max + 1]
time_XR = time_XR[id_min : id_max + 1]
output = {
"the_time": the_time,
"altitudes": e_here.altitudes,
"electrons": sol.y[0:length_one_density, :],
"O2m": sol.y[length_one_density : 2 * length_one_density, :],
"Xm": sol.y[2 * length_one_density : 3 * length_one_density, :],
"NOp": sol.y[3 * length_one_density : 4 * length_one_density, :],
"Yp": sol.y[4 * length_one_density : 5 * length_one_density, :],
"O2p": sol.y[5 * length_one_density : 6 * length_one_density, :],
"O4p": sol.y[6 * length_one_density : 7 * length_one_density, :],
"EUV": rad_here.EUV,
"SXR": SXR,
"HXR": HXR,
"time_XR": time_XR,
"tau_SXR": rad_here.tau_SXR,
"tau_HXR": rad_here.tau_HXR,
"tau_EUV": rad_here.tau_EUV,
"HF_AH_O": HF_abs_O,
"HF_AH_X": HF_abs_X,
"HF_Ec": HF_abs_Ec,
"f": f,
"time_HF": HF_time_array,
"date": la_date,
"glat": glat,
"glon": glon,
"sza": chi,
"alpha": alpha,
"HXR_bins": HXR_bins,
}
print("Done ! Everything is finished.")
return output