connor/astrokit
1
0 Fork 0
Code Issues 20 Pull Requests Actions Packages Projects Releases Wiki Activity

Added a basic Harris-Priester model (with one change for smoothness)

Browse Source
This commit is contained in:
connor committed 2026-10-02 14:25:50 -04:00
1 parent 36e96b0b65
commit b60dbacc67
12 files changed
+18412 -18
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

.githooks/pre-commit
+7 -2
View File
@@ -13,8 +13,11 @@ cargo clippy --all-targets --all-features -- -D warnings
echo "pre-commit: US spelling"
staged=$(git diff --cached --name-only --diff-filter=ACM -- '*.rs' '*.md' '*.toml')
if [ -n "$staged" ]; then
# "GFZ German Research Centre for Geosciences" is an organization's own name, carried
# verbatim because the space weather record's CC BY 4.0 license requires the
# attribution as written. Not ours to Americanize.
# shellcheck disable=SC2086
if git grep -wniE -e 'metre|metres|kilometre|kilometres' \
hits=$(git grep -wniE -e 'metre|metres|kilometre|kilometres' \
-e 'licence|defence|programme|catalogue|grey' \
-e 'manoeuvre|manoeuvres|manoeuvred|manoeuvring' \
-e 'normalise|normalised|normalises|normalising|normalisation' \
@@ -28,7 +31,9 @@ if [ -n "$staged" ]; then
-e 'initialise|initialised|initialisation|utilise|utilised' \
-e 'modelling|modelled|labelled|cancelled|travelled|fuelled' \
-e 'favour|favours|honour|honours|neighbour|neighbours' \
--cached -- $staged; then
--cached -- $staged | grep -v 'Research Centre for Geosciences' || true)
if [ -n "$hits" ]; then
echo "$hits"
echo
echo "pre-commit: British spellings found above. This project uses US English."
exit 1
examples/full.rs
+37 -14
View File
@@ -1,6 +1,6 @@
//! A full propagation: a spherical harmonic field, the Sun and Moon, sunlight pushing on
//! a faceted spacecraft and charging it, attitude carried alongside the orbit, and the
//! ground track that follows from both.
//! A full propagation: a spherical harmonic field, the Sun and Moon, sunlight pushing on a
//! faceted spacecraft and charging it, drag against a co-rotating atmosphere, attitude
//! carried alongside the orbit, and the ground track that follows from both.
//!
//! Run with `cargo run --example full`. Needs kernels — `scripts/fetch-kernels.sh`.
//!
@@ -10,6 +10,7 @@
use std::path::PathBuf;
use std::sync::Arc;
use astrokit::dynamics::Dynamics;
use astrokit::eclipse;
use astrokit::prelude::*;
use astrokit::{attitude, kernels};
@@ -71,6 +72,9 @@ fn main() -> Result<(), Box<dyn std::error::Error>> {
ThirdBody { body: Body::MOON },
ThirdBody { body: Body::SUN },
Sunshine,
Drag {
density: DensityModel::harris_priester(),
},
];
let field = GravityModel::Egm96.field(DEGREE, DEGREE)?;
@@ -99,8 +103,31 @@ fn main() -> Result<(), Box<dyn std::error::Error>> {
Full.state(*plain.run().states.last().unwrap()).position()
};
let point_mass = ending_with(&[TwoBody { body: Body::EARTH }]);
let without_third_bodies = ending_with(&forces[..1]);
let without_sunlight = ending_with(&forces[..3]);
// What each force is worth: the same arc with that one taken out. The field is the
// exception — removing it leaves nothing to orbit — so for that one the harmonics are
// swapped for a point mass instead.
let labels = [
"gravity harmonics",
"the moon",
"the sun",
"sunlight",
"drag",
];
let worth: Vec<f64> = (0..forces.len())
.map(|index| {
let mut list: Vec<Dynamics> = forces
.iter()
.enumerate()
.filter(|(at, _)| *at != index)
.map(|(_, force)| *force)
.collect();
if index == 0 {
list.insert(0, TwoBody { body: Body::EARTH });
}
(final_state.position() - ending_with(&list)).norm()
})
.collect();
let at_end = time::advance(epoch, arc.to_seconds());
let elements = final_state.keplerian(frame)?;
@@ -225,17 +252,13 @@ fn main() -> Result<(), Box<dyn std::error::Error>> {
}
);
println!(
"against a point mass {:.3} km",
"against a point mass {:.3} km",
(final_state.position() - point_mass).norm() / 1e3
);
println!(
"of that, the sun and moon {:.3} km",
(final_state.position() - without_third_bodies).norm() / 1e3
);
println!(
"of that, sunlight {:.3} km",
(final_state.position() - without_sunlight).norm() / 1e3
);
println!("what each force is worth over the day, by leaving it out:");
for (label, moved) in labels.iter().zip(&worth) {
println!(" {label:<18}{:>12.3} km", moved / 1e3);
}
// Dense output: any time in the arc, not only the ones a step landed on.
let midway = Full.state(solution.interpolate(arc.to_seconds() / 2.0));
scripts/fetch-space-weather.py
Executable
+276
View File
@@ -0,0 +1,276 @@
#!/usr/bin/env python3
"""Regenerates src/atmosphere/space_weather.rs, the drivers a density model reads.
Usage: scripts/fetch-space-weather.py [path to a local copy of the source file]
The record is emitted as Rust rather than as data to parse, the same way the gravity
coefficients are. Nothing at runtime reads a file or validates a field.
Unlike the gravity coefficients and the kernels, this record is *meant* to change: it gains
a row a day. So there is no digest to pin it to. What this records instead is the digest of
what was downloaded, with the date and the coverage, in the output's own header —
provenance for a file that moves rather than a lock on one that does not.
Source: GFZ Potsdam's combined Kp/ap/Ap/SN/F10.7 record. Chosen over Celestrak's
SW-All.csv for two reasons. It carries an explicit license where Celestrak states none, and
that matters for a file committed into an MIT crate. And Celestrak's runs decades past
today on its own predictions, which are not measurements; a propagation that silently used
them would not be reproducible.
The sunspot column is dropped. It is the one part of the GFZ file under CC BY-NC 4.0 rather
than CC BY 4.0, and no density model here wants it.
"""
import hashlib
import subprocess
import sys
from datetime import date, timedelta
from pathlib import Path
URL = "https://kp.gfz-potsdam.de/app/files/Kp_ap_Ap_SN_F107_since_1932.txt"
ATTRIBUTION = (
"GFZ German Research Centre for Geosciences and Geomagnetic Observatory Niemegk, "
"CC BY 4.0"
)
MISSING = -1.0
# Geomagnetic indices reach back to 1932 and the solar flux to February 1947, but the early
# flux record is sporadic — 105 missing days in 1947 and one hole of 44. From 1957 it is all
# but complete: 39 gaps in seventy years, none longer than five days, none longer than two
# since 2020. So that is where the retained record starts, which is also where the other
# published driver files start.
START = date(1957, 1, 1)
# A run of unmeasured days this long is refused rather than interpolated. Five is the worst
# the source holds from 1957 on, so this leaves room for it to change without quietly
# inventing a fortnight of weather.
LONGEST_FILL = 7
ROOT = Path(__file__).resolve().parent.parent
OUT = ROOT / "src" / "atmosphere" / "space_weather.rs"
# Modified Julian Day zero. The record is indexed by MJD because a UT day is not always
# 86400 seconds and MJD still increments by exactly one across a leap second.
MJD_EPOCH = date(1858, 11, 17)
# Field positions in the source, which is whitespace-aligned with 28 columns.
YEAR, MONTH, DAY = 0, 1, 2
AP_FIRST, AP_LAST = 15, 22
AP_DAILY = 23
F107_OBSERVED, F107_ADJUSTED = 25, 26
DEFINITIVE = 27
def source_text(argv):
if len(argv) > 1:
raw = Path(argv[1]).read_bytes()
print(f"reading {argv[1]}")
else:
print(f"downloading {URL}")
raw = subprocess.run(
["curl", "-sS", "--fail", "--max-time", "180", "-L", URL],
check=True,
capture_output=True,
).stdout
return raw.decode(), hashlib.sha256(raw).hexdigest()
def parse(text):
"""Daily rows from START onward, with the flux left as None where unmeasured."""
rows = []
for line in text.splitlines():
if line.startswith("#") or not line.strip():
continue
fields = line.split()
if len(fields) != 28:
raise SystemExit(f"expected 28 fields, found {len(fields)}: {line[:60]}")
when = date(int(fields[YEAR]), int(fields[MONTH]), int(fields[DAY]))
if when < START:
continue
def flux(at):
value = float(fields[at])
return None if value <= MISSING else value
rows.append(
{
"date": when,
"ap": [int(float(f)) for f in fields[AP_FIRST : AP_LAST + 1]],
"ap_daily": int(float(fields[AP_DAILY])),
"observed": flux(F107_OBSERVED),
"adjusted": flux(F107_ADJUSTED),
"definitive": int(fields[DEFINITIVE]),
"filled": False,
}
)
if not rows:
raise SystemExit(f"no rows at or after {START}")
# The record is addressed by day, so a missing date would shift everything after it.
for earlier, later in zip(rows, rows[1:]):
if later["date"] - earlier["date"] != timedelta(days=1):
raise SystemExit(f"missing day between {earlier['date']} and {later['date']}")
return rows
def trim(rows):
"""Drops any unmeasured days at either end.
START is a round date rather than a measurement, and the flux happens not to have been
taken on it. The far end matters more: an unmeasured day there would have to be
extrapolated rather than interpolated, and the whole point of refusing to propagate
past the record is not to do that.
"""
measured = [index for index, row in enumerate(rows) if row["observed"] is not None]
if not measured:
raise SystemExit("no measured solar flux at all")
return rows[measured[0] : measured[-1] + 1]
def fill(rows):
"""Interpolates the days the flux was not measured, marks them, counts them.
The gaps are isolated — one or two days, five at the very worst — and the flux moves a
few percent over that, so a straight line between the measurements either side is a
fair reconstruction. A long gap is not, and is refused.
"""
measured = [index for index, row in enumerate(rows) if row["observed"] is not None]
filled = 0
for before, after in zip(measured, measured[1:]):
span = after - before
if span == 1:
continue
if span - 1 >= LONGEST_FILL:
raise SystemExit(
f"{span - 1} unmeasured days from {rows[before + 1]['date']}, too many "
f"to interpolate"
)
for offset in range(1, span):
weight = offset / span
for key in ("observed", "adjusted"):
rows[before + offset][key] = (
rows[before][key] * (1.0 - weight) + rows[after][key] * weight
)
rows[before + offset]["filled"] = True
filled += 1
return filled
def main():
text, digest = source_text(sys.argv)
rows = trim(parse(text))
filled = fill(rows)
print(f"digest {digest}")
print(f"interpolated {filled} unmeasured days of {len(rows)}")
OUT.parent.mkdir(parents=True, exist_ok=True)
OUT.write_text(rust_source(rows, digest, date.today().isoformat(), filled))
print(
f"wrote {OUT.relative_to(ROOT)}: {len(rows)} days, "
f"{rows[0]['date']} to {rows[-1]['date']}, {OUT.stat().st_size / 1e6:.1f} MB"
)
def mjd(when):
return when.toordinal() - MJD_EPOCH.toordinal()
def wrapped(tokens, width=96, indent=" "):
"""The tokens as comma-separated lines, packed to `width`."""
lines, line = [], indent
for token in tokens:
piece = token + ","
if len(line) + len(piece) > width and line != indent:
lines.append(line)
line = indent
line += piece
if line != indent:
lines.append(line)
return "\n".join(lines)
def array(name, kind, doc, tokens, count):
return (
f"{doc}\n#[rustfmt::skip]\npub(super) static {name}: [{kind}; {count}] = [\n"
f"{wrapped(tokens)}\n];\n"
)
def rust_source(rows, digest, downloaded, filled):
"""The generated module: provenance in the docs, the record in six flat arrays."""
first, last = rows[0]["date"], rows[-1]["date"]
doc = f"""//! Solar and geomagnetic drivers for atmospheric density, one row per UT day.
//!
//! Generated by `scripts/fetch-space-weather.py`; do not edit by hand. Committed as Rust
//! rather than as data so that a clone needs no network, no data directory and no parser.
//!
//! Source: <{URL}>
//! downloaded {downloaded}, SHA-256 `{digest}`.
//!
//! Attribution: {ATTRIBUTION}.
//!
//! The digest records what was downloaded rather than pinning it: the source gains a row a
//! day, so there is nothing fixed to pin. Its sunspot column is dropped, being the one part
//! of it under a non-commercial license.
//!
//! Coverage: {first} to {last}, {len(rows)} days. Geomagnetic indices reach back to 1932
//! and the flux to 1947, but the early flux record is too sporadic to use: 105 missing days
//! in 1947 and one hole of 44. Of the days here, {filled} had no flux measurement and were
//! interpolated from their neighbors, never more than five in a row.
"""
parts = [
doc,
f"/// Modified Julian Day of the first row, {first}.\n"
f"pub(super) const FIRST_DAY: i64 = {mjd(first)};\n",
f"/// How many days the record holds.\npub(super) const DAYS: usize = {len(rows)};\n",
array(
"AP",
"[f64; 8]",
"/// The eight three-hour geomagnetic indices of each day, earliest first.",
["[" + ",".join(f"{v}." for v in row["ap"]) + "]" for row in rows],
"DAYS",
),
array(
"AP_DAILY",
"f64",
"/// The daily mean geomagnetic index.",
[f"{row['ap_daily']}." for row in rows],
"DAYS",
),
array(
"FLUX_OBSERVED",
"f64",
"/// Solar flux as measured, in solar flux units.",
[f"{row['observed']:.1f}" for row in rows],
"DAYS",
),
array(
"FLUX_ADJUSTED",
"f64",
"/// Solar flux scaled to one astronomical unit.",
[f"{row['adjusted']:.1f}" for row in rows],
"DAYS",
),
array(
"DEFINITIVE",
"bool",
"/// Whether the geomagnetic indices are settled rather than preliminary.",
["true" if row["definitive"] else "false" for row in rows],
"DAYS",
),
array(
"FLUX_MEASURED",
"bool",
"/// Whether the flux was measured that day rather than interpolated across a gap.",
["false" if row["filled"] else "true" for row in rows],
"DAYS",
),
]
return "\n".join(parts)
if __name__ == "__main__":
main()
src/atmosphere/drivers.rs
+525
View File
@@ -0,0 +1,525 @@
//! The solar and geomagnetic record a density model reads.
//!
//! Generated into `space_weather.rs` by
//! `scripts/fetch-space-weather.py` and compiled in, so a propagation needs no network, no
//! data directory and nothing to parse. Refreshing it is a deliberate act: a fetch inside a
//! propagation would mean the same program gave different answers on different days.
//!
//! The record is GFZ Potsdam's, used under CC BY 4.0, attributed to the
//! GFZ German Research Centre for Geosciences and Geomagnetic Observatory Niemegk. The
//! generated module carries the source URL, the download date and the digest of what was
//! downloaded.
//!
//! # The record ends, and so does the propagation
//!
//! Drivers are measurements. There are none for tomorrow, so an arc that reaches past the
//! record is refused rather than extrapolated. Predicting them is a future effort, recorded
//! in Issue #32.
//!
//! The usable end is earlier than the last row, because the 81-day solar flux average that
//! the model wants is *centered*: the last day that has one is forty days before the last
//! day that has a flux measurement. [`Drivers::covers`] accounts for that.
use core::fmt;
use std::sync::OnceLock;
use hifitime::Epoch;
use thiserror::Error;
use super::space_weather::{
AP, AP_DAILY, DAYS, DEFINITIVE, FIRST_DAY, FLUX_ADJUSTED, FLUX_MEASURED, FLUX_OBSERVED,
};
/// Days either side of a day that its centered flux average spans.
///
/// The average is over 81 days, so forty each way and the day itself.
pub const FLUX_AVERAGE_RADIUS: i64 = 40;
/// One day of the record.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct Day {
/// Solar flux as measured, in solar flux units.
pub flux_observed: f64,
/// Solar flux scaled to one astronomical unit.
pub flux_adjusted: f64,
/// Daily mean geomagnetic index.
pub ap_daily: f64,
/// The eight three-hour geomagnetic indices of the day, earliest first.
pub ap: [f64; 8],
/// Whether the geomagnetic indices are settled rather than preliminary.
pub definitive: bool,
/// Whether the flux was measured on this day rather than interpolated across a gap.
pub flux_measured: bool,
}
/// A day the record does not cover.
#[derive(Clone, Copy, Debug, Error, PartialEq, Eq)]
pub enum DriverError {
/// Before the record starts.
#[error("{} is before the driver record, which starts at {}", civil(*asked), civil(*first))]
BeforeRecord {
/// The day asked for.
asked: Epoch,
/// The first day the record holds.
first: Epoch,
},
/// After the record ends, which includes the future.
#[error(
"{} is past the driver record, which reaches {}; drivers are measurements and \
there are none for later, so refresh with scripts/fetch-space-weather.py or \
propagate inside the record",
civil(*asked),
civil(*last)
)]
PastRecord {
/// The day asked for.
asked: Epoch,
/// The last day usable, allowing for the centered flux average.
last: Epoch,
},
}
/// An epoch as the civil date it falls on, which is the resolution these errors are about.
fn civil(epoch: Epoch) -> String {
let (year, month, day, ..) = epoch.to_gregorian_utc();
format!("{year:04}-{month:02}-{day:02}")
}
/// The Modified Julian Day `epoch` falls on.
///
/// The record is one row per UT *day* and a UTC day is not always 86400 seconds, so the
/// index has to come from the calendar. MJD increments by exactly one across a leap second,
/// where a count of elapsed seconds would slip.
fn day_number(epoch: Epoch) -> i64 {
epoch.to_mjd_utc_days().floor() as i64
}
/// Midnight UTC on a Modified Julian Day.
fn epoch_of(day: i64) -> Epoch {
Epoch::from_mjd_utc(day as f64)
}
/// The record's row at `index`.
fn row(index: usize) -> Day {
Day {
flux_observed: FLUX_OBSERVED[index],
flux_adjusted: FLUX_ADJUSTED[index],
ap_daily: AP_DAILY[index],
ap: AP[index],
definitive: DEFINITIVE[index],
flux_measured: FLUX_MEASURED[index],
}
}
/// The whole record, indexed by day.
pub struct Drivers {
/// The 81-day centered flux average, one per day, `None` where the window runs off
/// either end of the record.
centered_flux: Vec<Option<f64>>,
}
impl fmt::Debug for Drivers {
/// The coverage rather than twenty-five thousand rows.
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
write!(
f,
"Drivers({} to {}, {DAYS} days)",
civil(self.first()),
civil(self.last_measured())
)
}
}
impl Drivers {
/// The compiled-in record.
///
/// The flux averages are accumulated on first use; the rows themselves are already
/// there, so there is nothing to read and nothing that can fail.
pub fn committed() -> &'static Self {
static RECORD: OnceLock<Drivers> = OnceLock::new();
RECORD.get_or_init(|| Self {
centered_flux: centered_averages(&FLUX_OBSERVED),
})
}
/// The first day the record holds.
pub fn first(&self) -> Epoch {
epoch_of(FIRST_DAY)
}
/// The last day the record holds a measurement for.
pub fn last_measured(&self) -> Epoch {
epoch_of(FIRST_DAY + DAYS as i64 - 1)
}
/// The last day a model can be evaluated on.
///
/// Forty days before [`Drivers::last_measured`], because the centered flux average
/// needs that much of a future to average over.
pub fn last_usable(&self) -> Epoch {
epoch_of(FIRST_DAY + DAYS as i64 - 1 - FLUX_AVERAGE_RADIUS)
}
/// How many days the record holds.
pub fn len(&self) -> usize {
DAYS
}
/// Whether the record holds no days at all.
pub fn is_empty(&self) -> bool {
DAYS == 0
}
/// Whether `epoch` falls on a day a model can be evaluated on.
pub fn covers(&self, epoch: Epoch) -> bool {
self.index_of(epoch).is_ok()
}
/// Checks that an arc from `start` to `end` is covered throughout.
///
/// Only the ends are examined, which is enough: the record is one contiguous run of
/// days, so an arc with both ends inside it has nothing missing in between. The order
/// of the two does not matter, since a propagation may run backwards.
///
/// # Errors
///
/// Naming the end that falls outside and the day the record reaches.
pub fn covers_arc(&self, start: Epoch, end: Epoch) -> Result<(), DriverError> {
self.index_of(start)?;
self.index_of(end)?;
Ok(())
}
/// The day `epoch` falls in.
///
/// # Errors
///
/// If the day is outside what the record can be evaluated on.
pub fn day(&self, epoch: Epoch) -> Result<Day, DriverError> {
Ok(row(self.index_of(epoch)?))
}
/// The 81-day centered average of the observed flux around `epoch`.
///
/// # Errors
///
/// If the day is outside what the record can be evaluated on.
pub fn centered_flux(&self, epoch: Epoch) -> Result<f64, DriverError> {
let index = self.index_of(epoch)?;
Ok(self.centered_flux[index].expect("inside the usable range"))
}
/// The observed flux of the day before `epoch`, which is what a model asks for.
///
/// # Errors
///
/// If the day is outside what the record can be evaluated on.
pub fn previous_flux(&self, epoch: Epoch) -> Result<f64, DriverError> {
let index = self.index_of(epoch)?;
Ok(FLUX_OBSERVED[index.saturating_sub(1)])
}
/// The seven geomagnetic indices a thermospheric model asks for.
///
/// In the order the model wants them: the daily mean, then the three-hour index for the
/// interval `epoch` falls in and for the three intervals before it, then the mean of
/// the eight intervals from 12 to 33 hours before, then the mean of the eight from 36
/// to 57 hours before. A storm's effect on density lags it by hours, which is what the
/// trailing means are for.
///
/// # Errors
///
/// If the day, or the two and a half days before it the means reach back over, is
/// outside what the record can be evaluated on.
pub fn ap_history(&self, epoch: Epoch) -> Result<[f64; 7], DriverError> {
let index = self.index_of(epoch)?;
let (_, _, _, hours, _, _, _) = epoch.to_gregorian_utc();
let interval = index as i64 * 8 + i64::from(hours) / 3;
// Three-hour intervals from the start of the record, so the lookback crosses days
// without any special case.
let at = |back: i64| -> f64 {
let wanted = interval - back;
if wanted < 0 {
// Before the record. The first days are 1957 and nothing propagates there,
// so holding the earliest value is enough not to produce a NaN.
return AP[0][0];
}
let (day, slot) = ((wanted / 8) as usize, (wanted % 8) as usize);
AP[day.min(DAYS - 1)][slot]
};
let mean = |first_back: i64| -> f64 {
(0..8).map(|offset| at(first_back + offset)).sum::<f64>() / 8.0
};
Ok([
AP_DAILY[index],
at(0),
at(1),
at(2),
at(3),
// Twelve hours back is four intervals, and the mean runs backwards from there.
mean(4),
// Thirty-six hours back is twelve.
mean(12),
])
}
fn index_of(&self, epoch: Epoch) -> Result<usize, DriverError> {
let asked_day = day_number(epoch);
if asked_day < FIRST_DAY {
return Err(DriverError::BeforeRecord {
asked: epoch,
first: self.first(),
});
}
let index = (asked_day - FIRST_DAY) as usize;
if index >= DAYS || self.centered_flux[index].is_none() {
return Err(DriverError::PastRecord {
asked: epoch,
last: self.last_usable(),
});
}
Ok(index)
}
}
/// The 81-day centered average of the observed flux for every day it can be had.
fn centered_averages(flux: &[f64]) -> Vec<Option<f64>> {
let radius = FLUX_AVERAGE_RADIUS as usize;
let window = 2 * radius + 1;
if flux.len() < window {
return vec![None; flux.len()];
}
// One running sum rather than a sum per day, which at eighty-one days apiece would be
// two million additions for the committed record.
let mut running: f64 = flux[..window].iter().sum();
let mut averages = vec![None; flux.len()];
averages[radius] = Some(running / window as f64);
for index in radius + 1..flux.len() - radius {
running += flux[index + radius] - flux[index - radius - 1];
averages[index] = Some(running / window as f64);
}
averages
}
#[cfg(test)]
mod tests {
use super::*;
use approx::assert_relative_eq;
fn at(year: i32, month: u8, day: u8) -> Epoch {
Epoch::from_gregorian_utc_hms(year, month, day, 12, 0, 0)
}
/// The day index is a Modified Julian Day, which is what the record is addressed by.
#[test]
fn test_the_day_index_is_a_modified_julian_day() {
// MJD is a fixed offset from the Unix day count; 1970-01-01 is MJD 40587.
assert_eq!(day_number(at(1970, 1, 1)), 40587);
// A leap year's end of February.
assert_eq!(day_number(at(2026, 3, 1)) - day_number(at(2026, 2, 28)), 1);
assert_eq!(day_number(at(2024, 3, 1)) - day_number(at(2024, 2, 28)), 2);
// And it round trips through the epoch it names, at any hour of the day.
for hour in [0, 1, 12, 23] {
let when = Epoch::from_gregorian_utc_hms(2003, 10, 29, hour, 30, 0);
assert_eq!(civil(epoch_of(day_number(when))), "2003-10-29");
}
}
/// A UT day is not always 86400 seconds, so the index has to come from the calendar.
/// Across the 2016 leap second, two successive noons are one day apart.
#[test]
fn test_the_index_is_a_calendar_day_not_a_count_of_seconds() {
let drivers = Drivers::committed();
let before = at(2016, 12, 31);
let after = at(2017, 1, 1);
// Really is 86401 seconds of elapsed time.
assert_relative_eq!(
(after.to_tai_seconds() - before.to_tai_seconds()),
86401.0,
epsilon = 1e-6
);
// And still adjacent days in the record.
assert_eq!(
drivers.index_of(after).unwrap() - drivers.index_of(before).unwrap(),
1
);
}
#[test]
fn test_the_committed_record_is_continuous() {
let drivers = Drivers::committed();
assert!(drivers.len() > 25_000, "only {} days", drivers.len());
assert_eq!(civil(drivers.first()), "1957-01-02");
// The span the record claims has to match the number of rows it holds, or some day
// between them is missing and every index after it is wrong.
let span = day_number(drivers.last_measured()) - day_number(drivers.first());
assert_eq!(span + 1, drivers.len() as i64);
}
/// Every value should be physically possible. A column swapped somewhere in the
/// pipeline would show up here rather than as a strange density much later.
#[test]
fn test_every_value_is_in_range() {
let drivers = Drivers::committed();
for index in 0..drivers.len() {
let day = row(index);
let when = civil(epoch_of(FIRST_DAY + index as i64));
assert!(
(40.0..1500.0).contains(&day.flux_observed),
"{when}: observed flux {}",
day.flux_observed
);
assert!(
(40.0..1500.0).contains(&day.flux_adjusted),
"{when}: adjusted flux {}",
day.flux_adjusted
);
// ap runs 0 to 400 by definition.
assert!(
(0.0..=400.0).contains(&day.ap_daily),
"{when}: Ap {}",
day.ap_daily
);
for value in day.ap {
assert!((0.0..=400.0).contains(&value), "{when}: ap {value}");
}
}
}
/// The running sum had better agree with adding the eighty-one days up.
#[test]
fn test_the_centered_average_is_the_average() {
let drivers = Drivers::committed();
let radius = FLUX_AVERAGE_RADIUS as usize;
for index in [radius, radius + 1, 5000, 20_000, drivers.len() - radius - 1] {
let direct: f64 = FLUX_OBSERVED[index - radius..=index + radius]
.iter()
.sum::<f64>()
/ (2 * radius + 1) as f64;
assert_relative_eq!(
drivers.centered_flux[index].unwrap(),
direct,
max_relative = 1e-12
);
}
// And it is absent exactly where the window runs off the end.
assert!(drivers.centered_flux[radius - 1].is_none());
assert!(drivers.centered_flux[drivers.len() - radius].is_none());
}
/// The centered average costs forty days off the usable end, which is the whole reason
/// `last_usable` exists separately from `last_measured`.
#[test]
fn test_the_usable_end_is_forty_days_inside_the_measured_one() {
let drivers = Drivers::committed();
let measured = drivers.last_measured();
let usable = drivers.last_usable();
assert_eq!(
day_number(measured) - day_number(usable),
FLUX_AVERAGE_RADIUS
);
assert!(drivers.covers(usable));
assert!(!drivers.covers(measured));
}
/// An arc is checked at its ends, which is enough for a contiguous record, and the
/// order of the ends does not matter because a propagation may run backwards.
#[test]
fn test_an_arc_is_covered_or_says_which_end_is_not() {
let drivers = Drivers::committed();
let last = drivers.last_usable();
let inside = at(2020, 6, 1);
let also_inside = at(2020, 6, 8);
assert!(drivers.covers_arc(inside, also_inside).is_ok());
// Backwards is the same arc.
assert!(drivers.covers_arc(also_inside, inside).is_ok());
// Ending exactly on the last usable day is still inside.
assert!(drivers.covers_arc(inside, last).is_ok());
// Reaching past the end, which is the case that matters: an arc that starts inside
// the record and runs into the future.
let over = drivers.covers_arc(last, at(2099, 1, 1)).unwrap_err();
assert!(matches!(over, DriverError::PastRecord { .. }), "{over:?}");
// And starting before it.
let under = drivers.covers_arc(at(1940, 1, 1), inside).unwrap_err();
assert!(
matches!(under, DriverError::BeforeRecord { .. }),
"{under:?}"
);
// A day past the usable end is already outside, forty days before the last
// measurement, which is the part that surprises.
assert!(drivers.covers_arc(inside, drivers.last_measured()).is_err());
}
#[test]
fn test_outside_the_record_says_which_side() {
let drivers = Drivers::committed();
let before = drivers.day(at(1950, 1, 1)).unwrap_err();
assert!(
matches!(before, DriverError::BeforeRecord { .. }),
"{before:?}"
);
// The future, which is the case that matters.
let ahead = drivers.day(at(2099, 1, 1)).unwrap_err();
assert!(matches!(ahead, DriverError::PastRecord { .. }), "{ahead:?}");
assert!(
ahead.to_string().contains("fetch-space-weather"),
"the message should name the refresh script: {ahead}"
);
}
/// The geomagnetic lookback crosses midnight without a special case, so a time early in
/// a day should reach into the day before rather than wrapping within its own.
#[test]
fn test_the_geomagnetic_lookback_crosses_days() {
let drivers = Drivers::committed();
let day = at(2000, 6, 15);
let index = drivers.index_of(day).unwrap();
// Noon falls in the fifth three-hour interval, so index 4.
let history = drivers.ap_history(day).unwrap();
assert_eq!(history[0], AP_DAILY[index]);
assert_eq!(history[1], AP[index][4]);
assert_eq!(history[2], AP[index][3]);
assert_eq!(history[3], AP[index][2]);
assert_eq!(history[4], AP[index][1]);
// At one in the morning the three before it are all in the previous day.
let early = Epoch::from_gregorian_utc_hms(2000, 6, 15, 1, 0, 0);
let history = drivers.ap_history(early).unwrap();
assert_eq!(history[1], AP[index][0]);
assert_eq!(history[2], AP[index - 1][7]);
assert_eq!(history[3], AP[index - 1][6]);
assert_eq!(history[4], AP[index - 1][5]);
// The trailing means are means of eight, so they sit inside the range of what they
// average rather than outside it.
for mean in [history[5], history[6]] {
assert!((0.0..=400.0).contains(&mean), "{mean}");
}
}
#[test]
fn test_the_previous_day_is_the_previous_day() {
let drivers = Drivers::committed();
let day = at(2003, 10, 29);
let index = drivers.index_of(day).unwrap();
assert_eq!(
drivers.previous_flux(day).unwrap(),
FLUX_OBSERVED[index - 1]
);
}
}
src/atmosphere/harris_priester.rs
+656
View File
@@ -0,0 +1,656 @@
//! Harris-Priester: a fixed density profile with a diurnal bulge.
//!
//! Two tabulated profiles, a minimum and a maximum, interpolated in altitude and blended
//! by where the spacecraft is relative to the bulge the Sun raises. The model reads no
//! space weather at all, being fitted to one level of solar activity. That makes it cheap
//! and predictable, and it is the reason to reach for something that does read the weather
//! when the actual state of the atmosphere matters.
//!
//! Harris and Priester (1962), as tabulated in Montenbruck and Gill, *Satellite Orbits*,
//! Table 3.8. The values were cross-checked against Orekit's copy of the same table.
//!
//! # The interpolation is not the published one
//!
//! Harris and Priester tabulate densities at fifty heights; the usual presentation then
//! interpolates exponentially within each interval. That form is continuous but its slope
//! is not, so the altitude derivative steps at every tabulated height — by a median of 4%
//! across the table, and by 29% at 120 and 130 km. Differentiating the drag force through
//! such a profile, which is what this crate's dual-number scalars are for, gives a
//! derivative that jumps as a spacecraft crosses a node.
//!
//! This implementation interpolates the logarithm of the density with a monotone cubic
//! instead (Fritsch-Carlson, usually called PCHIP). The tabulated values are still
//! reproduced exactly, the density still falls monotonically with altitude, and the
//! derivative is continuous everywhere.
//!
//! So the density between the tabulated heights is not the published model's. Against the
//! exponential interpolation it differs by:
//!
//! | 100-130 km | 130-200 km | 200-300 km | 300-500 km | 500-800 km | 800-1000 km |
//! |---|---|---|---|---|---|
//! | 13% | 2.6% | 0.27% | 0.25% | 0.17% | 0.47% |
//!
//! Over the band where the model is used for orbital drag that is a few tenths of a
//! percent. The bottom two rows are where the table is spaced 10 and 20 km and its scale
//! height steps by 40% between intervals, so neither interpolation carries much authority
//! there.
//!
//! Anyone reproducing a published Harris-Priester result, rather than wanting a cheap
//! smooth profile, should know that this is the difference.
use std::sync::OnceLock;
use nalgebra::Vector3;
use ordinary_diffeq::scalar::Scalar;
/// How many altitudes the profile tabulates.
const ROWS: usize = 50;
/// How many intervals they bound.
const INTERVALS: usize = ROWS - 1;
/// Altitude in m, then the minimum and maximum density in kg/m^3.
///
/// The two are equal below 130 km, where the bulge has not developed.
#[rustfmt::skip]
static PROFILE: [(f64, f64, f64); ROWS] = [
(100.0e3, 4.974e-07, 4.974e-07),
(120.0e3, 2.490e-08, 2.490e-08),
(130.0e3, 8.377e-09, 8.710e-09),
(140.0e3, 3.899e-09, 4.059e-09),
(150.0e3, 2.122e-09, 2.215e-09),
(160.0e3, 1.263e-09, 1.344e-09),
(170.0e3, 8.008e-10, 8.758e-10),
(180.0e3, 5.283e-10, 6.010e-10),
(190.0e3, 3.617e-10, 4.297e-10),
(200.0e3, 2.557e-10, 3.162e-10),
(210.0e3, 1.839e-10, 2.396e-10),
(220.0e3, 1.341e-10, 1.853e-10),
(230.0e3, 9.949e-11, 1.455e-10),
(240.0e3, 7.488e-11, 1.157e-10),
(250.0e3, 5.709e-11, 9.308e-11),
(260.0e3, 4.403e-11, 7.555e-11),
(270.0e3, 3.430e-11, 6.182e-11),
(280.0e3, 2.697e-11, 5.095e-11),
(290.0e3, 2.139e-11, 4.226e-11),
(300.0e3, 1.708e-11, 3.526e-11),
(320.0e3, 1.099e-11, 2.511e-11),
(340.0e3, 7.214e-12, 1.819e-11),
(360.0e3, 4.824e-12, 1.337e-11),
(380.0e3, 3.274e-12, 9.955e-12),
(400.0e3, 2.249e-12, 7.492e-12),
(420.0e3, 1.558e-12, 5.684e-12),
(440.0e3, 1.091e-12, 4.355e-12),
(460.0e3, 7.701e-13, 3.362e-12),
(480.0e3, 5.474e-13, 2.612e-12),
(500.0e3, 3.916e-13, 2.042e-12),
(520.0e3, 2.819e-13, 1.605e-12),
(540.0e3, 2.042e-13, 1.267e-12),
(560.0e3, 1.488e-13, 1.005e-12),
(580.0e3, 1.092e-13, 7.997e-13),
(600.0e3, 8.070e-14, 6.390e-13),
(620.0e3, 6.012e-14, 5.123e-13),
(640.0e3, 4.519e-14, 4.121e-13),
(660.0e3, 3.430e-14, 3.325e-13),
(680.0e3, 2.632e-14, 2.691e-13),
(700.0e3, 2.043e-14, 2.185e-13),
(720.0e3, 1.607e-14, 1.779e-13),
(740.0e3, 1.281e-14, 1.452e-13),
(760.0e3, 1.036e-14, 1.190e-13),
(780.0e3, 8.496e-15, 9.776e-14),
(800.0e3, 7.069e-15, 8.059e-14),
(840.0e3, 4.680e-15, 5.741e-14),
(880.0e3, 3.200e-15, 4.210e-14),
(920.0e3, 2.210e-15, 3.130e-14),
(960.0e3, 1.560e-15, 2.360e-14),
(1000.0e3, 1.150e-15, 1.810e-14),
];
/// How far east of the Sun the bulge sits, in radians.
///
/// The atmosphere takes time to heat, so the densest point trails the subsolar one by
/// around thirty degrees of rotation — two hours of local time.
pub const BULGE_LAG: f64 = 30.0 * core::f64::consts::PI / 180.0;
/// The exponent for an orbit of low inclination.
pub const LOW_INCLINATION_EXPONENT: f64 = 2.0;
/// The exponent for a polar orbit.
///
/// A polar orbit crosses the bulge quickly and at a steep angle, and the sharper falloff
/// this gives fits the observed drag better.
pub const POLAR_EXPONENT: f64 = 6.0;
/// The lowest altitude the profile covers, m.
pub fn lowest() -> f64 {
PROFILE[0].0
}
/// The highest altitude the profile covers, m.
pub fn highest() -> f64 {
PROFILE[PROFILE.len() - 1].0
}
/// The cubic that [`Interpolant`] fits to one interval, as coefficients in the normalized
/// coordinate `u = (altitude - lower) / span`.
///
/// Power basis rather than Hermite, so an evaluation is three multiply-adds and the
/// derivative another two, with nothing to recompute per call.
#[derive(Clone, Copy, Debug, Default)]
struct Cubic {
a: f64,
b: f64,
c: f64,
d: f64,
}
impl Cubic {
/// The value at `u`.
#[inline]
fn at<T: Scalar>(&self, u: T) -> T {
((u * self.d + self.c) * u + self.b) * u + self.a
}
/// The derivative with respect to `u`.
#[inline]
fn slope(&self, u: f64) -> f64 {
(u * 3.0 * self.d + 2.0 * self.c) * u + self.b
}
}
/// One profile's log density: a cubic per interval, and a straight line off each end.
struct Curve {
cubic: [Cubic; INTERVALS],
/// `ln(density)` and its altitude derivative at the first row, and at the last.
below: (f64, f64),
above: (f64, f64),
}
impl Curve {
/// `ln(density)` at `altitude`.
///
/// Outside the table the end slope is extended, which is an exponential atmosphere at
/// the end interval's scale height. Extending the cubic instead would turn over and
/// eventually go positive.
#[inline]
fn log_density<T: Scalar>(&self, altitude: T) -> T {
let height = altitude.re();
if height < PROFILE[0].0 {
let (value, slope) = self.below;
return (altitude - PROFILE[0].0) * slope + value;
}
if height > PROFILE[ROWS - 1].0 {
let (value, slope) = self.above;
return (altitude - PROFILE[ROWS - 1].0) * slope + value;
}
let index = interval(height);
let span = PROFILE[index + 1].0 - PROFILE[index].0;
self.cubic[index].at((altitude - PROFILE[index].0) / span)
}
/// `d(ln density)/d(altitude)`, in 1/m.
fn log_slope(&self, altitude: f64) -> f64 {
if altitude < PROFILE[0].0 {
return self.below.1;
}
if altitude > PROFILE[ROWS - 1].0 {
return self.above.1;
}
let index = interval(altitude);
let span = PROFILE[index + 1].0 - PROFILE[index].0;
self.cubic[index].slope((altitude - PROFILE[index].0) / span) / span
}
}
/// Both profiles, fitted once.
struct Interpolant {
minimum: Curve,
maximum: Curve,
}
/// The fit, built on first use and shared by every evaluation and every scalar.
///
/// The coefficients are real, and `T op f64` is available, so one table serves a plain
/// propagation and a differentiated one alike.
fn interpolant() -> &'static Interpolant {
static FIT: OnceLock<Interpolant> = OnceLock::new();
FIT.get_or_init(|| Interpolant {
minimum: fit(|row| row.1),
maximum: fit(|row| row.2),
})
}
/// Fits one profile, picking its density out of each row with `of`.
fn fit(of: fn(&(f64, f64, f64)) -> f64) -> Curve {
let height: Vec<f64> = PROFILE.iter().map(|row| row.0).collect();
let value: Vec<f64> = PROFILE.iter().map(|row| of(row).ln()).collect();
let slope = monotone_slopes(&height, &value);
let mut cubic = [Cubic::default(); INTERVALS];
for (index, cubic) in cubic.iter_mut().enumerate() {
let span = height[index + 1] - height[index];
let rise = value[index + 1] - value[index];
*cubic = Cubic {
a: value[index],
b: span * slope[index],
c: 3.0 * rise - span * (2.0 * slope[index] + slope[index + 1]),
d: -2.0 * rise + span * (slope[index] + slope[index + 1]),
};
}
Curve {
cubic,
below: (value[0], slope[0]),
above: (value[ROWS - 1], slope[ROWS - 1]),
}
}
/// Fritsch-Carlson slopes: the cubic through these points that cannot overshoot them.
///
/// Monotonicity is the property worth having here. A natural cubic spline would be
/// smoother still, C2 rather than C1, but nothing guarantees it will not put a rise into a
/// density profile that only falls.
fn monotone_slopes(height: &[f64], value: &[f64]) -> Vec<f64> {
let count = value.len();
let span: Vec<f64> = (0..count - 1).map(|i| height[i + 1] - height[i]).collect();
let secant: Vec<f64> = (0..count - 1)
.map(|i| (value[i + 1] - value[i]) / span[i])
.collect();
let mut slope = vec![0.0; count];
for i in 1..count - 1 {
// A sign change is a turning point, where the only non-overshooting slope is zero.
if secant[i - 1] * secant[i] <= 0.0 {
continue;
}
let (first, second) = (2.0 * span[i] + span[i - 1], span[i] + 2.0 * span[i - 1]);
slope[i] = (first + second) / (first / secant[i - 1] + second / secant[i]);
}
slope[0] = end_slope(span[0], span[1], secant[0], secant[1]);
slope[count - 1] = end_slope(
span[count - 2],
span[count - 3],
secant[count - 2],
secant[count - 3],
);
slope
}
/// The one-sided slope at an end, limited so the end interval stays monotone.
fn end_slope(near: f64, far: f64, secant: f64, next: f64) -> f64 {
let slope = ((2.0 * near + far) * secant - near * next) / (near + far);
if slope * secant <= 0.0 {
0.0
} else if secant * next <= 0.0 && slope.abs() > (3.0 * secant).abs() {
3.0 * secant
} else {
slope
}
}
/// The tabulated interval `altitude` falls in. Only valid inside the table.
fn interval(altitude: f64) -> usize {
match PROFILE.binary_search_by(|row| row.0.partial_cmp(&altitude).expect("no NaN")) {
// The top row opens no interval, so pair it with the one below.
Ok(index) => index.min(INTERVALS - 1),
Err(index) => index - 1,
}
}
/// Density at `altitude` above the ellipsoid, in kg/m^3.
///
/// `cosine_to_bulge` is the cosine of the angle from the bulge's apex, which
/// [`bulge_cosine`] computes, and `exponent` is 2 for an orbit of low inclination through
/// 6 for a polar one.
///
/// Outside the tabulated 100 to 1000 km the end interval's scale height is extended rather
/// than the density being cut to zero, which would put a step in the derivative at 1000 km.
/// Below 100 km the model has no business being used at all, and this only keeps the
/// arithmetic finite on the way down.
pub fn density<T: Scalar>(altitude: T, cosine_to_bulge: T, exponent: f64) -> T {
let fit = interpolant();
let minimum = fit.minimum.log_density(altitude).exp();
let maximum = fit.maximum.log_density(altitude).exp();
// The blend is cos^(n/2) of half the angle to the apex, written as a half-angle
// identity so there is no arc cosine and nothing to clamp.
let half = positive((cosine_to_bulge + 1.0) * 0.5);
minimum + (maximum - minimum) * half.powf(T::lift(exponent * 0.5))
}
/// The cosine of the angle between a position and the bulge's apex.
///
/// Both vectors are in the body-fixed frame, which is the frame that makes this easy: the
/// apex is the Sun direction turned east about the rotation axis, and in body-fixed
/// coordinates that axis is z.
pub fn bulge_cosine<T: Scalar>(position: Vector3<T>, toward_sun: Vector3<T>) -> T {
let (sine, cosine) = (T::lift(BULGE_LAG.sin()), T::lift(BULGE_LAG.cos()));
let apex = Vector3::new(
toward_sun[0] * cosine - toward_sun[1] * sine,
toward_sun[0] * sine + toward_sun[1] * cosine,
toward_sun[2],
);
position.normalize().dot(&apex.normalize())
}
/// The mean of the two profiles at `altitude`, in kg/m^3.
///
/// What the model gives halfway round from the bulge, and a reasonable single number for
/// the density at a height when the diurnal variation is not wanted.
pub fn mean_density(altitude: f64) -> f64 {
density(altitude, 0.0, LOW_INCLINATION_EXPONENT)
}
/// The scale height of the mean profile at `altitude`, in m.
///
/// Differentiated analytically rather than differenced. Unlike the published form, where
/// the scale height is one number per interval, this varies continuously with altitude.
pub fn scale_height(altitude: f64) -> f64 {
let fit = interpolant();
let minimum = fit.minimum.log_density(altitude).exp();
let maximum = fit.maximum.log_density(altitude).exp();
// The mean of two exponentials is not one, so differentiate the sum rather than
// reading a scale height off either profile.
let mean = 0.5 * (minimum + maximum);
let slope = 0.5
* (minimum * fit.minimum.log_slope(altitude) + maximum * fit.maximum.log_slope(altitude));
-mean / slope
}
fn positive<T: Scalar>(x: T) -> T {
if x.re() > 0.0 {
x
} else {
T::lift(0.0)
}
}
#[cfg(test)]
mod tests {
use super::*;
use approx::assert_relative_eq;
/// On a tabulated altitude the interpolation should give the table back.
#[test]
fn test_the_table_rows_are_reproduced() {
for (altitude, minimum, maximum) in PROFILE {
// Fully out of the bulge is the minimum, fully in it the maximum.
assert_relative_eq!(
density(altitude, -1.0, 2.0),
minimum,
max_relative = 1e-12,
epsilon = 1e-30
);
assert_relative_eq!(
density(altitude, 1.0, 2.0),
maximum,
max_relative = 1e-12,
epsilon = 1e-30
);
}
}
/// Density falls with altitude, everywhere, by a lot.
#[test]
fn test_density_falls_monotonically() {
let mut previous = f64::INFINITY;
for step in 0..=900 {
let altitude = 100.0e3 + 1.0e3 * step as f64;
let value = density(altitude, 0.0, 2.0);
assert!(
value < previous,
"{} km: {value:e} is not below {previous:e}",
altitude / 1e3
);
previous = value;
}
// Nine orders of magnitude across the table, which is what an atmosphere does.
let ratio = density(100.0e3, 0.0, 2.0) / density(1000.0e3, 0.0, 2.0);
assert!((1e7..1e10).contains(&ratio), "ratio {ratio:e}");
}
/// The interpolation is exponential, so the implied scale height between any two nearby
/// altitudes should be a physical one rather than an artifact.
#[test]
fn test_the_implied_scale_height_is_physical() {
for altitude in [150.0e3, 300.0e3, 500.0e3, 800.0e3] {
let step = 1.0e3;
let scale =
step / (density(altitude, 0.0, 2.0) / density(altitude + step, 0.0, 2.0)).ln();
assert!(
(10.0e3..200.0e3).contains(&scale),
"{} km: scale height {:.1} km",
altitude / 1e3,
scale / 1e3
);
}
}
/// The bulge is where the Sun has been, not where it is.
#[test]
fn test_the_bulge_trails_the_sun() {
let toward_sun = Vector3::new(1.0, 0.0, 0.0);
// Directly at the Sun is not the densest point.
let at_sun = bulge_cosine(toward_sun, toward_sun);
assert_relative_eq!(at_sun, BULGE_LAG.cos(), max_relative = 1e-12);
// The densest point is thirty degrees east of it.
let apex = Vector3::new(BULGE_LAG.cos(), BULGE_LAG.sin(), 0.0);
assert_relative_eq!(bulge_cosine(apex, toward_sun), 1.0, max_relative = 1e-12);
// And the opposite side of the world is the least dense.
assert_relative_eq!(bulge_cosine(-apex, toward_sun), -1.0, max_relative = 1e-12);
}
/// The diurnal variation is a real factor, and bigger higher up where the bulge is.
#[test]
fn test_the_bulge_matters_more_with_altitude() {
let at = |altitude: f64| density(altitude, 1.0, 2.0) / density(altitude, -1.0, 2.0);
// Below 130 km the two profiles are the same, so there is no variation at all.
assert_relative_eq!(at(120.0e3), 1.0, max_relative = 1e-12);
// And by 600 km it is a factor of several.
assert!(at(600.0e3) > 3.0, "only {:.2} at 600 km", at(600.0e3));
assert!(at(600.0e3) > at(300.0e3), "not growing with altitude");
}
/// A polar orbit's sharper exponent should never raise the density above what a
/// low-inclination one gives, and should fall away faster from the apex.
#[test]
fn test_a_sharper_exponent_narrows_the_bulge() {
let altitude = 500.0e3;
for cosine in [-0.5, 0.0, 0.5, 0.9] {
let low = density(altitude, cosine, LOW_INCLINATION_EXPONENT);
let polar = density(altitude, cosine, POLAR_EXPONENT);
assert!(
polar <= low * (1.0 + 1e-12),
"cos {cosine}: polar {polar:e} above low-inclination {low:e}"
);
}
// At the apex they agree, both being the maximum.
assert_relative_eq!(
density(altitude, 1.0, POLAR_EXPONENT),
density(altitude, 1.0, LOW_INCLINATION_EXPONENT),
max_relative = 1e-12
);
}
/// Outside the table it extends the end interval rather than stepping to zero, so the
/// derivative survives a propagation that strays out of range.
#[test]
fn test_outside_the_table_stays_finite_and_smooth() {
for altitude in [50.0e3, 99.0e3, 1001.0e3, 2000.0e3, 1.0e7] {
let value = density(altitude, 0.0, 2.0);
assert!(value.is_finite() && value > 0.0, "{altitude:e}: {value:e}");
}
// Still falling above the table.
assert!(density(1500.0e3, 0.0, 2.0) < density(1000.0e3, 0.0, 2.0));
// And rising below it.
assert!(density(90.0e3, 0.0, 2.0) > density(100.0e3, 0.0, 2.0));
// No step at the edges: approaching 1000 km from either side agrees.
let below = density(1000.0e3 - 1.0, 0.0, 2.0);
let above = density(1000.0e3 + 1.0, 0.0, 2.0);
assert_relative_eq!(below, above, max_relative = 1e-3);
}
#[cfg(feature = "autodiff")]
#[test]
fn test_the_density_differentiates() {
use num_dual::Dual64;
// Inside an interval the fit is a cubic, so a central difference is second order
// and the agreement is tight.
for altitude in [205.0e3, 450.0e3, 705.0e3] {
let dual = Dual64::from_re(altitude).derivative();
let got = density(dual, Dual64::from_re(0.3), 2.0);
let step = 1.0;
let numerical = (density(altitude + step, 0.3, 2.0)
- density(altitude - step, 0.3, 2.0))
/ (2.0 * step);
assert_relative_eq!(got.eps, numerical, max_relative = 1e-6);
// And it is negative, with a scale height of tens of kilometers.
let scale = -got.re / got.eps;
assert!(
(10.0e3..200.0e3).contains(&scale),
"{} km: scale height {:.1} km",
altitude / 1e3,
scale / 1e3
);
}
// And on a tabulated altitude, where the published form had no derivative at all.
// A central difference straddling a node is only first order, because the
// curvature steps there even though the slope does not, so this is necessarily
// looser than the interval case above.
for altitude in [200.0e3, 400.0e3, 700.0e3] {
let dual = Dual64::from_re(altitude).derivative();
let got = density(dual, Dual64::from_re(0.3), 2.0);
let step = 1.0;
let numerical = (density(altitude + step, 0.3, 2.0)
- density(altitude - step, 0.3, 2.0))
/ (2.0 * step);
assert_relative_eq!(got.eps, numerical, max_relative = 1e-3);
}
}
/// The published exponential interpolation, kept here to measure against.
fn published(altitude: f64, cosine_to_bulge: f64, exponent: f64) -> f64 {
let index = interval(altitude.clamp(PROFILE[0].0, PROFILE[ROWS - 1].0));
let (lower_altitude, lower_min, lower_max) = PROFILE[index];
let (upper_altitude, upper_min, upper_max) = PROFILE[index + 1];
let fraction = (altitude - lower_altitude) / (upper_altitude - lower_altitude);
let layer = |low: f64, high: f64| low * (high / low).powf(fraction);
let minimum = layer(lower_min, upper_min);
let maximum = layer(lower_max, upper_max);
let half = ((cosine_to_bulge + 1.0) * 0.5).max(0.0);
minimum + (maximum - minimum) * half.powf(exponent * 0.5)
}
/// How far apart the slopes either side of `node` are, probed over `step`.
fn slope_mismatch(of: impl Fn(f64) -> f64, node: f64, step: f64) -> f64 {
let below = (of(node) - of(node - step)) / step;
let above = (of(node + step) - of(node)) / step;
(above - below).abs() / below.abs()
}
/// The monotone fit has no corner at a tabulated altitude, which is the whole reason
/// it is not the published interpolation. The published form has one at every node.
#[test]
fn test_the_derivative_is_continuous_at_a_tabulated_altitude() {
let ours = |altitude| density(altitude, 0.3, 2.0);
let theirs = |altitude| published(altitude, 0.3, 2.0);
let mut cornered = 0;
for (node, _, _) in &PROFILE[1..ROWS - 1] {
let mismatch = slope_mismatch(ours, *node, 10.0);
assert!(
mismatch < 5e-3,
"corner at {} km: slopes differ by {mismatch:.2e}",
node / 1e3
);
if slope_mismatch(theirs, *node, 10.0) > 1e-2 {
cornered += 1;
}
}
// Most of the table is a corner in the published form; 120 km is the worst at 29%.
assert!(cornered > 20, "only {cornered} published corners");
assert!(slope_mismatch(theirs, 120.0e3, 10.0) > 0.25);
// What is left either side of a node is curvature, not a step, so it shrinks with
// the probe. A corner does not: the published mismatch is the same at any step.
let node = 120.0e3;
assert!(slope_mismatch(ours, node, 5.0) < 0.6 * slope_mismatch(ours, node, 10.0));
assert!(slope_mismatch(theirs, node, 5.0) > 0.9 * slope_mismatch(theirs, node, 10.0));
}
/// How far the monotone fit moves the density away from the published model. These are
/// the numbers quoted in the module documentation, so they are pinned here.
#[test]
fn test_the_fit_tracks_the_published_interpolation() {
// Band, and the deviation it is allowed. Tight where the model is used for orbital
// drag; loose at the bottom, where the table is spaced 10 and 20 km and its scale
// height steps by 40% between intervals.
let bands = [
(100.0e3, 130.0e3, 14.0e-2),
(130.0e3, 200.0e3, 3.0e-2),
(200.0e3, 300.0e3, 3.0e-3),
(300.0e3, 500.0e3, 3.0e-3),
(500.0e3, 800.0e3, 2.0e-3),
(800.0e3, 1000.0e3, 5.0e-3),
];
for (low, high, allowed) in bands {
let mut worst: f64 = 0.0;
for step in 0..=2000 {
let altitude = low + step as f64 * (high - low) / 2000.0;
let ours = density(altitude, 0.3, 2.0);
let ratio = ours / published(altitude, 0.3, 2.0) - 1.0;
worst = worst.max(ratio.abs());
}
assert!(
worst < allowed,
"{}-{} km: {:.3}% away from the published form, over {:.1}%",
low / 1e3,
high / 1e3,
worst * 100.0,
allowed * 100.0
);
}
}
/// `scale_height` is the mean profile's own log derivative, differentiated rather than
/// differenced, so it should agree with a difference of `mean_density`.
#[test]
fn test_the_scale_height_is_the_mean_profiles_slope() {
for altitude in [150.0e3, 205.0e3, 400.0e3, 455.0e3, 700.0e3, 900.0e3] {
let step = 1.0;
let slope =
(mean_density(altitude + step) - mean_density(altitude - step)) / (2.0 * step);
let differenced = -mean_density(altitude) / slope;
// A central difference of a cubic is second order, so this is tight away from
// a node and only first order on one.
assert_relative_eq!(scale_height(altitude), differenced, max_relative = 1e-3);
assert!(
(10.0e3..200.0e3).contains(&scale_height(altitude)),
"{} km: {:.1} km",
altitude / 1e3,
scale_height(altitude) / 1e3
);
}
}
/// Every tabulated row is still reproduced, so the fit interpolates the table rather
/// than smoothing it.
#[test]
fn test_the_fit_is_an_interpolation_not_a_smoothing() {
for (altitude, minimum, maximum) in PROFILE {
assert_relative_eq!(density(altitude, -1.0, 2.0), minimum, max_relative = 1e-12);
assert_relative_eq!(density(altitude, 1.0, 2.0), maximum, max_relative = 1e-12);
}
}
}
src/atmosphere/mod.rs
+124
View File
@@ -0,0 +1,124 @@
//! Atmospheric density, and the space weather it depends on.
pub mod drivers;
pub mod harris_priester;
mod space_weather;
pub use drivers::Drivers;
use hifitime::Epoch;
use nalgebra::Vector3;
use ordinary_diffeq::scalar::Scalar;
use serde::{Deserialize, Serialize};
use crate::body::Body;
use crate::cartesian::Cartesian;
use crate::frames::Frame;
/// Where and when a density model is being asked about.
///
/// Positions are in the body-fixed frame, which is the frame an atmosphere is at rest in
/// and the frame a diurnal bulge is easiest to place in.
#[derive(Clone, Copy, Debug)]
pub struct Conditions<T: Scalar> {
/// Position in the body-fixed frame, m.
pub position: Vector3<T>,
/// Geodetic altitude, m.
///
/// Carried rather than derived on demand because finding it is an iteration and more
/// than one model wants it.
pub altitude: T,
/// Unit vector toward the Sun, body-fixed.
///
/// Only meaningful when the model asked for it — see [`DensityModel::needs_sun`]. A
/// caller whose model does not is free to leave this zero rather than pay for an
/// ephemeris lookup.
pub toward_sun: Vector3<T>,
/// The body whose atmosphere it is.
pub body: Body,
/// When, for a model that reads a record of space weather.
pub epoch: Epoch,
}
impl<T: Scalar> Conditions<T> {
/// Conditions at a body-fixed position, working out the altitude.
pub fn new(position: Vector3<T>, toward_sun: Vector3<T>, body: Body, epoch: Epoch) -> Self {
let altitude = Cartesian::new(position, Vector3::zeros(), Frame::body_fixed(body))
.geodetic()
.expect("body-fixed by construction")
.height;
Self {
position,
altitude,
toward_sun,
body,
epoch,
}
}
}
/// Which model gives the density.
///
/// An enum rather than a trait, for the reason
/// [`SurfaceModel`](crate::spacecraft::SurfaceModel) is: the density has to come out in
/// whatever scalar is being propagated, so the method is generic, and a trait with a
/// generic method is not object-safe.
#[non_exhaustive]
#[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)]
pub enum DensityModel {
/// Harris-Priester: tabulated profiles with a diurnal bulge, and no space weather.
HarrisPriester {
/// 2 for an orbit of low inclination, 6 for a polar one. See
/// [`harris_priester::LOW_INCLINATION_EXPONENT`].
bulge_exponent: f64,
},
}
impl DensityModel {
/// Harris-Priester with the exponent for an orbit of low inclination.
pub fn harris_priester() -> Self {
Self::HarrisPriester {
bulge_exponent: harris_priester::LOW_INCLINATION_EXPONENT,
}
}
/// Harris-Priester with the exponent for a polar orbit.
pub fn polar_harris_priester() -> Self {
Self::HarrisPriester {
bulge_exponent: harris_priester::POLAR_EXPONENT,
}
}
/// Density, in kg/m^3.
pub fn density<T: Scalar>(&self, at: &Conditions<T>) -> T {
match self {
Self::HarrisPriester { bulge_exponent } => harris_priester::density(
at.altitude,
harris_priester::bulge_cosine(at.position, at.toward_sun),
*bulge_exponent,
),
}
}
/// Whether this model needs to know where the Sun is.
///
/// Asked before the ephemeris is looked up, so that a model with no diurnal term does
/// not pay for one.
pub fn needs_sun(&self) -> bool {
match self {
Self::HarrisPriester { .. } => true,
}
}
/// Whether this model reads a record of space weather, and so has a span of dates it
/// can be evaluated over.
///
/// Harris-Priester is fitted to one level of solar activity and reads nothing, so it
/// does not constrain the epoch. A model that does will, and a propagation reaching
/// past the record has to be refused rather than extrapolated.
pub fn needs_drivers(&self) -> bool {
match self {
Self::HarrisPriester { .. } => false,
}
}
}
src/atmosphere/space_weather.rs
+16355
View File
File diff suppressed because it is too large. Load diff
src/dynamics/mod.rs
+26 -2
View File
@@ -52,6 +52,7 @@ use ordinary_diffeq::scalar::Scalar;
use serde::{Deserialize, Serialize};
use std::ops::{Add, AddAssign};
use crate::atmosphere::DensityModel;
use crate::body::Body;
use crate::frames::Frame;
use crate::gravity::GravityModel;
@@ -327,6 +328,14 @@ pub enum Dynamics {
/// The perturbing body. Must not be the frame's center.
body: Body,
},
/// Drag against the atmosphere, which turns with the body.
///
/// Needs a [`Spacecraft`] for its surface. The relative wind comes out of the frame
/// rotation rather than from a wind model, so nothing else has to be supplied.
Drag {
/// Which model gives the density.
density: DensityModel,
},
/// What sunlight does: radiation pressure, and power if the layout is tracking it.
///
/// Needs a [`Spacecraft`] for its surface, and the eclipse geometry comes from the
@@ -356,6 +365,7 @@ impl Dynamics {
}
Self::ThirdBody { body } => Box::new(models::ThirdBody::new(body)),
Self::Sunshine => Box::new(models::Sunshine::new()),
Self::Drag { density } => Box::new(models::Drag::new(density)),
}
}
@@ -366,7 +376,9 @@ impl Dynamics {
pub fn needs_almanac(&self) -> bool {
match self {
Self::TwoBody { .. } => false,
Self::Gravity { .. } | Self::ThirdBody { .. } | Self::Sunshine => true,
Self::Gravity { .. } | Self::ThirdBody { .. } | Self::Sunshine | Self::Drag { .. } => {
true
}
}
}
@@ -387,7 +399,19 @@ impl Dynamics {
/// Whether this model reads the spacecraft's geometry, and so needs one.
pub fn needs_spacecraft(&self) -> bool {
matches!(self, Self::Sunshine)
matches!(self, Self::Sunshine | Self::Drag { .. })
}
/// The density model, if this spec has one.
///
/// For the pre-flight check: a model that reads space weather can only be evaluated
/// over the span the record covers, and a propagation reaching past that is worth
/// refusing before it starts rather than partway through.
pub fn density_model(&self) -> Option<DensityModel> {
match *self {
Self::Drag { density } => Some(density),
_ => None,
}
}
/// The perturbing body, if this spec is a third body.
src/dynamics/models/drag.rs
+393
View File
@@ -0,0 +1,393 @@
//! Atmospheric drag.
use nalgebra::Vector3;
use ordinary_diffeq::scalar::Scalar;
use super::super::{Context, Contribution, DynamicsModel};
use crate::atmosphere::{Conditions, DensityModel};
use crate::attitude;
use crate::body::Body;
use crate::frames::Frame;
use crate::state::{HasPosition, HasVelocity, State};
/// Drag against a co-rotating atmosphere.
///
/// `-0.5 rho C_d A |v| v`, with `v` the velocity relative to the air rather than to the
/// frame. That distinction is not small: the atmosphere turns with the body, so at the
/// equator the air is already moving at 465 m/s, which is six percent of orbital speed and
/// enters the force quadratically.
///
/// Taking the velocity in the body-fixed frame handles it exactly, because
/// [`Cartesian::to_frame`](crate::cartesian::Cartesian::to_frame)'s `R v + R_dot r` *is* the
/// relative velocity — the second term being the air's own motion. There is no separate
/// wind model and no angular rate to supply.
///
/// # What is not modeled
///
/// The force is along the relative wind, scaled by the surface's effective area. Molecules
/// really re-emit partly specularly and partly diffusely, which tilts the force off that
/// line and gives a torque; that needs a per-face force rather than a per-face area, and is
/// #30 along with the same change for radiation pressure.
///
/// Winds other than co-rotation are not modeled either, and neither is the drop in drag
/// coefficient with altitude as the flow goes from continuum to free-molecular.
pub struct Drag<T: Scalar> {
density: DensityModel,
sun: Frame,
scalar: core::marker::PhantomData<T>,
}
impl<T: Scalar> Drag<T> {
/// Drag with `density` giving the air.
pub fn new(density: DensityModel) -> Self {
Self {
density,
sun: Frame::inertial(Body::SUN),
scalar: core::marker::PhantomData,
}
}
/// The density model being used.
pub fn density(&self) -> DensityModel {
self.density
}
}
impl<T: Scalar, const D: usize, L: HasPosition<D> + HasVelocity<D>> DynamicsModel<T, D, L>
for Drag<T>
{
fn contribution(&self, t: T, y: &State<T, D, L>, context: &Context) -> Contribution<T> {
let almanac = context
.almanac()
.expect("drag needs an almanac; Propagation::run checks");
let spacecraft = context
.spacecraft()
.expect("drag needs a spacecraft to act on; Propagation::run checks");
let epoch = context.epoch_at(t);
let center = context.frame.center;
let (rotation, rate) = context
.frame
.rotation_to(Frame::body_fixed(center), epoch, almanac)
.expect("no rotation to the body-fixed frame; check the kernels cover the epoch");
let (rotation, rate) = (rotation.map(T::lift), rate.map(T::lift));
let position = y.position();
let fixed_position = rotation * position;
// The air is at rest in the rotating frame, so this is the relative wind.
let wind = rotation * y.velocity() + rate * position;
let speed = wind.norm();
// Only looked up when the model has a diurnal term to place.
let toward_sun = if self.density.needs_sun() {
let ephemeris = almanac
.translate_geometric(self.sun.anise(), Frame::inertial(center).anise(), epoch)
.expect("no ephemeris for the Sun; check the kernels cover the epoch");
let sun = Vector3::new(
T::lift(ephemeris.radius_km[0] * 1e3),
T::lift(ephemeris.radius_km[1] * 1e3),
T::lift(ephemeris.radius_km[2] * 1e3),
);
rotation * (sun - position).normalize()
} else {
Vector3::zeros()
};
let density =
self.density
.density(&Conditions::new(fixed_position, toward_sun, center, epoch));
// The flow arrives from the direction of travel, so the incident direction is
// against the relative wind, and a faceted surface wants it in its own frame.
let incident = -(rotation.transpose() * wind) / speed;
let incident = match L::try_attitude(&y.vector) {
Some(attitude) => attitude::to_dcm(attitude) * incident,
None => incident,
};
let area = spacecraft.surface.drag_coefficient_area(incident);
Contribution::from_force(rotation.transpose() * (-wind * (density * speed * area * 0.5)))
}
fn name(&self) -> &'static str {
"drag"
}
}
#[cfg(test)]
mod tests {
use super::*;
use crate::atmosphere::harris_priester;
use crate::cartesian::Cartesian;
use crate::dynamics::Dynamics;
use crate::kepler::{Keplerian, TranslationalLayout};
use crate::kernels;
use crate::propagation::{Accuracy, Propagation};
use crate::spacecraft::Spacecraft;
use crate::state::{Layout, Orbit};
use anise::almanac::Almanac;
use approx::assert_relative_eq;
use hifitime::{Duration, Epoch};
use std::sync::Arc;
const AREA: f64 = 10.0;
const DRAG_COEFFICIENT: f64 = 2.2;
fn epoch() -> Epoch {
Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0)
}
fn almanac() -> Option<Arc<Almanac>> {
let directory = kernels::test_directory()?;
Some(Arc::new(
kernels::load(&directory).expect("kernels are present"),
))
}
fn spacecraft() -> Spacecraft {
Spacecraft::cannonball(500.0, AREA, DRAG_COEFFICIENT, 1.3)
}
/// A prograde equatorial orbit, where the air's own motion is largest.
fn equatorial() -> Cartesian<f64> {
let radius = 6.778e6;
Cartesian::new(
Vector3::new(radius, 0.0, 0.0),
Vector3::new(0.0, (3.986004415e14 / radius).sqrt(), 0.0),
Frame::EARTH_J2000,
)
}
/// The force is against the wind, with the magnitude the definition gives — and the
/// wind is the air's, not the frame's. At the equator those differ by about 465 m/s,
/// which is six percent of orbital speed and enters quadratically.
#[test]
fn test_the_force_opposes_the_relative_wind_not_the_inertial_velocity() {
let Some(almanac) = almanac() else {
eprintln!("skipped: kernels not fetched, run scripts/fetch-kernels.sh");
return;
};
let state = equatorial();
let context = Context::new(epoch(), Frame::EARTH_J2000)
.with_spacecraft(spacecraft())
.with_almanac(almanac.clone());
// The relative velocity by an independent route: the public frame transform, whose
// `R v + R_dot r` is the same quantity the model forms internally.
let fixed = state
.to_frame(Frame::EARTH_FIXED, epoch(), &almanac)
.unwrap();
let relative_speed = fixed.velocity.norm();
// Which really is several hundred meters per second below the inertial speed.
let slip = state.velocity.norm() - relative_speed;
assert!(
(300.0..600.0).contains(&slip),
"the air should be moving: slip {slip:.1} m/s"
);
let force = Drag::<f64>::new(DensityModel::harris_priester())
.contribution(
0.0,
&Orbit
.with_position(state.position)
.with_velocity(state.velocity),
&context,
)
.force;
// Opposing the wind, taken back to the frame the force is reported in.
let against_wind = state
.to_frame(Frame::EARTH_FIXED, epoch(), &almanac)
.unwrap()
.velocity;
assert!(force.norm() > 0.0, "no drag at all");
// The density that would explain this force at each of the two speeds. Only one can
// be what the model used, and the readings differ by twelve percent here, so the
// plausible one identifies which. The relative reading has to be a thermosphere.
let implied = |speed: f64| force.norm() / (0.5 * AREA * DRAG_COEFFICIENT * speed * speed);
let at_relative = implied(relative_speed);
let at_inertial = implied(state.velocity.norm());
assert!(
at_relative / at_inertial > 1.1,
"the two readings are too close to tell apart"
);
assert_relative_eq!(
at_relative,
harris_priester::density(
fixed.geodetic().unwrap().height,
harris_priester::bulge_cosine(fixed.position, {
// The Sun direction the model would have used.
let sun = almanac
.translate_geometric(
Frame::inertial(Body::SUN).anise(),
Frame::EARTH_J2000.anise(),
epoch(),
)
.unwrap()
.radius_km
* 1e3;
let rotation = Frame::EARTH_J2000
.rotation_to(Frame::EARTH_FIXED, epoch(), &almanac)
.unwrap()
.0;
rotation * (sun - state.position).normalize()
}),
harris_priester::LOW_INCLINATION_EXPONENT,
),
max_relative = 1e-9
);
// And the direction is against the wind.
let direction = (Frame::EARTH_J2000
.rotation_to(Frame::EARTH_FIXED, epoch(), &almanac)
.unwrap()
.0
* force)
.normalize();
assert_relative_eq!(
direction.dot(&against_wind.normalize()),
-1.0,
max_relative = 1e-9
);
}
/// Drag takes energy out, so the orbit comes down. The rate is the classic
/// `2 pi (C_d A / m) rho a^2` per revolution, which for this orbit is tens of meters.
#[test]
fn test_drag_lowers_the_orbit() {
let Some(almanac) = almanac() else {
eprintln!("skipped: kernels not fetched, run scripts/fetch-kernels.sh");
return;
};
let orbit = Keplerian::new(
6.778e6,
1.0e-3,
51.6_f64.to_radians(),
0.4,
0.2,
0.0,
Frame::EARTH_J2000,
);
let arc = Duration::from_seconds(86400.0);
let mut propagation = Propagation::new(
epoch(),
Orbit.with_keplerian(orbit),
&[
Dynamics::TwoBody { body: Body::EARTH },
Dynamics::Drag {
density: DensityModel::harris_priester(),
},
],
)
.spacecraft(spacecraft())
.almanac(almanac.clone())
.duration(arc)
.accuracy(Accuracy::Medium);
let solution = propagation.run();
let end = Orbit
.state(*solution.states.last().unwrap())
.keplerian(Frame::EARTH_J2000)
.unwrap();
let fallen = orbit.semi_major_axis - end.semi_major_axis;
// Per-revolution decay times the revolutions in a day, as an order of magnitude.
let per_revolution = 2.0
* core::f64::consts::PI
* (AREA * DRAG_COEFFICIENT / 500.0)
* harris_priester::mean_density(400.0e3)
* orbit.semi_major_axis.powi(2);
let revolutions = arc.to_seconds() / orbit.period();
let expected = per_revolution * revolutions;
assert!(fallen > 0.0, "the orbit rose by {:.1} m", -fallen);
assert!(
(expected / 5.0..expected * 5.0).contains(&fallen),
"fell {fallen:.1} m against {expected:.1} m expected"
);
}
/// The density model changes the answer. With one model in the crate the knob is its
/// bulge exponent, which is what distinguishes a polar orbit's atmosphere from a low
/// inclination one.
#[test]
fn test_the_density_model_changes_the_answer() {
let Some(almanac) = almanac() else {
eprintln!("skipped: kernels not fetched, run scripts/fetch-kernels.sh");
return;
};
let orbit = Keplerian::new(6.778e6, 1.0e-3, 0.9, 0.4, 0.2, 0.0, Frame::EARTH_J2000);
let run = |density: DensityModel| {
let mut propagation = Propagation::new(
epoch(),
Orbit.with_keplerian(orbit),
&[
Dynamics::TwoBody { body: Body::EARTH },
Dynamics::Drag { density },
],
)
.spacecraft(spacecraft())
.almanac(almanac.clone())
.duration(Duration::from_seconds(86400.0))
.accuracy(Accuracy::Medium);
Orbit
.state(*propagation.run().states.last().unwrap())
.position()
};
let low = run(DensityModel::harris_priester());
let polar = run(DensityModel::polar_harris_priester());
let drag_free = {
let mut propagation = Propagation::new(
epoch(),
Orbit.with_keplerian(orbit),
&[Dynamics::TwoBody { body: Body::EARTH }],
)
.duration(Duration::from_seconds(86400.0))
.accuracy(Accuracy::Medium);
Orbit
.state(*propagation.run().states.last().unwrap())
.position()
};
// Drag moves the orbit by tens of kilometers over the day, either way.
let moved = (low - drag_free).norm();
assert!(moved > 10.0e3, "drag only moved it {moved:.1} m");
assert!((polar - drag_free).norm() > 10.0e3);
// The sharper polar bulge gives less density away from the apex, so it drags less
// and stays nearer the drag-free arc.
assert!(
(polar - drag_free).norm() < moved,
"the sharper bulge should drag less"
);
// And the two are far enough apart to be a usable knob: a good fraction of the
// drag itself, rather than a rounding difference on it.
let between = (low - polar).norm();
assert!(
between > 0.1 * moved,
"only {between:.1} m apart against {moved:.1} m of drag"
);
}
#[test]
#[should_panic(expected = "needs a spacecraft")]
fn test_drag_without_a_spacecraft_is_refused() {
let orbit = Keplerian::new(7.0e6, 0.0, 0.9, 0.0, 0.0, 0.0, Frame::EARTH_J2000);
let mut propagation = Propagation::new(
epoch(),
Orbit.with_keplerian(orbit),
&[Dynamics::Drag {
density: DensityModel::harris_priester(),
}],
)
.duration(Duration::from_seconds(60.0));
propagation.run();
}
}
src/dynamics/models/mod.rs
+2
View File
@@ -3,11 +3,13 @@
//! These are what [`Dynamics`](super::Dynamics) builds. A script composes the
//! enum; it only reaches in here to implement a model of its own.
pub mod drag;
pub mod gravity;
pub mod sunshine;
pub mod third_body;
pub mod two_body;
pub use drag::Drag;
pub use gravity::SphericalHarmonic;
pub use sunshine::Sunshine;
pub use third_body::ThirdBody;
src/lib.rs
+2
View File
@@ -9,6 +9,7 @@
#![cfg_attr(feature = "autodiff", doc = include_str!("../readme.md"))]
#![warn(missing_docs)]
pub mod atmosphere;
pub mod attitude;
pub mod body;
pub mod cartesian;
@@ -29,6 +30,7 @@ pub mod time;
/// The types you want in scope to set up and run a propagation.
pub mod prelude {
pub use crate::atmosphere::{Conditions, DensityModel};
pub use crate::body::Body;
pub use crate::cartesian::{Cartesian, FromCartesian};
pub use crate::dynamics::Dynamics::*;
src/propagation.rs
+9
View File
@@ -8,6 +8,7 @@ use nalgebra::SVector;
use ordinary_diffeq::prelude::*;
use serde::{Deserialize, Serialize};
use crate::atmosphere::Drivers;
use crate::attitude;
use crate::cartesian::Cartesian;
use crate::dynamics::{self, Context, Contribution, Dynamics, DynamicsModel};
@@ -472,6 +473,14 @@ where
}
}
}
if let Some(density) = spec.density_model() {
if density.needs_drivers() {
let end = crate::time::advance(self.epoch, duration.to_seconds());
if let Err(error) = Drivers::committed().covers_arc(self.epoch, end) {
panic!("{spec:?} reads space weather: {error}");
}
}
}
if let Some(body) = spec.third_body() {
if body.naif_id == self.frame.center.naif_id {
panic!(