ADCS.environment.igrf module

IGRF geomagnetic field model.

This is the default "igrf_numba" model of ADCS.environment.magnetic_field; ppigrf remains selectable as "igrf_ppigrf". It evaluates the same model from the same IAGA coefficients and agrees with ppigrf to ~1e-12 relative, but it is roughly 50-70x faster in batch and ~3000x faster per call (debug/debug_orbit/benchmark_igrf.py), which matters because ADCS.Orbital_State.get_b_eci() is on the simulation’s inner loop.

Why it was worth writing rather than tuning ppigrf

ppigrf is correct but its cost structure is wrong for this use:

  • igrf_gc calls read_shc on every invocation – reopening the .shc file and rebuilding two pandas DataFrames. There is no caching anywhere in the library. Evaluating one point therefore cost ~31 ms, essentially all of it re-reading a file that never changes.

  • Dates are deliberately kept out of broadcasting, so N coordinates and N dates produce an N x N result. Passing a trajectory produced N^2 field evaluations of which only the diagonal was kept.

  • The Legendre routine allocates (nmax+1)^2 arrays through a dict, then hstack``es 104 columns, and the field sum is built from several (N, 208) temporaries including a ``** over the whole array.

Here the coefficient file is parsed once and cached, the Legendre recursion keeps only the two previous degrees (so the per-sample scratch is 6 rows of 14 doubles rather than a 14x14 table), and the whole evaluation is a single numba serial kernel. numba is already a hard dependency of this package, so this adds no new install burden – unlike ESA’s eoxmagmod, which is not on PyPI and needs the NASA CDF library plus a Fortran build.

Coefficients

IGRF14.shc is the 14th-generation IGRF from IAGA Working Group V-MOD (https://doi.org/10.5281/zenodo.14012302), in the .shc layout used by ppigrf (MIT, Copyright (c) 2021 Karl M. Laundal). Valid 1900.0-2030.0.

ADCS.environment.igrf.coefficients_at(date)[source]

Gauss coefficients interpolated to a single date.

Parameters:

date (datetime.datetime) – Date at which to evaluate the model.

Returns:

(g, h), each of shape (NMAX + 1, NMAX + 1), Schmidt normalised.

Return type:

tuple

Example:

>>> from datetime import datetime
>>> g, h = coefficients_at(datetime(2024, 3, 1))
>>> g.shape
(14, 14)
ADCS.environment.igrf.igrf_gc(r, theta, phi, date, chunk=512)[source]

Geomagnetic field in geocentric spherical components.

Unlike ppigrf.igrf_gc, dates are broadcast elementwise against the coordinates rather than forming an outer product. Passing N positions and N dates returns N field vectors, not an N x N grid – which is what every caller in this package actually wanted, and what the np.diagonal workarounds in ADCS.orbits.orbit used to recover by hand.

Parameters:
  • r (float or numpy.ndarray) – Geocentric radius [km].

  • theta (float or numpy.ndarray) – Geocentric colatitude [deg].

  • phi (float or numpy.ndarray) – Geocentric longitude [deg], positive east.

  • date (datetime.datetime or sequence) – One date, or one date per sample.

  • chunk (int) – Samples per block of Legendre scratch reuse. Affects speed only, never the result.

Returns:

(Br, Btheta, Bphi) in nT, each with the broadcast shape of the coordinates. Btheta points south, Bphi east.

Return type:

tuple

Example:

>>> from datetime import datetime
>>> import numpy as np
>>> br, bt, bp = igrf_gc(6771.2, 45.0, 10.0, datetime(2024, 3, 1))
>>> bool(np.hypot(np.hypot(br, bt), bp) > 1e4)   # order 10^4 nT in LEO
True
ADCS.environment.igrf.NMAX = 13

Maximum spherical-harmonic degree of IGRF-14.

ADCS.environment.igrf.RE = 6371.2

Reference radius of the IGRF expansion [km]. This is the IGRF’s own constant, not an Earth radius from ADCS.helpers – changing it would silently rescale the field.