Modifying species and grids

An AuroraModel is built in two phases:

1. AuroraModel(...) is lightweight — it sets up the grids and reads the background atmosphere, but does not compute the expensive data (densities sampled on the grid, cross-sections, phase functions, scattering and cascading matrices).
2. initialize! does that heavy work. It is called automatically by run!(sim).

The gap between the two is an interception window: anything you change on the model before initialize!/run! is picked up when the derived data is built. This is how you can customize the atmosphere without needing to put your hands in the internals.

Changing the grids

The grid fields can be reassigned on an existing model. Doing so automatically marks the model as uninitialized, so the next run!(sim) rebuilds everything for the new grid — no manual re-initialization needed.

using AURORA
msis_file = find_msis_file()
iri_file  = find_iri_file()

model = AuroraModel([100, 300], 180:-90:0, 100, msis_file, iri_file)
initialize!(model)
model.initialized

true

model.energy_grid = EnergyGrid(500)   # also: model.altitude_grid, model.pitch_angle_grid
model.initialized                     # → false: derived data is now stale

false

initialize!(model)                    # rebuilds densities, cross-sections, … for the new grid
model.energy_grid.n

121

When a simulation already exists, just change the grid and call run!(sim) — it detects the change and rebuilds the model and its cache before solving.

Changing a species' density profile

Each species carries a density_profile: a callable mapping altitude (m) to density (m⁻³). initialize! samples it onto the altitude grid. Built-in options are MSISDensity and VectorDensity (your own altitude/density vectors, that will get pchip-interpolated in log-space).

For a custom analytic profile, wrap it in the @law macro. @law captures the law's source text so that it can be saved in physics_state.jld2. Bare anonymous functions are rejected as their source cannot be reconstructed when a saved model is reloaded.

# An analytic profile for N₂ (altitude in m → density in m⁻³)
model.species[:N2].density_profile = @law h -> 1e18 .* exp.(-(h .- 100e3) ./ 30e3)
initialize!(model)
model.species[:N2].density

368-element Vector{Float64}:
 1.0e18
 9.949876347809457e17
 9.899756425948795e17
 9.849642714327977e17
 9.799537679827698e17
 9.74944377611779e17
 9.69936344347815e17
 9.649299108622198e17
 9.599253184522903e17
 9.549228070241349e17
 ⋮
 2.969908081046221e15
 2.6998216467347665e15
 2.4476211445806335e15
 2.2127485658427435e15
 1.994615322253459e15
 1.7926051157286208e15
 1.6060770665088785e15
 1.434369072994877e15
 1.276801372059796e15

# Your own measured/downloaded profile, given as vectors:
model.species[:N2].density_profile = VectorDensity(altitude_m, density_m3)

# Or the default MSIS-backed profile, found/created by find_msis_file():
model.species[:N2].density_profile = MSISDensity(msis_file, :N2)

When a law needs parameters, use a functor

@law is for self-contained, closed-form laws: the expression may reference only its arguments and standard functions. A @law that closes over a local variable will be rejected (because the captured value can't be saved as source). When a law must carry parameters or tabulated data, define a small functor struct instead. Its fields are saved and can be reloaded with the model. Here is an example:

struct ExponentialDensity      # parameters live in fields → can round-trip through JLD2
    n0::Float64                # density at z_ref (m⁻³)
    z_ref::Float64             # reference altitude (m)
    H::Float64                 # scale height (m)
end
(d::ExponentialDensity)(z) = d.n0 .* exp.(-(z .- d.z_ref) ./ d.H)

model.species[:N2].density_profile = ExponentialDensity(1e18, 100e3, 30e3)

The same pattern applies to a parameterized cascading law or phase-function generator. (If you dig into the source code you might find out that this is how the default atomic-oxygen cascading law is built in order to carry some tabulated coefficients)

Overriding cross-sections, phase functions, or cascading

The same interception window lets you replace other per-species physics before initialize!. Set the generator (re-evaluated on every initialize!, so it tracks grid changes) rather than the materialized array where possible:

# Phase-function generator: (θ, E) -> (elastic, inelastic) matrices.
# Can be a named function, functor, or @law, but not an anonymous function.
model.species[:N2].phase_fcn_generator = my_custom_phase_function

# Cascading: ionization thresholds + a secondary-electron distribution law f(E_s, E_p)
law = @law (E_s, E_p) -> 1/(15.2^2 + E_s^2)
model.species[:O2].cascading_spec = AURORA.CascadingSpec("O2", [12.07, 16.1], law)

# Cross-sections can be pre-populated directly for a non-standard species (see below)

Adding, removing, or replacing species

Pass an explicit species tuple to the constructor. The defaults are N2Species/O2Species/OSpecies, which are helper functions accepting a density source (i.e. msis_file).

# Two species only:
model = AuroraModel(alt_lims, θ_lims, E_max, msis_file, iri_file;
                    species = (O2Species(msis_file), OSpecies(msis_file)))

A completely custom species needs its cascading law and a phase-function generator. Because the built-in cross-section library only knows N₂/O₂/O, pre-populate the cross-sections and excitation levels for a new gas in the interception window:

law  = @law (E_s, E_p) -> 1.0 / (12.0^2 + E_s^2)  # we are completely inventing here
spec = AURORA.CascadingSpec("Ar", [15.76, 27.63], law)
argon = NeutralSpecies(:Ar, MSISDensity(msis_file, :Ar);
                       cascading_spec = spec, phase_fcn_generator = phase_fcn_N2)

model = AuroraModel(alt_lims, θ_lims, E_max, msis_file, iri_file;
                    species = (N2Species(msis_file), O2Species(msis_file),
                               OSpecies(msis_file), argon))

model.species[:Ar].cross_sections    = my_sigma_matrix   # [n_levels × n_E]
model.species[:Ar].excitation_levels = my_levels_matrix  # [n_levels × 2]

run!(AuroraSimulation(model, flux, savedir; mode))