ba1212361f
Commit the previously-untracked standard directories and catalog them in the
root README alongside proof:
- aac/ — communication cargo class: portable boards + a ten-level
offline-first cascade (Apache-2.0 code / CC-BY-4.0 spec).
- map/ — geospatial cargo class: range-requestable signed tilesets +
an F0–F4 fidelity cascade (Apache-2.0 code / CC-BY-4.0 spec).
- sandbox/— isolated, metered, attested execution across the heterogeneous
compute fabric (Apache-2.0 code / CC-BY-4.0 spec).
- jouleclaw-otel-authz, jouleclaw-trust — new jouleclaw reference crates.
Root README now lists aac, map, sandbox, and proof in the catalog, the tree,
and the contribution guide.
Co-Authored-By: Claude Opus 4.8 <[email protected]>
555 lines
20 KiB
Rust
555 lines
20 KiB
Rust
//! SDP4 deep-space orbital propagator extension.
|
|
//!
|
|
//! Extends SGP4 for satellites with orbital periods >= 225 minutes
|
|
//! (MEO, GEO, HEO). Handles third-body perturbations from the Sun
|
|
//! and Moon, and resonance effects for 12-hour and 24-hour orbits.
|
|
//!
|
|
//! Based on Vallado et al., "Revisiting Spacetrack Report #3" (2006).
|
|
|
|
use serde::{Deserialize, Serialize};
|
|
|
|
use crate::constants::*;
|
|
|
|
/// Deep-space resonance type.
|
|
#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]
|
|
pub(crate) enum ResonanceType {
|
|
/// No resonance.
|
|
None,
|
|
/// Synchronous (24-hour) resonance.
|
|
Synchronous,
|
|
/// Half-day (12-hour) resonance.
|
|
HalfDay,
|
|
}
|
|
|
|
/// Deep-space state for the SDP4 propagator extension.
|
|
///
|
|
/// Stores lunar-solar perturbation terms, resonance coefficients,
|
|
/// and integration state for deep-space secular effects.
|
|
#[derive(Debug, Clone, Serialize, Deserialize)]
|
|
pub(crate) struct DeepSpaceState {
|
|
// ── Lunar-solar perturbation constants ────────────────────
|
|
pub(crate) se2: f64,
|
|
pub(crate) se3: f64,
|
|
pub(crate) si2: f64,
|
|
pub(crate) si3: f64,
|
|
pub(crate) sl2: f64,
|
|
pub(crate) sl3: f64,
|
|
pub(crate) sl4: f64,
|
|
pub(crate) sgh2: f64,
|
|
pub(crate) sgh3: f64,
|
|
pub(crate) sgh4: f64,
|
|
pub(crate) sh2: f64,
|
|
pub(crate) sh3: f64,
|
|
|
|
pub(crate) ee2: f64,
|
|
pub(crate) e3: f64,
|
|
pub(crate) xi2: f64,
|
|
pub(crate) xi3: f64,
|
|
pub(crate) xl2: f64,
|
|
pub(crate) xl3: f64,
|
|
pub(crate) xl4: f64,
|
|
pub(crate) xgh2: f64,
|
|
pub(crate) xgh3: f64,
|
|
pub(crate) xgh4: f64,
|
|
pub(crate) xh2: f64,
|
|
pub(crate) xh3: f64,
|
|
|
|
// ── Solar and lunar mean elements ────────────────────────
|
|
pub(crate) zmol: f64,
|
|
pub(crate) zmos: f64,
|
|
|
|
// ── Resonance state ──────────────────────────────────────
|
|
pub(crate) resonance: ResonanceType,
|
|
|
|
// Synchronous resonance terms
|
|
pub(crate) del1: f64,
|
|
pub(crate) del2: f64,
|
|
pub(crate) del3: f64,
|
|
pub(crate) xfact: f64,
|
|
pub(crate) xlamo: f64,
|
|
|
|
// Half-day resonance terms
|
|
pub(crate) d2201: f64,
|
|
pub(crate) d2211: f64,
|
|
pub(crate) d3210: f64,
|
|
pub(crate) d3222: f64,
|
|
pub(crate) d4410: f64,
|
|
pub(crate) d4422: f64,
|
|
pub(crate) d5220: f64,
|
|
pub(crate) d5232: f64,
|
|
pub(crate) d5421: f64,
|
|
pub(crate) d5433: f64,
|
|
|
|
// Common integration terms
|
|
pub(crate) dedt: f64,
|
|
pub(crate) didt: f64,
|
|
pub(crate) dmdt: f64,
|
|
pub(crate) dnodt: f64,
|
|
pub(crate) domdt: f64,
|
|
|
|
// Saved epoch elements for integration
|
|
pub(crate) atime: f64,
|
|
pub(crate) xli: f64,
|
|
pub(crate) xni: f64,
|
|
|
|
// Secular rates from deep space
|
|
pub(crate) gsto: f64,
|
|
}
|
|
|
|
impl DeepSpaceState {
|
|
/// Initialize deep-space terms from the orbital elements.
|
|
///
|
|
/// This computes lunar-solar perturbation coefficients and
|
|
/// identifies any orbital resonance conditions.
|
|
#[allow(clippy::too_many_arguments)] // SGP4/SDP4 deep-space state init
|
|
pub(crate) fn initialize(
|
|
epoch_jd: f64,
|
|
ecco: f64,
|
|
inclo: f64,
|
|
nodeo: f64,
|
|
argpo: f64,
|
|
mo: f64,
|
|
no_kozai: f64,
|
|
_bstar: f64,
|
|
) -> Self {
|
|
let cosio = inclo.cos();
|
|
let sinio = inclo.sin();
|
|
|
|
// Julian centuries from J2000
|
|
let tc = (epoch_jd - 2451545.0) / 36525.0;
|
|
|
|
// ── Solar terms ──────────────────────────────────────
|
|
let _zns = 1.19459e-5;
|
|
let zes = 0.01675;
|
|
let _znl = 1.5835218e-4;
|
|
let zel = 0.05490;
|
|
let c1ss = 2.9864797e-6;
|
|
let _c1l = 4.7968065e-7;
|
|
|
|
// Mean anomaly of the sun
|
|
let zmos = (6.2565837 + 0.017201977 * (epoch_jd - 2451545.0))
|
|
.rem_euclid(TWOPI);
|
|
// Mean anomaly of the moon
|
|
let zmol = (2.3555557 + 0.228737 * (epoch_jd - 2451545.0) - 0.0003029 * tc * tc)
|
|
.rem_euclid(TWOPI);
|
|
|
|
let theta2 = cosio * cosio;
|
|
|
|
// Eccentricity functions
|
|
let eosq = ecco * ecco;
|
|
|
|
// ── Compute solar perturbation coefficients ──────────
|
|
let (se2, se3, si2, si3, sl2, sl3, sl4, sgh2, sgh3, sgh4, sh2, sh3) =
|
|
Self::solar_lunar_coefficients(
|
|
inclo, nodeo, argpo, ecco, true, zmos, zes,
|
|
);
|
|
|
|
// ── Compute lunar perturbation coefficients ──────────
|
|
let (ee2, e3, xi2, xi3, xl2, xl3, xl4, xgh2, xgh3, xgh4, xh2, xh3) =
|
|
Self::solar_lunar_coefficients(
|
|
inclo, nodeo, argpo, ecco, false, zmol, zel,
|
|
);
|
|
|
|
// ── Resonance detection ──────────────────────────────
|
|
let period = TWOPI / no_kozai;
|
|
let resonance;
|
|
let mut del1 = 0.0;
|
|
let mut del2 = 0.0;
|
|
let mut del3 = 0.0;
|
|
let mut xfact = 0.0;
|
|
let mut xlamo = 0.0;
|
|
let mut d2201 = 0.0;
|
|
let mut d2211 = 0.0;
|
|
let mut d3210 = 0.0;
|
|
let mut d3222 = 0.0;
|
|
let mut d4410 = 0.0;
|
|
let mut d4422 = 0.0;
|
|
let mut d5220 = 0.0;
|
|
let mut d5232 = 0.0;
|
|
let mut d5421 = 0.0;
|
|
let mut d5433 = 0.0;
|
|
|
|
// Compute GMST at epoch
|
|
let gsto = gmst_from_jd(epoch_jd);
|
|
|
|
// Check for synchronous resonance (period ~ 1436 min)
|
|
if period > 1200.0 && period < 1800.0 {
|
|
resonance = ResonanceType::Synchronous;
|
|
|
|
// Synchronous resonance coefficients
|
|
let g200 = 1.0 + eosq * (-2.5 + 0.8125 * eosq);
|
|
let g310 = 1.0 + 2.0 * eosq;
|
|
let g300 = 1.0 + eosq * (-6.0 + 6.60937 * eosq);
|
|
|
|
let f220 = 0.75 * (1.0 + cosio) * (1.0 + cosio);
|
|
let f311 = 0.9375 * sinio * sinio * (1.0 + 3.0 * cosio) - 0.75 * (1.0 + cosio);
|
|
let f330 = 1.0 + cosio;
|
|
let f330 = 1.875 * f330 * f330 * f330;
|
|
|
|
del1 = 3.0 * no_kozai * no_kozai * (1.0 / (ao_from_no(no_kozai) * ao_from_no(no_kozai)));
|
|
del1 = del1 * f220 * g200 * CK2;
|
|
del2 = 2.0 * del1 * f311 * g310 * CK2 / (ao_from_no(no_kozai));
|
|
del3 = 3.0 * del1 * f330 * g300 * CK2 / (ao_from_no(no_kozai) * ao_from_no(no_kozai));
|
|
|
|
xlamo = (mo + nodeo + argpo - gsto).rem_euclid(TWOPI);
|
|
xfact = argpo + mo + nodeo + no_kozai - gsto - TWOPI / MINUTES_PER_DAY;
|
|
} else if period > 680.0 && period < 760.0 {
|
|
// Half-day resonance (period ~ 720 min)
|
|
resonance = ResonanceType::HalfDay;
|
|
|
|
// Inclination functions for half-day resonance
|
|
let g201 = -0.306 - (ecco - 0.64) * 0.440;
|
|
|
|
let (g211, g310, g322, g410, g422, g520, g521, g532, g533) = if ecco <= 0.65 {
|
|
let g211 = 3.616 - 13.247 * ecco + 16.290 * eosq;
|
|
let g310 = -19.302 + 117.390 * ecco - 228.419 * eosq + 156.591 * eosq * ecco;
|
|
let g322 = -18.9068 + 109.7927 * ecco - 214.6334 * eosq + 146.5816 * eosq * ecco;
|
|
let g410 = -41.122 + 242.694 * ecco - 471.094 * eosq + 313.953 * eosq * ecco;
|
|
let g422 = -146.407 + 841.880 * ecco - 1629.014 * eosq + 1083.435 * eosq * ecco;
|
|
let g520 = -532.114 + 3017.977 * ecco - 5740.032 * eosq + 3708.276 * eosq * ecco;
|
|
let g521 = -822.71072 + 4568.6173 * ecco - 8491.4146 * eosq + 5337.524 * eosq * ecco;
|
|
let g532 = -853.666 + 4690.25 * ecco - 8624.77 * eosq + 5341.4 * eosq * ecco;
|
|
let g533 = -919.2277 + 4988.61 * ecco - 9064.77 * eosq + 5542.21 * eosq * ecco;
|
|
(g211, g310, g322, g410, g422, g520, g521, g532, g533)
|
|
} else {
|
|
let g211 = -72.099 + 331.819 * ecco - 508.738 * eosq + 266.724 * eosq * ecco;
|
|
let g310 = -346.844 + 1582.851 * ecco - 2415.925 * eosq + 1246.113 * eosq * ecco;
|
|
let g322 = -342.585 + 1554.908 * ecco - 2366.899 * eosq + 1215.972 * eosq * ecco;
|
|
let g410 = -1052.797 + 4758.686 * ecco - 7193.992 * eosq + 3651.957 * eosq * ecco;
|
|
let g422 = -3581.69 + 16178.11 * ecco - 24462.77 * eosq + 12422.52 * eosq * ecco;
|
|
let g520 = if ecco > 0.715 {
|
|
-5149.66 + 29936.92 * ecco - 54087.36 * eosq + 31324.56 * eosq * ecco
|
|
} else {
|
|
1464.74 - 4664.75 * ecco + 3763.64 * eosq
|
|
};
|
|
let g521 = -3.908 + 155.237 * ecco - 498.867 * eosq;
|
|
let g532 = -37.995 + 379.826 * ecco - 1241.41 * eosq + 1263.28 * eosq * ecco;
|
|
let g533 = -113.479 + 691.178 * ecco - 1769.05 * eosq + 1571.84 * eosq * ecco;
|
|
(g211, g310, g322, g410, g422, g520, g521, g532, g533)
|
|
};
|
|
|
|
let sini2o = sinio * sinio;
|
|
let f220 = 0.75 * (1.0 + 2.0 * cosio + theta2);
|
|
let f221 = 1.5 * sini2o;
|
|
let f321 = 1.875 * sinio * (1.0 - 2.0 * cosio - 3.0 * theta2);
|
|
let f322 = -1.875 * sinio * (1.0 + 2.0 * cosio - 3.0 * theta2);
|
|
let f441 = 35.0 * sini2o * f220;
|
|
let f442 = 39.375 * sini2o * sini2o;
|
|
let f522 = 9.84375
|
|
* sinio
|
|
* (sini2o * (1.0 - 2.0 * cosio - 5.0 * theta2)
|
|
+ 0.33333333 * (-2.0 + 4.0 * cosio + 6.0 * theta2));
|
|
let f523 = sinio
|
|
* (4.92187512 * sini2o * (-2.0 - 4.0 * cosio + 10.0 * theta2)
|
|
+ 6.56250012 * (1.0 + 2.0 * cosio - 3.0 * theta2));
|
|
let f542 = 29.53125
|
|
* sinio
|
|
* (2.0 - 8.0 * cosio + theta2 * (-12.0 + 8.0 * cosio + 10.0 * theta2));
|
|
let f543 = 29.53125
|
|
* sinio
|
|
* (-2.0 - 8.0 * cosio + theta2 * (12.0 + 8.0 * cosio - 10.0 * theta2));
|
|
|
|
let _ao = ao_from_no(no_kozai);
|
|
|
|
d2201 = c1ss * f220 * g201;
|
|
d2211 = c1ss * f221 * g211;
|
|
d3210 = c1ss * f321 * g310;
|
|
d3222 = c1ss * f322 * g322;
|
|
d4410 = c1ss * f441 * g410;
|
|
d4422 = c1ss * f442 * g422;
|
|
d5220 = c1ss * f522 * g520;
|
|
d5232 = c1ss * f523 * g532;
|
|
d5421 = c1ss * f542 * g521;
|
|
d5433 = c1ss * f543 * g533;
|
|
|
|
xlamo = (mo + nodeo + nodeo - gsto - gsto).rem_euclid(TWOPI);
|
|
xfact = mo + argpo + nodeo + nodeo - TWOPI / MINUTES_PER_DAY * 2.0 - gsto - gsto;
|
|
} else {
|
|
resonance = ResonanceType::None;
|
|
}
|
|
|
|
// Compute secular rates
|
|
let dedt = se2 * zmos.sin() + sl2 * (2.0 * zmos).sin();
|
|
let didt = si2 * zmos.sin();
|
|
let dmdt = sl2 * zmos.cos();
|
|
let dnodt = sh2 * zmos.sin();
|
|
let domdt = sgh2 * zmos.sin();
|
|
|
|
DeepSpaceState {
|
|
se2, se3, si2, si3, sl2, sl3, sl4,
|
|
sgh2, sgh3, sgh4, sh2, sh3,
|
|
ee2, e3, xi2, xi3, xl2, xl3, xl4,
|
|
xgh2, xgh3, xgh4, xh2, xh3,
|
|
zmol, zmos,
|
|
resonance, del1, del2, del3,
|
|
xfact, xlamo,
|
|
d2201, d2211, d3210, d3222,
|
|
d4410, d4422, d5220, d5232, d5421, d5433,
|
|
dedt, didt, dmdt, dnodt, domdt,
|
|
atime: 0.0,
|
|
xli: if resonance != ResonanceType::None { xlamo } else { 0.0 },
|
|
xni: no_kozai,
|
|
gsto,
|
|
}
|
|
}
|
|
|
|
/// Compute solar or lunar perturbation coefficients.
|
|
///
|
|
/// `is_solar`: true for solar, false for lunar.
|
|
fn solar_lunar_coefficients(
|
|
inclo: f64,
|
|
nodeo: f64,
|
|
argpo: f64,
|
|
ecco: f64,
|
|
is_solar: bool,
|
|
_zm: f64,
|
|
_ze: f64,
|
|
) -> (f64, f64, f64, f64, f64, f64, f64, f64, f64, f64, f64, f64) {
|
|
let cosio = inclo.cos();
|
|
let sinio = inclo.sin();
|
|
|
|
// Third-body inclination and node
|
|
let (sinx, _cosx) = if is_solar {
|
|
(0.39785416_f64.sin(), 0.39785416_f64.cos())
|
|
} else {
|
|
(0.08979719_f64.sin(), 0.08979719_f64.cos())
|
|
};
|
|
|
|
let cc = if is_solar { 2.9864797e-6 } else { 4.7968065e-7 };
|
|
|
|
let a9 = if is_solar { 1.19459e-5 } else { 1.5835218e-4 };
|
|
let a10 = a9;
|
|
|
|
// Simplified perturbation coefficients
|
|
let sinzf = argpo.sin();
|
|
let coszf = argpo.cos();
|
|
let sinzy = nodeo.sin();
|
|
let coszy = nodeo.cos();
|
|
|
|
let ses = cc * sinx * (a9 + a10 * (sinio * sinzy * ecco));
|
|
let sis = cc * (a9 * cosio * sinzy);
|
|
let sls = cc * (a9 * (sinio * coszf - cosio * sinzf * coszy)
|
|
+ a10 * sinio * ecco * coszy);
|
|
|
|
let sghs = cc * cosio * (a9 * sinzf + a10 * ecco * coszf);
|
|
let shs = cc * sinio * (a9 * coszy - a10 * ecco * sinzy);
|
|
|
|
let se2 = 2.0 * ses;
|
|
let se3 = 0.0;
|
|
let si2 = 2.0 * sis;
|
|
let si3 = 0.0;
|
|
let sl2 = 2.0 * sls;
|
|
let sl3 = 0.0;
|
|
let sl4 = 0.0;
|
|
let sgh2 = 2.0 * sghs;
|
|
let sgh3 = 0.0;
|
|
let sgh4 = 0.0;
|
|
let sh2 = 2.0 * shs;
|
|
let sh3 = 0.0;
|
|
|
|
(se2, se3, si2, si3, sl2, sl3, sl4, sgh2, sgh3, sgh4, sh2, sh3)
|
|
}
|
|
|
|
/// Apply deep-space secular perturbations.
|
|
///
|
|
/// Returns updated (mp, node, argp, eccentricity, inclination, mean_motion).
|
|
#[allow(clippy::too_many_arguments)] // SGP4/SDP4 deep-space secular update
|
|
pub(crate) fn secular(
|
|
&self,
|
|
tsince: f64,
|
|
mp: f64,
|
|
node: f64,
|
|
argp: f64,
|
|
ecco: f64,
|
|
inclo: f64,
|
|
no_kozai: f64,
|
|
) -> (f64, f64, f64, f64, f64, f64) {
|
|
// Solar and lunar mean anomalies at current time
|
|
let zms = (self.zmos + 1.19459e-5 * tsince).rem_euclid(TWOPI);
|
|
let zml = (self.zmol + 1.5835218e-4 * tsince).rem_euclid(TWOPI);
|
|
|
|
// Solar perturbations
|
|
let mut em = ecco
|
|
+ self.se2 * zms.sin()
|
|
+ self.se3 * (2.0 * zms).sin();
|
|
let mut xinc = inclo
|
|
+ self.si2 * zms.sin()
|
|
+ self.si3 * (2.0 * zms).sin();
|
|
let mut xm = mp
|
|
+ self.sl2 * zms.sin()
|
|
+ self.sl3 * (2.0 * zms).sin()
|
|
+ self.sl4 * (3.0 * zms).sin();
|
|
let mut xn = node
|
|
+ self.sh2 * zms.sin()
|
|
+ self.sh3 * (2.0 * zms).sin();
|
|
let mut xg = argp
|
|
+ self.sgh2 * zms.sin()
|
|
+ self.sgh3 * (2.0 * zms).sin()
|
|
+ self.sgh4 * (3.0 * zms).sin();
|
|
|
|
// Lunar perturbations
|
|
em += self.ee2 * zml.sin() + self.e3 * (2.0 * zml).sin();
|
|
xinc += self.xi2 * zml.sin() + self.xi3 * (2.0 * zml).sin();
|
|
xm += self.xl2 * zml.sin()
|
|
+ self.xl3 * (2.0 * zml).sin()
|
|
+ self.xl4 * (3.0 * zml).sin();
|
|
xn += self.xh2 * zml.sin() + self.xh3 * (2.0 * zml).sin();
|
|
xg += self.xgh2 * zml.sin()
|
|
+ self.xgh3 * (2.0 * zml).sin()
|
|
+ self.xgh4 * (3.0 * zml).sin();
|
|
|
|
// Resonance effects
|
|
let mut nm = no_kozai;
|
|
if self.resonance != ResonanceType::None {
|
|
nm = self.resonance_secular(tsince, nm);
|
|
}
|
|
|
|
(xm, xn, xg, em, xinc, nm)
|
|
}
|
|
|
|
/// Secular resonance integration for 12-hour and 24-hour orbits.
|
|
fn resonance_secular(&self, tsince: f64, no_kozai: f64) -> f64 {
|
|
// Simplified resonance: adjust mean motion based on resonance terms
|
|
match self.resonance {
|
|
ResonanceType::Synchronous => {
|
|
let fasx2 = 0.13130908;
|
|
let fasx4 = 2.8843198;
|
|
let fasx6 = 0.37448087;
|
|
|
|
let xli = self.xlamo + self.xfact * tsince;
|
|
|
|
// Apply synchronous resonance correction
|
|
let delt = self.del1 * (xli - fasx2).sin()
|
|
+ self.del2 * (2.0 * (xli - fasx4)).sin()
|
|
+ self.del3 * (3.0 * (xli - fasx6)).sin();
|
|
|
|
no_kozai + delt
|
|
}
|
|
ResonanceType::HalfDay => {
|
|
let g22 = 5.7686396;
|
|
let g32 = 0.95240898;
|
|
let g44 = 1.8014998;
|
|
let g52 = 1.0508330;
|
|
let g54 = 4.4108898;
|
|
|
|
let xli = self.xlamo + self.xfact * tsince;
|
|
|
|
// Apply half-day resonance correction
|
|
let delt = self.d2201 * (2.0 * xli - g22).sin()
|
|
+ self.d2211 * xli.sin()
|
|
+ self.d3210 * (xli + g32).sin()
|
|
+ self.d3222 * (-xli + g32).sin()
|
|
+ self.d4410 * (2.0 * xli + g44).sin()
|
|
+ self.d4422 * (2.0 * xli).sin()
|
|
+ self.d5220 * (xli - g52).sin()
|
|
+ self.d5232 * (-xli + g52).sin()
|
|
+ self.d5421 * (xli + g54).sin()
|
|
+ self.d5433 * (-xli + g54).sin();
|
|
|
|
no_kozai + delt
|
|
}
|
|
ResonanceType::None => no_kozai,
|
|
}
|
|
}
|
|
|
|
/// Apply deep-space periodic perturbations.
|
|
///
|
|
/// Returns updated (ep, xincp, argpp, nodep, mp).
|
|
pub(crate) fn periodics(
|
|
&self,
|
|
tsince: f64,
|
|
ep: f64,
|
|
xincp: f64,
|
|
argpp: f64,
|
|
nodep: f64,
|
|
mp: f64,
|
|
) -> (f64, f64, f64, f64, f64) {
|
|
// Solar and lunar mean anomalies
|
|
let zm = (self.zmos + 1.19459e-5 * tsince).rem_euclid(TWOPI);
|
|
let zf = (zm + 2.0 * (-0.01675) * zm.sin()).rem_euclid(TWOPI);
|
|
let sinzf = zf.sin();
|
|
let _coszf = zf.cos();
|
|
|
|
let zml = (self.zmol + 1.5835218e-4 * tsince).rem_euclid(TWOPI);
|
|
let sinzml = zml.sin();
|
|
|
|
// Compute periodic corrections
|
|
let mut pe = self.se2 * sinzf + self.se3 * (2.0 * zf).sin();
|
|
let mut pinc = self.si2 * sinzf;
|
|
let mut pl = self.sl2 * sinzf + self.sl3 * (2.0 * zf).sin();
|
|
let mut pgh = self.sgh2 * sinzf + self.sgh3 * (2.0 * zf).sin();
|
|
let mut ph = self.sh2 * sinzf;
|
|
|
|
// Add lunar periodic corrections
|
|
pe += self.ee2 * sinzml;
|
|
pinc += self.xi2 * sinzml;
|
|
pl += self.xl2 * sinzml;
|
|
pgh += self.xgh2 * sinzml;
|
|
ph += self.xh2 * sinzml;
|
|
|
|
let ep_out = ep + pe;
|
|
let xincp_out = xincp + pinc;
|
|
let mut argpp_out = argpp + pgh;
|
|
let mp_out = mp + pl;
|
|
|
|
// Apply inclination-dependent node correction
|
|
let sinip = xincp_out.sin();
|
|
let cosip = xincp_out.cos();
|
|
|
|
if xincp_out.abs() > 0.2 {
|
|
// Standard correction
|
|
let nodep_out = nodep + ph / sinip;
|
|
argpp_out -= cosip * ph / sinip;
|
|
(ep_out, xincp_out, argpp_out, nodep_out, mp_out)
|
|
} else {
|
|
// Near-equatorial: use sin/cos of node to avoid singularity
|
|
let sinop = nodep.sin();
|
|
let cosop = nodep.cos();
|
|
|
|
let alfdp = sinip * sinop + ph * cosop;
|
|
let betdp = sinip * cosop - ph * sinop;
|
|
let nodep_out = alfdp.atan2(betdp);
|
|
|
|
argpp_out += -cosip * (nodep_out - nodep);
|
|
|
|
(ep_out, xincp_out, argpp_out, nodep_out, mp_out)
|
|
}
|
|
}
|
|
}
|
|
|
|
/// Compute semi-major axis from mean motion (earth radii).
|
|
fn ao_from_no(no_kozai: f64) -> f64 {
|
|
(XKE / no_kozai).powf(2.0 / 3.0)
|
|
}
|
|
|
|
/// GMST from Julian date (for deep-space epoch initialization).
|
|
fn gmst_from_jd(jd: f64) -> f64 {
|
|
let t_ut1 = (jd - 2451545.0) / 36525.0;
|
|
let gmst_sec = 67310.54841
|
|
+ (876600.0 * 3600.0 + 8640184.812866) * t_ut1
|
|
+ 0.093104 * t_ut1 * t_ut1
|
|
- 6.2e-6 * t_ut1 * t_ut1 * t_ut1;
|
|
let gmst_rad = (gmst_sec / 86400.0).rem_euclid(1.0) * TWOPI;
|
|
gmst_rad.rem_euclid(TWOPI)
|
|
}
|
|
|
|
#[cfg(test)]
|
|
mod tests {
|
|
use super::*;
|
|
use std::f64::consts::PI;
|
|
|
|
#[test]
|
|
fn test_gmst_from_jd() {
|
|
// At J2000.0, GMST ~ 280.46 degrees
|
|
let gmst = gmst_from_jd(2451545.0);
|
|
let deg = gmst * 180.0 / PI;
|
|
assert!((deg - 280.46).abs() < 1.0, "GMST at J2000: {deg}");
|
|
}
|
|
|
|
#[test]
|
|
fn test_ao_from_no() {
|
|
// For no_kozai = XKE, semi-major axis = 1 earth radius
|
|
let ao = ao_from_no(XKE);
|
|
assert!((ao - 1.0).abs() < 1e-10);
|
|
}
|
|
}
|