NRLMSISE-00, as a crate

Status: planned, not started Created: 2026-10-01 Why a separate crate: it is a self-contained empirical model with its own validation story, and nothing about it is specific to the library that will consume it.

What it is for

An atmospheric density model generic over its scalar type, so that a propagation using it can be differentiated as well as solved.

That is the whole reason to write one. Two Rust implementations already exist and are usable:

Crate	License	Validated against
`tobari`	MIT	`pymsis`, the official NRL Fortran wrapper
`satkit`	MIT OR Apache-2.0	GMAT drag cases

Both are f64. Depending on either would mean density could be computed but not differentiated, and that is not a small loss: density falls by a factor of e about every 50 km, so dropping its dependence on position makes a drag partial badly wrong rather than slightly wrong.

So the single acceptance criterion that justifies the work: the derivative of density with respect to altitude, taken through the model by a dual number, agrees with central differences and with -1/H from the scale height.

The validation target already exists

astrokit/tests/fixtures/nrlmsise00-cases.json — the seventeen standard test cases, inputs and outputs at seventeen significant digits. Copy it over.

Regenerable by astrokit/scripts/generate-msis-cases.sh, which fetches Dominik Brodowski's C translation of the NRL Fortran, compiles it with a driver that prints more digits than its own test does, runs it, and discards the source. What is committed is numbers, which are facts about the model rather than anyone's code. The reference's SHA-256 is in the fixture header.

The cases vary one input at a time from a baseline of day 172, 29000 s UT, 400 km, 60 N, 70 W, local solar time 16 h, F10.7 and its 81-day average 150, ap 4. Between them they reach 0 to 1000 km in altitude, the equator, local midnight, low and high flux, a storm, and the two cases that use the seven-element ap history rather than a daily value.

Make these numbers. Then widen with a grid.

Order of work

Coefficient tables, generated rather than transcribed. nrlmsise-00_data.c holds PT, PD (9 rows), PS, PDL (2), PTM, PDM (8), PTL (4), PMA (10), SAM and PAVGM — about 2700 numbers. Write a script that emits them as Rust statics, in the manner of astrokit/scripts/generate-gravity-field.py. Do it alongside the structure, since the shapes have to match what reads them. Transcribing by hand is the classic way to lose a week.
The small helpers, each with its own test. spline, splint, splini, densu, densm, ccor, ccor2, scalh, dnet, glatf, zeta. These are where a port quietly goes wrong and they are individually checkable, which the big ones are not.
globe7 and glob7s — the expansions in latitude, day of year, local time and flux. The largest pieces, and the ones the switches act on.
gts7, then gtd7 and gtd7d.
The seventeen cases.
A grid fixture: altitude by latitude by local time by flux by ap, from the same script, to catch what the seventeen happen not to exercise.

What genericity needs

Everything the model computes with is available on a scalar trait of the usual shape — exp, log, sqrt, tanh, cos, sin, powi, atan. Three rules, all learned the hard way in astrokit:

Every magnitude comparison branches on the real part. There are many: the altitude bands, the exponent clamps in ccor and dnet, the g0 and sumex guards. A dual number cannot be compared, and the branch is correct on the real part because the quantity being tested is a real one.
Coefficient tables and switch flags stay f64. T op f64 works, so lifting them would only cost memory and cache.
f64 on the left of an operator is not implemented. T * f64 is; f64 * T is not, and cannot be, because of the orphan rule. This has bitten twice in astrokit.

Traps, roughly by how expensive they are to find late

A consumer computing drag wants gtd7d, not gtd7. They differ only in the total mass density: gtd7d includes anomalous oxygen and gtd7 does not, which matters above about 500 km. The fixture gives all nine number densities, so both totals can be checked from it — gtd7d is 1.66e-24 * (4 He + 16 O + 28 N2 + 32 O2 + 40 Ar + H + 14 N + 16 anomalous O). Expose both and make the difference obvious in the names and the docs.
Fortran is 1-indexed and the C translation already shifted the arrays. Pick one source and stay with it. Mixing them is how an off-by-one gets in.
Units. The model returns cm^-3 and g/cm^3 unless its metric switch is set. Decide once whether the crate speaks SI at its boundary — it should — and convert in exactly one place. The fixture is in the model's own units and says so in its header.
switches[9] = -1 selects the ap-history mode, which is what a real propagation uses and what cases 15 and 16 exercise. The daily-ap mode is the simpler default and is not the one that matters.
Local solar time has to agree with the seconds of day and the longitude, as lst = sec/3600 + lon/15. The model does not check and gives quietly wrong answers when they disagree. Derive it inside rather than taking it as an input.
The 23 switches are a real part of the interface, not a detail: they turn individual terms on and off, and a -1 means something different from a 0. Model them as a type rather than an array of floats if it can be done without cost.

Checks beyond matching the reference

Density falls monotonically with altitude from 100 to 1000 km, and its logarithmic slope is a scale height in the tens of kilometers.
Day is denser than night at the same altitude; a storm is denser than quiet; high flux is denser than low.
It differentiates. The acceptance criterion above. Without this the crate has no reason to exist.

What the consumer already provides

astrokit::atmosphere::drivers reads a committed record of F10.7 and ap, and gives the previous day's observed flux, the 81-day centered average and the seven-element ap history in the order the model wants them. So the crate's interface should take those as plain numbers and own no data of its own.

One consequence worth knowing: the centered 81-day average needs forty days either side, so the usable end of any driver record is forty days before its last measurement.

Naming

nrlmsise00 is free on crates.io as of 2026-10-01 — confirm before publishing. MSIS is a registered trademark; tobari's module doc handles this well, claiming nominative fair use for identifying compatibility. The model itself is the work of Picone, Hedin and Drob at the US Naval Research Laboratory, believed to be in the public domain as a US Government work under 17 U.S.C. 105, though NRL attached no explicit license. Say so.

6.9 KiB Raw Permalink Blame History