Reference

AbstractGrid

Abstract supertype for all AURORA grids (altitude, energy, pitch-angle).

AltitudeGrid{FT, V<:AbstractVector{FT}} <: AbstractGrid

Altitude grid for the ionosphere simulation.

Fields

• h: altitude values (m)
• Δh: grid spacing (m)
• n: number of grid points
• bottom: bottom altitude (km)
• top: top altitude (km)

Constructor

AltitudeGrid(bottom, top; dz_max=25)

Create an altitude grid from bottom to top (in km), with maximum grid spacing dz_max (km).

AuroraModel{AG, EG, PAG, SD, IO, CS, FT, V}

Container for the grids, atmosphere, and collision data used by an AURORA simulation.

AuroraModel(altitude_lims, θ_lims, E_max, msis_file, iri_file, B_angle_to_zenith=0;
            verbose=true)

Load the atmosphere, the energy grid, and the collision data into an AuroraModel.

Calling

model = AuroraModel(altitude_lims, θ_lims, E_max, msis_file, iri_file, B_angle_to_zenith)

Inputs

• altitude_lims: the altitude limits, in km, for the bottom and top of the ionosphere in our simulation
• θ_lims: pitch-angle limits of the electron beams (e.g. 180:-10:0), where 180° corresponds to field aligned down, and 0° field aligned up. Vector [n_beam]
• E_max: upper limit for the energy grid (in eV)
• msis_file: path to the msis file to use
• iri_file: path to the iri file to use
• B_angle_to_zenith: angle between magnetic field and zenith (degrees)
• verbose=true: print progress messages when loading stuff

Returns

An AuroraModel with fields:

• altitude_grid::AltitudeGrid
• energy_grid::EnergyGrid
• pitch_angle_grid::PitchAngleGrid
• scattering::ScatteringData
• ionosphere::Ionosphere
• cross_sections::CrossSectionData
• B_angle_to_zenith::Real
• s_field::Vector

CrossSectionData{NT1<:NamedTuple, NT2<:NamedTuple}

Collision cross-sections and excitation thresholds for all neutral species.

EnergyGrid{FT, V<:AbstractVector{FT}} <: AbstractGrid

Energy grid for electron transport.

Fields

• E_edges: energy bin edges (eV), length n + 1
• E_centers: energy bin centers (eV), length n
• ΔE: energy bin widths (eV), length n
• n: number of energy bins
• E_max: requested maximum energy (eV)

Ionosphere{FT, V<:AbstractVector{FT}}

Atmospheric state containing neutral and electron densities and temperatures.

PitchAngleGrid{FT, V<:AbstractVector{FT}} <: AbstractGrid

Pitch-angle grid for electron beams.

ScatteringData{FT, A<:AbstractArray{FT}, M<:AbstractMatrix{FT}, V<:AbstractVector{FT}}

Pre-computed scattering probability matrices for beam-to-beam scattering.

TransportMatrices

Container for the matrices used in the electron transport equations.

TransportMatrices(n_altitude, n_angle, n_time, n_energy)

Construct an empty TransportMatrices container with zeros.

Crank_Nicolson_optimized!(Ie, t, model, v, matrices, iE, Ie_top, I0, cache; first_iteration = false)

Optimized Crank-Nicolson time-stepping scheme using direct nzval modification. This is an in-place version that modifies Ie directly to avoid allocations.

On first iteration, creates the sparsity pattern and mapping which are stored in cache. On subsequent iterations, only updates the nzval array directly.

Mathematical Background

The time-dependent electron transport equation is:

∂Ie/∂t + μ*v ∂Ie/∂z + A*Ie - ∫B*Ie'dΩ' - D*∇²Ie = Q

The Crank-Nicolson scheme uses implicit time-stepping with second-order accuracy:

(Ie^(n+1) - Ie^n)/Δt = [RHS^(n+1) + RHS^n] / 2

This leads to two matrices:

Mlhs * Ie^(n+1) = Mrhs * Ie^n + (Q^(n+1) + Q^n)/2
     ↑                  ↑
Implicit part      Explicit part

Where:

Mlhs = Ddt + μ*Ddz/2 + A/2 - B/2 - D*Ddiffusion/2
Mrhs = Ddt - μ*Ddz/2 - A/2 + B/2 + D*Ddiffusion/2

Both matrices have the same block structure as in steady-state:

┌─────────┬─────────┬─────────┐
│ Block   │ Block   │ Block   │  Each block is n_z × n_z
│ (1,1)   │ (1,2)   │ (1,3)   │
├─────────┼─────────┼─────────┤
│ Block   │ Block   │ Block   │  Off-diagonal: angular scattering
│ (2,1)   │ (2,2)   │ (2,3)   │
├─────────┼─────────┼─────────┤
│ Block   │ Block   │ Block   │  Diagonal: transport + diffusion
│ (3,1)   │ (3,2)   │ (3,3)   │
└─────────┴─────────┴─────────┘

Arguments

• Ie: pre-allocated output array [m⁻² s⁻¹], size (nz * nangle × n_t) to store results
• t: time grid [s]
• model: AuroraModel (s_field and pitch_angle_grid.μ_center are used)
• v: electron velocity [km/s]
• matrices::TransportMatrices: container with
  • A: electron loss rate [s⁻¹]
  • B: scattering matrix [s⁻¹], size (nz × nangle × n_angle)
  • D: pitch-angle diffusion coefficient [s⁻¹], size (n_angle,)
  • Q: source term [m⁻² s⁻¹] at each time step
  • Ddiffusion: spatial diffusion matrix (nz × nz)
• iE: current energy index
• Ie_top: boundary condition at top [m⁻² s⁻¹] at each time step
• I0: initial condition [m⁻² s⁻¹]
• cache: Cache object (must have fields for Mlhs, Mrhs, mappings, KLU, diff matrices)
• first_iteration: whether this is the first call

Ie_top_from_file(t, E, μ_center, filename; interpolation=:constant)

Load a time-dependent electron flux from a .mat file and resample it onto the simulation time grid t.

Arguments

• t: full simulation time grid (s). Range or Vector [n_t]
• E: energy grid (eV). Vector [n_E]
• μ_center: electron beams average pitch angle cosine. Vector [n_μ]
• filename: path to the .mat file containing the flux data

Keyword Arguments

• interpolation=: interpolation scheme used to resample the file's time grid onto the simulation time grid. Either:
  • :constant (default): each file value is held constant until the next sample.
  • :linear: linear interpolation between consecutive file samples.
  • :pchip: piecewise cubic Hermite interpolation.

Returns

• Ie_top: electron number flux (#e⁻/m²/s). Array [nμ, nt, n_E]

File format

The .mat file must contain:

• Ie_total: flux array of shape [n_μ, n_t_file, n_E]
• t_top: time grid (s) of the file, as a vector [n_t_file]. Optional when n_t_file == 1.

Time grid handling

The file's time grid (t_top) does not need to match the simulation time grid. It is treated as relative: t_top[1] is always aligned with t[1], regardless of the absolute values stored in the file. Only the duration and spacing of t_top matter. Any combination of time step sizes and grid lengths is supported:

• Different dt: the file is resampled onto the simulation grid via interpolation.
• File finer than simulation: the flux is point-sampled, which may miss rapid variations.
• File shorter than simulation: flux is set to zero for times beyond the file duration.
• File longer than simulation: the file is simply evaluated up to t[end].

Notes

• Upward-going beams (μ_center > 0) are always set to zero.
• The energy dimension of the file may be larger than length(E); only the first length(E) bins are used. It cannot be smaller.

Ie_top_modulated(IeE_tot, model, Beams, t, n_loop;
                 spectrum=:flat, E_min=0.0, E₀=nothing, ΔE=nothing,
                 modulation=:none, f=0.0, amplitude=1.0,
                 z_source=model.altitude_grid.h[end]/1e3, t_start=0.0, t_end=0.0)

Create a time-dependent electron flux distribution with configurable energy spectrum shape and temporal modulation.

Returns an electron flux distribution (in #e⁻/m²/s) such that when integrated over energy (at full modulation), the total energy flux IeE_tot is recovered.

Arguments

• IeE_tot: total energy flux (W/m²)
• model: AuroraModel describing the grids and atmosphere
• Beams: indices of the electron beams with precipitating flux
• t: time grid (s). Range or Vector [n_t]
• n_loop: number of loops (for repeated simulations)

Keyword Arguments

Energy spectrum shape

• spectrum=:flat: energy spectrum type. Either :flat or :gaussian
• E_min=0.0: minimum energy threshold (eV) - only used for spectrum=:flat
• E₀=nothing: characteristic/center energy (eV) - required for spectrum=:gaussian
• ΔE=nothing: energy width (eV) - required for spectrum=:gaussian

Temporal modulation

• modulation=:none: temporal modulation type. One of :none, :sinus, or :square
• f=0.0: modulation frequency (Hz) - used for :sinus and :square
• amplitude=1.0: modulation depth. 0 = constant flux, 1 = full on/off modulation. Values between 0 and 1 give partial modulation.

Time-dependent features

• z_source: altitude of the precipitation source (km). Default: top of ionosphere
• t_start=0.0: start time for smooth flux onset (s) - only used for modulation=:none
• t_end=0.0: end time for smooth flux onset (s) - only used for modulation=:none

Returns

• Ie_top: electron number flux (#e⁻/m²/s). Matrix [nbeams, nt, n_E]

Physics

The function creates an electron precipitation spectrum with:

1. Energy spectrum: Either flat (uniform in #e⁻/m²/s/eV above Emin) or Gaussian (centered at E₀ with width ΔE). Both are normalized so the total energy flux equals `IeEtot`.

2. Energy- and angle-dependent arrival times: Lower energy electrons travel slower, arriving later at the ionosphere top. Similarly, electrons with high pitch-angles travel longer paths. This creates natural dispersion.

3. Temporal modulation: The flux can be constant (:none), sinusoidally modulated (:sinus), or square-wave modulated (:square). The amplitude parameter controls the modulation depth.

Examples

Flat spectrum with smooth onset:

julia> msis_file = find_msis_file(verbose=false);

julia> iri_file = find_iri_file(verbose=false);

julia> model = AuroraModel((100, 600), 180:-10:0, 10e3, msis_file, iri_file; verbose=false);

julia> Ie = Ie_top_modulated(1e-2, model, 1:2, 0:0.01:1, 1;
                             spectrum=:flat, E_min=9000.0, t_start=0.0, t_end=0.1);

Gaussian spectrum with sinusoidal modulation at 10 Hz:

julia> msis_file = find_msis_file(verbose=false);

julia> iri_file = find_iri_file(verbose=false);

julia> model = AuroraModel((100, 600), 180:-10:0, 10e3, msis_file, iri_file; verbose=false);

julia> Ie = Ie_top_modulated(1e-2, model, 1:2, 0:0.001:0.5, 1;
                             spectrum=:gaussian, E₀=5000.0, ΔE=500.0,
                             modulation=:sinus, f=10.0, amplitude=1.0);

Ie_top_modulated(IeE_tot, model, Beams;
                 spectrum=:flat, E_min=0.0, E₀=nothing, ΔE=nothing)

Steady-state version of Ie_top_modulated that does not include time-dependent behavior.

This overload is designed for steady-state simulations and eliminates the need to specify time-related parameters (t, n_loop, modulation, f, amplitude, z_source, etc.). It internally calls the time-dependent version with minimal time grid (1:1:1) and n_loop=1.

Arguments

• IeE_tot: total energy flux (W/m²)
• model: AuroraModel describing the grids and atmosphere
• Beams: indices of the electron beams with precipitating flux

Keyword Arguments

• spectrum=:flat: energy spectrum type. Either :flat or :gaussian
• E_min=0.0: minimum energy threshold (eV) - only used for spectrum=:flat
• E₀=nothing: characteristic/center energy (eV) - required for spectrum=:gaussian
• ΔE=nothing: energy width (eV) - required for spectrum=:gaussian

Returns

• Ie_top: electron number flux (#e⁻/m²/s). Matrix [nbeams, 1, nE]

Examples

julia> msis_file = find_msis_file(verbose=false);

julia> iri_file = find_iri_file(verbose=false);

julia> model = AuroraModel((100, 600), 180:-10:0, 10e3, msis_file, iri_file; verbose=false);

julia> Ie = Ie_top_modulated(1e-2, model, 1:2; spectrum=:flat, E_min=9000.0);

Ie_with_LET(IeE_tot, E₀, model, Beams; low_energy_tail=true)

Return an electron spectra following a Maxwellian distribution with a low energy tail (LET)

This function is a corrected implementation of Meier/Strickland/Hecht/Christensen JGR 1989 (pages 13541-13552)

Arguments

• IeE_tot: total energy flux (W/m²)
• E₀: characteristic energy (eV)
• model: AuroraModel describing the grids and atmosphere
• Beams: indices of the electron beams with a precipitating flux

Keyword Arguments

• low_energy_tail=true: control the presence of a low energy tail

Returns:

• Ie_top: electron number flux (#e⁻/m²/s). Matrix [n_beams, 1, nE]

Important notes

This is a corrected version of the equations present in Meier et al. 1989 to match the results presented in Fig. 4 of their paper.
Changes were made to the factor b:

• no inverse

Examples:

Calling the function with flux only in the two first beams (0 to 20°) and an "isotropic" pitch-angle distribution.

julia> msis_file = find_msis_file(verbose=false);

julia> iri_file = find_iri_file(verbose=false);

julia> model = AuroraModel((100, 600), 180:-10:0, 10e3, msis_file, iri_file; verbose=false);

julia> Ie = AURORA.Ie_with_LET(1e-2, 1e3, model, 1:2);

Calling the function with flux only in the three first beams (0 to 30°) and a custom pitch-angle distribution (1/2 of the total flux in the first beam, 1/4 in the second beam and 1/4 in the third beam).

julia> msis_file = find_msis_file(verbose=false);

julia> iri_file = find_iri_file(verbose=false);

julia> model = AuroraModel((100, 600), 180:-10:0, 10e3, msis_file, iri_file; verbose=false);

julia> Ie = Ie_with_LET(1e-2, 1e3, model, 1:3);

_apply_modulation(t_shifted, modulation, f, amplitude, t_start, t_end)

Apply temporal modulation to the flux based on the shifted time grid.

Arguments

• t_shifted: time grid shifted for energy/angle-dependent delays
• modulation: :none, :sinus, or :square
• f: frequency for sinus/square modulation (Hz)
• amplitude: modulation depth (0 = constant, 1 = full on/off)
• t_start, t_end: smooth onset interval (only used for :none)

Returns

• Vector of modulation factors (0 to 1) for each time step

_flat_spectrum(IeE_tot_eV, E, dE, E_min)

Compute a flat (constant) differential number flux spectrum above Emin, normalized so that the total energy flux equals IeEtot_eV.

Returns a vector of differential number flux per eV (#e⁻/m²/s/eV) for each energy bin. The flux is zero below E_min.

_gaussian_spectrum(IeE_tot_eV, E, dE, E₀, ΔE)

Compute a Gaussian differential number flux spectrum centered at E₀ with width ΔE, normalized so that the total energy flux equals IeEtoteV.

The Gaussian shape is: Φ(E) ∝ exp(-(E - E₀)² / ΔE²)

Returns a vector of differential number flux per eV (#e⁻/m²/s/eV) for each energy bin.

_smooth_transition(x, x_start = 0.0, x_end = 1.0)

Create a smooth transition from 0 to 1 over the interval [x_start, x_end] using a C∞-smooth (infinitely differentiable) function.

• Returns 0 for x < x_start
• Returns 1 for x > x_end
• Smoothly transitions from 0 to 1 in between, with all derivatives continuous

Arguments

• x: The input value at which to evaluate the transition
• x_start: Start of the transition interval (default: 0.0)
• x_end: End of the transition interval (default: 1.0)

Returns

• A value between 0 and 1 representing the smooth transition

animate_Ie_in_time(directory_to_process; angles_to_plot=nothing, colorrange=nothing, ...)

Plot a heatmap of Ie over height and energy, and animate it in time. It will load the result files one by one. The animation will be saved as a .mp4 file under the directory_to_process.

Example

julia> directory_to_process = "Visions2/Alfven_475s";

# Using defaults for angles and colorrange:
julia> animate_Ie_in_time(directory_to_process)

# Or with custom angles and colorrange:
julia> angles_to_plot = [(180, 170)  (170, 150)  (150, 120)  (120, 100)  (100, 90);   # DOWN
                         (0, 10)     (10, 30)    (30, 60)    (60, 80)    (80, 90)];   # UP
julia> animate_Ie_in_time(directory_to_process; angles_to_plot, colorrange=(1e5, 1e9), plot_Ietop=true)

# Using nothing for empty panels:
julia> angles_to_plot = [(180, 90)  nothing;
                         (0, 45)    (45, 90)];
julia> animate_Ie_in_time(directory_to_process; angles_to_plot)

The angles_to_plot is a matrix of tuples, where each tuple defines a pitch-angle range from 0° to 180° (where 180° is field-aligned down and 0° is field-aligned up). A panel will be created for each matrix element at the corresponding row/column position. Angles > 90° are labeled as "DOWN", angles < 90° as "UP". Use nothing for empty panels.

The limits of angles_to_plot need to match existing limits of the beams used in the simulation. E.g. if θ_lims = 180:-10:0 was used in the simulation, (150, 120) will be fine as 150° and 120° exist as limits, but (155, 120) will not as 155° does not exist as a limit.

Arguments

• directory_to_process: directory containing the simulation results (absolute or relative path).

Keyword Arguments

• angles_to_plot = nothing: limits of the angles to plot as a matrix of tuples with angles in range 0-180°. Use nothing for empty panels. If the whole argument is nothing, uses the simulation grid with down-flux on the first row and up-flux on the second row.
• colorrange = nothing: limits for the colormap/colorbar as a tuple (min, max). If nothing, automatically computed as (maxvalue / 1e4, maxvalue) spanning 4 orders of magnitude.
• save_to_file = true: if true, saves the animation to a .mp4 file in the data directory.
• plot_Ietop = false: if true, also plots the precipitating Ie at the top of the ionosphere by loading it from the file Ie_top.mat.
• Ietop_angle_cone = [170, 180]: angle cone (in degrees) for the precipitating Ie plot.
• dt_steps: plot one frame every dt_steps timesteps.

beam_weight(θ_lims)

Return the beam weights of the pitch-angle beams delimited by θ_lims.

Calling

BeamW = beam_weight(θ_lims)

Inputs

• θ_lims : pitch-angle limits of all the beams, range or vector [n_beams + 1]

Outputs

• BeamW : solid angle of each pitch-angle beams (ster), vector [n_beams]

calculate_cascading_N2(E_grid, dE, lorentzian_width)

Calculate cascading transfer matrices for N₂ ionization.

Arguments

• E_grid: Energy grid for electrons (eV)
• dE: Energy grid step size (eV)
• lorentzian_width: Width of the Lorentzian distribution (eV)

Returns

• Q_transfer_matrix: Transfer matrix [nenergies, nenergies, n_thresholds]

Reference

Equation from Itikawa 1986 J. Phys. Chem. Ref. Data

calculate_cascading_N2_quadgk(E_grid, dE, lorentzian_width)

Calculate cascading transfer matrices for N₂ ionization using nested quadgk integrals.

This is an alternative implementation to calculate_cascading_N2 that uses two nested one-dimensional integrals instead of a hypercube transformation. The integration is performed over the physical domain directly:

• Outer integral: over degraded electron energy E_degraded
• Inner integral: over primary electron energy E_primary

Arguments

• E_grid: Energy grid for electrons (eV)
• dE: Energy grid step size (eV)
• lorentzian_width: Width of the Lorentzian distribution (eV)

Returns

• Q_transfer_matrix: Transfer matrix [nenergies, nenergies, n_thresholds]
• E_grid: Energy grid (returned for consistency)
• ionization_thresholds: Array of ionization threshold energies

Reference

Equation from Itikawa 1986 J. Phys. Chem. Ref. Data

calculate_cascading_O(E_grid, dE)

Calculate cascading transfer matrices for atomic O ionization.

Arguments

• E_grid: Energy grid for electrons (eV)
• dE: Energy grid step size (eV)

Returns

• Q_transfer_matrix: Transfer matrix [nenergies, nenergies, n_thresholds]

calculate_cascading_O2(E_grid, dE, lorentzian_width)

Calculate cascading transfer matrices for O₂ ionization.

Arguments

• E_grid: Energy grid for electrons (eV)
• dE: Energy grid step size (eV)
• lorentzian_width: Width of the Lorentzian distribution (eV)

Returns

• Q_transfer_matrix: Transfer matrix [nenergies, nenergies, n_thresholds]

calculate_cascading_O2_quadgk(E_grid, dE, lorentzian_width)

Calculate cascading transfer matrices for O₂ ionization using nested quadgk integrals.

This is an alternative implementation to calculate_cascading_O2 that uses two nested one-dimensional integrals instead of a hypercube transformation. The integration is performed over the physical domain directly:

• Outer integral: over degraded electron energy E_degraded
• Inner integral: over primary electron energy E_primary

Arguments

• E_grid: Energy grid for electrons (eV)
• dE: Energy grid step size (eV)
• lorentzian_width: Width of the Lorentzian distribution (eV)

Returns

• Q_transfer_matrix: Transfer matrix [nenergies, nenergies, n_thresholds]
• E_grid: Energy grid (returned for consistency)
• ionization_thresholds: Array of ionization threshold energies

calculate_cascading_O_quadgk(E_grid, dE)

Calculate cascading transfer matrices for atomic O ionization using nested quadgk integrals.

This is an alternative implementation to calculate_cascading_O that uses two nested one-dimensional integrals instead of a hypercube transformation. The integration is performed over the physical domain directly:

• Outer integral: over degraded electron energy E_degraded
• Inner integral: over primary electron energy E_primary

Arguments

• E_grid: Energy grid for electrons (eV)
• dE: Energy grid step size (eV)

Returns

• Q_transfer_matrix: Transfer matrix [nenergies, nenergies, n_thresholds]
• E_grid: Energy grid (returned for consistency)
• ionization_thresholds: Array of ionization threshold energies

calculate_density_from_Ie!(z, t_run, μ_lims, E_centers, v, Ie, n_e)

This function converts a particle flux Ie (#e⁻/m²/s) into a number density n_e (#e⁻/m³).

The particle flux Ie is defined along a magnetic field line and over an (Energy, pitch_angle)-grid. The number density n_e calculated is given along a magnetic field line and over an energy grid. That way, we have the density of electrons with a certain energy at a specific altitude and time.

Calling

calculate_density_from_Ie!(z, t_run, μ_lims, E_centers, v, Ie, n_e)

Inputs

• z: altitude (m), vector [nz]
• t_run: time (s), vector [nt]
• μ_lims: cosine of the pitch angle limits of the e- beams, vector [n_beam + 1]
• E_centers: center energy of the energy bins (eV), vector [nE]
• v: velocity corresponding to the E_centers (m/s), vector [nE]
• Ie: electron flux (#e⁻/m²/s), 3D array [n_beam * nz, nt, nE]
• n_e: electron density (#e⁻/m³), empty 3D array [nz, nt, nE]

calculate_e_transport(altitude_lims, θ_lims, E_max, B_angle_to_zenith, t_total, dt,
    msis_file, iri_file, savedir, INPUT_OPTIONS, CFL_number = 64;
    n_loop = nothing, max_memory_gb = 8, save_input_flux = true)

Run a time-dependent electron transport simulation and save the results to savedir.

Arguments

• altitude_lims: altitude range of the simulation (km). Tuple or Vector
• θ_lims: pitch-angle beam limits (degrees). Range or Vector [n_beams + 1]
• E_max: maximum energy of the simulation (eV)
• B_angle_to_zenith: angle between the magnetic field and the zenith (degrees)
• t_total: total simulation time (s)
• dt: output time step (s)
• msis_file: path to the MSIS neutral atmosphere file
• iri_file: path to the IRI ionosphere file
• savedir: path to the directory where results will be saved
• INPUT_OPTIONS: NamedTuple describing the incoming electron flux source
• CFL_number: CFL multiplier used to sub-sample the time grid (default: 64)

Keyword Arguments

• n_loop: number of loops to split the time grid into. Computed automatically from max_memory_gb if not provided
• max_memory_gb: memory budget used to compute n_loop automatically
• save_input_flux: if true, save the computed top-boundary flux Ie_top to Ie_incoming.mat in savedir before the simulation starts

calculate_e_transport_steady_state(altitude_lims, θ_lims, E_max, B_angle_to_zenith,
    msis_file, iri_file, savedir, INPUT_OPTIONS; save_input_flux = true)

Run a steady-state electron transport simulation and save the results to savedir.

Arguments

• altitude_lims: altitude range of the simulation (km). Tuple or Vector
• θ_lims: pitch-angle beam limits (degrees). Range or Vector [n_beams + 1]
• E_max: maximum energy of the simulation (eV)
• B_angle_to_zenith: angle between the magnetic field and the zenith (degrees)
• msis_file: path to the MSIS neutral atmosphere file
• iri_file: path to the IRI ionosphere file
• savedir: path to the directory where results will be saved
• INPUT_OPTIONS: NamedTuple describing the incoming electron flux source

Keyword Arguments

• save_input_flux: if true, save the computed top-boundary flux Ie_top to Ie_incoming.mat in savedir before the simulation starts

calculate_heating_rate(z, t, Ie_ztE_omni, E_centers, ne, Te)

Calculate the heating rate of thermal electrons by superthermal electrons through Coulomb collisions. The heating rate is the rate at which energy is transferred from superthermal electrons to thermal electrons.

Calling

heating_rate = calculate_heating_rate(z, t, Ie_ztE_omni, E_centers, ne, Te)

Inputs

• z: altitude (m). Vector [n_z]
• t: time (s). Vector [n_t]
• Ie_ztE_omni: omnidirectional electron flux (#e⁻/m²/s). 3D array [n_z, n_t, n_E]
• E_centers: energy bin centers (eV). Vector [n_E]
• ne: thermal electron density (m⁻³). Vector [n_z]
• Te: thermal electron temperature (K). Vector [n_z]

Output

• heating_rate: heating rate (eV/m³/s). 2D array [n_z, n_t]

calculate_iri_data(; year=2018, month=12, day=7, hour=11, minute=15,
                    lat=76, lon=5, height=85:1:700)

Calculate IRI-2020 ionospheric model data using Python interface.

This function calls the Python iri2020 package to compute ionospheric parameters including electron density, temperatures, and ion composition profiles. The data is returned as a matrix with a header row containing column names.

Keyword Arguments

• year::Int=2018: Year (defaults to Visions2 launch conditions)
• month::Int=12: Month (1-12)
• day::Int=7: Day of month (1-31)
• hour::Int=11: Hour in Universal Time (0-23)
• minute::Int=15: Minute (0-59)
• lat::Real=76: Geographic latitude in degrees North
• lon::Real=5: Geographic longitude in degrees East
• height::AbstractRange=85:1:700: Altitude range in km

Returns

• Tuple{Matrix, NamedTuple}:
  • Matrix with header row and IRI data (21 columns: height + 20 parameters)
  • NamedTuple with input parameters

Data Columns

The returned matrix contains the following columns:

1. height(km): Altitude
2. ne(m⁻³): Electron density
3. Tn(K): Neutral temperature
4. Ti(K): Ion temperature
5. Te(K): Electron temperature

6-12. nO⁺, nH⁺, nHe⁺, nO2⁺, nNO⁺, nCI, nN⁺: Ion densities (m⁻³) 13-18. NmF2, hmF2, NmF1, hmF1, NmE, hmE: Peak densities and heights

1. TEC: Total Electron Content
2. EqVertIonDrift: Equatorial vertical ion drift
3. foF2: F2 critical frequency

calculate_msis_data(; year=2018, month=12, day=7, hour=11, minute=15,
                     lat=76, lon=5, height=85:1:700)

Calculate NRLMSIS 2.1 atmospheric model data using Python interface.

This function calls the Python pymsis package to compute neutral atmosphere parameters including densities of various species (N₂, O₂, O, He, H, Ar, N, NO) and temperature profiles. The data is returned as a matrix with a header row containing column names.

Keyword Arguments

• year::Int=2018: Year (defaults to Visions2 launch conditions)
• month::Int=12: Month (1-12)
• day::Int=7: Day of month (1-31)
• hour::Int=11: Hour in Universal Time (0-23)
• minute::Int=15: Minute (0-59)
• lat::Real=76: Geographic latitude in degrees North
• lon::Real=5: Geographic longitude in degrees East
• height::AbstractRange=85:1:700: Altitude range in km

Returns

• Tuple{Matrix, NamedTuple}:
  • Matrix with header row and MSIS data (12 columns: height + 11 parameters)
  • NamedTuple with input parameters

Data Columns

The returned matrix contains the following columns:

1. height(km): Altitude
2. air(kg/m³): Total mass density
3. N₂(m⁻³): Molecular nitrogen density
4. O₂(m⁻³): Molecular oxygen density
5. O(m⁻³): Atomic oxygen density
6. He(m⁻³): Helium density
7. H(m⁻³): Atomic hydrogen density
8. Ar(m⁻³): Argon density
9. N(m⁻³): Atomic nitrogen density
10. anomalousO(m⁻³): Anomalous oxygen density
11. NO(m⁻³): Nitric oxide density
12. T(K): Temperature

Notes

• This function automatically downloads the SW-All.csv file (space weather data) if not already present in the pymsis package directory
• The geomagnetic activity parameter is set to -1 (use observational data)

calculate_n_loop(t, n_z, n_μ, n_E; max_memory_gb=8, verbose=true)

Calculate the optimal number of loops (n_loop) for the electron transport simulation based on memory constraints.

This function determines how many time-slices the simulation should be divided into to stay within the specified memory limit.

Arguments

• t: Time array after CFL refinement (from CFL_criteria)
• n_z::Int: Number of altitude grid points
• n_μ::Int: Number of pitch-angle beams
• n_E::Int: Number of energy grid points

Keyword Arguments

• max_memory_gb: Maximum memory to use (GB). Defaults to 8 GB.
• verbose::Bool=true: If true, print some information about the calculation

Returns

• n_loop::Int: Recommended number of loops

Example

t, CFL_factor = CFL_criteria(1.0, 0.001, h_atm, v_of_E(10000))
n_loop = calculate_n_loop(t, length(h_atm), n_μ, n_E)

calculate_scattered_flux!(result, B2B_inelastic, n, Ie_slice)

Calculate the flux of electrons after pitch-angle scattering by inelastic collisions.

Physics

The scattered flux at each altitude and angle is computed as: result[z, μ₁, t] = Σ_μ₂ n(z) x P(μ₁←μ₂) x Ie[z, μ₂, t]

where:

• n(z) is the neutral density at altitude z
• P(μ₁←μ₂) is the probability of scattering from pitch angle μ₂ to μ₁ (from B2B_inelastic)
• Ie[z, μ₂, t] is the incident electron flux before scattering

Implementation

This exploits the block-diagonal structure of the full scattering matrix to avoid storing and accessing a large sparse matrix. Instead, it computes the multiplication directly from the small scattering probability matrix and density profile.

Arguments

• result: Output array (nz x nμ, n_t) - scattered electron flux
• B2B_inelastic: Scattering probability matrix (nμ, nμ) - pitch angle redistribution
• n: Neutral density profile (n_z,) - altitude-dependent density [m⁻³]
• Ie_slice: Incident electron flux at the current energy (nz x nμ, n_t) - flux before scattering

calculate_scattering_matrices(θ_lims, n_direction = 720; verbose = true)

Calculate the scattering matrices for given pitch-angle limits θ_lims of the electron beams. Uses 720 directions by default for the start angle and the scattering angles, equivalent to (1/4)° steps.

Calling

P_scatter, Ω_subbeam_relative, θ₁ = calculate_scattering_matrices(θ_lims, n_direction)

Inputs

• θ_lims: pitch-angle limits of the electron beams (e.g. 180:-10:0), where 180° corresponds to field aligned down, and 0° field aligned up.
• n_direction: number of directions or sub-beams to use for the discretized calculations of the scattering matrices. Defaults to 720 when left empty.

Outputs

• P_scatter: probabilities for scattering in 3D from beam to beam. Matrix [n_direction x n_direction x n_beam]
• Ω_subbeam_relative: relative weight of each sub-beam within each beam, normalized so that summing along the sub-beams gives 1 for each beam. Matrix [n_beam x n_direction]
• θ₁: scattering angles used in the calculations. Vector [n_direction]

calculate_scattering_matrices_legacy(θ_lims, n_direction=720; verbose = true)

Calculate the scattering matrices for given pitch-angle limits θ_lims of the electron beams.

Calling

Pmu2mup, theta2beamW, BeamWeight_relative, θ₁ = calculate_scattering_matrices_legacy(θ_lims, n_direction)

Inputs

• θ_lims: pitch-angle limits of the electron beams (e.g. 180:-10:0), where 180° corresponds to field aligned down, and 0° field aligned up.
• n_direction: number of directions or sub-beams to use for the discretized calculations of the scattering matrices. Defaults to 720 when left empty.

Outputs

• Pmu2mup: probabilities for scattering in 3D from beam to beam. Matrix [n_direction x 2 * n_direction x n_beam]
• theta2beamW: weight of each sub-beam within each beam. Matrix [n_beam x n_direction]
• BeamWeight_relative: relative weight of each sub-beam within each beam. It is the same as theta2beamW but normalized so that summing along the sub-beams gives 1 for each beam. Matrix [n_beam x n_direction]
• θ₁: scattering angles used in the calculations. Vector [n_direction]

calculate_volume_excitation(z, t, Ie_ztE_omni, σ, n)

Calculate the volume-excitation-rate for an excitation of interest, produced by the electron flux Ie_ztE_omni that is summed over the beams (omnidirectional).

The excitation of interest is chosen through the cross-section σ given to the function. Note that the neutral density n should match the excitation of interest (e.g. use nN2 when calculating the volume-excitation-rate of the 4278Å optical emission).

Calling

Q = calculate_volume_excitation(z, t, Ie_ztE, σ, n)`

Inputs

• z: altitude (m). Vector [n_z]
• t: time (s). Vector [n_t]
• Ie_ztE_omni: omnidirectional electron flux (#e⁻/m²/s). 3D array [n_z, n_t, n_E]
• σ: excitation cross-section (m⁻²). Vector [n_E]
• n: density of exciteable atmospheric specie (m⁻³). Vector [n_z]

compute_F(Ie, μ_lims, v; vpar_edges=nothing) -> (F, vpar_edges, vpar_centers, Δvpar)

Reduce AURORA electron flux Ie to a 1D distribution F(v_parallel) while conserving total electron density.

If vpar_edges is nothing, a symmetric uniform grid is generated automatically with default_vpar_edges(v), spanning [-maximum(v), maximum(v)] and including an exact edge at v_parallel = 0.

compute_f(Ie, ΔE_J, ΔΩ, v) -> f

Convert AURORA electron flux Ie to full phase-space density f on the same [Nz, nμ, Nt, nE] indexing.

create_crank_nicolson_nzval_mappings(Mlhs, Mrhs, n_z, n_angle)

Create mappings for both Mlhs and Mrhs matrices. Returns (mappinglhs, mappingrhs).

create_crank_nicolson_sparsity_patterns(n_z, n_angle, μ, D, Ddiffusion)

Create the sparsity patterns for both Mlhs and Mrhs matrices. The structure is the same for both, only values differ.

Returns (Mlhs, Mrhs) with correct sparsity structure.

create_steady_state_nzval_mapping(Mlhs, n_z, n_angle)

Create a mapping from matrix block positions to nzval indices for efficient in-place updates.

Julia uses the CSC (Compressed Sparse Column) format for sparse matrices, with three arrays:

• colptr: Column pointers (which rows are in each column)
• rowval: Row indices of non-zero values
• nzval: The actual non-zero values

This function computes a mapping that tells us which index in nzval corresponds to each matrix element we want to update. With this mapping, we can directly modify nzval in place, avoiding expensive matrix reconstruction.

Returns a structured mapping where mapping[i1, i2][(:type, position)] gives the nzval index for the specified matrix element in block (i1, i2).

create_steady_state_sparsity_pattern(n_z, n_angle, μ, D, Ddiffusion)

Create the sparsity pattern (structure) of the steady-state LHS matrix once. This can be reused by only modifying nzval values, avoiding allocations.

Returns the sparse matrix with the correct structure.

default_vpar_edges(v) -> vpar_edges

Construct a symmetric, uniformly spaced v_parallel grid spanning [-maximum(v), maximum(v)] with an exact edge at v_parallel = 0.

downsampling_fluxes(directory_to_process, downsampling_factor)

This function extracts Ie from the simulation results in directory_to_process and downsample it in time. For example: if Ie is given with a time step of 1ms and we use a downsampling_factor of 10, this function will extract the values of Ie with a time step of 10ms. It will then save the results in a new subfolder calleddownsampled_10x, inside the directory_to_process.

Calling

downsampling_fluxes(directory_to_process, downsampling_factor)

Inputs

• directory_to_process: absolute or relative path to the simulation directory to process.
• downsampling_factor: downsampling factor for the time

Outputs

The downsampled electron fluxes Ie will be saved in a subfolder inside the directory_to_process.

excitation_4278(E)

Returns the electron-impact excitation cross-section for 4278Å from N2 ions.

Best fit to experiments (Tima Sergienko, personal communication).

Calling

σ = excitation_4278(E)

Input

• E: incoming electron energy (eV)

Output

• σ: emission cross-section (m²)

excitation_6730_N2(E)

Returns the electron-impact excitation cross-section for 6730Å bands from N2 for transitions 4-1 and 5-2.

Digitised and extrapolated from Lanchester et al. (2009), p. 2545. https://doi.org/10.5194/angeo-27-2881-2009

% parent N2 % products 0 (emission) % threshold 7.3532 (???) % units eV

Calling

σ = excitation_6730_N2(E)

Input

• E: incoming electron energy (eV)

Output

• σ: emission cross-section (m²)

excitation_7774_O(E)

Returns the electron-impact excitation cross-section for 7774Å from O. Produced by the transition groud-state –(different pathways)–> OI 3p 5P –(prompt emission of 7774Å)–> OI 3s 5S

Digitized and extrapolated from Julienne and Davis (1976), p. 1397. https://doi.org/10.1029/JA081i007p01397

% typed in by Mina Ashrafi July 2008 % parent O % products 0 (emission) % threshold 10.74 (Itikawa Corrected:BG-20191016) % units eV

Calling

σ = excitation_7774_O(E)

Input

• E: incoming electron energy (eV)

Output

• σ: emission cross-section (m²)

excitation_7774_O2(E)

Returns the electron-impact excitation cross-section for 7774Å from O2. Produced by the dissociative ionization/excitation of O2 –(different pathways)–> OI 3p 5P –(prompt emission of 7774Å)–> OI 3s 5S

From Erdman and Zipf (1987), p. 4540. https://doi.org/10.1063/1.453696

% typed in by Mykola/Nickolay Ivshenko August 2006 % parent O2 % products 0 (emission) % threshold 15.9 (???) Seems OK, 10.74 + 5.15 (O2-bond-energy) % units eV

Calling

σ = excitation_7774_O2(E)

Input

• E: incoming electron energy (eV)

Output

• σ: emission cross-section (m²)

excitation_8446_O(E)

Returns the electron-impact excitation cross-section for 8446Å from O. Produced by the transition ground-state –(different pathways)–> OI 3p 3P –(prompt emission of 8446Å)–> ??

From Itikawa and Ichimura (1990). https://doi.org/10.1063/1.555857

% parent O % products 0 (emission) % threshold 10.99 (Itikawa) Confirmed/BG-20191016 % units eV

Calling

σ = excitation_8446_O(E)

Input

• E: incoming electron energy (eV)

Output

• σ: emission cross-section (m²)

excitation_8446_O2(E)

Returns the electron-impact excitation cross-section for 8446Å from O2. Produced by the dissociative ionization/excitation of O2 –(different pathways)–> OI 3p 3P –(prompt emission of 8446Å)–> ??

From Schulman et al. (1985). Digitised from fig.7 by Daniel Whiter in March 2009. Value in paper is 2 +/- 15% @ 100eV. https://doi.org/10.1103/PhysRevA.32.2100

% parent O2 % products 0 (emission) % threshold 15.9 (???) - seems low, dissociation energy of O2 % is 5.16 eV and the energy-level of % O(3p3P) is 10.99 eV, so threshold has % to be 16.15 eV /BG 20180529

Calling

σ = excitation_8446_O2(E)

Input

• E: incoming electron energy (eV)

Output

• σ: emission cross-section (m²)

excitation_O1D(E)

Returns the electron-impact excitation cross-section for O1D (will emit 6300Å, not prompt).

From Itikawa and Ichimura (1990). https://doi.org/10.1063/1.555857

Calling

σ = excitation_O1D(E)

Input

• E: incoming electron energy (eV)

Output

• σ: emission cross-section (m²)

excitation_O1S(E)

Returns the electron-impact excitation cross-section for O1S (will emit 5577Å, not prompt).

From Itikawa and Ichimura (1990). https://doi.org/10.1063/1.555857

Calling

σ = excitation_O1S(E)

Input

• E: incoming electron energy (eV)

Output

• σ: emission cross-section (m²)

find_Ietop_file(path_to_directory)

Look for Ie_incoming file present in the directory given by path_to_directory. If several files are starting with the name "Ie_incoming", return an error. If only one file is found, return a string with the path to that file.

Calling

Ietop_file = find_Ietop_file(path_to_directory)

Inputs

• path_to_directory: path to a directory

Returns

• Ietop_file: path to the Ie_incoming file, in the form "pathtodirectory/Ieincoming*.mat"

find_cascading_file(E_grid, species_dir)

Search for a pre-computed cascading spectra file with matching energy grid.

Arguments

• E_grid: Energy grid to match
• species_dir: Directory containing cascading data files for the species

Returns

• (file_found, filepath): Tuple of boolean and filepath string

find_iri_file(; year=2018, month=12, day=7, hour=11, minute=15,
                lat=76, lon=5, height=85:1:700)

Find or create an IRI model data file for the specified conditions.

It first searches for an existing IRI file matching the given parameters. If no matching file is found, it calculates new IRI data using the Python iri2020 package and saves it to a file. The iri2020 package will compile and run some fortran code under the hood.

Keyword Arguments

• year::Int=2018: Year
• month::Int=12: Month (1-12)
• day::Int=7: Day of month (1-31)
• hour::Int=11: Hour in Universal Time (0-23)
• minute::Int=15: Minute (0-59)
• lat::Real=76: Geographic latitude in degrees North
• lon::Real=5: Geographic longitude in degrees East
• height::AbstractRange=85:1:700: Altitude range in km

Returns

• String: Full path to the IRI data file

Notes

• Default parameters correspond to the VISIONS-2 rocket launch conditions
• Files are stored in internal_data/data_electron/ directory

find_msis_file(; year=2018, month=12, day=7, hour=11, minute=15,
                lat=76, lon=5, height=85:1:700)

Find or create a MSIS model data file for the specified conditions.

It first searches for an existing MSIS file matching the given parameters. If no matching file is found, it calculates new MSIS data using the Python pymsis package and saves it to a file. The pymsis package will download, compile and run some fortran code under the hood.

Keyword Arguments

• year::Int=2018: Year
• month::Int=12: Month (1-12)
• day::Int=7: Day of month (1-31)
• hour::Int=11: Hour in Universal Time (0-23)
• minute::Int=15: Minute (0-59)
• lat::Real=76: Geographic latitude in degrees North
• lon::Real=5: Geographic longitude in degrees East
• height::AbstractRange=85:1:700: Altitude range in km

Returns

• String: Full path to the MSIS data file

Notes

• Default parameters correspond to the VISIONS-2 rocket launch conditions
• Files are stored in internal_data/data_neutrals/ directory

find_scattering_matrices(θ_lims, n_direction=720)

Look for scattering matrices that match the pitch-angle limits θ_lims and the number of direction/sub-beams n_direction. If a file is found, the scattering matrices are directly loaded. Otherwise, they are calculated and saved to a file.

Calling

P_scatter, Ω_subbeam_relative, θ₁ = find_scattering_matrices(θ_lims, n_direction)

Inputs

• θ_lims: pitch-angle limits of the electron beams (e.g. 180:-10:0), where 180° corresponds to field aligned down, and 0° field aligned up.
• n_direction: number of directions or sub-beams to use for the discretized calculations of the scattering matrices. Defaults to 720 when left empty.

Outputs

• P_scatter: probabilities for scattering in 3D from beam to beam. Matrix [n_direction x n_direction]
• Ω_subbeam_relative: relative weight of each sub-beam within each beam, normalized so that summing along the sub-beams gives 1 for each beam. Matrix [n_beam x n_direction]
• θ₁: scattering angles used in the calculations. Vector [n_direction]

get_cross_section(species_name, energy_grid)
get_cross_section(species_name, E_centers::AbstractVector)

Calculate the cross-section for a given species and their different energy states.

Calling

σ_N2 = get_cross_section("N2", energy_grid) σ_N2 = get_cross_section("N2", E_centers)

Inputs

• species_name: name of the species. String
• energy_grid: an EnergyGrid object, or
• E_centers: energy bin centers (eV). Vector [n_E]

Outputs

• σ_species: A matrix of cross-section values for each energy state, for the defined species

initialize_transport_matrices(model, t)

Create a TransportMatrices container initialized with zeros for A, B, D, Q and Ddiffusion.

interpolate_O_parameters(E_primary)

Interpolate energy-dependent parameters for atomic O ionization.

Arguments

• E_primary: Primary electron energy (eV)

Returns

• (A_factor, B_factor): Tuple of interpolated parameters

interpolate_iri_to_grid(iri_data, h_atm)

Interpolate IRI data to a custom altitude grid.

This function interpolates all altitude-dependent densities and temperatures from IRI to a new altitude grid, while preserving scalar parameters (peak densities/heights, TEC, etc.).

Arguments

• iri_data: NamedTuple from load_iri_data() or load_iri().data
• h_atm: Target altitude grid (m). Vector [nZ]

Returns

A NamedTuple with the same fields as iri_data, but with interpolated profiles:

• height_km: converted target altitude grid (km). Vector [nZ]
• ne: interpolated electron density (m⁻³). Vector [nZ]
• Tn: interpolated neutral temperature (K). Vector [nZ]
• Ti: interpolated ion temperature (K). Vector [nZ]
• Te: interpolated electron temperature (K). Vector [nZ]
• nO⁺, nH⁺, nHe⁺, etc.: interpolated ion densities (m⁻³). Vector [nZ]
• NmF2, hmF2, NmF1, etc.: preserved scalar parameters (unchanged from input)

interpolate_msis_to_grid(msis_data, h_atm)

Interpolate MSIS data to a custom altitude grid.

This function interpolates all altitude-dependent densities and temperature from MSIS to a new altitude grid. The interpolation is performed in log space for densities (which provides exponential extrapolation) and linear space for temperature.

Arguments

• msis_data: NamedTuple from load_msis_data() or load_msis().data
• h_atm: Target altitude grid (m). Vector [nZ]

Returns

A NamedTuple with the following fields (all vectors of length nZ):

• height_km: converted target altitude grid (km)
• N2: interpolated N₂ density (m⁻³)
• O2: interpolated O₂ density (m⁻³)
• O: interpolated O density (m⁻³)
• T: interpolated temperature (K)

Notes

• Densities are interpolated in log space for better exponential behavior
• Temperature is interpolated in linear space

interpolate_profile(data_values, data_altitude_km, target_altitude_m;
                   log_interpolation=true)

Interpolate a single profile from one altitude grid to another.

This is a helper function for interpolating individual data profiles. Interpolation can be performed in linear or logarithmic space.

Arguments

• data_values::Vector: The data to interpolate (e.g., density or temperature)
• data_altitude_km::Vector: Altitude grid of the input data (km)
• target_altitude_m::Vector: Target altitude grid for interpolation (m)

Keyword Arguments

• log_interpolation::Bool=true: If true, interpolation is done in log space (exponential extrapolation). Recommended for densities. Use false for temperatures.

Returns

• interpolated::Vector: Interpolated data on the target altitude grid

load_Ie(path)

Load a IeFlickering-NN.mat result file produced by AURORA.

Returns a named tuple with fields:

• Ie : [Nz, n_mu, Nt, nE] number flux [m^-2 s^-1], integrated over ΔE and ΔΩ
• E : energy bin left edges [eV], length nE
• μ_lims : cosine pitch-angle bin limits, length n_mu + 1
• t_run : time vector [s], length Nt
• h_atm : altitude grid [m], length Nz

load_cascading_matrices(filepath)

Load pre-computed cascading matrices from a file.

Arguments

• filepath: Path to the .mat file containing cascading data

Returns

• (Q_transfer_matrix, E_grid_for_Q, ionization_thresholds): Tuple of loaded data

load_cross_sections(energy_grid)
load_cross_sections(E_centers::AbstractVector)

Load the cross-sections of the neutrals species for their different energy states.

Calling

σ_neutrals = load_cross_sections(energy_grid) σ_neutrals = load_cross_sections(E_centers)

Inputs

• energy_grid: an EnergyGrid object, or
• E_centers: energy bin centers (eV). Vector [n_E]

Returns

• σ_neutrals: A named tuple containing the cross-sections (m⁻²) for N2, O2, and O.

load_electron_densities(iri_file, h_atm)

Load the electron density and temperature from an IRI file that was generated and saved using AURORA's IRI interface. Then interpolate the profiles over AURORA's altitude grid.

Calling

ne, Te = load_electron_densities(iri_file, h_atm)

Inputs

• iri_file: absolute path to the iri file to read ne and Te from. String
• h_atm: altitude (m). Vector [nZ]

Returns

• ne: e- density (m⁻³). Vector [nZ]
• Te: e- temperature (K). Vector [nZ]

See also

load_iri, interpolate_iri_to_grid

load_excitation_threshold()

Load the excitation thresholds or energy levels of the different states (vibrational, rotational, ionization, ...) of the neutrals species as specified in the XXlevels.dat files. The corresponding names of the states can be found in the XXlevels.name files. XX refers to N2, O2 or O.

Calling

collision_levels = load_excitation_threshold()

Returns

• collision_levels: A named tuple of matrices, namely (N2_levels, O2_levels, O_levels). The matrices have shape [n_levels x 2]. The first column contains the energy levels and the second column contains the number of secondaries associated to that level (is non-zero only for ionized states).

load_iri(iri_file)

Load an IRI file and return both parameters and data.

Arguments

• iri_file: Path to the IRI file to load

Returns

A NamedTuple with two fields:

• parameters: NamedTuple with (year, month, day, hour, minute, lat, lon, height)
• data: NamedTuple with all the IRI data columns

load_iri_data(iri_file)

Load IRI model data from a file into a structured NamedTuple.

Reads an IRI data file and organizes all columns into a NamedTuple with descriptive field names for easy access to densities, temperatures, and other ionospheric parameters.

Arguments

• iri_file::String: Path to the IRI data file

Returns

• NamedTuple: IRI data with the following fields (all vectors except where noted):
  • height_km: Altitude grid (km)
  • ne: Electron density (m⁻³)
  • Tn: Neutral temperature (K)
  • Ti: Ion temperature (K)
  • Te: Electron temperature (K)
  • nO⁺, nH⁺, nHe⁺, nO2⁺, nNO⁺, nCI, nN⁺: Ion densities (m⁻³)
  • NmF2, NmF1, NmE: Peak densities at F2, F1, E layers
  • hmF2, hmF1, hmE: Peak heights at F2, F1, E layers (km)
  • TEC: Total Electron Content
  • EqVertIonDrift: Equatorial vertical ion drift
  • foF2: F2 critical frequency

Notes

• All profile data are vectors with length equal to the iri altitude grid

load_msis(msis_file)

Load a MSIS file and return both parameters and data.

Arguments

• msis_file: Path to the MSIS file to load

Returns

A NamedTuple with two fields:

• parameters: NamedTuple with (year, month, day, hour, minute, lat, lon, height)
• data: NamedTuple with all the MSIS data columns

load_msis_data(msis_file)

Load MSIS model data from a file into a structured NamedTuple.

Reads a MSIS data file and organizes all columns into a NamedTuple with descriptive field names for easy access to densities, temperature, and other atmospheric parameters.

Arguments

• msis_file::String: Path to the MSIS data file

Returns

• NamedTuple: MSIS data with the following fields (all vectors):
  • height_km: Altitude grid (km)
  • air: Total mass density (kg/m³)
  • N2: Molecular nitrogen density (m⁻³)
  • O2: Molecular oxygen density (m⁻³)
  • O: Atomic oxygen density (m⁻³)
  • He: Helium density (m⁻³)
  • H: Atomic hydrogen density (m⁻³)
  • Ar: Argon density (m⁻³)
  • N: Atomic nitrogen density (m⁻³)
  • anomalousO: Anomalous oxygen density (m⁻³)
  • NO: Nitric oxide density (m⁻³)
  • T: Temperature (K)

Notes

• All data are vectors with length equal to the msis altitude grid

load_neutral_densities(msis_file, h_atm)

Load the neutral densities and temperature from a MSIS file that was generated and saved using AURORA's MSIS interface. Then interpolate the profiles over AURORA's altitude grid.

Upper boundary conditions are applied to smoothly transition the densities to zero.

Calling

n_neutrals, Tn = load_neutral_densities(msis_file, h_atm)

Inputs

• msis_file: absolute path to the msis file to read n_neutrals and Tn from. String
• h_atm: altitude (m). Vector [nZ]

Returns

• n_neutrals: neutral densities (m⁻³). Named tuple of vectors ([nZ], ..., [nZ])
• Tn: neutral temperature (K). Vector [nZ]

See also

load_msis, interpolate_msis_to_grid

load_parameters_iri(iri_file)

Load calculation parameters from an IRI data file header.

Reads the header section of an IRI file and extracts the input parameters that were used for the IRI model calculation.

Arguments

• iri_file::String: Path to the IRI data file

Returns

• NamedTuple: Parameters with fields:
  • year::Int: Year
  • month::Int: Month (1-12)
  • day::Int: Day (1-31)
  • hour::Int: Hour (0-23)
  • minute::Int: Minute (0-59)
  • lat::Real: Latitude (degrees North)
  • lon::Real: Longitude (degrees East)
  • height::AbstractRange: Altitude range (km)

load_parameters_msis(msis_file)

Load calculation parameters from a MSIS data file header.

Reads the header section of a MSIS file and extracts the input parameters that were used for the MSIS model calculation.

Arguments

• msis_file::String: Path to the MSIS data file

Returns

• NamedTuple: Parameters with fields:
  • year::Int: Year
  • month::Int: Month (1-12)
  • day::Int: Day (1-31)
  • hour::Int: Hour (0-23)
  • minute::Int: Minute (0-59)
  • lat::Real: Latitude (degrees North)
  • lon::Real: Longitude (degrees East)
  • height::AbstractRange: Altitude range (km)

load_scattering_matrices(θ_lims)

Load the scattering matrices for the given pitch-angle limits.

Calling

μ_lims, μ_center, scattering = load_scattering_matrices(θ_lims)

Inputs

• θ_lims: pitch angle limits of the e- beams (deg). Vector [n_beam + 1]

Outputs

• μ_lims: cosine of the pitch angle limits of the e- beams. Vector [n_beam + 1]
• μ_center: cosine of the pitch angle of the middle of the e- beams. Vector [n_beam]
• scattering: Tuple with several of the scattering informations, namely scattering = (P_scatter, Ω_subbeam_relative, Ω_beam)
  • P_scatter: probabilities for scattering in 3D from beam to beam. Matrix [n_direction x n_direction]
  • Ω_subbeam_relative: relative contribution from within each beam. Matrix [n_beam x n_direction]
  • Ω_beam: solid angle for each stream (ster). Vector [n_beam]
  • θ_scatter: scattering angles used in the calculations. Vector [n_direction]

loss_to_thermal_electrons(E, ne, Te)

Suprathermal electron energy loss function due to electron-electron collisions.

This function calculates the electron energy loss function due to electron-electron interaction. It uses the analytic form given for the energy-transfer rate from photoelectrons (or suprathermal electrons) to thermal electrons, given by Swartz and Nisbet (1971). The expression fits the classical formulation of Itikawa and Aono (1966) at low energies and gives a smooth transition to fit the quantum mechanical equation of Schunk and Hays (1971).

Arguments

• E::Real: Energy level [eV]. Scalar value.
• ne::Vector: Ambient electron concentration [/m³], length nZ.
• Te::Vector: Electron temperature [K], length nZ.

Returns

• Le::Vector: Electron energy loss function [eV/m], length nZ.

Notes

The paper by Swartz and Nisbet uses electron density in cm⁻³; here the constant is rescaled to use m⁻³ instead. We calculate the loss function dE/ds(E,ne,Te) directly and not as in Swartz and Nisbet dE/ds(E,ne,Te)/ne.

The loss is set to zero when the suprathermal electron energy E is below the thermal electron energy Ee = kB*Te/qₑ.

References

• Swartz, W. E., J. S. Nisbet, and A. E. S. Green (1971), Analytic expression for the energy transfer rate from photoelectrons to thermal electrons, J. Geophys. Res., 76(34), 8425-8426, doi: 10.1029/JA076i034p08425.
• Itikawa, Y., and O. Aono (1966), Energy change of a charged particle moving in a plasma, Phys. Fluids, 9, 1259-1261.
• Schunk, R. W., and P. B. Hays (1971), Photoelectron energy losses to thermal electrons, Planet. Space Sci., 19, 113-117.

make_Ie_top_file(directory_to_process)

Reads into a folder directory_to_process containing results from an AURORA.jl simulation and extracts the particle flux Ie (#e⁻/m²/s) at the top of the ionosphere (i.e. at the max altitude used in the simulation).

Calling

make_Ie_top_file(directory_to_process)

Inputs

• directory_to_process: absolute or relative path to the simulation directory to process.

make_altitude_grid(bottom_altitude, top_altitude; dz_max=25)

Create an altitude grid based on the altitude limits given as input. It uses constant steps of 150m for altitudes below 100km, and a non-linear grid above 100km.

Calling

h_atm = make_altitude_grid(bottom_altitude, top_altitude; max_dz=25)

Arguments

• bottom_altitude: altitude, in km, for the bottom of the simulation
• top_altitude: altitude, in km, for the top of the simulation

Keyword Arguments

• dz_max = 25: maximum step size, in km. Relevant for high altitudes where the numbers can get large.

Outputs

• h_atm: altitude (m). Vector [nZ]

make_column_excitation_file(directory_to_process)

Reads into a folder directory_to_process containing results from an AURORA.jl simulation, loads the volume excitation rates Q_XXXX (#excitation/m³/s) contained in the file Qzt_all_L.mat and integrate them in height, taking into account the finite speed of light.

The calculated colum-integrated excitation rates are saved to a file named I_lambda_of_t.mat. The column-integrated excitation rates are named "I_4278, I_6730, ...". They are all vectors in time (length n_t), and have units of (#excitation/m²/s).

Note that the function make_volume_excitation_file() needs to be run before this one, as we need the file Qzt_all_L.mat with the volume excitation rates.

Calling

make_current_file(directory_to_process)

Inputs

• directory_to_process: absolute or relative path to the simulation directory to process.

make_current_file(directory_to_process)

Reads into a folder directory_to_process containing results from an AURORA.jl simulation, loads the particle flux Ie (#e⁻/m²/s) and calculates the field-aligned current-density and field-aligned energy-flux for each height and through time.

The following variables are saved to a file named J.mat:

• J_up: Field-aligned current-density in the upward direction. 2D array [n_z, n_t]
• J_down: Field-aligned current-density in the downward direction. 2D array [n_z, n_t]
• IeE_up: Field-aligned energy-flux (eV/m²/s) in the upward direction. 2D array [n_z, n_t]
• IeE_down: Field-aligned energy-flux (eV/m²/s) in the downward direction. 2D array [n_z, n_t]

Calling

make_current_file(directory_to_process)

Inputs

• directory_to_process: absolute or relative path to the simulation directory to process.

make_density_file(directory_to_process)

This function reads into a folder directory_to_process containing results from an AURORA.jl simulation. It loads the particle flux Ie (#e⁻/m²/s), calculates the superthermal e- density n_e (#e⁻/m³) from it, and saves n_e into a new file "superthermal_e_density.mat".

The particle flux Ie is defined along a magnetic field line and over an (Energy, pitch_angle)-grid. The number density n_e calculated is given along a magnetic field line and over an energy grid. That way, we have the density of electrons with a certain energy at a specific altitude and time.

Calling

make_density_file(directory_to_process)

Inputs

• directory_to_process: absolute or relative path to the simulation directory to process

make_energy_grid(E_max)

Create an energy grid based on the maximum energy E_max given as input.

Calling

E_edges, E_centers, ΔE = make_energy_grid(E_max)

Inputs

• E_max: upper limit for the energy grid (in eV)

Outputs

• E_edges: energy bin edges (eV). Vector [nE + 1] (includes the last upper edge)
• E_centers: energy bin centers (eV). Vector [nE]
• ΔE: energy bin widths (eV). Vector [nE]

make_heating_rate_file(directory_to_process)

Reads into a folder directory_to_process containing results from an AURORA.jl simulation, loads the particle flux Ie (#e⁻/m²/s), and calculates the heating rate of thermal electrons by superthermal electrons.

The heating rate is the rate at which energy is transferred from superthermal electrons to thermal electrons through Coulomb collisions. It is saved as a function of altitude and time.

Calling

make_heating_rate_file(directory_to_process)

Inputs

• directory_to_process: absolute or relative path to the simulation directory to process.

make_psd_file(directory_to_process; compute=:both, vpar_edges=nothing, output_prefix="psd")

Reads into a folder directory_to_process containing results from an AURORA.jl simulation, loads particle flux files IeFlickering-NN.mat, converts them into phase-space density, and writes one output file per input as psd/psd-NN.mat.

Calling

make_psd_file(directory_to_process; compute=:both, vpar_edges=nothing, output_prefix="psd")

Inputs

• directory_to_process: absolute or relative path to the simulation directory to process.

Keyword Arguments

• compute: one of :f_only, :F_only, or :both.
• vpar_edges: custom v_parallel bin edges [m/s] used when computing F. If nothing, an automatic symmetric uniform interval grid is used, spanning [-maximum(v), maximum(v)] with an edge at v_parallel = 0.
• output_prefix: output filename prefix used as <output_prefix>-NN.mat.

make_psd_from_AURORA(path_to_file; compute=:f_only, vpar_edges=nothing)

Load an AURORA result file and compute full phase-space density f, the reduced distribution F(v_parallel), or both.

If vpar_edges is nothing and F is requested, the function uses an automatic symmetric uniform interval grid.

make_savedir(root_savedir, name_savedir; behavior = "default")

Return the path to the directory where the results will be saved. If the directory does not already exist, create it.

If the constructed savedir already exists and contains files starting with "IeFlickering-", a new directory is created to avoid accidental overwriting of results (e.g., savedir(1), savedir(2), etc.).

Calling

savedir = make_savedir(root_savedir, name_savedir)
savedir = make_savedir(root_savedir, name_savedir; behavior = "custom")

Arguments

• root_savedir::String: The root directory where the data will be saved. If empty or contains only spaces, it defaults to "backup".
• name_savedir::String: The name of the subdirectory to be created within root_savedir. If empty or contains only spaces, it defaults to the current date and time in the format "yyyymmdd-HHMM".
• behavior::String (optional): Determines how the full path is constructed.
  • "default": The path will be built starting under the data/ folder of the AURORA installation (i.e., AURORA_folder/data/root_savedir/name_savedir/, where AURORA_folder is the folder containing the AURORA code). This is the default behavior.
  • "custom": The path will be built as root_savedir/name_savedir/, with the argument root_savedir treated as an absolute or relative path. This allows for saving results in any location on the system. Useful if AURORA is installed as a dependency to some other project.

Returns

• savedir::String: The full path to the directory where the results will be saved.

make_volume_excitation_file(directory_to_process)

Reads into a folder directory_to_process containing results from an AURORA.jl simulation, loads the particle flux Ie (#e⁻/m²/s), and calculates the volume-excitation-rates. For prompt emissions, volume-excitation-rates correspond also to volume-emission-rates.

Calling

make_volume_excitation_file(directory_to_process)

Inputs

• directory_to_process: absolute or relative path to the simulation directory to process.

mu_avg(θ_lims)

Calculate the cosinus of the center of each pitch-angle beams delimited by θ_lims. This is for isotropically distributed fluxes within each beam, i.e the fluxes are weighted by sin(θ)

Calling

μ_center = mu_avg(θ_lims)

Inputs

• θ_lims : pitch-angle limits in degrees of all the beams, range or vector [n_beams + 1]

Outputs

• μ_center : cosine of the center of all the pitch-angle beams, vector [n_beams + 1]

n_neutrals(iono::Ionosphere)

Return neutral densities as a NamedTuple (; nN2, nO2, nO).

psd_grids(E_eV, μ_lims_cosine)

Compute energy, angle, and velocity grids needed by compute_f.

q2colem(t::Vector, z, Q, A = 1, τ = ones(length(z)))

Integrate the volume-excitation-rate (#exc/m³/s) to column-excitation-rate (#exc/m²/s).

Takes into account the time-delay between light emitted at different altitudes. Photons emitted at at altitude of 200km will arrive at the detector 100e3/3e8 = 0.333 ms later than electrons emitted at an altitude of 100km. This is a small time-shift, but it is close to the time-differences corresponding to the phase-shifts between auroral emissions varying at ~10Hz.

The einstein coefficient A and effective lifetime τ are optional (equal to one by default).

Calling

I = q2colem(t, z, Q, A, τ)

Inputs

• z: altitude (m). Vector [n_z]
• t: time (s). Vector [n_t]
• Q: volume-excitation-rate (#exc/m³/s) of the wavelength of interest. 2D array [n_z, n_t]
• A: einstein coefficient (s⁻¹). Scalar (Float or Int)
• τ: effective lifetime (s). Vector [n_z].

Output

• I: integrated column-excitation-rate (#exc/m²/s) of the wavelength of interest. Vector [n_t]

q2colem(t::Real, z, Q, A = 1, τ = ones(length(z)))

Same as above, except time is now a scalar (steady-state results). This is just a simple integration in height.

rename_if_exists(savefile)

This function takes a string as an input. If a file or folder with that name does not exist, it returns the same string back. But if the folder or file already exists, it appends a number between parenthesis to the name string.

For example, if the folder foo/ already exist and "foo" is given as input, the function will return a string "foo(1)" as an output. Similarly, if a file foo.txt already exists and "foo.txt" is given as input, the function will return a string "foo(1).txt". If the file foo(1).txt also already exist, the function will return a string "foo(2).txt", etc...

The function should support all types of extensions.

Calling

newsavefile = rename_if_exists(savefile)

restructure_streams_of_Ie(Ie, θ_lims, new_θ_lims)

Function that merges the streams of Ie that are given over θ_lims to fit the new_θ_lims of interest. It can be useful when wanting to merge some streams for plotting.

For example, if we have θ_lims = [180 160 140 120 100 90 80 60 40 20 0], and we want to plot with new_θ_lims = [(180, 160), (160, 120)], the function will keep the first stream as is and merge the streams (160°-140°) and (140°-120°) together into a new stream with limits (160°-120°).

Important: The limits in new_θ_lims need to match some existing limits in θ_lims. In the example above, new_θ_lims = [(180, 165)] would not have worked because 165° is not a limit that exists in θ_lims.

Entries in new_θ_lims can be nothing to leave a panel empty.

Calling

Ie_plot = restructure_streams_of_Ie(Ie, θ_lims, new_θ_lims)

Arguments

• Ie: array of electron flux with pitch-angle limits θ_lims. Of shape [n_z, n_μ, n_t, n_E].
• θ_lims: pitch-angle limits. Usually a vector or range.
• new_θ_lims: new pitch-angle limits. Given as an array of tuples with angles in the range 0-180° (where 180° is field-aligned down, 0° is field-aligned up). Use nothing for empty panels. For example:

julia> new_θ_lims = [(180, 170)  (170, 150)  (150, 120)  (120, 100)  (100, 90);  # DOWN
                     (0, 10)     (10, 30)    (30, 60)    (60, 80)    (80, 90)]   # UP

Returns

• Ie_plot: array of electron flux with the new pitch-angle limits new_θ_lims. Of shape [n_z, n_μ_new, n_t, n_E], where n_μ_new is the number of streams in new_θ_lims. The second dimension of Ie_plot is sorted such that the indices go along the first row of new_θ_lims, and then the second row. In our example with new_θ_lims from above, that would be $[1 2 3 4 5; 6 7 8 9 10]$.

save_cascading_matrices(Q_transfer_matrix, E_grid_for_Q, ionization_thresholds, species_dir, species_name)

Save calculated cascading matrices to a file.

Arguments

• Q_transfer_matrix: Transfer matrix to save
• E_grid_for_Q: Energy grid used for calculations
• ionization_thresholds: Ionization threshold energies
• species_dir: Directory to save the file
• species_name: Name of the species (for filename)

save_iri_data(iri_data, parameters)

Save IRI model data to a text file with metadata header.

Creates a formatted text file containing the IRI calculation parameters in the header followed by the data matrix. The filename is automatically generated based on the input parameters, and if a file with the same name exists, a unique name is created.

Arguments

• iri_data::Matrix: IRI data matrix from calculate_iri_data() (with header row)
• parameters::NamedTuple: Parameters used for IRI calculation, must contain:
  • year, month, day, hour, minute: Time specification
  • lat, lon: Location (degrees)
  • height: Altitude range (km)

Returns

• String: Full path to the created file

File Format

The file contains:

1. Header section with input parameters
2. Column headers
3. Data matrix (one row per altitude)

Filename Convention

iri_YYYYMMDD-HHMM_LATN-LONE.txt

Notes

• Files are saved to internal_data/data_electron/ directory
• Existing files are not overwritten; a suffix is added to the filename

save_msis_data(msis_data, parameters)

Save MSIS model data to a text file with metadata header.

Creates a formatted text file containing the MSIS calculation parameters in the header followed by the data matrix. The filename is automatically generated based on the input parameters, and if a file with the same name exists, a unique name is created.

Arguments

• msis_data::Matrix: MSIS data matrix from calculate_msis_data() (with header row)
• parameters::NamedTuple: Parameters used for MSIS calculation, must contain:
  • year, month, day, hour, minute: Time specification
  • lat, lon: Location (degrees)
  • height: Altitude range (km)

Returns

• String: Full path to the created file

File Format

The file contains:

1. Header section with input parameters
2. Column headers
3. Data matrix (one row per altitude)

Filename Convention

msis_YYYYMMDD-HHMM_LATN-LONE.txt

Notes

• Files are saved to internal_data/data_neutrals/ directory
• Existing files are not overwritten; a suffix is added to the filename

search_existing_iri_file(; year, month, day, hour, minute, lat, lon, height)

Search for an existing IRI data file matching the specified parameters.

This function scans the internal_data/data_electron/ directory for IRI files with matching time, location, and altitude grid parameters. It performs a quick pre-check on filenames before loading and comparing full parameters.

Keyword Arguments

• year::Int: Year
• month::Int: Month (1-12)
• day::Int: Day of month (1-31)
• hour::Int: Hour in Universal Time (0-23)
• minute::Int: Minute (0-59)
• lat::Real: Geographic latitude in degrees North
• lon::Real: Geographic longitude in degrees East
• height::AbstractRange: Altitude range in km

Returns

• Union{String, Nothing}: Full path to matching file, or nothing if not found

search_existing_msis_file(; year, month, day, hour, minute, lat, lon, height)

Search for an existing MSIS data file matching the specified parameters.

This function scans the internal_data/data_neutrals/ directory for MSIS files with matching time, location, and altitude grid parameters. It performs a quick pre-check on filenames before loading and comparing full parameters.

Keyword Arguments

• year::Int: Year
• month::Int: Month (1-12)
• day::Int: Day of month (1-31)
• hour::Int: Hour in Universal Time (0-23)
• minute::Int: Minute (0-59)
• lat::Real: Geographic latitude in degrees North
• lon::Real: Geographic longitude in degrees East
• height::AbstractRange: Altitude range in km

Returns

• Union{String, Nothing}: Full path to matching file, or nothing if not found

steady_state_scheme_optimized!(Ie, model, matrices, iE, Ie_top, cache; first_iteration = false)

Optimized steady-state scheme using direct nzval modification. This is an in-place version that modifies Ie directly to avoid allocations. This version avoids allocations by reusing the sparse matrix structure.

On first iteration, creates the sparsity pattern and mapping which are stored in cache. On subsequent iterations, only updates the nzval array directly.

Mathematical Background

The steady-state electron transport equation is:

μ ∂Ie/∂z + A*Ie - ∫B*Ie'dΩ' = Q

After spatial discretization, this becomes a linear system of coupled equations:

[μ*Ddz + A - B - D*Ddiffusion] * Ie = Q
         ↑
        Mlhs (the system matrix)

Where:

• Ddz = Ddz_Up or Ddz_Down (depending on sign of μ): spatial derivative operator
• A: electron loss rate matrix (diagonal)
• B: scattering operator matrix (couples different angles)
• D: diffusion coefficient (diagonal in angle space)
• Ddiffusion: second derivative operator for pitch-angle diffusion
• Q: source term (energy degradation and ionization)

The resulting sparse matrix Mlhs has a block structure:

┌─────────┬─────────┬─────────┐
│ Block   │ Block   │ Block   │  Each block is n_z × n_z
│ (1,1)   │ (1,2)   │ (1,3)   │  (n_z = number of altitudes)
├─────────┼─────────┼─────────┤
│ Block   │ Block   │ Block   │  Off-diagonal blocks (i1≠i2):
│ (2,1)   │ (2,2)   │ (2,3)   │  represent angular scattering (B matrix)
├─────────┼─────────┼─────────┤
│ Block   │ Block   │ Block   │  Diagonal blocks (i1=i2):
│ (3,1)   │ (3,2)   │ (3,3)   │  transport + loss + diffusion
└─────────┴─────────┴─────────┘

In the context of solving f(Ie) = Mlhs*Ie - Q = 0, the matrix Mlhs is the Jacobian:

Jacobian = ∂f/∂Ie = Mlhs

Arguments

• Ie: pre-allocated output array [m⁻² s⁻¹], size (nz * nangle) to store results
• model: AuroraModel (s_field and pitch_angle_grid.μ_center are used)
• matrices::TransportMatrices: container with
  • A: electron loss rate [s⁻¹]
  • B: scattering matrix [s⁻¹], size (nz × nangle × n_angle)
  • D: pitch-angle diffusion coefficient [s⁻¹], size (n_angle,)
  • Q: source term [m⁻² s⁻¹], size (nz × nangle × n_energy)
  • Ddiffusion: spatial diffusion matrix, size (nz × nz)
• iE: current energy index
• Ie_top: boundary condition at top [m⁻² s⁻¹]
• cache: Cache object storing Mlhs, mapping, KLU, and differentiation matrices
• first_iteration: whether this is the first call (creates structure) or subsequent (reuses structure)

update_crank_nicolson_matrices!(Mlhs, Mrhs, mapping_lhs, mapping_rhs,
                                A, B, D, Ddt, Ddiffusion, Ddz_Up, Ddz_Down, μ, z)

Update both Mlhs and Mrhs using pre-computed mappings. This is the fast path with zero allocations.

update_matrices!(matrices, model, phase_fcn_neutrals, iE, B2B_fragment)

Update the A and B matrices in place for a given energy level iE.

Arguments

• matrices::TransportMatrices: Container to update
• model: AuroraModel (grids + atmosphere + physics)
• phase_fcn_neutrals: Phase functions for all species
• iE: Current energy index
• B2B_fragment: Pre-computed beam-to-beam fragments

Returns

• B2B_inelastic_neutrals: Array of inelastic beam-to-beam matrices for cascading calculations

update_steady_state_matrix!(Mlhs, mapping, A, B, D, Ddiffusion, Ddz_Up, Ddz_Down, μ, z)

Update the sparse matrix values using pre-computed mapping. This avoids all allocations by directly modifying nzval.

v_of_E(E)

Calculate the velocity (in m/s) of an electron with energy E (in eV).

Calling

v = v_of_E(E)

Input

• E : energy in eV, can be a scalar, vector, range, ...

Output

• v : velocity in m/s

validate_θ_lims(θ_lims)

Validate that the pitch-angle limits θ_lims are correctly specified. Throws an ArgumentError if:

• θ_lims does not include 180° (field-aligned downward)
• θ_lims does not include 0° (field-aligned upward)
• θ_lims is not in descending order

Inputs

• θ_lims: pitch-angle limits of the electron beams (e.g. 180:-10:0)