Atmospheric drag and NRLMSISE-00 #11
Assignees
connor (Connor Johnstone)
Clear assignees
connor (Connor Johnstone)
No Assignees
Notifications
Due Date
No due date set.
Dependencies
No dependencies set.
Reference: connor/astrokit#11
Reference in new issue
Block a user
Drag behind a density model, with the models themselves split by where they live.
Harris-Priester lives here. A fixed pair of profiles with a diurnal bulge, reading no space weather at all. Cheap, smooth, and wrong in a predictable direction, which is what makes it a useful partner to a model that does read the weather: the two give two density realizations of the same trajectory.
NRLMSISE-00 is a separate crate, by the owner's decision, built and validated on its own. Both existing Rust implementations (
tobari,satkit, both MIT-compatible) aref64-only, and density falls by a factor of e every 50 km, so depending on one would make a drag partial badly wrong rather than slightly wrong. The plan for that port is at~/docs/projects/nrlmsise00/plan.md; the validation target is committed here astests/fixtures/nrlmsise00-cases.jsonand regenerable byscripts/generate-msis-cases.sh. This crate will take it as an optional dependency behind a feature, as a variant of the density enum.What is left in this issue: the
DensityModelenum, theDragdynamics with its force and its relative-wind handling, and refusing a propagation whose arc reaches past the driver record.M1.
Plan for porting NRLMSISE-00, generic over the scalar
Decided after checking what exists:
tobari(MIT, clean-room, validated againstpymsis)and
satkit(MIT/Apache) both implement it, and both aref64-only. Depending on eitherwould mean a propagation with drag could be solved but not differentiated, and that is not a
small loss — density changes by a factor of e every ~50 km, so dropping its dependence on
position makes a drag partial badly wrong rather than slightly wrong.
The validation target exists already
tests/fixtures/nrlmsise00-cases.jsonis committed: the seventeen standard test cases, withinputs and outputs at seventeen significant digits, regenerable by
scripts/generate-msis-cases.sh. That script fetches Brodowski's C translation of the NRLFortran, compiles it with a driver that prints more digits than its own test does, runs it
and throws the source away. What is committed is numbers, which are facts about the model.
The reference's digest 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 — reaching 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 the daily value.
So the port has a target before a line of it is written. Make these numbers.
Order of work, bottom up
nrlmsise-00_data.cholds PT,PD (9 rows), PS, PDL (2), PTM, PDM (8), PTL (4), PMA (10), SAM and PAVGM — about 2700
numbers. A script in the manner of
scripts/generate-gravity-field.pyshould emit themas Rust statics. Do this in the same sitting as the structure, since the shapes have to
match what reads them.
spline,splint,splini,densu,densm,ccor,ccor2,scalh,dnet,glatf,zeta. These are where a portquietly goes wrong, and they are individually checkable.
globe7andglob7s, the expansions in latitude, day of year, local time and flux.The largest pieces and the ones the switches act on.
gts7, thengtd7andgtd7d.the same script, to catch anything the seventeen happen not to exercise.
What genericity needs
Everything the model computes with is available on
T: Scalar—exp,log,sqrt,tanh,cos,sin,powi,atan. Three rules:.re(). There are many: altitude bands, theexponent clamps in
ccoranddnet, theg0/sumexguards.f64.T op f64already works, so lifting them would onlycost memory. Same for the 24 switch flags.
f64on the left of an operator. It is not implemented;T * f64is. This hasbitten twice already, in
third_bodyand ingravity.Traps, in rough order of how expensive they are to find late
gtd7dfor drag, notgtd7. They differ only in the total mass density:gtd7dincludes anomalous oxygen and
gtd7does not. Above about 500 km that matters. Thefixture gives all nine number densities, so both totals can be checked from it —
gtd7dis1.66e-24 * (4 He + 16 O + 28 N2 + 32 O2 + 40 Ar + H + 14 N + 16 anomalous O).from the Fortran would have to. Pick one source and stay with it.
is SI, so convert at the boundary and nowhere else. The fixture is in the model's own
units and says so.
switches[9] = -1is the ap-history mode, which is what a real propagation uses andwhat cases 15 and 16 exercise. The daily-ap mode is the simpler default and is not the
one that matters.
(
lst = sec/3600 + lon/15). The model does not check, and gives quietly wrong answers ifthey disagree. Derive it rather than taking it as an input.
atmosphere::driversgives the previous day's observedflux, the 81-day centered average and the seven-element ap history in the order the model
wants. Note that the centered average costs forty days off the usable end of the record.
Checks beyond matching the reference
slope is the scale height, which should come out in the tens of kilometers.
differences, and against
-1/Hfrom the scale height. This is the whole reason forporting it rather than depending on one of the existing crates, so it is the test that
justifies the work.
Size: large. L is right.
The plan comment above is superseded: MSIS is a separate crate and the plan now lives at
~/docs/projects/nrlmsise00/plan.md, which is fuller. Left here for the history rather than deleted.Landed so far:
atmosphere::drivers— the committed F10.7 and ap record, parsed on first use, with the previous day's flux, the 81-day centered average and the seven-element ap history. Ten tests. Note the centered average costs forty days off the usable end of the record.atmosphere::harris_priester— the fifty-row profile, exponential in altitude, with the diurnal bulge thirty degrees east of the Sun. Nine tests, including that it differentiates and that the implied scale height comes out in the tens of kilometers.One property worth knowing about Harris-Priester: the interpolation is exponential inside a tabulated interval and only continuous across one, so the derivative steps at each of the fifty tabulated heights. That is the published form and smoothing it would make it a different model, but it means a propagation crossing those heights meets a small corner each time. There is a test recording it. Found because an autodiff test disagreed with central differences by 4%, which turned out to be my sampling two altitudes that are exactly table nodes.
Drag is in, with two density models. Ready for review; MSIS still to come as its own crate.
Dynamics::Drag { density }, withDensityModelan enum for the same reasonSurfaceModelis one: the density has to come out in the propagated scalar, so the method is generic and a trait with a generic method is not object-safe. An external MSIS crate joins as a variant behind a feature.The relative wind comes out of the frame rotation, with no wind model and no angular rate to supply:
Cartesian::to_frame'sR v + R_dot ris the velocity relative to the co-rotating air, the second term being the air's own motion. That is not a small correction — at the equator the air is already doing 465 m/s, six percent of orbital speed, and it enters the force quadratically. There is a test that pins it down by asking which of the two speeds explains the force: the readings differ by twelve percent, and only the relative one gives a density that is a thermosphere.Exponential, one layer, added on the owner's suggestion for something easily differentiable. Harris-Priester's interpolation is exponential inside a tabulated interval and only continuous across one, so its slope steps fifty times between 100 and 1000 km; the single layer has no corners at all, which makes it the model to reach for when a derivative is the point.Exponential::fitted_to(altitude)takes both the density and the scale height from the Harris-Priester mean there, so no second table had to be sourced.Fitted at an orbit's own altitude the two agree to a few tens of meters over a day against kilometers of drag; fitted elsewhere they diverge, which is the knob for running one trajectory under two density realizations.
In
examples/full.rs, with each force now priced by leaving it out rather than by nested slices — which had silently stopped meaning what their labels said once drag was added. Over a day at 400 km:That is the textbook hierarchy for the altitude, and the Sun and Moon at 35 m and 6 m agree with the 45 m the two together gave in #9. Drag cost 7 steps a day out of 1420, so the Harris-Priester corners turn out not to matter.
Left: refusing an arc that reaches past the driver record, which nothing yet needs —
DensityModel::needs_driversis false for both of these and the hook is there for when a model returns true.The driver coverage check is wired in, so #11 is complete apart from NRLMSISE-00 arriving as its own crate.
Drivers::covers_arc(start, end)examines the two ends, which is enough for a record that is one contiguous run of days, and does not care which end is which because a propagation may run backwards.Propagation::runcalls it for any spec whose density model reads space weather, before integrating, and the error names the day the record reaches and the script that refreshes it.Neither model here reads the weather, so the branch cannot fire yet. The logic it guards is tested directly: an arc inside the record, backwards, ending exactly on the last usable day, reaching past it, starting before it, and — the one that catches people — ending on the last measured day, which is already outside, the centered 81-day flux average costing forty days off the end.
211 unit and 6 doc tests, clippy, docs, benchmarks, all three examples and the hook clean.
The part of this that lives here is done, in
b60dbac. Leaving open for NRLMSISE-00.What landed:
Dynamics::Drag { density }andatmosphere::DensityModelas an enum, for the sameobject-safety reason as
SurfaceModel: the density has to come out in whatever scalar isbeing propagated, so the method is generic and a trait with a generic method is not
object-safe.
moves the orbit 15.0 km over a day against 59.5 km of total drag.
Cartesian::to_frame'sR v + R_dot ris alreadythe velocity relative to the co-rotating air. At the equator the air does 465 m/s, six
percent of orbital speed, entering the force quadratically, so getting this wrong would
be a twelve percent error in the force rather than a rounding one. The test identifies
which speed the model used by asking which one implies a density that is a thermosphere.
generated into
src/atmosphere/space_weather.rsas Rust rather than parsed at runtime.GFZ over Celestrak because it states a license, and because Celestrak's file runs to 2041
on its own predictions, which are not measurements.
for the 81-day centered flux average costing forty days off the usable end — so an arc
ending on the last measured day is already outside.
Measured over a day at 400 km, each force priced by leaving it out: harmonics 1028.8 km,
drag 35.4 km, the Sun 0.035 km, the Moon 0.006 km, sunlight 0.006 km.
The Harris-Priester interpolation here is deliberately not the published one. The usual
presentation interpolates exponentially within each tabulated interval, which is continuous
but steps the altitude derivative at every one of the fifty tabulated heights — a median of
4%, and 29% at 120 and 130 km. That cannot be differentiated across, which matters because
the drag partial is wanted exactly where a spacecraft crosses those heights. This
implementation fits a monotone cubic to the logarithm of the density instead. The table is
still reproduced exactly and the profile still falls monotonically. Measured deviation from
the published form: 0.27% over 200-300 km, 0.25% over 300-500, 0.17% over 500-800, 0.47%
over 800-1000, and 13% over 100-130 km where the table is spaced 10 and 20 km and nothing
orbits. It is documented on the module under its own heading so nobody mistakes it for
published Harris-Priester.
Per call it is also cheaper than the published form, 16.2 ns against 19.9 ns, the cubic
beating two
powf. End to end that is invisible: at degree 70 the gravity field is 6.76 usper derivative evaluation and the density is 0.2% of it.
An
Exponentialsingle-layer model was added alongside, on the reasoning thatHarris-Priester had corners. Once it did not, that reason was gone and the model was
removed.
What is left, and why this stays open: NRLMSISE-00. That is now entirely the separate
crate's responsibility, and the validation fixture and its generator have moved out of this
repository into that project along with the plan. This crate's remaining part is small —
take it as an optional dependency behind a feature and add a
DensityModelvariant. Thehook for that already exists:
DensityModel::needs_driversand the arc-coverage check arewritten and tested, and simply return false for every model here, because neither reads
space weather.
Note the issue body above is now stale where it says the fixture is committed here.
The NRLMSISE-00 half is done, in
74595fa, which completes this issue. The drag modelitself landed earlier in
b60dbac, described in the comment above.It lives in its own crate, as decided:
nrlmsise00, now 0.2.0 on crates.io, genericover its scalar so density can be differentiated. This crate takes it as an optional
dependency behind a
msisfeature and adds aDensityModel::Nrlmsise00variant. Thefeature turns on that crate's
fastpath, which trades the reference's exact bits forroughly half the time: over its 2880-case grid the largest value it moves moves by 3.9e-14
of itself, against a model whose own agreement with the atmosphere is tens of percent.
src/atmosphere/msis.rsis only the join — kilometers and degrees and a Modified JulianDay on that side, meters and radians and an
Epochon this one — with the model behind aOnceLockso its coefficient tables are built once rather than per derivative evaluation.Two things that changed here as a consequence
This crate's own driver record is gone.
drivers.rs, the 1.4 MB generatedspace_weather.rsandscripts/fetch-space-weather.pyare deleted, and the arc check inPropagation::runnow asks the atmosphere crate's record instead. The two were the sameGFZ download a day apart and that one is a superset — it has the centered flux average,
the same Modified Julian Day indexing and a complete driver set per day. With the feature
off nothing reads space weather, so nothing is lost.
Conditionscarries a wholeGeodeticrather than a bare altitude, because this modelwants latitude and longitude and recomputing them would mean a second iteration of the
geodetic conversion. It dropped its separate
bodyfield, whichGeodeticalready had.The part that was not obvious
This crate depends on
ordinary-diffeqand onnrlmsise00, and each defined its ownScalartrait — structurally near-identical, implemented for exactly the same types, andmutually unsatisfiable, so a caller generic over one could not call code generic over the
other.
The resolution was that
num-dualalready carries the shared vocabulary.DualNum<Primitive = f64>gives the arithmetic, the transcendentals andFrom<f64>;DualStruct<Real = f64>givesre(); and both are implemented for plainf64as well asevery dual type. So both crates now inherit those rather than redeclaring them —
ordinary-diffeq0.8.0 andnrlmsise000.2.0 — and a type satisfying one satisfies theother. This crate needed no source change for it:
reis inherited rather than redeclared,so every call site kept compiling, and
liftstays as a provided method sincenum-dualhas no such name.
Worth recording for anyone tempted by the alternative: a newtype adapter, or a fourth crate
defining the shared trait. One was written before the check that
num-dualalready had it.Measured, in this crate rather than quoted
fastConditions::new, the geodetic conversionThe whole drag path is therefore about 79 ns with the fixed profile and 1.0 us with
NRLMSISE-00, against a degree-70 gravity field at 6.76 us per derivative evaluation.
Profiling this is what produced
4ad8bc4: the geodetic conversion was 181 ns, eight timesthe density model it exists to feed, because its iteration carried the latitude and so paid
an
atan2per pass that the next pass undid with asin. Carrying the sine and cosineinstead took it to 66.8 ns with values unchanged to 1.1e-16 rad.
Validated: a 5.7x density swing at 400 km between the 2009 solar minimum and the 2014
maximum, less sensitivity at 200 km, three orders of magnitude from 200 to 800 km, and
Harris-Priester landing within a factor of a few at solar maximum, which is the activity it
is fitted for. An arc reaching past the measured record is refused before integrating.
Known limitations
above a fitted height by a test on the altitude alone, so the density steps: eight places
between 72.5 and 1000 km, none above 500. In the band that matters here the steps are
5.3e-5 of the total at 300 km, 3.5e-5 at 450, 2.3e-6 at 240 — far below anything a step
controller notices, but genuine discontinuities rather than corners, and documented on
the variant. Smoothing them would make it a different model.
exactly the differentiability reason. See the earlier comment.
repository; it moved to the atmosphere crate along with its generator.