astrodynamics 0.14.0

//! TLE pass prediction over a ground station.
//!
//! This module owns the legacy Orbis pass orchestration: coarse elevation
//! sampling, horizon-crossing bisection, and peak-elevation search. The
//! low-level propagation and frame transforms stay delegated to the core SGP4
//! and frames modules.

use crate::constants::time::{MICROSECONDS_PER_DAY_I64, SECONDS_PER_DAY, SECONDS_PER_DAY_I64};
use crate::constants::units::MICROSECONDS_PER_SECOND_I64;
use crate::frames::transforms::{
    gcrs_to_topocentric_compute, teme_to_gcrs_compute, GeodeticStationKm, TemeStateKm,
};
use crate::sgp4::{ElementSet, JulianDate, OpsMode, Prediction, Satellite};
use crate::time::scales::TimeScales;

const UNIX_EPOCH_JDN: i64 = 2_440_588;

const BISECT_ITERATIONS: usize = 20;
const GOLDEN_ITERATIONS: usize = 30;
const GOLDEN_RESPHI: f64 = 0.381_966_011_250_105_1;

/// UTC instant represented as unix microseconds.
#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord)]
pub struct UtcInstant {
    unix_microseconds: i64,
}

impl UtcInstant {
    /// Construct from unix microseconds.
    pub fn from_unix_microseconds(unix_microseconds: i64) -> Self {
        Self { unix_microseconds }
    }

    /// Construct from UTC calendar fields.
    pub fn from_utc(
        year: i32,
        month: i32,
        day: i32,
        hour: i32,
        minute: i32,
        second: i32,
        microsecond: i32,
    ) -> Option<Self> {
        if !(1..=12).contains(&month)
            || !(1..=31).contains(&day)
            || !(0..=23).contains(&hour)
            || !(0..=59).contains(&minute)
            || !(0..=60).contains(&second)
            || !(0..=999_999).contains(&microsecond)
        {
            return None;
        }

        let days = julian_day_number(year, month, day) - UNIX_EPOCH_JDN;
        let seconds_of_day = hour as i64 * 3600 + minute as i64 * 60 + second as i64;
        Some(Self {
            unix_microseconds: days * MICROSECONDS_PER_DAY_I64
                + seconds_of_day * MICROSECONDS_PER_SECOND_I64
                + microsecond as i64,
        })
    }

    /// Unix microseconds.
    pub fn unix_microseconds(self) -> i64 {
        self.unix_microseconds
    }

    fn add_microseconds(self, delta: i64) -> Self {
        Self {
            unix_microseconds: self.unix_microseconds + delta,
        }
    }

    fn diff_microseconds(self, earlier: Self) -> i64 {
        self.unix_microseconds - earlier.unix_microseconds
    }

    fn diff_seconds(self, earlier: Self) -> i64 {
        self.diff_microseconds(earlier) / MICROSECONDS_PER_SECOND_I64
    }

    fn components(self) -> UtcComponents {
        let seconds = div_floor(self.unix_microseconds, MICROSECONDS_PER_SECOND_I64);
        let microsecond = rem_floor(self.unix_microseconds, MICROSECONDS_PER_SECOND_I64);
        let days = div_floor(seconds, SECONDS_PER_DAY_I64);
        let second_of_day = seconds - days * SECONDS_PER_DAY_I64;
        let (year, month, day) = civil_from_days(days);

        UtcComponents {
            year,
            month,
            day,
            hour: (second_of_day / 3600) as i32,
            minute: ((second_of_day % 3600) / 60) as i32,
            second: (second_of_day % 60) as i32,
            microsecond: microsecond as i32,
        }
    }

    fn time_scales(self) -> TimeScales {
        let c = self.components();
        TimeScales::from_utc(
            c.year,
            c.month,
            c.day,
            c.hour,
            c.minute,
            c.second as f64 + c.microsecond as f64 / 1_000_000.0,
        )
    }

    fn sgp4_julian_date(self) -> JulianDate {
        let c = self.components();
        let jdn = julian_day_number(c.year, c.month, c.day);
        let jd_midnight = jdn as f64 - 0.5;
        let frac = (c.hour as f64) / 24.0
            + (c.minute as f64) / 1440.0
            + (c.second as f64) / SECONDS_PER_DAY
            + (c.microsecond as f64) / MICROSECONDS_PER_DAY_I64 as f64;
        JulianDate(jd_midnight, frac)
    }
}

#[derive(Debug, Clone, Copy)]
struct UtcComponents {
    year: i32,
    month: i32,
    day: i32,
    hour: i32,
    minute: i32,
    second: i32,
    microsecond: i32,
}

/// Ground-station geodetic coordinates.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct GroundStation {
    pub latitude_deg: f64,
    pub longitude_deg: f64,
    pub altitude_m: f64,
}

/// Pass-prediction options.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct PassPredictionOptions {
    pub min_elevation_deg: f64,
    pub step_seconds: i64,
}

impl Default for PassPredictionOptions {
    fn default() -> Self {
        Self {
            min_elevation_deg: 0.0,
            step_seconds: 60,
        }
    }
}

/// Predicted visible pass.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct PredictedPass {
    pub rise: UtcInstant,
    pub set: UtcInstant,
    pub max_elevation_deg: f64,
    pub max_elevation_time: UtcInstant,
}

/// Topocentric look angle from a ground station to a TLE satellite.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct LookAngle {
    pub azimuth_deg: f64,
    pub elevation_deg: f64,
    pub range_km: f64,
}

/// One member of a TLE-backed constellation.
#[derive(Debug, Clone, PartialEq)]
pub struct ConstellationMember {
    pub catalog_number: String,
    pub elements: ElementSet,
}

/// One satellite visible from a ground station at an instant.
#[derive(Debug, Clone, PartialEq)]
pub struct VisibleSatellite {
    pub catalog_number: String,
    pub azimuth_deg: f64,
    pub elevation_deg: f64,
    pub range_km: f64,
    pub position_km: [f64; 3],
}

/// Error while computing a TLE look angle.
#[derive(Debug, Clone, PartialEq)]
pub enum LookAngleError {
    Init(crate::sgp4::Error),
    Propagate(crate::sgp4::Error),
}

/// Propagate a pre-parsed SGP4 element set and compute its topocentric look angle.
pub fn look_angle(
    elements: &ElementSet,
    ground_station: GroundStation,
    datetime: UtcInstant,
) -> Result<LookAngle, LookAngleError> {
    let satellite = Satellite::from_elements_with_opsmode(elements, OpsMode::Afspc)
        .map_err(LookAngleError::Init)?;
    let pred = satellite
        .propagate_jd(datetime.sgp4_julian_date())
        .map_err(LookAngleError::Propagate)?;
    Ok(look_angle_from_teme_prediction(
        &pred,
        datetime,
        ground_station,
    ))
}

/// Find constellation members above an elevation threshold at one instant.
///
/// Invalid element sets or per-satellite propagation failures are skipped,
/// matching the legacy Orbis `visible_from/4` behavior through
/// `propagate_all/2`.
pub fn visible_from_constellation(
    members: &[ConstellationMember],
    ground_station: GroundStation,
    datetime: UtcInstant,
    min_elevation_deg: f64,
) -> Vec<VisibleSatellite> {
    let mut visible = Vec::new();

    for member in members {
        let satellite =
            match Satellite::from_elements_with_opsmode(&member.elements, OpsMode::Afspc) {
                Ok(satellite) => satellite,
                Err(_) => continue,
            };
        let pred = match satellite.propagate_jd(datetime.sgp4_julian_date()) {
            Ok(pred) => pred,
            Err(_) => continue,
        };
        let look = look_angle_from_teme_prediction(&pred, datetime, ground_station);

        if look.elevation_deg >= min_elevation_deg {
            visible.push(VisibleSatellite {
                catalog_number: member.catalog_number.clone(),
                azimuth_deg: look.azimuth_deg,
                elevation_deg: look.elevation_deg,
                range_km: look.range_km,
                position_km: pred.position,
            });
        }
    }

    visible.sort_by(|a, b| {
        b.elevation_deg
            .partial_cmp(&a.elevation_deg)
            .unwrap_or(std::cmp::Ordering::Equal)
    });
    visible
}

/// Predict visible passes for a pre-parsed SGP4 element set.
///
/// Invalid element sets or per-sample propagation failures are treated as
/// below-horizon samples, matching the legacy Orbis public behavior.
pub fn predict_passes(
    elements: &ElementSet,
    ground_station: GroundStation,
    start_time: UtcInstant,
    end_time: UtcInstant,
    options: PassPredictionOptions,
) -> Vec<PredictedPass> {
    if options.step_seconds <= 0 {
        return Vec::new();
    }

    let satellite = match Satellite::from_elements_with_opsmode(elements, OpsMode::Afspc) {
        Ok(satellite) => satellite,
        Err(_) => return Vec::new(),
    };

    let samples = coarse_scan(
        &satellite,
        ground_station,
        start_time,
        end_time,
        options.step_seconds,
    );

    extract_passes(&samples, &satellite, ground_station)
        .into_iter()
        .filter(|pass| pass.max_elevation_deg >= options.min_elevation_deg)
        .collect()
}

fn coarse_scan(
    satellite: &Satellite,
    ground_station: GroundStation,
    start_time: UtcInstant,
    end_time: UtcInstant,
    step_seconds: i64,
) -> Vec<(UtcInstant, f64)> {
    let total_seconds = end_time.diff_seconds(start_time);
    let num_steps = (total_seconds / step_seconds).max(0);

    (0..=num_steps)
        .map(|i| {
            let dt = start_time.add_microseconds(i * step_seconds * MICROSECONDS_PER_SECOND_I64);
            (dt, elevation_at(satellite, dt, ground_station))
        })
        .collect()
}

fn extract_passes(
    samples: &[(UtcInstant, f64)],
    satellite: &Satellite,
    ground_station: GroundStation,
) -> Vec<PredictedPass> {
    let mut rise_time = match samples.first() {
        Some((dt, el)) if *el >= 0.0 => Some(*dt),
        _ => None,
    };
    let mut passes = Vec::new();

    for pair in samples.windows(2) {
        let (dt_a, el_a) = pair[0];
        let (dt_b, el_b) = pair[1];

        if rise_time.is_none() && el_a < 0.0 && el_b >= 0.0 {
            rise_time = Some(bisect_crossing(satellite, ground_station, dt_a, dt_b));
        } else if let Some(rise) = rise_time {
            if el_a >= 0.0 && el_b < 0.0 {
                let set = bisect_crossing(satellite, ground_station, dt_a, dt_b);
                passes.push(build_pass(satellite, ground_station, rise, set));
                rise_time = None;
            }
        }
    }

    passes
}

fn bisect_crossing(
    satellite: &Satellite,
    ground_station: GroundStation,
    dt_low: UtcInstant,
    dt_high: UtcInstant,
) -> UtcInstant {
    let mut lo = dt_low;
    let mut hi = dt_high;
    let mut el_lo = elevation_at(satellite, lo, ground_station);

    for _ in 0..BISECT_ITERATIONS {
        let mid = midpoint_instant(lo, hi);
        let el_mid = elevation_at(satellite, mid, ground_station);

        if same_sign(el_lo, el_mid) {
            lo = mid;
            el_lo = el_mid;
        } else {
            hi = mid;
        }
    }

    midpoint_instant(lo, hi)
}

fn same_sign(a: f64, b: f64) -> bool {
    (a >= 0.0 && b >= 0.0) || (a < 0.0 && b < 0.0)
}

fn midpoint_instant(a: UtcInstant, b: UtcInstant) -> UtcInstant {
    a.add_microseconds(b.diff_microseconds(a) / 2)
}

fn build_pass(
    satellite: &Satellite,
    ground_station: GroundStation,
    rise: UtcInstant,
    set: UtcInstant,
) -> PredictedPass {
    let (max_elevation_deg, max_elevation_time) =
        find_max_elevation(satellite, ground_station, rise, set);

    PredictedPass {
        rise,
        set,
        max_elevation_deg,
        max_elevation_time,
    }
}

fn find_max_elevation(
    satellite: &Satellite,
    ground_station: GroundStation,
    rise: UtcInstant,
    set: UtcInstant,
) -> (f64, UtcInstant) {
    let total_us = set.diff_microseconds(rise);
    let mut a = 0_i64;
    let mut b = total_us;

    for _ in 0..GOLDEN_ITERATIONS {
        let span = b - a;
        let x1 = (a as f64 + GOLDEN_RESPHI * span as f64).round() as i64;
        let x2 = (b as f64 - GOLDEN_RESPHI * span as f64).round() as i64;

        let dt1 = rise.add_microseconds(x1);
        let dt2 = rise.add_microseconds(x2);
        let el1 = elevation_at(satellite, dt1, ground_station);
        let el2 = elevation_at(satellite, dt2, ground_station);

        if el1 > el2 {
            b = x2;
        } else {
            a = x1;
        }
    }

    let best_us = (a + b) / 2;
    let best_dt = rise.add_microseconds(best_us);
    let best_el = elevation_at(satellite, best_dt, ground_station);
    (best_el, best_dt)
}

fn elevation_at(satellite: &Satellite, datetime: UtcInstant, ground_station: GroundStation) -> f64 {
    let pred = match satellite.propagate_jd(datetime.sgp4_julian_date()) {
        Ok(pred) => pred,
        Err(_) => return -90.0,
    };
    look_angle_from_teme_prediction(&pred, datetime, ground_station).elevation_deg
}

fn look_angle_from_teme_prediction(
    pred: &Prediction,
    datetime: UtcInstant,
    ground_station: GroundStation,
) -> LookAngle {
    let ts = datetime.time_scales();
    let (gcrs_position, _) = teme_to_gcrs_compute(
        &TemeStateKm {
            position_km: [pred.position[0], pred.position[1], pred.position[2]],
            velocity_km_s: [pred.velocity[0], pred.velocity[1], pred.velocity[2]],
        },
        &ts,
        false,
    );
    let (azimuth, elevation, range) = gcrs_to_topocentric_compute(
        [gcrs_position.0, gcrs_position.1, gcrs_position.2],
        &GeodeticStationKm {
            latitude_deg: ground_station.latitude_deg,
            longitude_deg: ground_station.longitude_deg,
            altitude_km: ground_station.altitude_m / 1000.0,
        },
        &ts,
        false,
    );
    LookAngle {
        azimuth_deg: azimuth,
        elevation_deg: elevation,
        range_km: range,
    }
}

fn julian_day_number(year: i32, month: i32, day: i32) -> i64 {
    let a = (14 - month) / 12;
    let y = year + 4800 - a;
    let m = month + 12 * a - 3;
    (day + (153 * m + 2) / 5 + 365 * y + y / 4 - y / 100 + y / 400 - 32045) as i64
}

fn civil_from_days(days_since_unix_epoch: i64) -> (i32, i32, i32) {
    let z = days_since_unix_epoch + 719_468;
    let era = if z >= 0 { z } else { z - 146_096 } / 146_097;
    let doe = z - era * 146_097;
    let yoe = (doe - doe / 1460 + doe / 36_524 - doe / 146_096) / 365;
    let y = yoe + era * 400;
    let doy = doe - (365 * yoe + yoe / 4 - yoe / 100);
    let mp = (5 * doy + 2) / 153;
    let d = doy - (153 * mp + 2) / 5 + 1;
    let m = mp + if mp < 10 { 3 } else { -9 };
    let year = y + if m <= 2 { 1 } else { 0 };
    (year as i32, m as i32, d as i32)
}

fn div_floor(a: i64, b: i64) -> i64 {
    let q = a / b;
    let r = a % b;
    if r != 0 && (r > 0) != (b > 0) {
        q - 1
    } else {
        q
    }
}

fn rem_floor(a: i64, b: i64) -> i64 {
    a - div_floor(a, b) * b
}

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

    fn iss_2024_12_19_elements() -> ElementSet {
        ElementSet {
            epoch_year_two_digit: 24,
            epoch_days: 354.52609954,
            bstar: 0.000_370_420_000_000_000_05,
            mean_motion_dot: 0.00020888,
            mean_motion_double_dot: 0.0,
            eccentricity: 0.0006955,
            argument_of_perigee_deg: 37.7614,
            inclination_deg: 51.6393,
            mean_anomaly_deg: 87.9783,
            mean_motion_rev_per_day: 15.49970085,
            right_ascension_deg: 213.2584,
            catalog_number: 0,
        }
    }

    fn iss_2024_01_01_elements() -> ElementSet {
        ElementSet {
            epoch_year_two_digit: 24,
            epoch_days: 1.5,
            bstar: 0.000_102_70,
            mean_motion_dot: 0.000_167_17,
            mean_motion_double_dot: 0.0,
            eccentricity: 0.000_264_4,
            argument_of_perigee_deg: 250.3037,
            inclination_deg: 51.6400,
            mean_anomaly_deg: 109.7782,
            mean_motion_rev_per_day: 15.49560812,
            right_ascension_deg: 208.8657,
            catalog_number: 25_544,
        }
    }

    fn iss_fixture_elements() -> ElementSet {
        ElementSet {
            epoch_year_two_digit: 26,
            epoch_days: 95.55331950,
            bstar: 0.000_164_20,
            mean_motion_dot: 0.000_085_43,
            mean_motion_double_dot: 0.0,
            eccentricity: 0.000_635_1,
            argument_of_perigee_deg: 274.8255,
            inclination_deg: 51.6328,
            mean_anomaly_deg: 85.2008,
            mean_motion_rev_per_day: 15.4878698,
            right_ascension_deg: 299.5432,
            catalog_number: 25_544,
        }
    }

    fn css_fixture_elements() -> ElementSet {
        ElementSet {
            epoch_year_two_digit: 26,
            epoch_days: 95.32454765,
            bstar: 0.000_372_23,
            mean_motion_dot: 0.000_331_73,
            mean_motion_double_dot: 0.0,
            eccentricity: 0.000_355_7,
            argument_of_perigee_deg: 129.2727,
            inclination_deg: 41.4682,
            mean_anomaly_deg: 230.8429,
            mean_motion_rev_per_day: 15.6194274,
            right_ascension_deg: 45.9319,
            catalog_number: 48_274,
        }
    }

    fn fregat_fixture_elements() -> ElementSet {
        ElementSet {
            epoch_year_two_digit: 26,
            epoch_days: 95.51225242,
            bstar: 0.012_423,
            mean_motion_dot: 0.000_085_41,
            mean_motion_double_dot: 0.0,
            eccentricity: 0.095_504_7,
            argument_of_perigee_deg: 120.8974,
            inclination_deg: 51.6426,
            mean_anomaly_deg: 248.9327,
            mean_motion_rev_per_day: 12.40936816,
            right_ascension_deg: 220.2066,
            catalog_number: 49_271,
        }
    }

    #[test]
    fn utc_instant_round_trips_calendar_fields() {
        let instant = UtcInstant::from_utc(2024, 12, 19, 7, 3, 11, 825_435).unwrap();
        assert_eq!(instant.unix_microseconds(), 1_734_591_791_825_435);
        let c = instant.components();
        assert_eq!(
            (
                c.year,
                c.month,
                c.day,
                c.hour,
                c.minute,
                c.second,
                c.microsecond
            ),
            (2024, 12, 19, 7, 3, 11, 825_435)
        );
    }

    #[test]
    fn iss_london_pass_matches_legacy_orbis_bits() {
        let start = UtcInstant::from_utc(2024, 12, 19, 0, 0, 0, 0).unwrap();
        let end = UtcInstant::from_utc(2024, 12, 19, 12, 0, 0, 0).unwrap();
        let station = GroundStation {
            latitude_deg: 51.5074,
            longitude_deg: -0.1278,
            altitude_m: 11.0,
        };

        let passes = predict_passes(
            &iss_2024_12_19_elements(),
            station,
            start,
            end,
            PassPredictionOptions::default(),
        );

        assert_eq!(passes.len(), 1);
        let pass = passes[0];
        assert_eq!(pass.rise.unix_microseconds(), 1_734_604_991_825_435);
        assert_eq!(pass.set.unix_microseconds(), 1_734_605_533_400_371);
        assert_eq!(
            pass.max_elevation_time.unix_microseconds(),
            1_734_605_261_892_583
        );
        assert_eq!(pass.max_elevation_deg.to_bits(), 0x4029_1832_84c1_525f);

        let high = predict_passes(
            &iss_2024_12_19_elements(),
            station,
            start,
            end,
            PassPredictionOptions {
                min_elevation_deg: 30.0,
                step_seconds: 60,
            },
        );
        assert!(high.is_empty());
    }

    #[test]
    fn iss_london_look_angle_matches_legacy_orbis_bits() {
        let datetime = UtcInstant::from_utc(2024, 1, 1, 12, 0, 0, 0).unwrap();
        let station = GroundStation {
            latitude_deg: 51.5,
            longitude_deg: -0.1,
            altitude_m: 11.0,
        };

        let look = look_angle(&iss_2024_01_01_elements(), station, datetime).unwrap();

        assert_eq!(look.azimuth_deg.to_bits(), 0x406f_f4aa_a5f4_2254);
        assert_eq!(look.elevation_deg.to_bits(), 0xc042_8a29_691f_1ca2);
        assert_eq!(look.range_km.to_bits(), 0x40c0_4e5e_046d_c53b);
    }

    #[test]
    fn constellation_visible_from_matches_legacy_orbis_bits() {
        let datetime = UtcInstant::from_utc(2026, 4, 5, 13, 16, 46, 804_800).unwrap();
        let station = GroundStation {
            latitude_deg: 51.5074,
            longitude_deg: -0.1278,
            altitude_m: 11.0,
        };
        let members = vec![
            ConstellationMember {
                catalog_number: "25544".to_string(),
                elements: iss_fixture_elements(),
            },
            ConstellationMember {
                catalog_number: "48274".to_string(),
                elements: css_fixture_elements(),
            },
            ConstellationMember {
                catalog_number: "49271".to_string(),
                elements: fregat_fixture_elements(),
            },
        ];

        let visible = visible_from_constellation(&members, station, datetime, -90.0);

        assert_eq!(
            visible
                .iter()
                .map(|sat| sat.catalog_number.as_str())
                .collect::<Vec<_>>(),
            ["49271", "25544", "48274"]
        );

        assert_eq!(visible[0].elevation_deg.to_bits(), 0xc03d_2a2f_bd8b_c9ba);
        assert_eq!(visible[0].azimuth_deg.to_bits(), 0x4063_ad5e_7f91_98bd);
        assert_eq!(visible[0].range_km.to_bits(), 0x40c0_7b37_546a_e871);
        assert_eq!(
            visible[0]
                .position_km
                .iter()
                .map(|value| value.to_bits())
                .collect::<Vec<_>>(),
            [
                0x40b0_1a88_998c_fb44,
                0x40b7_988c_0568_1811,
                0xc0a3_68d5_167e_3886,
            ]
        );

        assert_eq!(visible[1].elevation_deg.to_bits(), 0xc046_1d42_86e2_51f9);
        assert_eq!(visible[1].azimuth_deg.to_bits(), 0x4071_0d2c_4f97_dcbb);
        assert_eq!(visible[1].range_km.to_bits(), 0x40c2_8587_4aa0_b9c8);

        assert_eq!(visible[2].elevation_deg.to_bits(), 0xc051_0316_042a_4a06);
        assert_eq!(visible[2].azimuth_deg.to_bits(), 0x4073_5174_638d_55d0);
        assert_eq!(visible[2].range_km.to_bits(), 0x40c7_e8e0_e370_5793);

        assert!(visible_from_constellation(&members, station, datetime, -20.0).is_empty());
    }
}