# Here the different clases are defined
from pyfiri.firi import firi2018
from pysolar.solar import get_altitude
import datetime as dt
from datetime import datetime
import msise00
import numpy as np
import pytz
import netCDF4
import pickle
import lir_achem.compute_source as cs
[docs]
class ion_species:
"""Each species is declared as a class, with two attributes: their altitudes (in km) and number densities (in cm-3)
Only altitudes is required for the initialisation, as all others are initialised to zeros
:param altitudes: Altitudes in the D-region (in km)
:param density: Number density (in cm-3)
:param z: Number charge of the ion
"""
def __init__(self, altitudes, z):
"""All densities are initialised at zeros"""
self.altitudes = altitudes
self.density = np.zeros(np.shape(altitudes))
self.z = z
[docs]
class electrons:
"""Each species is declared as a class, with two attributes: their altitudes (in km) and number densities (in cm-3)
The electrons are considered separately as they are initialized from the FIRI model
Class attributes:
:param altitudes: Altitudes in the D-region (in km)
:param densities: Number densities (in cm-3) in the D-region
:param temperature: Temperature (in K), assumed to be the same as the neutrals
Parameters for the initialisation:
:param glat: Latitude (in °) of the point of interest
:param glon: Longitude (in °) of the point of interest
:param time_here: Datetime, time of the initialisation. It must contain the information on the timezone in order to compute the sza
:param f107: f10.7 index of the day in question
:param n_here: Neutrals class instance
"""
def __init__(self, glat, glon, time_here, f107, n_here):
"""
:param glat: Latitude (in °) of the point of interest
:param glon: Longitude (in °) of the point of interest
:param time_here: Datetime, time of the initialisation. It must contain the information on the timezone in order to compute the sza
:param f107: f10.7 index of the day in question
:param n_here: Neutrals class instance
"""
# Make time_here aware of the timezone
time_here_aware = pytz.utc.localize(time_here)
# Compute solar zenith angle
sza = float(90) - get_altitude(glat, glon, time_here_aware)
# Compute day of year
doy = time_here_aware.timetuple().tm_yday
# Initialise the electron densities
model = firi2018()
profile = model.get_profile(doy=doy, chi=sza, lat=glat, f10_7=f107)
self.densities = profile.values.squeeze() * 1e-6 # Cm-3
self.altitudes = profile["alt"].values
id_min = np.argmin(np.abs(self.altitudes[0] - n_here.altitudes))
id_max = np.argmin(np.abs(self.altitudes[-1] - n_here.altitudes))
self.temperature = n_here.Tn[id_min : id_max + 1]
[docs]
class neutrals:
"""The neutral species are declared apart from the others. Because they dominate in the D-region, we consider that they do not interact in the ion dynamics other than as a reservoir for the chemical reactions
All neutral densities are initialised from the MSISE-00 model.
Parameters for the initialisation:
:param time_here: Time of interest (datetime, must be timezone blind)
:param altitudes: Altitude array (in km)
:param glat: Latitude of interest (in °)
:param glon: Longitude of interest (in °)
Class attributes:
:param glat: Latitude (in °) of the point of interest
:param glon: Longitude (in °) of the point of interest
:param altitudes: Altitudes in the D-region (in km)
:param N2: Number densities (in cm-3) of N2 in the D-region
:param O2: Number densities (in cm-3) of O2 in the D-region
:param O: Number densities (in cm-3) of O in the D-region
:param He: Number densities (in cm-3) of He in the D-region
:param Ar: Number densities (in cm-3) of Ar in the D-region
:param H: Number densities (in cm-3) of H in the D-region
:param N: Number densities (in cm-3) of N in the D-region
:param AO: Number densities (in cm-3) of anomalous O in the D-region
:param NO: Number densities (in cm-3) of NO in the D-region
:param H2O: Number densities (in cm-3) of H2O in the D-region
:param Tn: Temperature (in K) in the D-region
:param Total: Total mass density (in kg.cm-3)
:param f107: Value of f107. It is not an attribute of the neutral densities themselves, but it is easier to store it here
:param M: Addition of O2 and N2 (here because we use it a lot)
:param O3: Number density (in cm-3) of ozone in the D-region
"""
def __init__(
self, time_here, altitudes, glat, glon, alpha=1, manual_ini=False, ini_file=None
):
"""Initialisation thorught the MSISE00 model
:param time_here: Time of interest (datetime, must be timezone blind)
:param altitudes: Altitude array (in km)
:param glat: Latitude of interest (in °)
:param glon: Longitude of interest (in °)
:param alpha: Float, number to multiply the NO density from Mitra 1966 to match FIRI electron densities at 75 km
:param manual_ini: Boolean (Default: False). If True, initialise the neutrals class from a precomputed pickle file given by 'ini_file'
:param ini_file: Path to the initialisation file, if manual_ini=True.
"""
if manual_ini: # Initialisation from pre-computed profiles
if ini_file is not None:
self.init_manually(ini_file)
else: # The initialisation file has not been specified
print("!!!!!!! The neutral initialisation file is missing !!!!!!!!")
print("Initialisation of neutrals as usual")
manual_ini = False
if not manual_ini: # Initialisation from scratch
self.glat = glat
self.glon = glon
atmos = msise00.run(time_here, altitudes, glat, glon)
# Variables based on temperature
Tn = atmos["Tn"].squeeze()
Tn = np.array(Tn.values)
self.Tn = Tn
# Neutral densities (cm^-3)
self.N2 = np.array(atmos["N2"].squeeze() * 1e-6)
self.O2 = np.array(atmos["O2"].squeeze() * 1e-6)
self.O = np.array(atmos["O"].squeeze() * 1e-6)
self.He = np.array(atmos["He"].squeeze() * 1e-6)
self.Ar = np.array(atmos["Ar"].squeeze() * 1e-6)
self.H = np.array(atmos["H"].squeeze() * 1e-6)
self.N = np.array(atmos["N"].squeeze() * 1e-6)
self.Total = np.array(atmos["Total"].squeeze() * 1e-6)
self.AO = np.array(atmos["AnomalousO"].squeeze() * 1e-6)
self.altitudes = atmos.alt_km.values
# Extend O to lower altitudes, assuming it is constant below 73 km (MSISE limit)
# See for example graph p282 of Brasseur & Solomon, 2005
self.O[0:73] = self.O[73] * np.ones(np.shape(self.O[0:73]))
# NO density (from Mitra, 1966)
self.NO = (1e-2 * np.exp(-3300 / self.Tn) * self.O2 + 5e-7 * self.O) * alpha
# Value of f107 from msise
self.f107 = atmos.f107
# From Reid, 1971
# We assume a mixing ratio of 3 ppm
self.H2O = (
3e-6
/ (1 + 3e-6)
* (self.N2 + self.O2 + self.O + self.He + self.Ar + self.H + self.N)
)
# From Levine, 1985 (p176)
self.M = self.N2 + self.O2
k2 = 6e-34 * (self.Tn / 300) ** (-2.3)
J3 = 9100e-6
J4 = 1100e-6
JO3 = J3 + J4
self.O3 = k2 * self.M * self.O2 * self.O / JO3
[docs]
def init_manually(self, file):
"""Initialise manually the neutrals from a pickle file previously computed
:param file: Path to the pickle file"""
# Open the file
with open(file, "rb") as f:
input = pickle.load(f)
# Distribute everything into the correct attributes
self.glat = input["glat"]
self.glon = input["glon"]
self.altitudes = input["altitudes"]
self.N2 = input["N2"]
self.O2 = input["O2"]
self.O = input["O"]
self.He = input["He"]
self.Ar = input["Ar"]
self.H = input["H"]
self.N = input["N"]
self.AO = input["AO"]
self.NO = input["NO"]
self.H2O = input["H2O"]
self.Tn = input["Tn"]
self.Total = input["Total"]
self.f107 = input["f107"]
self.M = input["M"]
self.O3 = input["O3"]
[docs]
class radiation:
"""EUV, hard X-rays and soft X-rays from the Sun. The absorbed fluxes are given by
abs_flux = flux*np.exp(-tau)
where tau depends on the wavelength
:param HXR: Hard X-ray flux at the time demanded
:param SXR: Soft X-ray flux at the time demanded
:param EUV: EUV flux at the time demanded
:param EUV_times: Time array for EUV data (from GOES daily average file)
:param EUV_array: Lyman-alpha data from GOES
:param XR_times: Time array for the XR from GOES
:param SXR_array: SXR data from GOES
:param HXR_array: HXR data from GOES
:param tau_HXR: Absorbed hard X-ray flux
:param tau_SXR: Absorbed soft X-ray flux
:param tau_EUV: Absorbed EUV flux
:param Ch: Chapman function for all altitudes
:param H: Atmosphere scale height (in km)
:param altitudes_D: Altitudes (in km) in the D-region
:param bin_0: Integrated HXR flux between 0.05 and 0.1 nm
:param tau_HXR_0: Tau_HXR for bin_0
:param bin_1: Integrated HXR flux between 0.1 and 0.15 nm
:param tau_HXR_1: Tau_HXR for bin_1
:param bin_2: Integrated HXR flux between 0.15 and 0.2 nm
:param tau_HXR_2: Tau_HXR for bin_2
:param bin_3: Integrated HXR flux between 0.2 and 0.25 nm
:param tau_HXR_3: Tau_HXR for bin_3
:param bin_4: Integrated HXR flux between 0.25 and 0.3 nm
:param tau_HXR_4: Tau_HXR for bin_4
:param bin_5: Integrated HXR flux between 0.3 and 0.35 nm
:param tau_HXR_5: Tau_HXR for bin_5
:param bin_6: Integrated HXR flux between 0.35 and 0.4 nm
:param tau_HXR_6: Tau_HXR for bin_6
:param Phi_HXR_bins: Phi_HXR in the shape (bins * altitudes). It is computed in 'cs.compute_photon_flux' and called in 'chemistry_mr_eq'
"""
def __init__(self, today, file_EUV, file_XR, altitudes_D, neutrals_here, chi):
"""
Relies on past measurements from GOES data
:param today: Datetime, time of interest. Must be timezone unaware
:param file_EUV: Path to the EUV file (daily EUV average)
:param file_XR: Path to the XR data (1s measurements)
:param altitudes_D: Altitudes (in km) in the D-region
:param neutrals_here: Neutrals class instance
:param chi: Solar zenith angle (in °)"""
# Load the files
X_rays = netCDF4.Dataset(file_XR, "r")
EUV = netCDF4.Dataset(file_EUV, "r")
# Get X-ray flux
time_XR = np.arange(0, 24, 1 / 3600)
id = np.argmin(
np.abs(today.hour + today.minute / 60 + today.second / 3600 - time_XR)
)
self.SXR = X_rays.variables["xrsb_flux"][id]
self.HXR = X_rays.variables["xrsa_flux"][id]
# Bin HXR flux
self.bin_HXR_flux()
# Get the daily average for EUV flux
# EUV GOES files times are given in s since 2000_01_01 12:00 UTC
ref_EUV_date = dt.datetime(2000, 1, 1, 12, 0, 0, 0) # Defined in the EUV file
nb_seconds = (today - ref_EUV_date).total_seconds()
EUV_times = EUV.variables["time"][:]
# Might not be the most elegant method
id = np.argmin(np.abs(nb_seconds - EUV_times))
self.EUV = EUV.variables["irr_1216"][id] # W/m2
# Compute absorbed fluxes
self.altitudes_D = altitudes_D
tau_EUV, tau_SXR, tau_HXR, tau_HXR_bins, Ch, H = cs.compute_tau(
self, neutrals_here, chi
)
self.tau_EUV = tau_EUV
self.tau_SXR = tau_SXR
self.tau_HXR = tau_HXR
self.tau_HXR_0 = tau_HXR_bins[0, :]
self.tau_HXR_1 = tau_HXR_bins[1, :]
self.tau_HXR_2 = tau_HXR_bins[2, :]
self.tau_HXR_3 = tau_HXR_bins[3, :]
self.tau_HXR_4 = tau_HXR_bins[4, :]
self.tau_HXR_5 = tau_HXR_bins[5, :]
self.tau_HXR_6 = tau_HXR_bins[6, :]
# Restric Ch and H to altitudes in radiation class
id_min = np.argmin(np.abs(self.altitudes_D[0] - neutrals_here.altitudes))
id_max = np.argmin(np.abs(self.altitudes_D[-1] - neutrals_here.altitudes))
Ch = Ch[id_min : id_max + 1]
H = H[id_min : id_max + 1]
# Keep track of H and Ch
self.Ch = Ch
self.H = H
# Keep in mind the file contents to avoid opening them again
self.EUV_array = EUV.variables["irr_1216"][:]
self.EUV_times = EUV.variables["time"][:]
self.XR_times = time_XR
self.SXR_array = X_rays.variables["xrsb_flux"][:]
self.HXR_array = X_rays.variables["xrsa_flux"][:]
[docs]
def get_flux_now(self, time_s, today):
"""Updates the XR flux (the EUV is a daily average)
:param time_s: Time since the start (in s)
:param today: Datetime of the start of the computation"""
id = np.argmin(
np.abs(
today.hour
+ today.minute / 60
+ (today.second + time_s) / 3600
- self.XR_times
)
)
self.SXR = self.SXR_array[id]
self.HXR = self.HXR_array[id]
# Update HXR bins
self.bin_HXR_flux()
[docs]
def update_chi(self, today, n_here):
"""Updates the solar zenith angle and recomputes taus, Ch and H
:param today: Datetime, time of interest. Must be timezone unaware
:param glat/glon: geographic latitude and longitude
:param n_here: Neutrals class instance
"""
# Make time_here aware of the timezone
time_here_aware = pytz.utc.localize(today)
# Compute solar zenith angle
chi = float(90) - get_altitude(n_here.glat, n_here.glon, time_here_aware)
# Recompute taus
tau_EUV, tau_SXR, tau_HXR, tau_HXR_bins, Ch, H = cs.compute_tau(
self, n_here, chi
)
# Keep track of taus
self.tau_EUV = tau_EUV
self.tau_SXR = tau_SXR
self.tau_HXR = tau_HXR
self.tau_HXR_0 = tau_HXR_bins[0, :]
self.tau_HXR_1 = tau_HXR_bins[1, :]
self.tau_HXR_2 = tau_HXR_bins[2, :]
self.tau_HXR_3 = tau_HXR_bins[3, :]
self.tau_HXR_4 = tau_HXR_bins[4, :]
self.tau_HXR_5 = tau_HXR_bins[5, :]
self.tau_HXR_6 = tau_HXR_bins[6, :]
# Restric Ch and H to altitudes in radiation class
id_min = np.argmin(np.abs(self.altitudes_D[0] - n_here.altitudes))
id_max = np.argmin(np.abs(self.altitudes_D[-1] - n_here.altitudes))
Ch = Ch[id_min : id_max + 1]
H = H[id_min : id_max + 1]
# Keep track of H and Ch
self.Ch = Ch
self.H = H
[docs]
def bin_HXR_flux(self):
"""Bins GOES HXR flux (0.05 - 0.4 nm) into bins with 0.05 nm length, following Siskind et al., 2022
Does not return anything, but puts the integrated flux in each bin as attributes of the radiation class.
For example, bin_0 contains the integrated flux between 0.05 and 0.1 nm.
It is considered that the spectrum is linear between 0.05 and 0.15 nm, and increases by a factor of 100.
It similarly increases by a factor 5 between 0.15 and 0.4 nm (See Fig. 1, Siskind et al., 2022)
If M is the flux at 0.05 nm, the total integral between 0.05 nm and 0.4 nm is given by:
A = M (0.5*0.1nm*99 + 0.5*0.25nm*400 + 0.25nm*99) = M * B"""
# Computation of B
B = (
0.5 * 0.1 * 99 + 0.5 * 0.25 * 400 + 0.25 * 99
) # We do not really care about nm here, [M] = Flux/nm
# Computation of M
M = self.HXR / B
# Bins 0 and 1 (fast linear increase, ax + b)
a = 99 * M / (0.15 - 0.05)
b = M - 0.05 / (0.15 - 0.05) * 99 * M
self.bin_0 = (0.5 * a * (0.1) ** 2 + b * 0.1) - (
0.5 * a * (0.05) ** 2 + b * 0.05
)
self.bin_1 = (0.5 * a * (0.15) ** 2 + b * 0.15) - (
0.5 * a * (0.1) ** 2 + b * 0.1
)
# Bins 2 to 6 (slower linear increase, ax + b)
a = (500 - 100) * M / (0.4 - 0.15)
b = 100 * M - 0.15 / (0.4 - 0.15) * 400 * M
self.bin_2 = (0.5 * a * (0.2) ** 2 + b * 0.2) - (
0.5 * a * (0.15) ** 2 + b * 0.15
)
self.bin_3 = (0.5 * a * (0.25) ** 2 + b * 0.25) - (
0.5 * a * (0.2) ** 2 + b * 0.2
)
self.bin_4 = (0.5 * a * (0.3) ** 2 + b * 0.3) - (
0.5 * a * (0.25) ** 2 + b * 0.25
)
self.bin_5 = (0.5 * a * (0.35) ** 2 + b * 0.35) - (
0.5 * a * (0.3) ** 2 + b * 0.3
)
self.bin_6 = (0.5 * a * (0.4) ** 2 + b * 0.4) - (
0.5 * a * (0.35) ** 2 + b * 0.35
)