kshana 0.37.0

// SPDX-License-Identifier: AGPL-3.0-only
//! Launch-window & ascent geometry: the two-body relations a mission analyst uses
//! at the trade-study tier to turn a launch-site latitude and a target orbit into
//! launch azimuth(s), the minimum reachable inclination, the Earth-rotation
//! velocity bonus, plane-change (dogleg) Δv, and the number of daily launch
//! opportunities.
//!
//! HONEST SCOPE (MODELLED): spherical-Earth two-body geometry. The launch azimuth
//! is the geometric relation `sin(Az) = cos(i)/cos(lat)`; it does **not** apply the
//! launch-velocity-triangle correction for the rotating-Earth boost (that boost is
//! reported separately as `site_rotation_speed` / the eastward Δv saving, for the
//! analyst to fold in). No ascent trajectory, drag-loss, or steering model — this
//! is pre-Phase-A geometry, not a trajectory optimiser.

use crate::forces::EARTH_ROTATION_RATE;
use crate::orbit::{MU_EARTH, R_EARTH_EQUATORIAL_M};
use crate::portable_math::PortableFloat;
use serde::Deserialize;

/// Circular orbital speed (m/s) at geometric `altitude_m` above the equatorial
/// radius: `sqrt(mu / r)`.
pub fn circular_velocity(altitude_m: f64) -> f64 {
    (MU_EARTH / (R_EARTH_EQUATORIAL_M + altitude_m)).sqrt()
}

/// The minimum orbital inclination directly reachable from a site at geodetic
/// `lat_rad` without a plane-change dogleg: `|lat|`.
pub fn min_inclination(lat_rad: f64) -> f64 {
    lat_rad.abs()
}

/// Eastward surface speed (m/s) of a launch site at `lat_rad` from Earth's
/// rotation: `omega · R_eq · cos(lat)` — the free velocity a posigrade (easterly)
/// launch starts with (≈465 m/s at the equator).
pub fn site_rotation_speed(lat_rad: f64) -> f64 {
    EARTH_ROTATION_RATE * R_EARTH_EQUATORIAL_M * lat_rad.pcos()
}

/// Surface speed (m/s) of an Earth-fixed site at geodetic `lat_rad`, `lon_rad` on the
/// spherical Earth (radius `R_eq`), at TT Julian Date `jd_tt`, rotating about the true pole:
/// `ω_eff · R_eq · |ẑ_CIP × r̂|` with `ω_eff = Ω (1 − LOD / 86400 s)` and the Celestial
/// Intermediate Pole placed in the terrestrial frame at `(x_p, −y_p, 1)` from the IERS
/// polar motion. `eop` supplies `x_p`, `y_p` and the excess length of day `LOD`. Unlike
/// [`site_rotation_speed`], which rotates about the terrestrial z axis at the nominal rate,
/// this carries the polar-motion tilt (amplified by `tan(lat)`) and the day-length change.
/// The precession-nutation rate of the pole in space (about 1e-11 rad/s) is left out.
pub fn site_rotation_speed_at(
    lat_rad: f64,
    lon_rad: f64,
    jd_tt: f64,
    eop: &crate::eop::EopSeries,
) -> f64 {
    let (_, xp, yp) = eop.frame_args_tt(jd_tt);
    let lod_s = eop.lod_ms_tt(jd_tt) * 1e-3;
    let omega = EARTH_ROTATION_RATE * (1.0 - lod_s / 86_400.0);
    let n = (xp * xp + yp * yp + 1.0).sqrt();
    let pole = [xp / n, -yp / n, 1.0 / n];
    let site = [
        lat_rad.pcos() * lon_rad.pcos(),
        lat_rad.pcos() * lon_rad.psin(),
        lat_rad.psin(),
    ];
    let d = pole[0] * site[0] + pole[1] * site[1] + pole[2] * site[2];
    omega * R_EARTH_EQUATORIAL_M * (1.0 - d * d).max(0.0).sqrt()
}

/// Plane-change Δv (m/s) to rotate a velocity of magnitude `v_orbit` through
/// `delta_i_rad`: the vector relation `2·v·sin(Δi/2)`.
pub fn plane_change_dv(v_orbit: f64, delta_i_rad: f64) -> f64 {
    2.0 * v_orbit * (delta_i_rad.abs() / 2.0).psin()
}

/// Geometric launch azimuth(s) (rad clockwise from north, in `[0, 2π)`) to reach
/// orbital `inclination_rad` from a site at geodetic `lat_rad`, two-body:
/// `sin(Az) = cos(i)/cos(lat)`. Returns `(ascending, descending)` azimuths
/// (the descending pass is `π − ascending`). Errors when `i < |lat|` (the target
/// inclination is unreachable without a dogleg) or `i > π − |lat|`.
pub fn launch_azimuth(lat_rad: f64, inclination_rad: f64) -> Result<(f64, f64), String> {
    let cl = lat_rad.pcos();
    if cl.abs() < 1e-12 {
        return Err("launch site at the pole has no defined azimuth".to_string());
    }
    let s = inclination_rad.pcos() / cl;
    if !(-1.0..=1.0).contains(&s) {
        return Err(format!(
            "inclination {:.2}° unreachable from latitude {:.2}° without a dogleg \
             (need |lat| ≤ i ≤ 180−|lat|)",
            inclination_rad.to_degrees(),
            lat_rad.to_degrees()
        ));
    }
    let asc = s.pasin(); // in [-π/2, π/2]; for prograde i this is the NE/SE azimuth
    let asc = asc.rem_euclid(std::f64::consts::TAU);
    let desc = (std::f64::consts::PI - s.pasin()).rem_euclid(std::f64::consts::TAU);
    Ok((asc, desc))
}

/// Number of daily launch opportunities to a given inclination from a site at
/// `lat_rad`: two (an ascending and a descending node pass) in general, one when
/// the target inclination equals the site latitude (the orbit plane is tangent to
/// the site's latitude circle), and none when it is unreachable.
pub fn daily_launch_opportunities(lat_rad: f64, inclination_rad: f64) -> u8 {
    let l = lat_rad.abs();
    let i = inclination_rad;
    if i < l - 1e-9 || i > std::f64::consts::PI - l + 1e-9 {
        0
    } else if (i - l).abs() < 1e-9 || (i - (std::f64::consts::PI - l)).abs() < 1e-9 {
        1
    } else {
        2
    }
}

fn lw_default_lat() -> f64 {
    28.5
}
fn lw_default_inc() -> f64 {
    51.6
}
fn lw_default_alt() -> f64 {
    400.0
}

/// Unit and provenance class for every numeric field the `launch-window` report emits.
/// The closed-form rows are the two-body relations documented on the free functions
/// above, each checkable by hand from the site latitude and the target orbit.
const UNITS: &[crate::field_schema::FieldUnit] = {
    use crate::field_schema::{FieldUnit, ProvenanceClass::*};
    &[
        FieldUnit {
            path: "site_lat_deg",
            unit: "deg",
            provenance: Input,
            definition: "launch-site geodetic latitude",
        },
        FieldUnit {
            path: "target_inclination_deg",
            unit: "deg",
            provenance: Input,
            definition: "target orbital inclination",
        },
        FieldUnit {
            path: "altitude_km",
            unit: "km",
            provenance: Input,
            definition: "target circular-orbit geometric altitude above the equatorial radius",
        },
        FieldUnit {
            path: "min_inclination_deg",
            unit: "deg",
            provenance: ClosedForm,
            definition: "|site latitude| — the lowest inclination a direct launch reaches",
        },
        FieldUnit {
            path: "circular_velocity_m_s",
            unit: "m/s",
            provenance: ClosedForm,
            definition: "sqrt(mu_Earth / (R_eq + altitude))",
        },
        FieldUnit {
            path: "site_rotation_speed_m_s",
            unit: "m/s",
            provenance: ClosedForm,
            definition: "omega_Earth * R_eq * cos(lat): the eastward speed a posigrade \
                         launch starts with",
        },
        FieldUnit {
            path: "site_lon_deg",
            unit: "deg",
            provenance: Input,
            definition: "launch-site longitude, east; used only by the true-pole site speed",
        },
        FieldUnit {
            path: "site_rotation_speed_true_pole_m_s",
            unit: "m/s",
            provenance: Computed,
            definition: "Omega (1 - LOD/86400 s) * R_eq * |z_CIP x r_site| at the scenario epoch, \
                         the pole and length of day from the supplied IERS finals2000A rows; \
                         null when no epoch is given",
        },
        FieldUnit {
            path: "daily_opportunities",
            unit: "count",
            provenance: Computed,
            definition: "ascending plus descending node passes per day: 2 in general, 1 when \
                         the target inclination equals the site latitude, 0 when unreachable",
        },
        FieldUnit {
            path: "launch_azimuth_deg.ascending",
            unit: "deg",
            provenance: ClosedForm,
            definition: "clockwise from north, from sin(Az) = cos(i)/cos(lat); the whole \
                         object is null when the inclination is unreachable directly",
        },
        FieldUnit {
            path: "launch_azimuth_deg.descending",
            unit: "deg",
            provenance: ClosedForm,
            definition: "the descending-pass mirror, 180 deg - ascending",
        },
        FieldUnit {
            path: "dogleg_plane_change_dv_m_s",
            unit: "m/s",
            provenance: ClosedForm,
            definition: "2 * v_circular * sin(|lat - i| / 2); null whenever a direct \
                         azimuth exists",
        },
    ]
};

/// The `launch-window` scenario: launch azimuth(s), minimum inclination, circular
/// velocity, Earth-rotation bonus, dogleg Δv (if the target is below the site
/// latitude) and daily opportunities for a site latitude + target orbit.
#[derive(Deserialize)]
pub struct LaunchWindowScenario {
    /// Launch-site geodetic latitude (deg).
    #[serde(default = "lw_default_lat")]
    pub site_lat_deg: f64,
    /// Target orbital inclination (deg).
    #[serde(default = "lw_default_inc")]
    pub target_inclination_deg: f64,
    /// Target circular-orbit altitude (km).
    #[serde(default = "lw_default_alt")]
    pub altitude_km: f64,
    /// Launch-site longitude (deg east), used only for the true-pole site speed.
    #[serde(default)]
    pub site_lon_deg: f64,
    /// Launch epoch, UTC, `YYYY-MM-DDTHH:MM:SS`. With [`Self::eop_finals2000a`] it selects
    /// the IERS polar motion and length of day for `site_rotation_speed_true_pole_m_s`.
    #[serde(default)]
    pub epoch: Option<String>,
    /// The body of a real IERS `finals2000A` file covering [`Self::epoch`] (the `--eop`
    /// CLI path inlines it). Required together with `epoch`.
    #[serde(default)]
    pub eop_finals2000a: Option<String>,
}

impl LaunchWindowScenario {
    /// The site speed about the true pole at the scenario epoch ([`site_rotation_speed_at`]),
    /// or `None` when no epoch is given. Errors when only one of `epoch` and
    /// `eop_finals2000a` is given, when the file has no readable row, or when the epoch lies
    /// outside the file's span (the series is never extrapolated here).
    fn true_pole_site_speed(&self, lat_rad: f64) -> Result<Option<f64>, String> {
        let (epoch, body) = match (&self.epoch, &self.eop_finals2000a) {
            (None, None) => return Ok(None),
            (Some(e), Some(b)) => (e, b),
            _ => {
                return Err("the true-pole site speed needs both `epoch` and \
                            `eop_finals2000a` (a real IERS finals2000A body)"
                    .to_string())
            }
        };
        if !(-180.0..=360.0).contains(&self.site_lon_deg) {
            return Err("site_lon_deg must be in [-180, 360]".to_string());
        }
        let jd_utc = crate::leo_pass::parse_epoch_jd(epoch)?;
        let records = crate::eop::parse_all(body);
        let (Some(first), Some(last)) = (records.first(), records.last()) else {
            return Err("eop_finals2000a contained no readable IERS finals2000A rows".to_string());
        };
        let mjd_utc = jd_utc - crate::timescales::MJD_OFFSET;
        if mjd_utc < first.mjd || mjd_utc > last.mjd {
            return Err(format!(
                "epoch {epoch} (MJD {mjd_utc:.3}) is outside the EOP rows supplied \
                 (MJD {} to {})",
                first.mjd, last.mjd
            ));
        }
        let eop = crate::eop::EopSeries::new(records);
        let jd_tt = crate::timescales::utc_to_tt(jd_utc);
        Ok(Some(site_rotation_speed_at(
            lat_rad,
            self.site_lon_deg.to_radians(),
            jd_tt,
            &eop,
        )))
    }
}

impl LaunchWindowScenario {
    /// Run the scenario, returning `(json, summary)`.
    pub fn run_json(&self) -> Result<(String, String), String> {
        if !(-90.0..=90.0).contains(&self.site_lat_deg) {
            return Err("site_lat_deg must be in [-90, 90]".to_string());
        }
        if !(0.0..=180.0).contains(&self.target_inclination_deg) {
            return Err("target_inclination_deg must be in [0, 180]".to_string());
        }
        if !self.altitude_km.is_finite() || self.altitude_km <= 0.0 {
            return Err("altitude_km must be finite and positive".to_string());
        }
        let lat = self.site_lat_deg.to_radians();
        let inc = self.target_inclination_deg.to_radians();
        let alt_m = self.altitude_km * 1000.0;
        let v_orbit = circular_velocity(alt_m);
        let i_min = min_inclination(lat).to_degrees();
        let opportunities = daily_launch_opportunities(lat, inc);
        let true_pole_speed = self.true_pole_site_speed(lat)?;

        // Direct azimuths when reachable; otherwise the dogleg Δv from i_min.
        // A direct azimuth exists exactly when inc ≥ |lat| (cos i / cos lat ∈ [-1, 1]),
        // so the `Ok` arm never needs a dogleg — the dogleg is the `Err` (inc < |lat|) arm.
        let (azimuths, dogleg_dv) = match launch_azimuth(lat, inc) {
            Ok((asc, desc)) => (Some((asc.to_degrees(), desc.to_degrees())), None),
            Err(_) => (None, Some(plane_change_dv(v_orbit, (lat - inc).abs()))),
        };

        let json = serde_json::json!({
            "kind": "launch-window",
            "label": "MODELLED — two-body spherical-Earth launch geometry; azimuth \
                      is the geometric sin(Az)=cos(i)/cos(lat) relation (no rotating-Earth \
                      velocity-triangle correction, no ascent/drag-loss model)",
            "units": crate::field_schema::units_block(UNITS),
            "site_lat_deg": self.site_lat_deg,
            "target_inclination_deg": self.target_inclination_deg,
            "altitude_km": self.altitude_km,
            "min_inclination_deg": i_min,
            "circular_velocity_m_s": v_orbit,
            "site_rotation_speed_m_s": site_rotation_speed(lat),
            "site_lon_deg": self.site_lon_deg,
            "epoch": self.epoch,
            "site_rotation_speed_true_pole_m_s": true_pole_speed,
            "daily_opportunities": opportunities,
            "launch_azimuth_deg": azimuths.map(|(a, d)| serde_json::json!({"ascending": a, "descending": d})),
            "dogleg_plane_change_dv_m_s": dogleg_dv,
        });
        let summary = match azimuths {
            Some((asc, _)) => format!(
                "launch-window: lat {:.1}° -> i {:.1}°: Az {:.1}° asc, v_circ {:.0} m/s, \
                 +{:.0} m/s Earth-rotation bonus, {} opportunities/day (MODELLED)",
                self.site_lat_deg,
                self.target_inclination_deg,
                asc,
                v_orbit,
                site_rotation_speed(lat),
                opportunities
            ),
            None => format!(
                "launch-window: i {:.1}° < lat {:.1}° -> direct launch impossible; \
                 dogleg Δv {:.0} m/s (MODELLED)",
                self.target_inclination_deg,
                self.site_lat_deg,
                dogleg_dv.unwrap_or(0.0)
            ),
        };
        let json = serde_json::to_string_pretty(&json).map_err(|e| e.to_string())?;
        Ok((json, summary))
    }
}

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

    #[test]
    fn due_east_launch_reaches_inclination_equal_to_latitude() {
        // From any site, a due-east (Az=90°) launch reaches i = latitude.
        for lat_deg in [0.0_f64, 28.5, 51.6] {
            let lat = lat_deg.to_radians();
            let (asc, _desc) = launch_azimuth(lat, lat).unwrap();
            assert!(
                (asc.to_degrees() - 90.0).abs() < 1e-6,
                "lat {lat_deg}: due-east azimuth"
            );
        }
    }

    #[test]
    fn ksc_to_iss_inclination_is_the_textbook_45_degrees() {
        // KSC (28.5°N) to the ISS inclination (51.6°): the classic ~45° azimuth.
        let (asc, desc) = launch_azimuth(28.5_f64.to_radians(), 51.6_f64.to_radians()).unwrap();
        assert!(
            (asc.to_degrees() - 44.98).abs() < 0.1,
            "asc {:.2}",
            asc.to_degrees()
        );
        // The descending opportunity is the south-east mirror (180 − asc).
        assert!((desc.to_degrees() - (180.0 - asc.to_degrees())).abs() < 1e-9);
    }

    #[test]
    fn polar_launch_is_due_north_and_south() {
        let (asc, desc) =
            launch_azimuth(28.5_f64.to_radians(), std::f64::consts::FRAC_PI_2).unwrap();
        assert!(asc.to_degrees().abs() < 1e-6 || (asc.to_degrees() - 360.0).abs() < 1e-6);
        assert!((desc.to_degrees() - 180.0).abs() < 1e-6);
    }

    #[test]
    fn inclination_below_latitude_is_unreachable_directly() {
        // Can't reach i=10° directly from 28.5°N — needs a dogleg.
        assert!(launch_azimuth(28.5_f64.to_radians(), 10.0_f64.to_radians()).is_err());
        assert_eq!(
            daily_launch_opportunities(28.5_f64.to_radians(), 10.0_f64.to_radians()),
            0
        );
    }

    #[test]
    fn earth_rotation_bonus_is_465_m_s_at_the_equator() {
        assert!(
            (site_rotation_speed(0.0) - 465.1).abs() < 1.0,
            "equator bonus"
        );
        // Falls off as cos(lat); ~0 at the pole.
        assert!(site_rotation_speed(std::f64::consts::FRAC_PI_2) < 1e-6);
        assert!(site_rotation_speed(60.0_f64.to_radians()) < site_rotation_speed(0.0));
    }

    #[test]
    fn plane_change_dv_matches_the_vector_triangle() {
        // 10° plane change at ~7.7 km/s ≈ 1.34 km/s.
        let v = 7700.0;
        let dv = plane_change_dv(v, 10.0_f64.to_radians());
        assert!((dv - 1342.0).abs() < 5.0, "plane-change dv {dv}");
        // A full reversal (180°) costs 2v; zero change costs nothing.
        assert!((plane_change_dv(v, std::f64::consts::PI) - 2.0 * v).abs() < 1e-6);
        assert_eq!(plane_change_dv(v, 0.0), 0.0);
    }

    #[test]
    fn circular_velocity_is_about_7_67_km_s_in_leo() {
        let v = circular_velocity(400_000.0);
        assert!((7600.0..7700.0).contains(&v), "LEO circular velocity {v}");
    }

    #[test]
    fn daily_opportunities_are_two_when_strictly_between_the_bounds() {
        let lat = 28.5_f64.to_radians();
        assert_eq!(daily_launch_opportunities(lat, 51.6_f64.to_radians()), 2);
        assert_eq!(daily_launch_opportunities(lat, lat), 1); // tangent case
    }

    #[test]
    fn scenario_runs_reproducibly_and_is_modelled() {
        let scn = LaunchWindowScenario {
            site_lat_deg: 28.5,
            target_inclination_deg: 51.6,
            altitude_km: 400.0,
            site_lon_deg: 0.0,
            epoch: None,
            eop_finals2000a: None,
        };
        let (j1, _s) = scn.run_json().unwrap();
        let (j2, _s) = scn.run_json().unwrap();
        assert_eq!(j1, j2, "scenario must be reproducible");
        let v: serde_json::Value = serde_json::from_str(&j1).unwrap();
        assert_eq!(v["kind"], "launch-window");
        assert!(v["label"].as_str().unwrap().contains("MODELLED"));
        assert!(!j1.contains("VALIDATED"));
        assert!((v["launch_azimuth_deg"]["ascending"].as_f64().unwrap() - 44.98).abs() < 0.1);
        assert_eq!(v["daily_opportunities"], 2);
    }

    #[test]
    fn scenario_reports_a_dogleg_when_inclination_is_below_latitude() {
        let scn = LaunchWindowScenario {
            site_lat_deg: 51.6,
            target_inclination_deg: 28.5,
            altitude_km: 400.0,
            site_lon_deg: 0.0,
            epoch: None,
            eop_finals2000a: None,
        };
        let (j, _s) = scn.run_json().unwrap();
        let v: serde_json::Value = serde_json::from_str(&j).unwrap();
        assert!(
            v["launch_azimuth_deg"].is_null(),
            "no direct azimuth below latitude"
        );
        assert!(v["dogleg_plane_change_dv_m_s"].as_f64().unwrap() > 0.0);
    }

    // Two real IERS finals2000A rows (MJD 59579 and 59580, flag I with the Bulletin B block).
    const EOP_ROWS: &str = "211231 59579.00 I  0.056257 0.000030  0.275943 0.000035  I-0.1104179 0.0000019  0.1927 0.0016  I     0.073    0.060    -0.273    0.299  0.056304  0.275973 -0.1104355     0.040    -0.287  \n22 1 1 59580.00 I  0.054644 0.000026  0.276986 0.000032  I-0.1104988 0.0000023 -0.0267 0.0022  I     0.095    0.060    -0.250    0.299  0.054574  0.276983 -0.1105197     0.059    -0.259  \n";

    #[test]
    fn scenario_reports_the_true_pole_site_speed_at_an_epoch() {
        let scn = |epoch: Option<&str>, eop: Option<&str>| LaunchWindowScenario {
            site_lat_deg: 62.9,
            target_inclination_deg: 90.0,
            altitude_km: 400.0,
            site_lon_deg: 10.0,
            epoch: epoch.map(str::to_string),
            eop_finals2000a: eop.map(str::to_string),
        };
        // No epoch: the field is present and null.
        let (j, _) = scn(None, None).run_json().unwrap();
        let v: serde_json::Value = serde_json::from_str(&j).unwrap();
        assert!(v["site_rotation_speed_true_pole_m_s"].is_null());
        // With an epoch and EOP rows: exactly site_rotation_speed_at on the same inputs.
        let (j, _) = scn(Some("2022-01-01T00:00:00"), Some(EOP_ROWS))
            .run_json()
            .unwrap();
        let v: serde_json::Value = serde_json::from_str(&j).unwrap();
        let eop = crate::eop::EopSeries::from_finals2000a(EOP_ROWS);
        let jd_tt =
            crate::timescales::utc_to_tt(crate::timescales::julian_date(2022, 1, 1, 0, 0, 0.0));
        let want =
            site_rotation_speed_at(62.9_f64.to_radians(), 10.0_f64.to_radians(), jd_tt, &eop);
        assert_eq!(
            v["site_rotation_speed_true_pole_m_s"].as_f64().unwrap(),
            want
        );
        assert!((want / site_rotation_speed(62.9_f64.to_radians()) - 1.0).abs() < 1e-5);
        // Half the inputs, or an epoch outside the rows: an error, never a clamp.
        assert!(scn(Some("2022-01-01T00:00:00"), None).run_json().is_err());
        assert!(scn(None, Some(EOP_ROWS)).run_json().is_err());
        assert!(scn(Some("2022-01-05T00:00:00"), Some(EOP_ROWS))
            .run_json()
            .is_err());
    }
}