Reference

AbstractGrid

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

AbstractMode

Abstract supertype for simulation modes.

Concrete subtypes:

• SteadyStateMode: solves each time step independently as a steady-state problem
• TimeDependentMode: time-stepping with a Crank-Nicolson scheme

The mode controls both the numerical scheme used at each energy step and how time parameters (if any) are resolved into an internal time configuration.

AbstractModulation

Abstract supertype for all temporal modulation types.

Concrete subtypes: ConstantModulation, SinusoidalFlickering, SquareFlickering, SmoothOnset.

AbstractSpectrum

Abstract supertype for all electron energy spectrum shapes.

Concrete subtypes: FlatSpectrum, GaussianSpectrum, MaxwellianSpectrum, FileSpectrum.

AbstractTimeConfig

Abstract supertype for resolved time configurations used by simulations.

Concrete subtypes:

• SingleStepConfig: explicit single-step configuration for steady-state simulations
• UniformTimeGrid: simple uniform grid for multi-step steady-state simulations
• RefinedTimeGrid: CFL-refined grid for time-dependent simulations

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, SP, 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;
            species = (N2Species(msis_file), O2Species(msis_file), OSpecies(msis_file)))

Build a lightweight AuroraModel. Only the grids and ionosphere are set up; the expensive computation (scattering matrices, species densities, cross sections, cascading data) is deferred to initialize!(model), which is called automatically by run!(sim).

Inputs

• altitude_lims: altitude limits (km) for the bottom and top of the ionosphere
• θ_lims: pitch-angle beam limits (e.g. 180:-10:0), 180° = field-aligned down
• E_max: upper energy-grid limit (eV)
• msis_file: path to the MSIS atmosphere file
• iri_file: path to the IRI ionosphere file
• B_angle_to_zenith: angle between the magnetic field and zenith (degrees, default 0)

Keyword Arguments

• species: tuple (or any iterable) of NeutralSpecies instances that the model should simulate. Defaults to the standard atmosphere: N₂, O₂, O from msis_file. Pass a custom tuple to add, remove, or replace species:

  # Two species example:
  model = AuroraModel(...; species = (O2Species(msis_file), OSpecies(msis_file)))
  
  # Four species example, with a custom 4th gas
  # Laws (density, cascading, phase) must be @law-wrapped, a functor, or a named function to
  # ensure reproducibility when saved to physics_state.jld2:
  custom_sp = NeutralSpecies(:MyGas, @law(z -> 1e18 .* exp.(-z ./ 30e3));
                             cascading_spec = my_spec, phase_fcn_generator = phase_fcn_N2)
  model = AuroraModel(...; species = (N2Species(msis_file), O2Species(msis_file),
                                      OSpecies(msis_file),  custom_sp))
  # Interception window: pre-populate custom cross sections before run!/initialize!
  model.species[end].cross_sections    = my_sigma_matrix   # [n_levels × n_E]
  model.species[end].excitation_levels = my_levels_matrix  # [n_levels × 2]
  run!(sim)

Returns

An uninitialized AuroraModel. Call initialize!(model) or run!(sim) to complete setup.

AuroraOutputManager(savedir; overwrite=false, compress=false)

Configuration struct that controls where and how simulation output is written.

Arguments

• savedir: directory where output files will be written (created by run! if absent). May be a relative or absolute path. If empty (or only whitespace), output goes to backup/<yyyymmdd-HHMM>/ in the current working directory.

Keyword arguments

• overwrite::Bool=false: if false (the default), run! errors when simulation_data.nc already exists in savedir; set to true to allow overwriting.

• compress: zlib compression level for all NetCDF variables.

  • false or 0 (default): no compression
  • true: level 4
  • 1–9: exact deflate level

  Note about compression: from experience, zlib saves only a few percent of file size while slowing writes and reads severalfold. Use with caution.

Output layout

savedir/
├── config.toml               # scalar simulation parameters
├── inputs/
│   ├── atmosphere.nc         # neutral/ionosphere density profiles
│   └── physics_state.jld2   # full AuroraModel (species + physics matrices)
├── simulation_data.nc        # time-dependent flux output (appended per loop)
└── analysis/                 # derived quantities (written by make_* functions)

Examples

# Full control:
out = AuroraOutputManager("my_run"; compress=true)   # deflate level 4
out = AuroraOutputManager("my_run"; compress=6)      # deflate level 6
sim = AuroraSimulation(model, flux, out; mode=TimeDependentMode(duration=0.5, dt=0.001))

# Convenience — pass a plain String and defaults apply:
sim = AuroraSimulation(model, flux, "my_run"; mode=SteadyStateMode())

AuroraSimulation{M<:AuroraModel, F<:InputFlux, S<:AbstractMode}

A complete simulation configuration, bundling the physical model, input flux, mode strategy, output configuration, resolved time configuration, and lazy simulation cache.

Use initialize! to allocate the workspace explicitly, or call run! directly to auto-initialize and execute the simulation.

Examples

# Single-step steady-state
sim = AuroraSimulation(model, flux, "my_run"; mode=SteadyStateMode())

# Multi-step steady-state
sim = AuroraSimulation(model, flux, "my_run"; mode=SteadyStateMode(duration=0.5, dt=0.01))

# Time-dependent
sim = AuroraSimulation(model, flux, "my_run";
                       mode=TimeDependentMode(duration=0.5, dt=0.001, CFL_number=128))

# Full output control via AuroraOutputManager
out = AuroraOutputManager("my_run"; overwrite=true, compress=false)
sim = AuroraSimulation(model, flux, out; mode=SteadyStateMode())

BlockIndices

Pre-computed nzval indices for a single (i1, i2) block of the transport matrix.

For off-diagonal blocks (i1 ≠ i2) only diag is populated (scattering coupling), and the boundary fields are zero.

For diagonal blocks (i1 = i2) the tridiagonal entries and boundary condition indices are stored. bc_last_sub is zero for downward streams (no sub-diagonal in the last-row boundary condition).

CachePolicy(; force_recompute=false, save_cache=true,
            cache_root=default_cache_root())

Control how internal scattering and cascading caches are reused.

• force_recompute: ignore compatible cache files already present on disk
• save_cache: skip writing newly computed caches when set to false
• cache_root: parent directory that contains the e_cascading/ and e_scattering/ cache subdirectories. By default this is pkgdir(AURORA, "internal_data").

ColumnExcitationResult

Result of the column-integrated excitation rate computation. Contains the column-integrated emission intensities as a function of time, as saved in analysis/column_excitation.nc.

Fields

• I_4278, I_6730, I_7774, I_8446, I_O1D, I_O1S: column-integrated intensities (photons/m²/s) as Vector{Float64} (n_t)
• I_7774_O, I_7774_O2: species-resolved contributions to I_7774
• I_8446_O, I_8446_O2: species-resolved contributions to I_8446
• t: time grid (s) as Vector{Float64}

ConstantModulation()

No temporal modulation — the flux is constant in time.

Used as the default modulation for steady-state simulations and for file-based input where the time dependence is already contained in the file.

DegradedCascadingIntegrand(E_primary_bin_min, E_primary_bin_max, threshold, secondary_law)

Callable integrand used by hcubature to integrate the degraded-primary transfer matrix.

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)

ExprLaw

Serialization-friendly wrapper around a user-supplied physical law. This can be a density profile, a cascading secondary-distribution law, a phase-function generator, etc. It stores the law's source as a string so the model can round-trip through JLD2: only the string is written to disk, and the callable is rebuilt by parsing and evaluating that string on load.

Build one with the @law macro:

density         = @law z -> 1e18 .* exp.(-(z .- 100e3) ./ 30e3)
secondary_e_law = @law (E_s, E_p) -> 1 / (11.4^2 + E_s^2)

To ensure reproducibility, laws must be self-contained: closed-form expressions referencing only their arguments and names available in AURORA/Base. A law that needs to carry e.g. parameters should be constructed as a functor struct instead. To ensure this, trying to build a @law that captures local variables will throw an error.

FileSpectrum(filename; interpolation=:constant)

Energy spectrum loaded from a .mat file.

The file contains the complete flux distribution (possibly time-dependent). When used inside an InputFlux, this spectrum type requires ConstantModulation because the file already contains the full temporal behavior.

Arguments

• filename: path to the .mat file containing the flux data

Keyword Arguments

• interpolation=:constant: interpolation scheme for resampling the file's time grid. Either :constant, :linear, or :pchip.

File format

The .mat file must contain:

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

Examples

FileSpectrum("path/to/flux.mat")
FileSpectrum("path/to/flux.mat"; interpolation=:linear)

FlatSpectrum(IeE_tot; E_min=0.0)

Flat (constant) differential number flux spectrum above a minimum energy.

The flux is uniform in #e⁻/m²/s/eV for energies above E_min, and zero below. Normalized so that the integrated field-aligned (vertical) energy flux equals IeE_tot.

Arguments

• IeE_tot: total field-aligned (vertical) energy flux (W/m²), i.e. the energy flux crossing the horizontal top boundary.

Keyword Arguments

• E_min=0.0: minimum energy threshold (eV). Flux is zero below this energy.

Examples

FlatSpectrum(1e-2)                  # flat from lowest energy
FlatSpectrum(1e-2; E_min=2900.0)    # flat above 2900 eV

GaussianSpectrum(IeE_tot, E₀, ΔE)

Gaussian energy spectrum centered at E₀ with width ΔE.

The spectral shape is Φ(E) ∝ exp(-(E - E₀)² / ΔE²), normalized so that the integrated field-aligned (vertical) energy flux equals IeE_tot.

Arguments

• IeE_tot: total field-aligned (vertical) energy flux (W/m²), i.e. the energy flux crossing the horizontal top boundary.
• E₀: center energy (eV)
• ΔE: energy width (eV)

Examples

GaussianSpectrum(1e-2, 5000.0, 500.0)

IeTopResult

Electron flux at the top altitude of the model, as saved in analysis/Ie_top.nc by make_Ie_top_file. It contains all beams — both downward (precipitation) and upward (backscattered) — and is derived from the simulation output, so it is distinct from the Ie_input boundary condition stored in simulation_data.nc.

Fields

• Ietop: flux array of size (nbeams x nt x n_E), units: m⁻² s⁻¹
• t: time grid (s), relative so always start at 0
• E_edges: energy bin edges (eV), length n_E + 1
• E_centers: energy bin centres (eV), length n_E
• ΔE: energy bin widths (eV), length n_E
• μ_lims: cosine-of-pitch-angle limits, length n_beams + 1

InputFlux{S<:AbstractSpectrum, M<:AbstractModulation}

A unified electron precipitation input, combining an energy spectrum shape with a temporal modulation and beam selection.

Fields

• spectrum: the energy spectrum shape (any AbstractSpectrum)
• modulation: the temporal modulation (any AbstractModulation)
• beams: indices of the electron beams with precipitating flux
• z_source: altitude of the precipitation source (km)

Constructors

# Steady-state (constant modulation is the default)
InputFlux(FlatSpectrum(1e-2; E_min=2900); beams=1:2)
InputFlux(MaxwellianSpectrum(1e-2, 1000); beams=1:2)

# Time-dependent with modulation
InputFlux(FlatSpectrum(1e-2; E_min=100), SinusoidalFlickering(5.0);
          beams=1, z_source=3000)

# From file (modulation must be and is by default ConstantModulation)
InputFlux(FileSpectrum("path.mat"))

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

Electron background: electron density and temperature. Neutral species densities are owned by the individual NeutralSpecies objects.

MSISDensity(msis_file, species)

Density profile backed by an MSIS file. Callable on an altitude grid (m), returns the density vector (m⁻³) for the requested species (:N2, :O2, or :O).

Example

profile = MSISDensity(find_msis_file(), :N2)
n = profile(altitude_grid.h)   # Vector{Float64}, length == length(altitude_grid.h)

MaxwellianSpectrum(IeE_tot, E₀; low_energy_tail=true)

Maxwellian energy spectrum with optional low-energy tail (LET).

Based on the corrected implementation of Meier/Strickland/Hecht/Christensen JGR 1989 (pages 13541-13552).

Arguments

• IeE_tot: total field-aligned (vertical) energy flux (W/m²), i.e. the energy flux crossing the horizontal top boundary.
• E₀: characteristic energy (eV)

Keyword Arguments

• low_energy_tail=true: include a low-energy tail following Meier et al. 1989

Examples

MaxwellianSpectrum(1e-2, 1000.0)
MaxwellianSpectrum(1e-2, 1000.0; low_energy_tail=false)

NeutralSpecies{S<:CascadingSpec, C<:SpeciesCascadingCache, P}

All per-species data needed to advance the transport equation through one neutral species.

Fields

• name::Symbol: short identifier (e.g. :N2, :O2, :O)
• density_profile: callable h_atm (m) → density (m⁻³) used to (re)sample density. Can be an MSISDensity, a VectorDensity, or any callable. Untyped so it can be replaced freely before calling initialize!(model).
• density::Vector{Float64}: density profile sampled on the model altitude grid (m⁻³). Empty until initialize!(model) is called.
• cross_sections::Matrix{Float64}: collision cross sections, shape [n_levels × n_E] (m²). Auto-loaded from the built-in library when empty; pre-populate before initialize!(model) to supply custom data for a non-standard species. Empty until initialize!(model) is called.
• excitation_levels::Matrix{Float64}: excitation thresholds and secondary counts, shape [n_levels × 2]. First column = threshold energy (eV), second column = number of secondaries produced (non-zero only for ionizing channels). Auto-loaded from the built-in library when empty; pre-populate before initialize!(model) to supply custom data for a non-standard species. Empty until initialize!(model) is called.
• phase_fcn_generator: callable (θ, E) -> (phaseE, phaseI) used to (re)build phase_fcn whenever the pitch-angle or energy grid changes. Untyped so it can be replaced freely before calling initialize!(model).
• phase_fcn::P: tuple (phaseE, phaseI) of [n_θ × n_E] matrices materialized from the generator on the model's scattering θ grid and energy centers. Holds 0×0 placeholder matrices until initialize!(model) is called.
• cascading_spec::S: ionization thresholds and secondary-distribution law used to build the cascading transfer matrices
• cascading_data::C: cascading transfer matrices (populated lazily by load_or_compute_cascading!)

NeutralSpecies(name, density_profile; cascading_spec, phase_fcn_generator)

Build a lightweight NeutralSpecies. All grid-dependent fields (density, cross_sections, excitation_levels, phase_fcn) start empty and are populated by initialize!(model).

OperatorDiagonals

Dense vectors extracted once from the sparse finite-difference operators Ddz_Up, Ddz_Down and Ddiffusion. Storing them avoids repeated CSC look-ups during the value-update phase.

Fields with suffix _diag hold the main diagonal, _sub the sub-diagonal (offset −1), and _super the super-diagonal (offset +1).

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

Pitch-angle grid for electron beams.

RefinedTimeGrid

Resolved time configuration for time-dependent simulations.

The coarse save grid (t_save) is constructed from dt and duration, then the fine internal grid (t) is obtained by subdividing each save interval until the CFL condition is satisfied. The end result is a grid with CFL_factor sub-steps per save interval.

Loops are partitioned by save intervals using a balanced partition (loop sizes differ by at most one), so every loop holds between fld(n_save, n_loop) and cld(n_save, n_loop) save intervals and none is ever empty (for n_loop ≤ n_save).

See loop_save_count, loop_save_start, loop_internal_count, loop_internal_start for helpers to index into these grids per loop.

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

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

DegradedCascadingIntegrand(E_primary_bin_min, E_primary_bin_max, threshold, secondary_law)

Callable integrand used by hcubature to integrate the secondary transfer matrix.

SimulationResult

The raw electron-flux output of a simulation, as loaded from simulation_data.nc by load_results.

Fields

• Ie : electron number flux [n_z, n_μ, n_t, n_E] (m⁻² s⁻¹). See the note on Output & data: this is a number flux, already integrated over each energy bin and beam solid angle.
• t : time axis [n_t] (s)
• h_atm : altitude grid [n_z] (m)
• E_centers : energy bin centres [n_E] (eV)
• E_edges : energy bin edges [n_E+1] (eV)
• ΔE : energy bin widths [n_E] (eV) (= diff(E_edges))
• μ_lims : pitch-angle cosine bin boundaries [n_μ+1]
• beam_weights: solid-angle beam weights [n_μ] (sr)
• savedir : directory the result was loaded from, or nothing

SingleStepConfig()

Explicit single-step time configuration.

Used for steady-state simulations without temporal evolution.

SinusoidalFlickering(f; amplitude=1.0)

Sinusoidal temporal modulation of the electron flux.

The flux oscillates between (1 - amplitude) and 1 of the base flux, using the pattern 1 - cos²(πft). With amplitude=1.0, the flux modulates fully between 0 and 1.

Arguments

• f: modulation frequency (Hz)

Keyword Arguments

• amplitude=1.0: modulation depth. 0 = constant flux, 1 = full on/off modulation.

Examples

SinusoidalFlickering(10.0)                  # 10 Hz, full modulation
SinusoidalFlickering(5.0; amplitude=0.5)    # 5 Hz, partial modulation

SmoothOnset(t_start, t_end)

Smooth flux onset using a C∞-smooth transition function.

The flux transitions smoothly from 0 to 1 over the interval [t_start, t_end]. Before t_start, the flux is 0; after t_end, the flux is 1.

When t_start == t_end, this acts as a step function (Heaviside) at t_start.

Arguments

• t_start: start time of the transition (s)
• t_end: end time of the transition (s)

Examples

SmoothOnset(0.0, 0.1)    # smooth ramp from t=0 to t=0.1 s
SmoothOnset(0.0, 0.0)    # step function at t=0

SquareFlickering(f; amplitude=1.0)

Square-wave temporal modulation of the electron flux.

The flux alternates between (1 - amplitude) and 1 of the base flux at the given frequency. With amplitude=1.0, the flux alternates fully between 0 and 1.

Arguments

• f: modulation frequency (Hz)

Keyword Arguments

• amplitude=1.0: modulation depth. 0 = constant flux, 1 = full on/off modulation.

Examples

SquareFlickering(10.0)                  # 10 Hz, full modulation
SquareFlickering(5.0; amplitude=0.5)    # 5 Hz, partial modulation

SteadyState()
SteadyState(; duration, dt)

Convenience alias for SteadyStateMode.

SteadyStateMode()
SteadyStateMode(; duration, dt)

Solve each time step independently as a steady-state (time-independent) transport problem.

When called without arguments, a single steady-state solve is performed.

When called with duration and dt, a UniformTimeGrid is built and each point is solved independently. This enables, e.g., a flickering input where each time step is solved in steady-state.

Examples

# Single-step steady-state
sim = AuroraSimulation(model, flux, savedir; mode=SteadyStateMode())

# Multi-step steady-state (each time point solved independently)
sim = AuroraSimulation(model, flux, savedir; mode=SteadyStateMode(duration=0.5, dt=0.01))

TimeDependent(; duration, dt, CFL_number=64, max_memory_gb=8.0, n_loop=nothing)

Convenience alias for TimeDependentMode.

TimeDependentMode(; duration, dt, CFL_number=64, max_memory_gb=8.0, n_loop=nothing)

Solve the transport equation in a time-dependent manner with a Crank-Nicolson scheme.

Keyword Arguments

• duration: total simulation time (s)
• dt: time step for saving data (s). The internal time step may be finer to satisfy the CFL condition.
• CFL_number = 64: CFL stability factor. Crank-Nicolson is unconditionally stable, so large values are acceptable (it will affect the accuracy though).
• max_memory_gb = 8.0: memory budget (GB) used to auto-split the simulation into loops.
• n_loop = nothing: manually specify the number of loops. Overrides max_memory_gb.

Examples

sim = AuroraSimulation(model, flux, savedir;
                       mode=TimeDependentMode(duration=0.5, dt=0.001, CFL_number=128))

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.

The scratch buffers that depend on the number of species and on the scattering grid (B2B_inelastic_neutrals, phase_fcn_e, phase_fcn_i) are created empty here and sized later in initialize_transport_matrices.

UniformTimeGrid(duration, dt)

A simple uniform time grid for multi-step steady-state simulations.

Each time point is solved independently — no CFL refinement or loop splitting is needed.

VectorDensity(h, n)

Density profile defined by user-supplied altitude (h, m) and density (n, m⁻³) vectors. Callable on any altitude grid (m); evaluates via PCHIP interpolation in log-space, consistent with AURORA's MSIS interpolation convention.

Example

profile = VectorDensity(h_msis_m, n_N2)
n = profile(altitude_grid.h)

VolumeExcitationResult

Result of the volume excitation rate computation. Contains the volume excitation/ionization rates on the altitude x time grid, as saved in analysis/volume_excitation.nc.

Fields

• Q4278, Q6730, Q7774, Q8446, QO1D, QO1S: volume excitation rates (photons/m³/s) as Matrix{Float64} (nz x nt)
• Q7774_O, Q7774_O2: species-resolved contributions to Q7774
• Q8446_O, Q8446_O2: species-resolved contributions to Q8446
• QOi, QO2i, QN2i: ionization rates (/m³/s) as Matrix{Float64} (nz x nt)
• h_atm: altitude grid (m) as Vector{Float64}
• t: time grid (s) as Vector{Float64}
• savedir: directory from which the result was loaded, or nothing

Crank_Nicolson!(Ie, t, model, v, matrices, iE, Ie_top, I0, cache)

Solve the time-dependent electron transport equation for energy level iE using the Crank-Nicolson implicit scheme.

On the first call the sparse matrix structures, nzval index arrays, and operator diagonals are computed and stored in cache. On subsequent calls only the numerical values in Mlhs.nzval / Mrhs.nzval are updated (zero allocations on the hot path).

Mathematical Background

The time-dependent electron transport equation is:

∂Ie/∂t + μ·v ∂Ie/∂z + A·Ie − ∫B·Ie′ dΩ′ − D·∇²Ie = Q

Crank-Nicolson gives second-order accuracy in time:

Mlhs · Ie^(n+1)  =  Mrhs · Ie^n  +  (Q^(n+1) + Q^n)/2

with

Mlhs = 1/(v·Δt) + μ·Ddz/2 + A/2 − B/2 − D·Ddiffusion
Mrhs = 1/(v·Δt) − μ·Ddz/2 − A/2 + B/2 + D·Ddiffusion

Both matrices share the same block structure as the steady-state system:

┌─────────┬─────────┬─────────┐
│ 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 (n_z * n_angle, n_t)
• 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, B, D, Q, Ddiffusion
• iE: current energy index
• Ie_top: boundary condition at top [m⁻² s⁻¹] at each time step
• I0: initial condition [m⁻² s⁻¹]
• cache: SolverCache storing Mlhs, Mrhs, indices, op_diags, KLU

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.

N2Species(density_profile)
N2Species(msis_file::AbstractString)

Construct the default N₂ species with cross sections from the built-in e_N2_* library, phase function from phase_fcn_N2, and cascading described by DefaultCascadingSpecN2.

density_profile can be an MSISDensity, a VectorDensity, or any callable h_atm (m) → density (m⁻³). Passing an MSIS file path string is a shorthand for MSISDensity(msis_file, :N2). Grid-dependent fields are populated later by initialize!(model).

O2Species(density_profile)
O2Species(msis_file::AbstractString)

Default O₂ species, analogous to N2Species.

OSpecies(density_profile)
OSpecies(msis_file::AbstractString)

Default O species, analogous to N2Species.

animate_Ie_in_time(directory_to_process;
                   angles_to_plot=nothing, colorrange=nothing, save_to_file=true,
                   plot_input=false, input_angle_cone=[170, 180], dt_steps=1,
                   framerate=30, max_bytes=512*1024^2)

Plot a heatmap of Ie over height and energy, and animate it in time. The animation is saved as a .mp4 file under the directory_to_process if save_to_file = true.

Requires a Makie backend (e.g. using CairoMakie or using GLMakie).

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 an animation.mp4 file (or animation_with_precipitation.mp4 when plot_input = true) inside directory_to_process.
• plot_input = false: if true, also plots the precipitating electron flux at the top of the ionosphere by loading it from simulation_data.nc.
• input_angle_cone = [170, 180]: pitch-angle cone (degrees) used to select and sum beams for the precipitation overlay. Only used when plot_input = true.
• dt_steps = 1: plot one frame every dt_steps timesteps. Increase to speed up rendering.
• framerate = 30: framerate of the animation in frames per second.
• max_bytes = 512 * 1024^2: per-chunk memory budget for streaming the flux from simulation_data.nc.

Example

# Use default angles and colorrange:
julia> animate_Ie_in_time(directory_to_process)

# Use 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_input=true)

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

Notes

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.

append_chunk_nc!(ds, Ie_chunk, t_chunk, sim)

Reshape Ie_chunk from [n_z·n_μ, n_t, n_E] to [n_z, n_μ, n_t, n_E] and append it — together with t_chunk — to the unlimited time dimension of the open NCDataset ds. Calls sync after writing so output is crash-safe.

apply_modulation(mod::AbstractModulation, t_shifted)

Apply temporal modulation to a time grid, returning a vector of modulation factors (values between 0 and 1) for each time step.

t_shifted is the time grid shifted for energy- and angle-dependent electron travel times, important in case a z_source was specified in the InputFlux.

apply_modulation(::ConstantModulation, t_shifted)

Returns ones — no temporal modulation. The flux is simply turned on at t=0.

apply_modulation(mod::SinusoidalFlickering, t_shifted)

Sinusoidal modulation: oscillates between (1 - amplitude) and 1.

apply_modulation(mod::SmoothOnset, t_shifted)

Smooth onset using a C∞-smooth transition function over [t_start, t_end].

apply_modulation(mod::SquareFlickering, t_shifted)

Square-wave modulation: alternates between (1 - amplitude) and 1.

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]

build_spatial_operators(z; half_weight=false)

Construct the upwind/downwind spatial differentiation matrices Ddz_Up and Ddz_Down from the altitude grid z.

When half_weight=true (used by Crank-Nicolson), the coefficients are halved so that the matrices represent Ddz/2 directly.

calculate_cascading_matrices(spec, E_edges; verbose=true)

Calculate the energy-degradation transfer matrix for a species defined by its CascadingSpec. The matrices can later be used to directly get the degraded primary and secondary electron distributions for any given primary energy and ionization threshold.

The outer loop structure is identical for all species. The only species-specific ingredients are

• spec.secondary_law(E_secondary, E_primary) -> Float64, which describes how secondary electrons distribute in energy given a primary electron at the current integration energy.
• spec.ionization_thresholds, which define the ionization thresholds for the species and thus the number of transfer matrices to calculate.

Arguments

• spec::CascadingSpec contains species name, ionization thresholds and secondary distribution law
• E_edges: Energy grid edges to match (eV)

Returns

• (primary_transfer_matrix, secondary_transfer_matrix, E_edges, ionization_thresholds): degraded-primary transfer matrix [nE, nE, nthresholds], secondary transfer matrix [nE, nE, nthresholds], energy grid edges, and ionization thresholds

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(n_z, n_μ, n_E, n_t, N_neutrals, CFL_factor, n_save;
                 max_memory_gb=8, verbose=true) → Int

Choose the number of loops (n_loop) for a time-dependent simulation so that the peak memory stays within max_memory_gb.

Using the split returned by estimate_simulation_memory, the peak memory of an n_loop-loop run is ≈ fixed_gb + scaling_gb / n_loop. The smallest n_loop satisfying fixed_gb + scaling_gb / n_loop ≤ max_memory_gb is therefore

n_loop = ceil(scaling_gb / (max_memory_gb - fixed_gb))

The result is capped at n_save, because a loop cannot hold fewer than one save interval. If even one save interval per loop (or the fixed allocations alone) exceeds the budget, the maximum split n_loop = n_save is returned together with a warning.

Arguments

• n_z, n_μ, n_E: altitude, pitch-angle and energy grid sizes
• n_t: number of internal time steps over the whole simulation (full CFL-refined grid)
• N_neutrals: number of neutral species
• CFL_factor: internal sub-steps per save interval
• n_save: total number of save intervals

Keyword Arguments

• max_memory_gb: memory budget (GB). Defaults to 8 GB.
• verbose::Bool=true: if true, print a short summary of the calculation

Returns

• n_loop::Int: recommended number of loops

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_omni, σ, 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]

check_n_loop(n_loop, n_z, n_μ, n_E, n_t, N_neutrals, CFL_factor; verbose=true)

Verify that running with n_loop loops fits in the detected RAM. Errors if the estimated peak memory exceeds the available RAM (always), and warns if it exceeds half of it (only when verbose).

Collect nzval indices for a given diagonal offset within the interior rows 2:(n_z-1). offset is 0 for main diagonal, +1 for super, -1 for sub.

compute_F!(F, F_local, Ie, μ_lims, v, vpar_edges, Δvpar) -> F

In-place core of compute_F. Accumulates the reduction into the work buffer F_local (shape [Nz, Nvpar, Nt], zeroed on entry) and writes the [Nvpar, Nz, Nt] result into F. Both buffers are reused across time-chunks when streaming; vpar_edges/Δvpar are fixed up front so every chunk maps onto the same bins.

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!(f, Ie, factor) -> f

In-place core of compute_f: write the phase-space density into f (shape matching Ie) using the precomputed factor from psd_f_factor. Folding the division into factor turns the hot loop into a single multiply, and reusing f avoids a full-chunk allocation per time-chunk.

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.

compute_flux(flux::InputFlux, model::AuroraModel, t)

Compute the full 3D electron flux array for a time-dependent simulation.

Evaluates the energy spectrum, distributes it over beams, and applies temporal modulation (including energy- and angle-dependent electron travel-time delays if a z_source was specified in the InputFlux).

The spectrum is spread over the selected downward-going beams proportionally to their solid angle, and normalized so that the field-aligned (vertical) energy flux crossing the horizontal top boundary equals IeE_tot, independent of the beam selection.

Arguments

• flux: an InputFlux describing the precipitation
• model: an AuroraModel with grids and atmosphere/ionosphere
• t: simulation time grid (s). Range or Vector

Returns

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

compute_flux(flux::InputFlux, model::AuroraModel)

Compute the electron flux array for a steady-state simulation.

This method is for steady-state runs: it evaluates the energy spectrum and distributes it over beams without any time dependence or travel-time delays. As in the time-dependent method, the spectrum is normalized so that the field-aligned (vertical) energy flux equals IeE_tot, independent of the beam selection.

Only ConstantModulation is allowed for steady-state. Other modulation types will raise an error.

Arguments

• flux: an InputFlux describing the precipitation
• model: an AuroraModel with grids and atmosphere

Returns

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

create_simulation_nc(sim) → NCDataset

Create <savedir>/simulation_data.nc with all dimensions and static variables. Returns the open dataset; the caller is responsible for calling close(ds).

The input boundary flux is written immediately into the file (before any solver loops run) as a separate Ie_input variable on its own fixed time_input dimension. The unlimited time dimension is left empty at creation time and is populated by subsequent append_chunk_nc! calls.

create_transport_sparsity_pattern(n_z, n_angle, μ, D, Ddiffusion; include_rhs=false)

Build the sparse matrix with the correct sparsity structure for the transport system. The returned matrix contains placeholder values; the actual physics values are filled in later by the update functions.

When include_rhs=true a second matrix Mrhs (same structure but without boundary rows) is also returned – used by Crank-Nicolson.

Returns Mlhs or (Mlhs, Mrhs).

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.

estimate_simulation_memory(n_z, n_μ, n_E, n_t, N_neutrals, CFL_factor)
    → (; fixed_gb, scaling_gb, total_gb)

Estimate the simulation's peak memory footprint (in decimal GB) by summing the actual working arrays allocated in build_simulation_cache.

The footprint is split into two parts:

• scaling_gb: arrays laid out along the internal time grid — cache.Ie, matrices.Q, the energy-degradation working arrays, and the subsampled cache.Ie_save buffer. Evaluated here at the full internal length n_t; splitting the run into n_loop loops shrinks this contribution to ≈ scaling_gb / n_loop.
• fixed_gb: arrays whose size does not depend on the loop split — most importantly cache.Ie_top, which always stores the entire input flux over the full time grid, plus cache.I0, the transport/scattering matrices and the sparse solver factorisation.

total_gb = fixed_gb + scaling_gb is the single-loop (n_loop = 1) footprint, so the peak memory of an n_loop-loop run is ≈ fixed_gb + scaling_gb / n_loop.

evaluate_spectrum(spec::AbstractSpectrum, model::AuroraModel)

Evaluate the energy spectrum on the model's energy grid, returning a vector of differential number flux per eV (#e⁻/m²/s/eV) of length n_E.

This function does not distribute over beams or apply temporal modulation. Use compute_flux for the full 3D flux array.

evaluate_spectrum(spec::FlatSpectrum, model::AuroraModel)

Evaluate a flat spectrum on the model's energy grid.

evaluate_spectrum(spec::GaussianSpectrum, model::AuroraModel)

Evaluate a Gaussian spectrum on the model's energy grid.

evaluate_spectrum(spec::MaxwellianSpectrum, model::AuroraModel)

Evaluate a Maxwellian spectrum (with optional low-energy tail) on the model's energy grid.

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²)

extract_nzval_indices(M, n_z, n_angle)

Walk the CSC structure of sparse matrix M and return a Matrix{BlockIndices} (size n_angle × n_angle) that maps each block's tridiagonal positions to their index in M.nzval.

This is computed once after sparsity-pattern creation, then used on every energy step to write values directly into nzval without any Dict look-up.

extract_operator_diagonals(Ddz_Up, Ddz_Down, Ddiffusion)

Extract the tridiagonal entries of the three finite-difference operators into dense vectors for fast element-wise access.

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

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.

foreach_Ie_time_chunk(f, sim_dir; max_bytes = 512 * 1024^2)

Stream the electron flux Ie from simulation_data.nc in sim_dir over contiguous time-chunks, calling f(Ie_chunk, t_range) for each. Ie_chunk::Array{Float64,4} has shape [n_z, n_μ, length(t_range), n_E] and t_range is the corresponding range of time indices.

The chunk length is the largest that keeps one Ie_chunk under max_bytes, so peak memory is bounded regardless of the total run length. On disk Ie is chunked one time-slice per chunk, so each slab read is contiguous.

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.

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

get_level_names(species_name)

Return the names of the excited/ionized states for a given species as a Vector{String}.

Example

get_level_names("N2")  # → ["_elastic", "_rot0_2", ..., "_ionization"]

initialize!(model::AuroraModel; verbose=true, policy=CachePolicy())

Perform all heavy setup for model:

1. Compute s_field and ionosphere from the current altitude grid.
2. Compute (or load from cache) the scattering matrices.
3. For each species: sample the density profile, load cross sections and excitation levels (skipped if already pre-populated), build phase functions, and load/compute cascading transfer matrices.

Called internally by initialize!(sim) and/or run!(sim).

initialize!(sim::AuroraSimulation;
            force_recompute=false,
            save_cache=true,
            cache_root=default_cache_root())

Allocate or re-allocate the working cache for sim.

Keywords

• force_recompute: ignore compatible on-disk cascading caches and rebuild them
• save_cache: if false, skip writing newly computed cascading caches to disk
• cache_root: parent directory that contains the cascading and scattering cache folders

initialize_transport_matrices(model, t)

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

input_flux_at_save_cadence(sim) → (t_top, Ie_top_3D)

Return the input boundary flux subsampled to the save cadence as [n_μ, n_t, n_E], together with the corresponding time axis (in seconds).

• Single-step steady-state: one slice at t = 0.
• Multi-step steady-state: all steps at the uniform time grid.
• Time-dependent: subsampled from the internal grid to the save cadence via 1 : CFL_factor : end.

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

is_multi_step(mode::SteadyStateMode)

Return true if the mode is configured for multi-step steady-state.

load_Ie_top(directory::String)

Load the electron flux at the top altitude of the model from analysis/Ie_top.nc in directory. Contains all beams (downward precipitation and upward backscatter). Returns an IeTopResult.

load_column_excitation(sim::AuroraSimulation)

Load the column-integrated excitation intensities from the simulation's save directory.

load_column_excitation(directory::String)

Load the column-integrated excitation intensities from analysis/column_excitation.nc in the given directory. Returns a ColumnExcitationResult.

load_coordinates(sim_dir) → NamedTuple

Load everything in SimulationResult except the bulky Ie array: the coordinate vectors (t, h_atm, E_centers, E_edges, ΔE, μ_lims, beam_weights), the Ie dimensions (n_z, n_μ, n_t, n_E), and savedir. Cheap to call.

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_or_compute_cascading!(cache::SpeciesCascadingCache, energy_grid; verbose=true, policy=CachePolicy())

Populate the cascading transfer matrices for one neutral species.

If a compatible (same grid, same version_AURORA) JLD2 cache file is found on disk, matrices are loaded from there. Otherwise they are computed from scratch (and optionally saved to disk depending on CachePolicy options).

load_or_compute_scattering(θ_lims, n_direction=720; verbose=true, policy=CachePolicy())

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

Calling

P_scatter, Ω_subbeam_relative, θ₁ = load_or_compute_scattering(θ_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.
• policy.cache_root: parent directory that contains the e_scattering/ cache subtree

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]

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_results(sim_dir::AbstractString; tidx=:, zidx=:, μidx=:, eidx=:, max_bytes=2*1024^3)
load_results(sim::AuroraSimulation; kwargs...)

Load the raw electron-flux output from simulation_data.nc in sim_dir (or in the simulation's save directory). Returns a SimulationResult.

The whole Ie array is read into memory eagerly. For large time-dependent runs this can be many GB; use the tidx/zidx/μidx/eidx selectors to load only a sub-box, or stream the data with foreach_Ie_time_chunk. As a safety net load_results refuses to allocate more than max_bytes (default 2 GiB) and reports how to narrow the load; pass max_bytes=Inf to disable the guard.

Selectors

Each selector is a Colon (:, the default — load everything) or a contiguous integer range (e.g. tidx = 1:100). The matching coordinate vectors and edge-derived quantities (E_edges, ΔE, μ_lims) are sliced consistently. Non-contiguous indexing is unsupported; read the NCDataset directly if you need it.

Example

res = load_results("my_run")              # everything
res = load_results("my_run"; tidx=1:100)  # first 100 time slices only
res.Ie                                     # flux [n_z, n_μ, n_t, n_E]
res.t                                      # time axis (s)

load_volume_excitation(sim::AuroraSimulation)

Load the volume excitation rates from the simulation's save directory.

load_volume_excitation(directory::String)

Load the volume excitation rates from analysis/volume_excitation.nc in the given directory. Returns a VolumeExcitationResult.

loop_internal_count(time::RefinedTimeGrid, i_loop) → Int

Number of internal time steps (columns of cache.Ie) for loop i_loop, including the shared boundary point that initialises from I0.

loop_internal_start(time::RefinedTimeGrid, i_loop) → Int

Return index in time.t (the full internal grid) of the boundary point that starts loop i_loop.

loop_save_count(time::RefinedTimeGrid, i_loop) → Int

Number of save intervals covered by loop i_loop under a balanced partition: the first rem(n_save, n_loop) loops carry cld(n_save, n_loop) intervals and the rest carry fld(n_save, n_loop). Loop sizes differ by at most one and no loop is ever empty (for n_loop ≤ n_save). The counts sum to time.n_save.

loop_save_start(time::RefinedTimeGrid, i_loop) → Int

Return index in time.t_save of the first (boundary) save point of loop i_loop.

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)

Read simulation_data.nc from directory_to_process, extract the electron flux at the top of the ionosphere (highest altitude), and write results to analysis/Ie_top.nc.

Calling

make_Ie_top_file(directory_to_process)

Inputs

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

make_Ie_top_file(sim::AuroraSimulation)

Convenience wrapper that calls make_Ie_top_file on sim.output.savedir.

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)

Read analysis/volume_excitation.nc from directory_to_process, integrate volume-excitation-rates in altitude (accounting for finite photon travel time), and write results to analysis/column_excitation.nc.

Returns a ColumnExcitationResult.

make_column_excitation_file(sim::AuroraSimulation)

Convenience wrapper that calls make_column_excitation_file on sim.output.savedir.

make_current_file(directory_to_process)

Read simulation_data.nc from directory_to_process, compute field-aligned current density and energy flux, and write results to analysis/currents.nc.

Calling

make_current_file(directory_to_process)

Keyword arguments

• max_bytes: per-chunk memory budget for streaming the flux (default 512 MiB).

make_current_file(sim::AuroraSimulation)

Convenience wrapper that calls make_current_file on sim.output.savedir.

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]
• E_centers: energy bin centers (eV). Vector [nE]
• ΔE: energy bin widths (eV). Vector [nE]

make_heating_rate_file(directory_to_process)

Read simulation_data.nc and inputs/atmosphere.nc from directory_to_process, compute the thermal-electron heating rate by superthermal electrons, and write results to analysis/heating_rate.nc.

Calling

make_heating_rate_file(directory_to_process)

Inputs

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

Keyword arguments

• max_bytes: per-chunk memory budget for streaming the flux (default 512 MiB).

make_heating_rate_file(sim::AuroraSimulation)

Convenience wrapper that calls make_heating_rate_file on sim.output.savedir.

make_psd_file(directory_to_process; compute=:both, vpar_edges=nothing, max_bytes=512*1024^2, compress=false, show_progress=false)

Read simulation_data.nc from directory_to_process, convert electron flux to phase-space density, and write results to analysis/psd.nc.

The flux is streamed over time-chunks and the phase-space density is written to disk chunk by chunk, so peak memory stays bounded even for large runs.

Keyword Arguments

• compute: one of :f_only, :F_only, or :both.
• vpar_edges: custom v_parallel bin edges [m/s]. If nothing, an automatic symmetric uniform grid is used spanning [-maximum(v), maximum(v)].
• max_bytes: per-chunk memory budget for streaming the flux (default 512 MiB).
• compress: zlib compression level for the f and F variables, with the same semantics as in AuroraOutputManager: false/0 (default, no compression), true (level 4), or an exact level 1–9.
• show_progress: show a ProgressMeter progress bar while streaming chunks (default false).

make_psd_file(sim::AuroraSimulation; kwargs...)

Convenience wrapper that calls make_psd_file on sim.output.savedir.

make_volume_excitation_file(directory_to_process)

Read simulation_data.nc and inputs/atmosphere.nc from directory_to_process, compute volume-excitation-rates for all tracked optical emissions and ionizations, and write results to analysis/volume_excitation.nc.

Returns a VolumeExcitationResult.

Keyword arguments

• max_bytes: per-chunk memory budget for streaming the flux (default 512 MiB).

make_volume_excitation_file(sim::AuroraSimulation)

Convenience wrapper that calls make_volume_excitation_file on sim.output.savedir.

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]

plot_column_excitation!(ax, data::ColumnExcitationResult;
                        wavelengths=[:I_4278, :I_6730, :I_7774, :I_8446, :I_O1D, :I_O1S],
                        kwargs...)

Plot column-integrated emission intensities onto an existing Axis.

Requires a Makie backend (e.g. using CairoMakie or using GLMakie).

Keyword Arguments

• wavelengths: vector of symbols selecting which emission lines to plot. Available values are :I_4278, :I_6730, :I_7774, :I_8446, :I_O1D, :I_O1S. Defaults to all six. O(¹D) and O(¹S) lines are plotted with a dashed linestyle.
• kwargs: additional keyword arguments are forwarded to each underlying lines! call.

plot_excitation!(ax, data::VolumeExcitationResult;
                 field=:total, time_index=nothing, colorrange=nothing,
                 colormap=:viridis, show_contours=false, show_height_of_max=true,
                 kwargs...)

Plot a volume excitation or ionization rate onto an existing Axis.

Requires a Makie backend (e.g. using CairoMakie or using GLMakie).

• If data contains one time step, or if time_index is given, plots an altitude profile (returns a Lines plot).
• Otherwise, plots a time-altitude heatmap (returns a Heatmap).

Keyword Arguments

• field = :total: which rate to plot. Can be any field of VolumeExcitationResult: :Q4278, :Q6730, :Q7774, :Q8446, :QO1D, :QO1S, :QOi, :QO2i, :QN2i, or :total which sums the three ionization rates (:QN2i + :QO2i + :QOi).
• time_index = nothing: if given, plot a single altitude profile at that time index. If nothing and the data has more than one time step, a heatmap is plotted.
• colorrange = nothing: (min, max) limits for the colorscale. Defaults to (max_value / 1e4, max_value) spanning 4 orders of magnitude.
• colormap = :viridis: Makie colormap for the heatmap.
• show_contours = false: if true, overlays black contour lines on the heatmap (heatmap mode only).
• show_height_of_max = true: if true, overlays a dashed red line at the altitude of peak rate at each time step (heatmap mode only).
• kwargs: additional keyword arguments are forwarded to the underlying heatmap! or lines! call.

plot_input(sim::AuroraSimulation)

Plot the input electron flux for a simulation.

For time-dependent simulations, produces a heatmap of flux vs energy and time for each active beam. For steady-state simulations, produces a line plot of flux vs energy.

Requires a Makie backend (e.g. using CairoMakie or using GLMakie).

Examples

using CairoMakie
sim = AuroraSimulation(model, flux, savedir;
                       mode=TimeDependentMode(duration=0.5, dt=0.001))
fig = plot_input(sim)

plot_input!(ax, data::IeTopResult;
            beams=1, colorrange=nothing, colormap=:inferno, kwargs...)

Plot the input electron flux stored in an IeTopResult onto an existing Axis.

Requires a Makie backend (e.g. using CairoMakie or using GLMakie).

Keyword Arguments

• beams = 1: pitch-angle beam index, or vector of indices, to include. Selected beams are summed and normalised by their combined solid angle, then converted to differential energy flux (#e⁻/m²/s/eV/ster).
• colorrange = nothing: (min, max) limits for the colorscale. Defaults to (max_value / 1e4, max_value) spanning 4 orders of magnitude.
• colormap = :inferno: Makie colormap for the heatmap.
• kwargs...: additional keyword arguments are forwarded to the underlying heatmap! call.

plot_model(model::AuroraModel; panels=[:all])

Plot the model setup: atmosphere, energy grid, cross-sections, phase functions, and scattering data.

Returns a Dict{Symbol, Figure} with one figure per panel. Use the panels keyword to select which panels to plot.

Requires a Makie backend (e.g. using CairoMakie or using GLMakie).

Available panels

:atmosphere, :energy_levels, :energy_grid, :cross_sections, :phase_functions, :scattering, or :all (default) to plot everything.

Examples

using CairoMakie
figs = plot_model(model)
figs = plot_model(model; panels=[:atmosphere, :cross_sections])

psd_f_factor(ΔE_J, ΔΩ, v) -> factor[nμ, nE]

Per-(beam, energy) conversion factor m_e / (ΔE_J * v^2 * ΔΩ) used by compute_f!. It depends only on the grids, so when streaming it is built once and reused across time-chunks.

psd_grids(E_centers, ΔE, μ_lims_cosine)

Compute velocity and pitch-angle grids needed by compute_f and compute_F.

• E_centers is the energy bin centers (eV). Vector [nE]
• ΔE is the energy bin widths (eV). Vector [nE]

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]

read_atmosphere_nc(sim_dir) → NamedTuple

Read inputs/atmosphere.nc from sim_dir. Returns a named tuple with fields h_atm, ne, Te, and one field per species (e.g. nN2, nO2, nO).

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)

resolve_deflatelevel(compress) -> Int

Convert a user-facing compress value into a zlib deflate level: false → 0, true → 4, an integer 0–9 → itself.

resolve_savedir(savedir) -> String

Normalise a user-supplied savedir. A non-empty path (relative or absolute) is returned unchanged; an empty or whitespace-only path falls back to backup/<yyyymmdd-HHMM> in the current working directory.

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]$.

run!(sim::AuroraSimulation)

Execute the simulation described by sim.

Automatically calls initialize! when needed, then dispatches to the appropriate execution path based on the selected mode.

Examples

sim = AuroraSimulation(model, flux, savedir;
                       mode=TimeDependentMode(duration=0.5, dt=0.001, CFL_number=128))
run!(sim)

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

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

steady_state_scheme!(Ie, model, matrices, iE, Ie_top, cache)

Solve the steady-state electron transport equation for energy level iE.

On the first call the sparse matrix structure, nzval index arrays, and operator diagonals are computed and stored in cache. On subsequent calls only the numerical values in Mlhs.nzval are updated (zero allocations on the hot path).

Mathematical Background

The steady-state electron transport equation is:

μ ∂Ie/∂z + A·Ie − ∫B·Ie′ dΩ′ = Q

After spatial discretization this becomes the linear system Mlhs · Ie = Q with:

Mlhs = μ·Ddz + diag(A) − B − D·Ddiffusion

The matrix has a block structure indexed by pitch-angle pairs (i1, i2):

┌─────────┬─────────┬─────────┐
│ 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 (i1≠i2): scattering (−B)
│ (2,1)   │ (2,2)   │ (2,3)   │
├─────────┼─────────┼─────────┤
│ Block   │ Block   │ Block   │  Diagonal (i1=i2): transport + loss + diffusion
│ (3,1)   │ (3,2)   │ (3,3)   │
└─────────┴─────────┴─────────┘

Arguments

• Ie: pre-allocated output array [m⁻² s⁻¹], size n_z * n_angle
• model: AuroraModel (s_field and pitch_angle_grid.μ_center are used)
• matrices::TransportMatrices: container with A, B, D, Q, Ddiffusion
• iE: current energy index
• Ie_top: boundary condition at top [m⁻² s⁻¹]
• cache: SolverCache storing Mlhs, indices_lhs, op_diags, KLU

update_crank_nicolson_matrices!(Mlhs, Mrhs, idx_lhs, idx_rhs,
                                A, B, D, ddt, op, μ, n_z)

Fill both Mlhs and Mrhs with the Crank-Nicolson operator values using the pre-computed BlockIndices arrays and dense OperatorDiagonals.

ddt is the scalar 1/(v·Δt) (constant for all altitudes).

The physics formulas (per stream direction) are:

Mlhs =  ddt  +  μ·Ddz  +  A/2  −  B/2  −  D·Ddiffusion
Mrhs =  ddt  −  μ·Ddz  −  A/2  +  B/2  +  D·Ddiffusion

where Ddz already contains the /2 factor (built with half_weight=true).

update_matrices!(matrices, model, 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)
• 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, indices, A, B, D, op, μ, n_z)

Fill the sparse matrix Mlhs with the steady-state transport operator values using the pre-computed indices (a Matrix{BlockIndices}) and dense operator diagonals op (OperatorDiagonals).

The physics formula for the system matrix is (per stream direction):

Mlhs = μ * Ddz  +  diag(A)  −  B  −  D * Ddiffusion

where Ddz = Ddz_Down for downward (μ < 0) or Ddz_Up for upward (μ > 0).

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

@law expr

Wrap a law expr (e.g. z -> exp(-z)) into an ExprLaw, capturing its source text so it can be saved and reloaded verbatim. Use this instead of a bare anonymous function for any density profile, cascading law, or phase-function generator.