17 KiB
nrlmsise00
The NRLMSISE-00 empirical atmosphere model, generic over its scalar type.
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 byScalarcarries them.DualNum,DualStructandSignedare re-exported from thescalarmodule 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 onf64when 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)), wherew * xis shared across the whole evaluation andpis one of the coefficient table's own numbers. Writtencos(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
powis five calls.powis 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.