sgp 1.0.0

// Copyright (c) [2025] [Qin Chen <qchen2015@hotmail.com>]
// [SGP] is licensed under Mulan PSL v2.
// You can use this software according to the terms and conditions of the Mulan PSL v2.
// You may obtain a copy of Mulan PSL v2 at:
//          http://license.coscl.org.cn/MulanPSL2
// THIS SOFTWARE IS PROVIDED ON AN "AS IS" BASIS, WITHOUT WARRANTIES OF ANY KIND,
// EITHER EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO NON-INFRINGEMENT,
// MERCHANTABILITY OR FIT FOR A PARTICULAR PURPOSE.
// See the Mulan PSL v2 for more details.

use std::f64::consts::PI;

// Common helpers used across the SGP ports.
//
// This module exposes small utility types and functions that are shared by the
// propagator and the transform code. The goal is to keep thin, well-documented
// helpers here (TLE parsing, time conversions, simple math helpers) so other
// modules can import them from a single place.

/// A simple 3-vector of f64 used for position/velocity results.
///
/// This lightweight container is intentionally public and Copy so it can be
/// returned by value from the propagation routines such as `sgp8` and `sgp4`.
#[derive(Debug, Clone, Copy)]
pub struct Vector3 {
    pub x: f64,
    pub y: f64,
    pub z: f64,
}

/// Satellite input parameters required by the SGP ports.
///
/// Fields match the MATLAB `satdata` structure used by the reference SGP
/// implementations. Angles are stored in radians where appropriate and mean
/// motion is stored in radians per minute.
#[derive(Debug, Clone, Copy)]
pub struct SatData {
    pub xno: f64,
    pub xincl: f64,
    pub eo: f64,
    pub bstar: f64,
    pub omegao: f64,
    pub xnodeo: f64,
    pub xmo: f64,
}

pub fn satdata_from_tle(
    mean_motion_rev_per_day: f64,
    incl_deg: f64,
    ecc: f64,
    bstar: f64,
    argp_deg: f64,
    raan_deg: f64,
    m_deg: f64,
) -> SatData {
    let twopi = 2.0 * PI;
    let minutes_per_day = 1440.0;
    SatData {
        xno: mean_motion_rev_per_day * twopi / minutes_per_day,
        xincl: incl_deg.to_radians(),
        eo: ecc,
        bstar,
        omegao: argp_deg.to_radians(),
        xnodeo: raan_deg.to_radians(),
        xmo: m_deg.to_radians(),
    }
}

/// Create `SatData` from typical TLE numeric fields.
///
/// - `mean_motion_rev_per_day` : mean motion in revolutions per day (as in a TLE)
/// - `incl_deg`, `argp_deg`, `raan_deg`, `m_deg` : angles in degrees
/// - `ecc`, `bstar` : unitless eccentricity and B* drag coefficient
///
/// The function converts degrees->radians and rev/day->rad/min to match the
/// internal units used by the propagator.

/// Parse two TLE lines (standard 69-char lines) into `SatData`.
pub fn parse_tle_lines(line1: &str, line2: &str) -> Option<SatData> {
    if line2.len() < 63 || line1.len() < 61 {
        return None;
    }
    let incl = line2.get(8..16)?.trim().parse::<f64>().ok()?;
    let raan = line2.get(17..25)?.trim().parse::<f64>().ok()?;
    let ecc_str = line2.get(26..33)?.trim();
    let ecc = if ecc_str.starts_with('.') {
        ecc_str
            .trim_start_matches('.')
            .parse::<f64>()
            .ok()? / 10f64.powi(ecc_str.len() as i32)
    } else {
        ecc_str.parse::<f64>().ok().map(|v| v / 1e7)?
    };
    let omega = line2.get(34..42)?.trim().parse::<f64>().ok()?;
    let m = line2.get(43..51)?.trim().parse::<f64>().ok()?;
    let no = line2.get(52..63)?.trim().parse::<f64>().ok()?;
    let bstar_base = line1.get(53..59)?.trim();
    let bstar_exp = line1.get(59..61)?.trim();
    let bstar_val = bstar_base.parse::<f64>().ok().unwrap_or(0.0)
        * 1e-5
        * 10f64.powf(bstar_exp.parse::<f64>().unwrap_or(0.0));

    Some(satdata_from_tle(no, incl, ecc, bstar_val, omega, raan, m))
}

/// Parse a TLE pair (line1, line2) and return `SatData`.
///
/// The parser uses fixed-column slices consistent with standard TLE format:
/// - line1 columns 54..59 and 60..61 contain the B* mantissa/exponent used
///   to reconstruct the B* coefficient.
/// - line2 contains inclination, RAAN, eccentricity (packed), arg of perigee,
///   mean anomaly and mean motion. The function returns `None` on malformed
/// input.

pub fn fmod2p(mut x: f64) -> f64 {
    let two_pi = 2.0 * PI;
    x %= two_pi;
    if x < 0.0 {
        x += two_pi;
    }
    x
}

pub fn actan(sinf: f64, cosf: f64) -> f64 {
    sinf.atan2(cosf)
}

/// Wrap angle into the range [0, 2*pi).
///
/// This is a thin helper used to mimic MATLAB's `fmod2p` behaviour in the
/// reference codebase.


/// Extract epoch from a TLE line 1 (columns 19:32, MATLAB-style) and return
/// the Modified Julian Day (MJD) for the epoch at UTC midnight+fraction.
///
/// Returns `None` if the TLE line is malformed or parsing fails.
pub fn tle_epoch_mjd(line1: &str) -> Option<f64> {
    // TLE epoch format: yyddd.dddddd (columns 19..32, 1-based). Use safe `get`.
    let epoch_str = line1.get(18..32)?.trim();
    if epoch_str.len() < 3 {
        return None;
    }
    let yy = epoch_str.get(0..2)?.trim().parse::<i32>().ok()?;
    let doy = epoch_str.get(2..)?.trim().parse::<f64>().ok()?;
    let year = if yy < 57 { 2000 + yy } else { 1900 + yy };
    let (mon, day, hr, minute, sec) = days2mdh(year, doy);
    Some(mjday(year, mon, day, hr, minute, sec))
}

/// Parse TLE epoch (yyddd.dddddd) from line 1 and return MJD (UTC)
///
/// The TLE epoch format stores year as two digits and day-of-year with
/// fractional days. This helper converts that representation into a
/// Modified Julian Day (floating point, including fractional day).
///
/// Example:
///
/// ```rust
/// let l1 = "1 00000U 00000A   23123.12345678  .00000000  00000-0  00000-0 0  9991";
/// let mjd = sgp4rs::tle_epoch_mjd(l1).unwrap();
/// ```
///
/// Returns `None` if parsing fails.

// Mjday (moved here from transforms.rs)
pub fn mjday(year: i32, month: i32, day: i32, hour: i32, minute: i32, sec: f64) -> f64 {
    let mut year = year;
    let mut month = month;
    if month <= 2 {
        month += 12;
        year -= 1;
    }
    let b = if (10000 * year + 100 * month + day) <= 15821004 {
        -2 + ((year + 4716) as f64 / 4.0).floor() as i32 - 1179
    } else {
        (year / 400) - (year / 100) + (year / 4)
    };
    let mjd_midnight = 365 * year - 679004 + b + ((30.6001 * (month as f64 + 1.0)).floor() as i32);
    let frac_of_day = (hour as f64 + minute as f64 / 60.0 + sec / 3600.0) / 24.0;
    mjd_midnight as f64 + frac_of_day
}
/// Convert fractional day-of-year into month, day, hour, minute and second.
///
/// - `year` : full year (e.g., 2023)
/// - `days` : day-of-year, with fractional part expressing the time of day
///
/// Returns `(month, day, hour, minute, second)` which can be fed to `mjday`.
pub fn days2mdh(year: i32, days: f64) -> (i32, i32, i32, i32, f64) {
    let mut lmonth = [31i32; 12];
    lmonth.fill(31);
    lmonth[1] = 28;
    lmonth[3] = 30;
    lmonth[5] = 30;
    lmonth[8] = 30;
    lmonth[10] = 30;
    let dayofyr = days.floor() as i32;
    if (year - 1900) % 4 == 0 {
        lmonth[1] = 29;
    }
    let mut i = 0usize;
    let mut inttemp = 0i32;
    while dayofyr > inttemp + lmonth[i] && i < 11 {
        inttemp += lmonth[i];
        i += 1;
    }
    let mon = (i + 1) as i32;
    let day = dayofyr - inttemp;
    let temp = (days - dayofyr as f64) * 24.0;
    let hr = temp.floor() as i32;
    let temp2 = (temp - hr as f64) * 60.0;
    let minute = temp2.floor() as i32;
    let sec = (temp2 - minute as f64) * 60.0;
    (mon, day, hr, minute, sec)
}

/// Compute several time-scale differences used by the transforms.
///
/// Returns a 5-tuple:
/// `(UT1-TAI, UTC->GPS, UT1->GPS, TT-UTC, GPS-UTC)` in seconds where
/// appropriate. The helper uses the constant offsets `TT-TAI = 32.184 s` and
/// `GPS-TAI = -19 s` (GPS is 19 s behind TAI by definition).
pub fn timediff(ut1_utc: f64, tai_utc: f64) -> (f64, f64, f64, f64, f64) {
    let tt_tai = 32.184;
    let gps_tai = -19.0;
    let _tt_gps = tt_tai - gps_tai;
    let _tai_gps = -gps_tai;
    let ut1_tai = ut1_utc - tai_utc;
    let utc_tai = -tai_utc;
    let utc_gps = utc_tai - gps_tai;
    let ut1_gps = ut1_tai - gps_tai;
    let tt_utc = tt_tai - utc_tai;
    let gps_utc = gps_tai - utc_tai;
    (ut1_tai, utc_gps, ut1_gps, tt_utc, gps_utc)
}