Added an example

Cargo.toml

@@ -45,3 +45,10 @@ serde_json = { version = "1.0", features = ["arbitrary_precision"] }
 [[bench]]
 name = "density"
 harness = false
+# Measures the dual-number scalars and the driver lookup, so it needs both.
+required-features = ["autodiff", "drivers"]
+
+[[example]]
+name = "density"
+# Reads the compiled-in driver record and seeds a dual number, so it needs both.
+required-features = ["drivers", "autodiff"]

examples/density.rs

@@ -0,0 +1,156 @@
+//! The atmosphere a satellite flies through, and how fast it thins with altitude.
+//!
+//! Run with `cargo run --example density`.
+//!
+//! Three things, which between them are the whole of what this crate is for: a vertical
+//! profile of the model, the same altitude across the solar cycle from the compiled-in
+//! driver record, and the altitude derivative of the density — which comes out of the
+//! model itself rather than out of a finite difference, because the model is generic over
+//! its scalar type.
+
+use nrlmsise00::{drivers, Model, Output, Scalar};
+use num_dual::Dual64;
+
+/// The species of [`Output::densities`], in its order. Index five is the total mass rather
+/// than a number density, so it is not a species and is skipped.
+const SPECIES: [&str; 9] = ["He", "O", "N2", "O2", "Ar", "", "H", "N", "O*"];
+
+/// The species with the most atoms or molecules per cubic meter, and its share of them.
+fn dominant<T: Scalar>(output: &Output<T>) -> (&'static str, f64) {
+    let mut total = 0.0;
+    let mut best = ("", 0.0);
+    for (index, (density, name)) in output.densities().iter().zip(SPECIES).enumerate() {
+        if index == 5 {
+            continue;
+        }
+        let density = density.re();
+        total += density;
+        if density > best.1 {
+            best = (name, density);
+        }
+    }
+    (best.0, if total > 0.0 { best.1 / total } else { 0.0 })
+}
+
+fn main() {
+    let model = Model::new();
+
+    // The day of the Gannon storm, the largest geomagnetic storm in twenty years. The
+    // drivers for it come out of the record compiled into the crate, so the date is the
+    // only thing that has to be supplied.
+    let when = drivers::mjd(2024, 5, 11, 12, 0, 0.0);
+    let (latitude, longitude) = (40.0, -105.0);
+    let place = format!("{latitude:.0} N, {:.0} W", -longitude);
+    let set = drivers::Drivers::embedded()
+        .at(when)
+        .expect("2024 is inside the driver record");
+
+    println!("NRLMSISE-00 at {place}, 2024-05-11 12:00 UTC, during the Gannon storm");
+    println!(
+        "  F10.7 {:.1}, its 81-day average {:.1}, daily Ap {:.0}",
+        set.f107,
+        set.f107_average,
+        match set.magnetic {
+            nrlmsise00::MagneticActivity::History { ap } => ap[0],
+            nrlmsise00::MagneticActivity::Daily { ap } => ap,
+            nrlmsise00::MagneticActivity::Off { daily_ap } => daily_ap,
+        }
+    );
+
+    // A vertical profile. The derivative of the density with respect to altitude comes
+    // from seeding the altitude as a dual number, so the scale height below is the model's
+    // own slope rather than a difference of two evaluations.
+    println!();
+    println!("  altitude      density   temperature   scale height   most of it");
+    println!("        km       kg/m^3             K             km");
+    for altitude in [0.0, 20.0, 50.0, 90.0, 150.0, 250.0, 400.0, 600.0, 1000.0] {
+        let output = model
+            .gtd7d_at(
+                when,
+                Dual64::from(altitude).derivative(),
+                Dual64::lift(latitude),
+                Dual64::lift(longitude),
+            )
+            .expect("inside the driver record");
+
+        // rho = rho0 * exp(-z/H), so H = -rho / (d rho / dz).
+        let scale_height = -output.total_mass.re / output.total_mass.eps;
+        let (species, share) = dominant(&output);
+        println!(
+            "  {altitude:>8.0}   {:>10.3e}   {:>11.1}   {scale_height:>12.1}   \
+             {species:<3} {:>3.0}%",
+            output.total_mass.re,
+            output.temperature.re,
+            share * 100.0
+        );
+    }
+
+    // The same altitude through two solar cycles. Nothing changes but the date, and the
+    // drivers that the date selects.
+    println!();
+    println!("  The thermosphere at 400 km follows the Sun:");
+    for (year, what) in [
+        (2001, "cycle 23 maximum"),
+        (2008, "the deepest minimum of the space age"),
+        (2014, "cycle 24 maximum, a weaker one"),
+        (2020, "cycle 25 minimum"),
+        (2024, "cycle 25 maximum"),
+    ] {
+        let midsummer = drivers::mjd(year, 6, 21, 12, 0, 0.0);
+        let output = model
+            .gtd7d_at(midsummer, 400.0, latitude, longitude)
+            .expect("inside the driver record");
+        let day = drivers::Drivers::embedded()
+            .day(midsummer)
+            .expect("inside the driver record");
+        println!(
+            "  {year}-06-21   {:>10.3e} kg/m^3   F10.7 {:>5.1}   {what}",
+            output.total_mass, day.flux
+        );
+    }
+
+    // What a drag calculation actually asks for, and a check that the derivative is the
+    // one a finite difference converges to.
+    println!();
+    println!("  What a drag partial needs at 400 km:");
+    let seeded = model
+        .gtd7d_at(
+            when,
+            Dual64::from(400.0).derivative(),
+            Dual64::lift(latitude),
+            Dual64::lift(longitude),
+        )
+        .expect("inside the driver record");
+    let at = |altitude: f64| {
+        model
+            .gtd7d_at(when, altitude, latitude, longitude)
+            .expect("inside the driver record")
+            .total_mass
+    };
+    // A thousandth of the scale height. The difference's own truncation error grows as
+    // the square of the step, so at this size it is down around a part in five million,
+    // and what is left over below is that rather than anything the model did.
+    let step = 0.07;
+    let difference = (at(400.0 + step) - at(400.0 - step)) / (2.0 * step);
+    println!(
+        "    density                {:>12.5e} kg/m^3",
+        seeded.total_mass.re
+    );
+    println!(
+        "    from the model         {:>12.5e} kg/m^3 per km",
+        seeded.total_mass.eps
+    );
+    println!(
+        "    central difference     {:>12.5e} kg/m^3 per km",
+        difference
+    );
+    println!(
+        "    they agree to          {:>12.1e} relative",
+        (seeded.total_mass.eps - difference).abs() / difference.abs()
+    );
+    println!("    and the residual there is the difference's error, not the model's");
+    println!(
+        "    scale height           {:>12.1} km",
+        -seeded.total_mass.re / seeded.total_mass.eps
+    );
+}

hooks/pre-push

@@ -19,4 +19,7 @@ check() {
 
 check "the tests" cargo test --workspace
 check "the tests without default features" cargo test --no-default-features
+# Every target, not only the tests: an example or a benchmark that needs a feature it did
+# not declare compiles fine by default and fails only here.
+check "every target without default features" cargo build --no-default-features --all-targets
 check "the docs" cargo doc --no-deps --workspace

readme.md

@@ -58,6 +58,10 @@ assert!((scale_height - 55.8).abs() < 0.1);
 # fn main() {}
 ```
 
+`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.
+
 ## `gtd7` or `gtd7d`
 
 They differ in one number. `gtd7`'s total mass density sums helium, atomic oxygen,