satkit 0.22.0

//! Orbit Propagation Settings

use crate::earthgravity::GravityModel;
use crate::orbitprop::tides::TideModel;
use crate::orbitprop::Precomputed;
use crate::TimeLike;

use super::error::{Error, Result};

/// Choice of ODE integrator for orbit propagation
#[derive(Debug, Clone, Copy, PartialEq, Eq, serde::Serialize, serde::Deserialize, Default)]
pub enum Integrator {
    /// Verner 9(8) with 8th-degree dense output, 21 stages (16 + 5 for the
    /// interpolant; default). With `enable_interp = false` the propagator
    /// runs the 16-stage [`Integrator::RKV98NoInterp`] tableau instead —
    /// same order and error control, 24% fewer force evaluations per step.
    #[default]
    RKV98,
    /// Verner 9(8) without interpolation, 16 stages. Selecting it explicitly
    /// is equivalent to `RKV98` with `enable_interp = false`.
    RKV98NoInterp,
    /// Verner 8(7) with 7th-degree dense output, 17 stages (13 + 4 for the interpolant)
    RKV87,
    /// Verner 6(5) with 6th-degree dense output, 10 stages
    RKV65,
    /// Tsitouras 5(4) with FSAL, 7 stages
    RKTS54,
    /// RODAS4 — L-stable Rosenbrock 4(3), 6 stages. For stiff problems (re-entry, low perigee).
    /// Does not support state transition matrix propagation or dense output interpolation.
    RODAS4,
    /// Gauss-Jackson 8 — 8th-order fixed-step multistep predictor-corrector
    /// specialised for 2nd-order ODEs (Berry & Healy 2004). The dominant
    /// integrator in high-precision astrodynamics codes (GMAT, STK, ODTK).
    /// Typically uses 3-10× fewer force evaluations than a Runge-Kutta method
    /// of comparable accuracy on smooth orbit propagation problems.
    ///
    /// Uses a fixed step size set via [`PropSettings::gj_step_seconds`].
    /// Supports dense output (quintic Hermite interpolation between steps
    /// when `enable_interp` is set) but not state transition matrix
    /// propagation (C=7). Not recommended for highly eccentric orbits or
    /// integration across discontinuities (eclipse boundaries, maneuvers).
    GaussJackson8,
}

impl std::fmt::Display for Integrator {
    fn fmt(&self, f: &mut std::fmt::Formatter) -> std::fmt::Result {
        match self {
            Self::RKV98 => write!(f, "RKV98 (9th order, 21 stages)"),
            Self::RKV98NoInterp => write!(f, "RKV98NoInterp (9th order, 16 stages)"),
            Self::RKV87 => write!(f, "RKV87 (8th order, 17 stages)"),
            Self::RKV65 => write!(f, "RKV65 (6th order, 10 stages)"),
            Self::RKTS54 => write!(f, "RKTS54 (5th order, 7 stages, FSAL)"),
            Self::RODAS4 => write!(f, "RODAS4 (4th order, 6 stages, L-stable)"),
            Self::GaussJackson8 => write!(f, "Gauss-Jackson 8 (8th order, fixed-step multistep)"),
        }
    }
}

/// Propagation settings
///
/// These include
///
/// * `gravity_degree` - maximum degree of spherical harmonic gravity model.  Default is 4
/// * `gravity_order` - maximum order of spherical harmonic gravity model.  Default is same as `gravity_degree`.
///   Must be ≤ `gravity_degree`.
/// * `gravity_model` - gravity model to use.  Default is EGM96.  Options: EGM96, JGM3, JGM2, ITUGrace16
/// * `abs_error` - the maximum absolute error for the infinity norm of the state in Runge-Kutta integrator.  Default is 1e-8
/// * `rel_error` - the maximum relative error for the infinity norm of the state in Runge-Kutta integrator.  Default is 1e-8
/// * `use_spaceweather` -  Do we use space weather when computing the atmospheric density.  Default is true
/// * `use_sun_gravity` - Do we include sun third-body gravitational perturbation.  Default is true
/// * `use_moon_gravity` - Do we include moon third-body gravitational perturbation.  Default is true
/// * `tide_model` - Solid Earth tide model. Default is `TideModel::SolidStep1`
///   (IERS 2010 §6.2.1 frequency-independent Love-number response).
/// * `use_relativistic_correction` - Include the Schwarzschild
///   post-Newtonian acceleration (IERS 2010 §10.3 Eq. 10.12). Default is true.
///   Its position effect depends on the orbit, propagation arc, and fitted
///   parameters; computational cost is negligible.
/// * `enable_interp` - Do we enable interpolation of the state between begin and end times.  Default is true.
///   Setting it false skips storing dense output and, with the default `RKV98`, switches to the
///   16-stage tableau (24% fewer force evaluations per step).
/// * `integrator` - which Runge-Kutta integrator to use.  Default is RKV98
/// * `max_steps` - maximum number of integrator steps before the propagator
///   aborts with [`numeris::ode::OdeError::MaxStepsExceeded`] (adaptive
///   solvers) or a Gauss-Jackson max-steps error. Default is 1_000_000,
///   which covers very long propagation arcs (e.g., ~700 days of GJ8 at
///   60 s step) with plenty of headroom. Lower if you want a tighter
///   runaway-propagation safeguard.
/// * `initial_step_secs` - first step (seconds) the adaptive integrators
///   attempt. Default `None`: derived from the initial state, the tolerances
///   and the integrator order (about 170 s for RKV98 at 1e-9 in LEO). Set it
///   to warm-start from a previous arc's [`PropagationResult::next_step_secs`].
///
#[derive(Debug, Clone, serde::Serialize, serde::Deserialize)]
pub struct PropSettings {
    pub gravity_degree: u16,
    pub gravity_order: u16,
    pub gravity_model: GravityModel,
    pub abs_error: f64,
    pub rel_error: f64,
    pub use_spaceweather: bool,
    pub use_sun_gravity: bool,
    pub use_moon_gravity: bool,
    pub tide_model: TideModel,
    pub use_relativistic_correction: bool,
    /// Store dense output so [`PropagationResult::interp`] works between
    /// the begin and end times. Default `true`. When `false`, no dense
    /// output is stored and [`Integrator::RKV98`] runs its 16-stage
    /// no-interpolant tableau (same order and error control, 24% fewer
    /// force evaluations per step).
    pub enable_interp: bool,
    pub integrator: Integrator,
    /// Fixed step size (seconds) used by [`Integrator::GaussJackson8`].
    /// Ignored by adaptive integrators. Typical values: 30-120 s for LEO,
    /// 60-300 s for MEO, 300-600 s for GEO. Default: 60 s.
    pub gj_step_seconds: f64,
    /// Maximum number of integrator steps before the propagator aborts
    /// with a max-steps error. Applies to all integrators (adaptive
    /// Runge-Kutta / Rosenbrock via [`numeris::ode::AdaptiveSettings`] and
    /// Gauss-Jackson 8 via its own settings). Default: 1_000_000.
    pub max_steps: usize,
    /// Fail a propagation whose span (plus integrator padding) extends past
    /// the end of the loaded Earth-orientation-parameter (EOP) table with
    /// [`Error::EopCoverage`](super::Error::EopCoverage), instead of holding
    /// the last EOP row constant with a one-time warning. Off by default —
    /// propagating into next year is a legitimate use and constant
    /// extrapolation is a reasonable best effort (polar motion and UT1−UTC
    /// drift by ~0.1″ / ~10 ms over a few months, i.e. metres at LEO) —
    /// but set it for precision work and refresh the data files
    /// (`satkit::utils::update_datafiles()`) when it trips. See
    /// [`crate::earth_orientation_params::coverage`]. Default: `false`.
    #[serde(default)]
    pub require_eop_coverage: bool,
    /// First step (seconds) the adaptive integrators attempt. `None` (the
    /// default) derives it from the initial state and the tolerances as
    /// `1.5 · |r|/|v| · tol^(1/(p+1))`, with `p` the integrator's order and
    /// `tol = rel_error + abs_error/|r|` — within a factor of ~2.5 of the
    /// settled stride across the Runge-Kutta integrators from 1e-6 to 1e-12
    /// (about 170 s for RKV98 at 1e-9 in LEO), which the step controller
    /// closes within a step or two. The
    /// integrator's own starting-step heuristic is not used: it is scale
    /// sensitive and, for an orbit in metres and seconds, starts several
    /// orders of magnitude below the working step and spends a dozen steps
    /// (half the derivative evaluations of a one-hour arc at 1e-9
    /// tolerance) growing into it.
    ///
    /// Set it explicitly to warm-start a follow-on arc from the previous
    /// arc's [`PropagationResult::next_step_secs`], which continues at full
    /// stride, or to override the state-derived default. It is a magnitude:
    /// backward propagation applies the sign, and a value longer than the
    /// arc is clamped to it. Ignored by [`Integrator::GaussJackson8`] (fixed
    /// step). Zero or non-finite values make `propagate` fail with
    /// [`numeris::ode::OdeError::InvalidInitialStep`].
    #[serde(default)]
    pub initial_step_secs: Option<f64>,
    /// Regenerable ephemeris/EOP cache; excluded from serialization (a
    /// deserialized `PropSettings` recomputes it lazily as needed).
    #[serde(skip)]
    pub precomputed: Option<Precomputed>,
}

impl Default for PropSettings {
    fn default() -> Self {
        Self {
            gravity_degree: 4,
            gravity_order: 4,
            gravity_model: GravityModel::EGM96,
            abs_error: 1e-8,
            rel_error: 1e-8,
            use_spaceweather: true,
            use_sun_gravity: true,
            use_moon_gravity: true,
            tide_model: TideModel::default(),
            use_relativistic_correction: true,
            enable_interp: true,
            integrator: Integrator::default(),
            gj_step_seconds: 60.0,
            max_steps: 1_000_000,
            require_eop_coverage: false,
            initial_step_secs: None,
            precomputed: None,
        }
    }
}

impl PropSettings {
    /// Set gravity degree and order, with validation
    ///
    /// # Arguments
    /// * `degree` - Maximum degree of spherical harmonic gravity model
    /// * `order` - Maximum order (must be ≤ degree)
    ///
    /// # Errors
    /// Returns error if order > degree, or if degree exceeds
    /// [`MAX_GRAVITY_DEGREE`](crate::earthgravity::MAX_GRAVITY_DEGREE) (40).
    pub fn set_gravity(&mut self, degree: u16, order: u16) -> Result<()> {
        Self::check_gravity(degree, order)?;
        self.gravity_degree = degree;
        self.gravity_order = order;
        Ok(())
    }

    /// Validate the gravity degree / order currently held by these settings.
    ///
    /// [`propagate`](super::propagate) calls this on entry so that settings
    /// built by struct literal (bypassing [`set_gravity`](Self::set_gravity))
    /// are checked too: a degree above
    /// [`MAX_GRAVITY_DEGREE`](crate::earthgravity::MAX_GRAVITY_DEGREE) used to
    /// be silently evaluated at 40, and an order above the degree silently
    /// clamped.
    pub fn validate_gravity(&self) -> Result<()> {
        Self::check_gravity(self.gravity_degree, self.gravity_order)
    }

    fn check_gravity(degree: u16, order: u16) -> Result<()> {
        use crate::earthgravity::MAX_GRAVITY_DEGREE;
        if degree > MAX_GRAVITY_DEGREE {
            return Err(Error::InvalidGravityDegree {
                degree,
                max: MAX_GRAVITY_DEGREE,
            });
        }
        if order > degree {
            return Err(Error::InvalidGravityOrder { order, degree });
        }
        Ok(())
    }

    /// Compute the required padding (in seconds) beyond the nominal
    /// propagation interval when building a [`Precomputed`] interp table
    /// for these settings.
    ///
    /// Most integrators only evaluate the force within the nominal
    /// interval and need the default padding. The
    /// [`Integrator::GaussJackson8`] integrator, however, runs a
    /// symmetric ±4·h_gj startup around the starting epoch — requiring
    /// the interp table to cover times up to `4·gj_step_seconds` outside
    /// the interval on the startup side. A small safety margin is added
    /// to guard against floating-point round-off.
    pub fn required_precompute_padding(&self) -> f64 {
        match self.integrator {
            Integrator::GaussJackson8 => {
                // 4 startup steps backward from epoch + safety margin
                (4.0 * self.gj_step_seconds.abs() + 60.0)
                    .max(crate::orbitprop::precomputed::DEFAULT_PADDING_SECS)
            }
            Integrator::RKV98
            | Integrator::RKV98NoInterp
            | Integrator::RKV87
            | Integrator::RKV65
            | Integrator::RKTS54
            | Integrator::RODAS4 => crate::orbitprop::precomputed::DEFAULT_PADDING_SECS,
        }
    }

    /// Precompute terms between begin and end instants
    ///
    /// # Arguments
    /// * `begin` - Begin instant
    /// * `end` - End instant
    ///
    /// Pre-computes inertial to earth-fixed rotation vector (used for Earth gravity calculation),
    /// sun, and moon positions between the begin and end instants.  These are used in the
    /// force model when propagating orbits
    ///
    /// Pre-computing these terms means the settings can be used for multiple propagations
    /// between the same begin and end instants without needing to recompute these terms each time.
    /// (significant speedup when propagating many orbits over the same time span)
    ///
    /// The precomputed range is automatically padded to accommodate the
    /// selected integrator. For [`Integrator::GaussJackson8`] the padding
    /// is extended to `4·gj_step_seconds + 60 s` on each end to cover the
    /// backward RK4 startup stencil; for other integrators a default
    /// 240 s padding is used.
    ///
    /// # Errors
    /// Returns error if precomputation fails
    ///
    /// # Example
    /// ```
    /// use satkit::orbitprop::PropSettings;
    /// use satkit::Instant;
    ///
    /// let begin = Instant::now();
    /// let end = begin + satkit::Duration::from_hours(1.0);
    /// let mut props = PropSettings::default();
    /// props.precompute_terms(&begin, &end).unwrap();
    ///
    /// ```
    pub fn precompute_terms<T: TimeLike>(&mut self, begin: &T, end: &T) -> Result<()> {
        let padding = self.required_precompute_padding();
        self.precomputed = Some(Precomputed::new_padded(begin, end, 60.0, padding)?);
        Ok(())
    }

    pub fn precompute_terms_with_step<T: TimeLike>(
        &mut self,
        begin: &T,
        end: &T,
        step_secs: f64,
    ) -> Result<()> {
        let padding = self.required_precompute_padding();
        self.precomputed = Some(Precomputed::new_padded(begin, end, step_secs, padding)?);
        Ok(())
    }
}

impl std::fmt::Display for PropSettings {
    fn fmt(&self, f: &mut std::fmt::Formatter) -> std::fmt::Result {
        write!(
            f,
            r#"Orbit Propagation Settings
            Gravity Degree: {},
            Gravity Order: {},
            Gravity Model: {},
            Max Abs Error: {:e},
            Max Rel Error: {:e},
            Space Weather: {},
            Sun Gravity: {},
            Moon Gravity: {},
            Tide Model: {:?},
            GR (IERS 10.12): {},
            Interpolation: {},
            Integrator: {},
            Max Steps: {},
            Require EOP Coverage: {},
            Initial Step: {},
            {}"#,
            self.gravity_degree,
            self.gravity_order,
            self.gravity_model,
            self.abs_error,
            self.rel_error,
            self.use_spaceweather,
            self.use_sun_gravity,
            self.use_moon_gravity,
            self.tide_model,
            self.use_relativistic_correction,
            self.enable_interp,
            self.integrator,
            self.max_steps,
            self.require_eop_coverage,
            self.initial_step_secs.map_or_else(
                || "auto (state, tolerance, order)".to_string(),
                |h| format!("{h} s")
            ),
            self.precomputed.as_ref().map_or_else(
                || "No Precomputed".to_string(),
                |p| format!("Precomputed: {} to {}", p.begin, p.end)
            )
        )
    }
}

#[cfg(test)]
mod test {
    use super::*;

    #[test]
    fn testdisplay() {
        let props = PropSettings::default();
        println!("props = {}", props);
    }

    #[test]
    fn set_gravity_rejects_degree_above_max() {
        let mut s = PropSettings::default();
        assert!(s.set_gravity(40, 40).is_ok());
        assert_eq!((s.gravity_degree, s.gravity_order), (40, 40));
        assert!(matches!(
            s.set_gravity(41, 41),
            Err(Error::InvalidGravityDegree {
                degree: 41,
                max: 40
            })
        ));
        assert!(matches!(
            s.set_gravity(360, 360),
            Err(Error::InvalidGravityDegree { degree: 360, .. })
        ));
        assert!(matches!(
            s.set_gravity(20, 21),
            Err(Error::InvalidGravityOrder {
                order: 21,
                degree: 20
            })
        ));
        // A rejected call leaves the settings untouched.
        assert_eq!((s.gravity_degree, s.gravity_order), (40, 40));
    }

    #[test]
    fn validate_gravity_catches_struct_literal() {
        let s = PropSettings {
            gravity_degree: 41,
            gravity_order: 41,
            ..Default::default()
        };
        assert!(matches!(
            s.validate_gravity(),
            Err(Error::InvalidGravityDegree { .. })
        ));
        let s = PropSettings {
            gravity_degree: 8,
            gravity_order: 9,
            ..Default::default()
        };
        assert!(matches!(
            s.validate_gravity(),
            Err(Error::InvalidGravityOrder { .. })
        ));
        assert!(PropSettings::default().validate_gravity().is_ok());
    }
}