Evaluated at f64 it gives density and temperature. Evaluated at a forward-mode dual number it gives density and its exact derivatives, which is the reason this crate exists: density falls by a factor of e about every fifty kilometers, so a drag partial that drops its dependence on position is badly wrong rather than slightly wrong. Two other Rust implementations exist and both are f64 only.

use nrlmsise00::{Input, MagneticActivity, Model};

let model = Model::new();
// The scalar is named, because a call built only from literals has nothing else to fix it.
let output = model.gtd7d::<f64>(&Input::consistent(
    172,                                   // day of year
    29000.0,                               // seconds of universal time
    400.0,                                 // altitude, km
    60.0,                                  // latitude, degrees north
    -70.0,                                 // longitude, degrees east
    150.0,                                 // 81-day average F10.7, centered
    150.0,                                 // previous day's F10.7
    MagneticActivity::Daily { ap: 4.0 },
));

// SI at the boundary: kg/m^3 and m^-3. Local solar time comes out at 3.4 h here, since
// `consistent` derives it from the universal time and the longitude.
assert!((output.total_mass - 2.400e-12).abs() < 0.001e-12);
assert!((output.temperature - 1098.2).abs() < 0.1);

And the derivative, through the model rather than around it:

use num_dual::{Dual64, DualNum};
use nrlmsise00::{Input, MagneticActivity, Model, Scalar};

let model = Model::new();
let altitude = Dual64::from(400.0).derivative();
let output = model.gtd7d(&Input::consistent(
    172,
    Dual64::lift(29000.0),
    altitude,
    Dual64::lift(60.0),
    Dual64::lift(-70.0),
    Dual64::lift(150.0),
    Dual64::lift(150.0),
    MagneticActivity::Daily { ap: Dual64::lift(4.0) },
));

// d(rho)/dz, in kg/m^3 per km, and the scale height it implies.
let scale_height = -output.total_mass.re / output.total_mass.eps;
assert!((scale_height - 55.8).abs() < 0.1);

cargo run --example density prints a vertical profile, the same altitude across two solar cycles, and the altitude derivative a drag calculation needs, checked against a central difference.

The scalar type is `num-dual`'s

Scalar is a marker over [num_dual::DualNum] with one method of its own, lift, and it is blanket-implemented:

pub trait Scalar: DualNum<Primitive = f64> + DualStruct<Real = f64> + Copy {
    fn lift(x: f64) -> Self { Self::from(x) }
}

impl<T: DualNum<Primitive = f64> + DualStruct<Real = f64> + Copy> Scalar for T {}

So every type that implements num-dual's trait at f64 precision is a scalar here — f64, each of its dual numbers, nestings of them, and anything a caller writes. There is nothing to implement and nothing to enumerate, which is the point: an earlier version declared exp, ln, powf and the rest itself and implemented them for a fixed list of types, so a caller whose scalar was a dual type the list had missed was stopped by the orphan rule with no way out, and a caller generic over another crate's scalar trait of the same shape was adapting between two names for one thing.

Two things follow for a caller:

Generic code imports only Scalar. The other methods are its supertraits', and a type parameter bounded by Scalar carries them. DualNum, DualStruct and Signed are re-exported from the scalar module for code that calls them on a concrete dual number.
A call built only from float literals has to name its scalar, as model.gtd7d::<f64>(..). Inference will not settle a literal on f64 when the blanket bound admits every dual number as well. Code that gets its position from somewhere typed never meets this.

num-dual is therefore not optional, and this crate takes it with default-features = false: that drops nalgebra, which only the vector dual types need and which this crate never names. A caller who wants DualSVec64 enables it from their own manifest, and the blanket implementation covers those types without this crate knowing about them.

`gtd7` or `gtd7d`

They differ in one number. gtd7's total mass density sums helium, atomic oxygen, molecular nitrogen, molecular oxygen, argon, hydrogen and atomic nitrogen. gtd7d also includes anomalous oxygen, which is hot atomic oxygen and is 18.6% of the mass at 1000 km in the model's own standard case. A drag calculation wants gtd7d. The individual number densities are identical either way.

How closely it matches the reference

By default, exactly. Not approximately and not to a tolerance: every number this crate returns at f64 is the same bits the reference C returns, over every case there is to compare against.

What	Cases	Values compared
The standard test cases	17	187
A wide input grid, fifty altitudes against sixty scenarios under eight switch settings	2880	31680
The individual helpers, expansions and profiles	819	several thousand

The comparisons are on the bits rather than against a tolerance because an approximate one would not notice a term dropped from an expansion whose contribution happens to be small where it is sampled. The fixtures hold the reference's own output at seventeen significant digits and are regenerated by the scripts in scripts/, which fetch Dominik Brodowski's C translation of the NRL Fortran, run it, and discard the source: what is committed is numbers, which are facts about the model rather than anyone's code.

With the fast feature, no. That is what the feature is for, and the section on it below says what it costs.

At a dual type the value agrees to about 1e-14 rather than exactly, because num-dual computes a division as a reciprocal and a multiply. At f64 every method is a thin forward to the std function of the same name, which is why that path is the exact one. See the scalar module docs.

Where it deliberately differs from the reference

The reference keeps its working state in file-scope statics, so two of its answers depend on what was asked before them. Both are fixed here, and both are fixed by computing a value the reference reads without having written rather than by changing any arithmetic:

At exactly 32.5 km, it splines through temperature nodes its own call never filled.
With switch 15 off, below 72.5 km, it scales every species from a mixed nitrogen density its own call never computed.

These nodes and that density do not depend on altitude, so filling them reproduces what a preceding call would have left, which is why the grid fixture still matches at 32.5 km in sixty scenarios. Nothing else departs from the reference's control flow.

The drivers do almost nothing below 80 km

The model runs from the ground up — that is what the E in MSISE is, where MSIS-86 started at 85 km — but the solar and magnetic inputs barely reach the bottom of it, and the reference says so in its own input notes:

f107, f107A, and ap effects are neither large nor well established below 80 km and these parameters should be set to 150., 150., and 4. respectively.

Measured, against the drivers of the May 2024 Gannon storm, which was the largest geomagnetic storm in twenty years:

Altitude	Change in total mass density
0, 10, 20, 50 km	+0.770%, the same figure at every one
70 km	+0.753%
80 km	+5.2%
90 km	-15.9%
200 km	+71%
400 km	+225%

The figure below 62.5 km is identical at every altitude because of how the model is built there: every species is the fully mixed molecular nitrogen density times a fixed ratio, and that density is scaled from a single value computed at 72.5 km. So the drivers move the whole lower profile by one factor and change its shape not at all. Above 62.5 km the mixing ratio starts to vary and the figure drifts.

Below 80 km, then, pass 150, 150 and 4 as the reference asks — DriverSet::LOWER_ATMOSPHERE is exactly those three. Between 80 and 120 km is where the drivers do something the reference calls not well established, and the sign of it is not what you would guess: a severe storm lowers the density at 90 km by a sixth.

Model::gtd7d_at does not substitute them for you down there, deliberately. Switching drivers at a threshold would put a 5.2% step in the density at exactly 80 km, which is thirty-five times the largest discontinuity the model itself has and would be this crate's invention rather than the model's; and it would make the convenience entry point disagree with gtd7d given the same drivers. So it hands you the measurements at every altitude and says so, and below 80 km you pass the fixed set instead.

It is not continuous in altitude

Worth knowing before integrating through it. Each species' mixed correction is dropped above a fitted altitude by a test on the altitude alone, and the model's three vertical regions are joined by switches rather than smoothly, so the density steps. Sweeping 72.5 to 1000 km at ten meters and comparing each interval against the model's own derivative finds eight places where it does, and no others:

Altitude	Total mass density	What happens there
72.5 km	1.0e-3	the join between the model's two halves
123.435 km	1.5e-3	`za`, the join between the Bates profile and the spline below it
160 km	7.1e-5	molecular nitrogen's mixed correction is dropped
240 km	2.3e-6	argon's
250 km	1.9e-5	molecular oxygen's
300 km	5.3e-5	atomic oxygen's, and the lower thermosphere temperature variation
428.475 km	6.9e-7	atomic nitrogen's turbopause blend saturates, a step of about exp(-10)
450 km	3.5e-5	atomic nitrogen's mixed correction is dropped

Individual species move more than the total does — up to 3e-3 — because the species that steps is usually not the one carrying the mass. Above 500 km there are no jumps at all.

The temperature is continuous at all of these but za, where the exospheric temperature steps by nearly two fifths: 1027.3 K to 1426.4 K in the standard case. That is the reference's documented behavior rather than an artifact, its header saying the exospheric temperature "is set to global average for altitudes below 120 km", and it is what makes 123.435 km the largest discontinuity in the model above 72.5 km.

What is not implemented

ghp7, the reference's pressure-level entry point, which answers at a given pressure rather than a given altitude. It is a secant iteration wrapped around gtd7, so a caller who wants it can write it against this crate; nothing in the model needs it.

Speed, and the `fast` feature

	400 km, ap history	400 km, daily ap	20 km, ap history
default, bit-exact	1.44 us	1.09 us	3.52 us
with `fast`	0.83 us	0.68 us	2.37 us
the reference C at `-O2`, for scale	1.97 us		4.24 us

One core of an Intel Core Ultra 9 285K. Low altitude costs more because the mesosphere adds two more spline layers on top of a thermosphere evaluation.

Two thirds of the default build's time is inside the maths library: one evaluation at 400 km makes 274 calls into it, 119 of them cos. So what makes this model fast is making fewer of those calls rather than tightening the arithmetic between them, and fast takes that count from 274 to 98.

fast is not a faithful reproduction of NRLMSISE-00. It computes the same model by different arithmetic, in two places:

Phase-shifted cosines by the angle-sum identity. Almost every cosine in the expansions is cos(w * (x - p)), where w * x is shared across the whole evaluation and p is one of the coefficient table's own numbers. Written cos(wx)cos(wp) + sin(wx)sin(wp), with each phase's sine and cosine worked out once when the model is built, a hundred cosine calls become six. The identity is exact in real arithmetic and is two roundings rather than one in floating point.
Constant powers by repeated squaring. The magnetic activity function asks for the 3rd, 4th, 8th, 12th and 19th powers of one number. Seven multiplications give all five, where pow is five calls. pow is correctly rounded; a product of products is not.

What that costs, measured over the 2880-case grid rather than argued: 11778 of 31680 values move, and the largest moves by 3.9e-14 of itself, which is 303 units in the last place. For scale, the model's own agreement with the real atmosphere is a few tens of percent, and a dual-number evaluation of even the exact path differs from f64 by around 1e-14 for unrelated reasons. The test suite holds the default to the reference's exact bits and fast to 1e-13 relative, so both are checked on every run.

Take fast if you want this model and want it quickly. Leave it off if you want the reference's answer, which is why the crate exists.

Units and switches

The crate answers in SI by default — m^-3, kg/m^3 and K — using the model's own metric switch, so the conversion happens in one place and is the reference's own. Units::Native gives the cm^-3 and g/cm^3 the published cases are quoted in.

The twenty-three switches are named rather than indexed, and each has three states: the reference's 2 keeps a term's cross terms while dropping the term itself. The reference's switch 9 does two jobs, one of which is to select between a daily ap value and a seven-element history; that is a statement about the data, so it lives on the input as MagneticActivity and the mode cannot disagree with the values.

Drivers

The solar and geomagnetic values the model reads are measurements, and this crate does not own them: f107, f107_average and the ap history are plain numbers on the input, so a caller with its own record passes them straight in and pays nothing for the convenience it is not using.

The drivers feature, on by default, is for callers who would rather not. It compiles in GFZ Potsdam's record — 25475 days, 1957 to the present, 535 KB — and looks a date up in it:

# #[cfg(feature = "drivers")]
# fn main() {
use nrlmsise00::{drivers, Model};

let model = Model::new();
let when = drivers::mjd(2015, 3, 15, 8, 3, 20.0);   // or hifitime's to_mjd_utc_days()
let output = model
    .gtd7d_at::<f64>(when, 400.0, 60.0, -70.0)
    .expect("inside the record");

assert!((output.total_mass - 2.3e-12).abs() < 0.5e-12);
# }
# #[cfg(not(feature = "drivers"))]
# fn main() {}

That picks the previous day's flux, the 81-day centered average and the seven-element magnetic history, and derives the day of the year, the time of day and a consistent local solar time. Drivers::embedded() exposes the record itself if you want the parts.

Nothing is read or allocated until something asks for a day. Measured: a lookup is 70 ns against an evaluation's 1.76 us, the first call is 278 us to accumulate the flux averages, and a caller supplying its own drivers touches none of it. Dates are Modified Julian Days and the module needs no time library; drivers::mjd builds one from a civil date.

The record is measurements, so a date past it is refused rather than extrapolated, and the usable range is forty days shorter than the record at both ends, because the 81-day flux average is centered. Refreshing it is a deliberate act — scripts/fetch-drivers.py — not something that happens inside an evaluation: a program that fetched its own weather would answer differently on different days, and nothing here should acquire a network failure mode. drivers::provenance() returns the source, license, digest and coverage of whatever is compiled in.

The model is not this crate's work

NRLMSISE-00 is the work of Mike Picone, Alan Hedin and Doug Drob at the US Naval Research Laboratory, published as Picone et al. (2002), NRLMSISE-00 empirical model of the atmosphere: Statistical comparisons and scientific issues, J. Geophys. Res. 107(A12). It is believed to be in the public domain as a work of the US Government under 17 U.S.C. 105; NRL attached no explicit license to it. The coefficients compiled into this crate are the model's own.

MSIS is a registered trademark of the Government of the United States of America, as represented by the Secretary of the Navy. It is used here only to identify which model this implements, which is nominative fair use; this crate is not endorsed by or affiliated with NRL or the US Navy.

Later versions of the model, NRLMSIS 2.0 and 2.1, are not here and cannot be on these terms. NRL licenses them for academic and non-commercial use only; a translation into another language is a derivative work under those terms, which must be delivered to NRL and published as open source; and the methods are covered by US Patent 10,641,925, which runs to 2039. A permissively licensed port of 2.x is not something anyone can give you.

License

MIT, for this implementation. The full text is in LICENSE, along with the attribution and trademark notes above.

readme.md

nrlmsise00

The scalar type is num-dual's

gtd7 or gtd7d