empyrean 0.7.0

//! Ephemeris generation (predicted RA/Dec for observers).

use crate::context::Context;
use crate::coordinate::Frame;
use crate::error::{Error, Result};
use crate::observers::{Observer, obs_code_from_bytes};
use crate::orbit::Orbit;
use crate::propagate::{ForceModelTier, PropagationConfig, UncertaintyMethod};
use std::ffi::CStr;

/// One row of predicted astrometry.
#[derive(Debug, Clone, PartialEq)]
pub struct EphemerisEntry {
    /// Orbit identifier.
    pub orbit_id: String,
    /// Epoch.
    pub epoch: crate::Epoch,
    /// Right ascension (degrees).
    pub ra_deg: f64,
    /// Declination (degrees).
    pub dec_deg: f64,
    /// Topocentric distance (AU).
    pub rho_au: f64,
    /// Radial velocity (AU/day).
    pub vrho_au_day: f64,
    /// RA rate (deg/day).
    pub vra_deg_day: f64,
    /// Dec rate (deg/day).
    pub vdec_deg_day: f64,
    /// One-way light time (days). NaN if unavailable.
    pub light_time_days: f64,
    /// Phase angle (degrees). NaN if unavailable.
    pub phase_angle_deg: f64,
    /// Solar elongation (degrees). NaN if unavailable.
    pub elongation_deg: f64,
    /// Heliocentric distance (AU). NaN if unavailable.
    pub heliocentric_distance_au: f64,
    /// Predicted apparent magnitude. NaN if unavailable.
    pub mag: f64,
    /// Magnitude uncertainty. NaN if unavailable.
    pub mag_sigma: f64,
    /// Topocentric zenith angle (degrees). NaN if unavailable.
    pub zenith_angle_deg: f64,
    /// Topocentric azimuth, East of North (degrees). NaN if unavailable.
    pub azimuth_deg: f64,
    /// Local hour angle (degrees). NaN if unavailable.
    pub hour_angle_deg: f64,
    /// Angular separation from the Moon (degrees). NaN if unavailable.
    pub lunar_elongation_deg: f64,
    /// Position angle of motion, East of North (degrees). NaN if unavailable.
    pub position_angle_deg: f64,
    /// Total apparent sky-plane rate of motion (degrees/day). NaN if unavailable.
    pub sky_rate_deg_day: f64,
    /// MPC observatory code.
    pub obs_code: String,
}

impl EphemerisEntry {
    pub(crate) fn from_ffi(e: &empyrean_sys::EmpyreanEphemerisEntry) -> Self {
        let orbit_id = if e.orbit_id.is_null() {
            String::new()
        } else {
            unsafe { CStr::from_ptr(e.orbit_id).to_string_lossy().into_owned() }
        };
        Self {
            orbit_id,
            epoch: crate::Epoch::from_mjd_tdb(e.epoch_mjd_tdb),
            ra_deg: e.ra_deg,
            dec_deg: e.dec_deg,
            rho_au: e.rho_au,
            vrho_au_day: e.vrho_au_day,
            vra_deg_day: e.vra_deg_day,
            vdec_deg_day: e.vdec_deg_day,
            light_time_days: e.light_time_days,
            phase_angle_deg: e.phase_angle_deg,
            elongation_deg: e.elongation_deg,
            heliocentric_distance_au: e.heliocentric_distance_au,
            mag: e.mag,
            mag_sigma: e.mag_sigma,
            zenith_angle_deg: e.zenith_angle_deg,
            azimuth_deg: e.azimuth_deg,
            hour_angle_deg: e.hour_angle_deg,
            lunar_elongation_deg: e.lunar_elongation_deg,
            position_angle_deg: e.position_angle_deg,
            sky_rate_deg_day: e.sky_rate_deg_day,
            obs_code: obs_code_from_bytes(&e.obs_code),
        }
    }
}

/// Observation-sensitivity row: the partial derivatives of the sky-plane
/// observable w.r.t. the input state, for one `(orbit, observer, epoch)`.
/// Produced when the ephemeris uncertainty method traced the STM.
#[derive(Debug, Clone, PartialEq)]
pub struct ObservationSensitivity {
    /// Orbit identifier.
    pub orbit_id: String,
    /// Object identifier, when present.
    pub object_id: Option<String>,
    /// MPC observatory code.
    pub obs_code: String,
    /// Observation epoch (MJD TDB).
    pub epoch_mjd_tdb: f64,
    /// Solved-parameter dimension: 6 (state) or 9 (state + non-grav).
    pub n_params: u8,
    /// Jacobian ∂(observable)/∂(input), row-major `[6][n_params]` flattened.
    /// Empty when this epoch carries no Jacobian.
    pub jacobian: Vec<f64>,
    /// Hessian ∂²(observable)/∂(input)², row-major `[6][n_params][n_params]`
    /// flattened. Empty unless a second-order method ran.
    pub hessian: Vec<f64>,
    /// Frame of the input axis (Frame enum as int).
    pub frame: i32,
    /// Origin of the input axis (NAIF id).
    pub origin: i32,
}

impl ObservationSensitivity {
    pub(crate) fn from_ffi(e: &empyrean_sys::EmpyreanObservationSensitivity) -> Self {
        let orbit_id = if e.orbit_id.is_null() {
            String::new()
        } else {
            unsafe { CStr::from_ptr(e.orbit_id).to_string_lossy().into_owned() }
        };
        let object_id = if e.object_id.is_null() {
            None
        } else {
            Some(unsafe { CStr::from_ptr(e.object_id).to_string_lossy().into_owned() })
        };
        let jacobian = if e.jacobian.is_null() {
            Vec::new()
        } else {
            unsafe { std::slice::from_raw_parts(e.jacobian, e.jacobian_len).to_vec() }
        };
        let hessian = if e.hessian.is_null() {
            Vec::new()
        } else {
            unsafe { std::slice::from_raw_parts(e.hessian, e.hessian_len).to_vec() }
        };
        Self {
            orbit_id,
            object_id,
            obs_code: obs_code_from_bytes(&e.obs_code),
            epoch_mjd_tdb: e.epoch_mjd_tdb,
            n_params: e.n_params,
            jacobian,
            hessian,
            frame: e.frame,
            origin: e.origin,
        }
    }
}

/// Result of [`Context::generate_ephemeris`]: the per-`(orbit, observer,
/// epoch)` ephemeris entries plus the observation-sensitivity rows (empty
/// unless the uncertainty method traced the STM).
#[derive(Debug, Clone, PartialEq)]
pub struct EphemerisResult {
    /// Ephemeris entries (RA/Dec + diagnostics), one per observation.
    pub entries: Vec<EphemerisEntry>,
    /// Observation-sensitivity rows (Jacobian/Hessian). Empty on the
    /// f64-only path.
    pub sensitivity: Vec<ObservationSensitivity>,
}

/// Ephemeris-generation configuration.
///
/// Mirrors the upstream `EphemerisConfig`: drives the inner
/// propagation that brings each orbit to its observation epoch, plus
/// the light-time iteration loop and the diagnostics toggle that
/// gates phase-angle / elongation / heliocentric-distance / magnitude
/// computation. Use [`EphemerisConfig::default`] for sensible
/// production settings.
#[derive(Debug, Clone, PartialEq)]
pub struct EphemerisConfig {
    /// Inner propagation configuration applied while bringing each
    /// orbit to its observation epoch. Sets force model, uncertainty
    /// method, and output frame for the trajectory.
    pub propagation: PropagationConfig,
    /// Maximum iterations for light-time convergence. 0 selects the
    /// upstream default (3).
    pub max_light_time_iterations: usize,
    /// Tolerance (days) for light-time convergence. 0.0 selects the
    /// upstream default (1e-10).
    pub light_time_tolerance_days: f64,
    /// Whether to compute phase-angle, elongation, heliocentric
    /// distance, and apparent magnitude. Skip these when only RA/Dec
    /// are needed (DC inner loop) for a small speedup.
    pub compute_diagnostics: bool,
}

impl Default for EphemerisConfig {
    fn default() -> Self {
        Self {
            propagation: PropagationConfig::default(),
            max_light_time_iterations: 0,
            light_time_tolerance_days: 0.0,
            compute_diagnostics: true,
        }
    }
}

impl EphemerisConfig {
    /// Build the C-ABI representation. Returns the FFI struct plus a
    /// keepalive that owns the raw arrays the FFI struct points into.
    pub(crate) fn to_ffi_with(
        &self,
    ) -> (
        empyrean_sys::EmpyreanEphemerisConfig,
        crate::propagate::PropConfigKeep,
    ) {
        let (prop_ffi, keep) = self.propagation.to_ffi_with();
        let cfg = empyrean_sys::EmpyreanEphemerisConfig {
            propagation: prop_ffi,
            max_light_time_iterations: self.max_light_time_iterations,
            light_time_tolerance_days: self.light_time_tolerance_days,
            compute_diagnostics: u8::from(self.compute_diagnostics),
        };
        (cfg, keep)
    }

    /// Convenience builder: a config carrying just the requested force
    /// model, defaults for everything else.
    ///
    /// The internal propagation runs in **EclipticJ2000** because
    /// villeneuve's ephemeris pipeline assumes that integration frame
    /// when it converts the propagated state to ICRF for the
    /// observer-relative geometry. The user-facing RA/Dec output is
    /// still in ICRF — only the integration frame is overridden here.
    pub fn with_force_model(force_model: ForceModelTier) -> Self {
        Self {
            propagation: PropagationConfig {
                force_model,
                frame: Frame::EclipticJ2000,
                ..PropagationConfig::default()
            },
            ..Self::default()
        }
    }
}

impl Context {
    /// Generate predicted ephemeris for orbits as seen by observers.
    ///
    /// Returns `num_orbits * num_observers` entries (orbit-major order).
    pub fn generate_ephemeris(
        &self,
        orbits: &[Orbit],
        observers: &[Observer],
        config: &EphemerisConfig,
    ) -> Result<EphemerisResult> {
        let _ = (Frame::ICRF, UncertaintyMethod::FirstOrder); // suppress unused-import in default-config branch
        let mut _orbit_keep: Vec<crate::orbit::OrbitFfiKeep> = Vec::with_capacity(orbits.len());
        let ffi_orbits: Vec<_> = orbits
            .iter()
            .map(|o| {
                let (ffi, keep) = o.to_ffi_with_keep()?;
                _orbit_keep.push(keep);
                Ok(ffi)
            })
            .collect::<Result<Vec<_>>>()?;
        let ffi_observers: Vec<empyrean_sys::EmpyreanObserver> = observers
            .iter()
            .map(|o| {
                let mut code_bytes = [0u8; 4];
                let src = o.obs_code.as_bytes();
                let n = src.len().min(3);
                code_bytes[..n].copy_from_slice(&src[..n]);
                Ok(empyrean_sys::EmpyreanObserver {
                    obs_code: code_bytes,
                    epoch_mjd_tdb: o.epoch.mjd_tdb()?,
                    x: o.position[0],
                    y: o.position[1],
                    z: o.position[2],
                    vx: o.velocity[0],
                    vy: o.velocity[1],
                    vz: o.velocity[2],
                    observing_night: o.observing_night,
                })
            })
            .collect::<Result<Vec<_>>>()?;

        let mut result = empyrean_sys::EmpyreanEphemerisResult {
            entries: std::ptr::null_mut(),
            num_entries: 0,
            sensitivity: std::ptr::null_mut(),
            num_sensitivity: 0,
        };
        let (ffi_config, _config_keep) = config.to_ffi_with();
        let code = unsafe {
            empyrean_sys::empyrean_generate_ephemeris(
                self.as_raw(),
                ffi_orbits.as_ptr(),
                ffi_orbits.len(),
                ffi_observers.as_ptr(),
                ffi_observers.len(),
                &ffi_config,
                &mut result,
            )
        };
        if code != 0 {
            return Err(Error::capture(code));
        }
        let entries = unsafe {
            std::slice::from_raw_parts(result.entries, result.num_entries)
                .iter()
                .map(EphemerisEntry::from_ffi)
                .collect()
        };
        let sensitivity = if result.sensitivity.is_null() {
            Vec::new()
        } else {
            unsafe {
                std::slice::from_raw_parts(result.sensitivity, result.num_sensitivity)
                    .iter()
                    .map(ObservationSensitivity::from_ffi)
                    .collect()
            }
        };
        unsafe { empyrean_sys::empyrean_ephemeris_result_free(&mut result) };
        Ok(EphemerisResult {
            entries,
            sensitivity,
        })
    }
}