cheatah
Source BigBrain
Class

space::irbem::Igrf

cheatah-space v0.1.0-alpha — Biome Standard 0.6.5-alpha

The IGRF internal geomagnetic field at one instant, ready to evaluate.

The date is bound once, at construction, by at — which is where the two things that can fail live: an out-of-range year, and the choice between interpolating between tabulated epochs and extrapolating with the published secular variation. evaluate is therefore total, branch-free and allocation-free, and an ephemeris of a million points pays the coefficient interpolation once instead of a million times. at is constexpr, so a model for a fixed epoch costs nothing at run time at all.

Template parameters
NMAX

the truncation degree, 1..13. 13 is IGRF-14 in full; 10 is what IRBEM uses internally. See the file-level note on degree.

P

the precision policy (SoundPrecision). P::integrand is what the Legendre and longitude recursions are carried in; P::accum — which the concept pins to double — is what the ~200 harmonic terms are summed in. That split is exactly the one docs/ERROR_BUDGET.md §3 argues for: one fp32 evaluation costs ~6e-8 and is free against the discretization floor, but a naive fp32 sum of 10³ terms costs ~1.2e-4 and lands on it. Roundoff in an evaluation does not grow with the size of the problem; roundoff in a reduction does.

Functions

fn static std::optional< Igrf > at(double decimal_year) #

The model at a decimal year.

Between tabulated epochs the coefficients are linear in time — the IGRF's own defined interpolation, not an approximation to something smoother. From latest_epoch_year to latest_year they are the last main field plus the published secular variation times the elapsed years. At a tabulated epoch the result is bit-identical to the table.

Parameters
decimal_year

the epoch, e.g. 2007.5 for mid-2007.

Returns

the model, or std::nullopt when decimal_year lies outside [earliest_year, latest_year] or is NaN. The empty case is never quietly clamped: a 2035 ephemeris is a question this model cannot answer, and saying so is the answer.

Complexity

O(NMAX²) — one pass over the 105 coefficient slots.

Allocation

none.

fn double year() const #

The epoch this model was built for.

Returns

the decimal year passed to at.

Complexity

O(1).

Allocation

none.

fn integrand g(int n, int m) const #

A time-interpolated g Gauss coefficient.

Parameters
n

the degree.

m

the order.

Returns

gⁿₘ in nT, or exactly zero when (n, m) is outside this truncation — the coefficients above degree, and every n < 1, are structurally absent from the model rather than merely unknown.

Complexity

O(1).

Allocation

none.

fn integrand h(int n, int m) const #

A time-interpolated h Gauss coefficient.

Parameters
n

the degree.

m

the order.

Returns

hⁿₘ in nT, or exactly zero outside the truncation. hⁿ₀ is zero by definition — sin 0φ is identically zero, so the coefficient has no meaning.

Complexity

O(1).

Allocation

none.

fn evaluate · 2 overloads
FieldVector< Frame::GEO > evaluate(const Position< Frame::GEO > &p) const#
FieldVector< Frame::GEO > evaluate(const Position< Frame::SPH > &p) const#

The internal field at a geocentric Cartesian position — the hot kernel.

Writing u = z/r, x̂ = x/r, ŷ = y/r, Rₙ = r^-(n+2), Gⁿₘ = gⁿₘcₘ + hⁿₘsₘ, and summing over the model, T0 = Σ (n+1)Rₙ Gⁿₘ Aⁿₘ T1 = Σ Rₙ Gⁿₘ A'ⁿₘ T2 = Σ m Rₙ Gⁿₘ Aⁿₘ U = Σ m Rₙ Aⁿₘ (gⁿₘcₘ₋₁ + hⁿₘsₘ₋₁) V = Σ m Rₙ Aⁿₘ (hⁿₘcₘ₋₁ - gⁿₘsₘ₋₁)

the Cartesian field is exactly Bx = x̂(T0 + u·T1 + T2) - U By = ŷ(T0 + u·T1 + T2) - V Bz = u·T0 - (1-u²)T1 + u·T2

The 1/sin θ of the eastward component and the sinᵐθ of the Legendre functions have both cancelled analytically — cₘcos φ + sₘsin φ = sin θ·cₘ₋₁ is the identity that does it — so the pole needs no test and no epsilon. The m = 0 term reads a padded zero slot for cₘ₋₁ and is multiplied by m, so it contributes exactly nothing without a branch.

Parameters
p

the position, geocentric Cartesian, in units of reference_radius_km. Inside the Earth is arithmetically fine and physically meaningless — the series diverges below the core; callers that care enforce r >= 1 themselves. r == 0 is the one input with no answer and yields infinities rather than a silent number.

Returns

the field in nT, in the same geographic frame.

Complexity

O(NMAX²) — 105 coefficient slots at degree 13, ~1900 flops, one sqrt, no transcendentals, no branches in the summation.

Allocation

none — three fixed std::arrays on the stack, 2 KiB at degree 13.

fn constexpr Igrf()=default #

Only at builds one, because only at can decide whether the date is answerable.

Constants & variables

var static std::size_t kSlots #

Flat coefficient slots at this truncation; the packing matches tables::igrf14_g, so truncating the model is a prefix copy rather than a re-index.

var static auto kNorm #

The compile-time normalisation table — 224 square roots that never run.

var static int degree #

The truncation degree this instantiation evaluates.

var static double reference_radius_km #

The IGRF reference radius a, in kilometres — and hence the Earth radius that Position's Cartesian components are counted in, throughout this module.

var static double earliest_year #

The first tabulated epoch. Before this the model is not defined and at reports so.

var static double latest_epoch_year #

The last tabulated main-field epoch; beyond it the published secular variation is used.

var static double latest_year #

The end of the secular-variation prediction.

Past this the model is not defined; at reports so rather than extrapolating a five-year linear prediction indefinitely.

var std::array< integrand, kSlots > g_ #

g coefficients in nT at year_, triangular packing.

var std::array< integrand, kSlots > h_ #

h coefficients in nT at year_, triangular packing.

var double year_ #

The epoch the coefficients were interpolated to.

Types

type typename P::integrand integrand #

The type the recursions and the stored coefficients are carried in.

type typename P::accum accum #

The type the harmonic sums are accumulated in — double, by SoundPrecision.