From e9d1d3771a266391e5626e2eebcceb11d62533bc Mon Sep 17 00:00:00 2001 From: Connor Johnstone Date: Fri, 2 Oct 2026 18:38:27 -0400 Subject: [PATCH] Decided not to step the density at 150km, but made it explicit --- examples/density.rs | 14 ++++++-- readme.md | 39 +++++++++++++++++++++++ src/drivers.rs | 28 ++++++++++++++++ src/model.rs | 11 +++++++ tests/drivers.rs | 78 +++++++++++++++++++++++++++++++++++++++++++++ 5 files changed, 167 insertions(+), 3 deletions(-) diff --git a/examples/density.rs b/examples/density.rs index 95da39c..ac09ab7 100644 --- a/examples/density.rs +++ b/examples/density.rs @@ -57,9 +57,17 @@ fn main() { } ); - // 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. + // A vertical profile, from the ground: the model reaches the surface, which is what + // distinguishes MSISE from the MSIS-86 that preceded it. + // + // The storm's drivers are passed at every altitude, which below 80 km is not what the + // reference asks for — it says to use 150, 150 and 4 there, the effects being neither + // large nor well established that low. Large they are not: it is 0.77% at the bottom + // of this table. The readme has the measurements. + // + // 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"); diff --git a/readme.md b/readme.md index 65b5ca9..e7eb2f7 100644 --- a/readme.md +++ b/readme.md @@ -110,6 +110,45 @@ These nodes and that density do not depend on altitude, so filling them reproduc 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 diff --git a/src/drivers.rs b/src/drivers.rs index bc5d80c..d3af36b 100644 --- a/src/drivers.rs +++ b/src/drivers.rs @@ -111,6 +111,29 @@ pub struct DriverSet { pub magnetic: MagneticActivity, } +impl DriverSet { + /// What the reference asks for below 80 km, where the real values do not belong. + /// + /// Its input notes say that the flux and magnetic 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 May 2024 Gannon storm, the largest geomagnetic + /// storm in twenty years, using the day's real drivers instead of these moves the + /// total mass density by 0.77% at 50 km and by the same 0.77% at 20, 10 and 0 km — + /// the same figure at each, because below 62.5 km every species is the fully mixed + /// nitrogen density times a fixed ratio, so the drivers scale the whole profile and + /// change its shape not at all. + /// + /// [`Drivers::at`] and [`crate::Model::gtd7d_at`] do **not** substitute these for you + /// at low altitude, deliberately: the switch 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 it would be ours rather than the model's. Below 80 km, pass these. + pub const LOWER_ATMOSPHERE: Self = Self { + f107: 150.0, + f107_average: 150.0, + magnetic: MagneticActivity::Daily { ap: 4.0 }, + }; +} + /// A day the record cannot answer for. #[derive(Clone, Copy, Debug, PartialEq)] pub enum DriverError { @@ -304,6 +327,11 @@ impl Drivers { /// model's own order, reaching back fifty-seven hours from the three-hour interval the /// moment falls in. /// + /// These are the measurements, whatever altitude they are going to be used at. Below + /// 80 km the reference asks for fixed values instead; see + /// [`DriverSet::LOWER_ATMOSPHERE`] for what they are and why this does not quietly + /// use them. + /// /// # Errors /// /// If the day is outside what the record can be evaluated on. diff --git a/src/model.rs b/src/model.rs index 3161bec..14b30ab 100644 --- a/src/model.rs +++ b/src/model.rs @@ -172,6 +172,17 @@ impl Model { /// `f64` rather than generic because a table lookup has no derivative to carry; /// position stays generic, so a drag partial still comes out. /// + /// # Below 80 km it gives the record's drivers, which is not what the reference asks + /// + /// The reference wants the flux and magnetic inputs fixed at 150, 150 and 4 below + /// 80 km, the effects being neither large nor well established that low. This passes + /// the day's measurements at every altitude instead, for two reasons: substituting at + /// a threshold would step the density by 5.2% at exactly 80 km, which is thirty-five + /// times the largest discontinuity the model itself has, and it would make this + /// disagree with [`Model::gtd7d`] given the same drivers. The difference it leaves is + /// 0.77% at 50 km and below. For work down there, assemble the input with + /// [`crate::drivers::DriverSet::LOWER_ATMOSPHERE`]. + /// /// ``` /// use nrlmsise00::{drivers, Model}; /// diff --git a/tests/drivers.rs b/tests/drivers.rs index 06959c2..1b089d8 100644 --- a/tests/drivers.rs +++ b/tests/drivers.rs @@ -338,3 +338,81 @@ fn test_the_provenance_is_compiled_in() { assert!(note.contains("sha256"), "no digest"); assert!(note.contains("fetch-drivers.py"), "no way to regenerate it"); } + +/// The drivers barely reach below 80 km, which is why the reference asks for fixed values +/// there, and why the convenience entry point does not quietly substitute them. +#[test] +fn test_the_drivers_hardly_matter_below_80_km() { + use nrlmsise00::drivers::DriverSet; + use nrlmsise00::Input; + + let model = Model::new(); + // The Gannon storm, which is as far from quiet as the record gets. + let when = drivers::mjd(2024, 5, 11, 12, 0, 0.0); + let measured = Drivers::embedded().at(when).expect("inside the record"); + let seconds = drivers::seconds_of_day(when); + + let density = |set: DriverSet, altitude_km: f64| { + model + .gtd7d(&Input { + day_of_year: drivers::day_of_year(when), + seconds_of_day: seconds, + altitude_km, + latitude_deg: 40.0, + longitude_deg: -105.0, + local_solar_time_h: seconds / 3600.0 - 105.0 / 15.0, + f107_average: set.f107_average, + f107: set.f107, + magnetic: set.magnetic, + }) + .total_mass + }; + + assert_eq!(DriverSet::LOWER_ATMOSPHERE.f107, 150.0); + assert_eq!(DriverSet::LOWER_ATMOSPHERE.f107_average, 150.0); + assert_eq!( + DriverSet::LOWER_ATMOSPHERE.magnetic, + MagneticActivity::Daily { ap: 4.0 } + ); + + // Below 62.5 km the drivers scale the whole profile by one factor, so the storm's + // effect is the same figure at every altitude down there and it is under a percent. + let mut shares = Vec::new(); + for altitude in [0.0, 10.0, 20.0, 50.0, 60.0] { + let quiet = density(DriverSet::LOWER_ATMOSPHERE, altitude); + let stormy = density(measured, altitude); + let share = (stormy - quiet) / quiet; + assert!( + (0.005..0.01).contains(&share), + "at {altitude} km the storm moves the density by {:.3}%", + share * 100.0 + ); + shares.push(share); + } + for share in &shares { + assert!( + (share - shares[0]).abs() < 1e-12, + "the factor should be the same at every altitude below 62.5 km: {shares:?}" + ); + } + + // Where they do matter, they matter a great deal, which is the other half of the point. + for (altitude, least) in [(200.0, 0.5), (400.0, 2.0)] { + let quiet = density(DriverSet::LOWER_ATMOSPHERE, altitude); + let stormy = density(measured, altitude); + assert!( + (stormy - quiet) / quiet > least, + "at {altitude} km the storm only moves the density by {:.1}%", + (stormy - quiet) / quiet * 100.0 + ); + } + + // And the entry point that looks the drivers up does not substitute: it agrees with + // the measured set at low altitude, not with the fixed one. + let looked_up = model + .gtd7d_at(when, 50.0, 40.0, -105.0) + .expect("inside the record") + .total_mass; + assert_eq!(looked_up, density(measured, 50.0)); + assert_ne!(looked_up, density(DriverSet::LOWER_ATMOSPHERE, 50.0)); +}