satkit 0.16.2

//!
//! JPL Solar System Ephemerides
//!
//! # Introduction
//!
//! This module provides high-precision ephemerides
//! for bodies in the solar system, calculated from the
//! Jet Propulsion Laboratory (JPL) ephemeris data files
//!
//! # Links
//!
//! Ephemerides files can be found at:
//! <https://ssd.jpl.nasa.gov/ftp/eph/planets//>
//!
//!
//! # Notes
//!
//! For little-endian systems, download from the "Linux" subdirectory
//! For big-endian systems, download from the "SunOS" subdirectory
//!

use crate::solarsystem::SolarSystem;

use crate::utils::{datadir, download_if_not_exist};

use std::sync::OnceLock;

use crate::mathtypes::*;
use crate::{Instant, TimeLike, TimeScale};

use anyhow::{bail, Result};

impl TryFrom<i32> for SolarSystem {
    type Error = ();
    fn try_from(v: i32) -> Result<Self, Self::Error> {
        match v {
            x if x == Self::Mercury as i32 => Ok(Self::Mercury),
            x if x == Self::Venus as i32 => Ok(Self::Venus),
            x if x == Self::EMB as i32 => Ok(Self::EMB),
            x if x == Self::Mars as i32 => Ok(Self::Mars),
            x if x == Self::Jupiter as i32 => Ok(Self::Jupiter),
            x if x == Self::Saturn as i32 => Ok(Self::Saturn),
            x if x == Self::Uranus as i32 => Ok(Self::Uranus),
            x if x == Self::Neptune as i32 => Ok(Self::Neptune),
            x if x == Self::Pluto as i32 => Ok(Self::Pluto),
            x if x == Self::Moon as i32 => Ok(Self::Moon),
            x if x == Self::Sun as i32 => Ok(Self::Sun),
            _ => Err(()),
        }
    }
}

/// JPL Ephemeris Structure
///
/// Included ephemerides and solar system constants loaded from the
/// JPL ephemerides file.
///
/// Also includes functions to compute heliocentric and geocentric
/// positions and velocities of solar system bodies as a function
/// of time
#[derive(Debug)]
struct JPLEphem {
    /// Version of ephemeris code
    _de_version: i32,
    /// Julian date of start of ephemerides database
    jd_start: f64,
    /// Julian date of end of ephemerides database
    jd_stop: f64,
    /// Step size in Julian date
    jd_step: f64,
    /// Length of 1 astronomical unit, km
    _au: f64,
    /// Earth/Moon Ratio
    emrat: f64,

    // Offset lookup table
    ipt: [[usize; 3]; 15],
    consts: std::collections::HashMap<String, f64>,
    cheby: DMatrix<f64>,
}

/// Pre-computed parameters for Chebyshev polynomial evaluation
struct ChebySetup {
    t_seg: f64,
    offset0: usize,
    int_num: usize,
    nsubint: usize,
    ncoeff: usize,
}

macro_rules! dispatch_ncoeff {
    ($self:expr, $method:ident, $setup:expr) => {
        match $setup.ncoeff {
            6 => $self.$method::<6>($setup),
            7 => $self.$method::<7>($setup),
            8 => $self.$method::<8>($setup),
            10 => $self.$method::<10>($setup),
            11 => $self.$method::<11>($setup),
            12 => $self.$method::<12>($setup),
            13 => $self.$method::<13>($setup),
            14 => $self.$method::<14>($setup),
            _ => bail!("Invalid body"),
        }
    };
}

fn jplephem_singleton() -> &'static Result<JPLEphem> {
    static INSTANCE: OnceLock<Result<JPLEphem>> = OnceLock::new();
    INSTANCE.get_or_init(|| JPLEphem::from_file("linux_p1550p2650.440"))
}

impl JPLEphem {
    fn consts(&self, s: &str) -> Option<&f64> {
        self.consts.get(s)
    }

    /// Compute Chebyshev setup parameters for a given body and time
    fn cheby_setup(&self, body: SolarSystem, tm: &Instant) -> Result<ChebySetup> {
        let tt = tm.as_jd_with_scale(TimeScale::TT);
        if self.jd_start > tt || self.jd_stop < tt {
            bail!("Invalid Julian date: {}", tt);
        }

        let t_int = (tt - self.jd_start) / self.jd_step;
        let int_num = t_int.floor() as usize;
        let bidx = body as usize;

        let ncoeff = self.ipt[bidx][1];
        let nsubint = self.ipt[bidx][2];

        let t_int_2 = (t_int - int_num as f64) * nsubint as f64;
        let sub_int_num = t_int_2.floor() as usize;
        let t_seg = 2.0f64.mul_add(t_int_2 - sub_int_num as f64, -1.0);

        let offset0 = self.ipt[bidx][0] - 1 + sub_int_num * ncoeff * 3;

        Ok(ChebySetup {
            t_seg,
            offset0,
            int_num,
            nsubint,
            ncoeff,
        })
    }

    /// Construct a JPL Ephemerides object from the provided binary data file
    ///
    ///
    /// # Return
    ///
    /// * Object holding all of the JPL ephemerides in memory that can be
    ///   queried to find solar system body position in heliocentric or
    ///   geocentric coordinate system as function of time
    ///
    /// # Example
    ///
    /// ```
    ///
    /// use satkit::jplephem;
    /// use satkit::SolarSystem;
    /// use satkit::Instant;
    ///
    /// // Construct time: March 1 2021 12:00pm UTC
    /// let t = Instant::from_datetime(2021, 3, 1, 12, 0, 0.0).unwrap();
    ///
    /// // Find geocentric moon position at this time in the GCRF frame
    /// let p = jplephem::geocentric_pos(SolarSystem::Moon, &t).unwrap();
    /// println!("p = {}", p);
    /// ```
    ///
    fn from_file(fname: &str) -> Result<Self> {
        use std::collections::HashMap;
        use std::path::PathBuf;

        // Dimensions of ephemeris for given index
        const fn dimension(idx: usize) -> usize {
            if idx == 11 {
                2
            } else if idx == 14 {
                1
            } else {
                3
            }
        }

        // Open the file
        let path = datadir().unwrap_or_else(|_| PathBuf::from(".")).join(fname);
        if !path.is_file() {
            println!("Downloading JPL Ephemeris file.  File size is approx. 100MB");
        }
        download_if_not_exist(&path, None)?;

        // Read in bytes
        let raw = std::fs::read(path)?;
        let title: &str = std::str::from_utf8(&raw[0..84])?;

        // Get version
        let de_version: i32 = title[26..29].parse()?;

        let jd_start = f64::from_le_bytes(raw[2652..2660].try_into()?);
        let jd_stop: f64 = f64::from_le_bytes(raw[2660..2668].try_into()?);
        let jd_step: f64 = f64::from_le_bytes(raw[2668..2676].try_into()?);
        let n_con: i32 = i32::from_le_bytes(raw[2676..2680].try_into()?);
        let au: f64 = f64::from_le_bytes(raw[2680..2688].try_into()?);
        let emrat: f64 = f64::from_le_bytes(raw[2688..2696].try_into()?);

        // Get table
        let ipt: [[usize; 3]; 15] = {
            let mut ipt: [[usize; 3]; 15] = [[0, 0, 0]; 15];
            let mut idx = 2696;
            #[allow(clippy::needless_range_loop)]
            for ix in 0..15 {
                for iy in 0..3 {
                    ipt[ix][iy] = u32::from_le_bytes(raw[idx..(idx + 4)].try_into()?) as usize;
                    idx += 4;
                }
            }

            ipt[12][0] = ipt[12][1];
            ipt[12][1] = ipt[12][2];
            ipt[12][2] = ipt[13][0];

            if de_version > 430 && n_con != 400 {
                if n_con > 400 {
                    let idx = ((n_con - 400) * 6) as usize;
                    ipt[13][0] = u32::from_le_bytes(raw[idx..(idx + 4)].try_into()?) as usize;
                } else {
                    ipt[13][0] = 1_usize;
                }
            }

            // Check for garbage data not populated in earlier files
            if ipt[13][0] != (ipt[12][0] + ipt[12][1] * ipt[12][2] * 3)
                || ipt[14][0] != (ipt[13][0] + ipt[13][1] * ipt[13][2] * 3)
            {
                ipt.iter_mut().skip(13).for_each(|x| {
                    x[0] = 0;
                    x[1] = 0;
                    x[2] = 0;
                });
            }
            ipt
        };

        // Kernel size
        let kernel_size: usize = {
            let mut ks: usize = 4;
            ipt.iter().enumerate().for_each(|(ix, _)| {
                ks += 2 * ipt[ix][1] * ipt[ix][2] * dimension(ix);
            });

            ks
        };

        Ok(Self {
            _de_version: de_version,
            jd_start,
            jd_stop,
            jd_step,
            _au: au,
            emrat,
            ipt,
            consts: {
                let mut hm = HashMap::new();

                // Read in constants defined in file
                for ix in 0..n_con {
                    let sidx: usize = kernel_size * 4 + ix as usize * 8;
                    let eidx: usize = sidx + 8;
                    let val: f64 = f64::from_le_bytes(raw[sidx..eidx].try_into()?);

                    let stridx: usize = if ix >= 400 {
                        (84 * 3 + 400 * 6 + 5 * 8 + 41 * 4 + ix * 6) as usize
                    } else {
                        (84 * 3 + ix * 6) as usize
                    };
                    let s = String::from_utf8(raw[stridx..(stridx + 6)].to_vec())?;

                    hm.insert(String::from(s.trim()), val);
                }
                hm
            },
            cheby: {
                // we are going to do this unsafe since I can't find a
                // fast way to do it otherwise
                let ncoeff: usize = (kernel_size / 2) as usize;
                let nrecords = ((jd_stop - jd_start) / jd_step) as usize;
                let record_size = (kernel_size * 4) as usize;
                let mut v: DMatrix<f64> = DMatrix::zeros(ncoeff, nrecords);

                if raw.len() < record_size * 2 + ncoeff * nrecords * 8 {
                    bail!("Invalid record size for cheby data");
                }

                unsafe {
                    std::ptr::copy_nonoverlapping(
                        raw.as_ptr().add(record_size * 2) as *const f64,
                        v.as_mut_slice().as_mut_ptr(),
                        ncoeff * nrecords,
                    );
                }
                v
            },
        })
    }

    fn body_pos_optimized<const N: usize>(&self, setup: &ChebySetup) -> Result<Vector3> {
        let mut t = Vector::<N>::zeros();
        t[0] = 1.0;
        t[1] = setup.t_seg;
        for j in 2..N {
            t[j] = (2.0 * setup.t_seg).mul_add(t[j - 1], -t[j - 2]);
        }

        let mut pos = Vector3::zeros();
        for ix in 0..3 {
            let mut sum = 0.0;
            for k in 0..N {
                sum += self.cheby[(setup.offset0 + N * ix + k, setup.int_num)] * t[k];
            }
            pos[ix] = sum;
        }

        Ok(pos * 1.0e3)
    }

    /// Return the position of the given body in the Barycentric
    /// coordinate system (origin is solarsystem barycenter)
    ///
    /// # Inputs
    ///
    ///  * `body` - the solar system body for which to return position
    ///  * `tm` - The time at which to return position
    ///
    /// # Return
    ///
    ///  * 3-vector of cartesian Heliocentric position in meters
    ///
    ///
    /// # Notes:
    ///  * Positions for all bodies are natively relative to solar system barycenter,
    ///    with exception of moon, which is computed in Geocentric system
    ///  * EMB (2) is the Earth-Moon barycenter
    ///  * The sun position is relative to the solar system barycenter
    ///    (it will be close to origin)
    fn barycentric_pos(&self, body: SolarSystem, tm: &Instant) -> Result<Vector3> {
        let setup = self.cheby_setup(body, tm)?;
        dispatch_ncoeff!(self, body_pos_optimized, &setup)
    }
    /// Return the position & velocity the given body in the barycentric coordinate system
    /// (origin is solar system barycenter)
    ///
    /// # Arguments
    ///  * `body` - the solar system body for which to return position
    ///  * `tm` - The time at which to return position
    ///
    /// # Return
    ///  * 2-element tuple with following values:
    ///    * 3-vector of cartesian Heliocentric position in meters
    ///    * 3-vector of cartesian Heliocentric velocity in meters / second
    ///
    ///
    /// # Notes:
    ///  * Positions for all bodies are natively relative to solar system barycenter,
    ///    with exception of moon, which is computed in Geocentric system
    ///  * EMB (2) is the Earth-Moon barycenter
    ///  * The sun position is relative to the solar system barycenter
    ///    (it will be close to origin)
    fn barycentric_state(&self, body: SolarSystem, tm: &Instant) -> Result<(Vector3, Vector3)> {
        let setup = self.cheby_setup(body, tm)?;
        dispatch_ncoeff!(self, body_state_optimized, &setup)
    }

    fn body_state_optimized<const N: usize>(
        &self,
        setup: &ChebySetup,
    ) -> Result<(Vector3, Vector3)> {
        let mut t = Vector::<N>::zeros();
        let mut v = Vector::<N>::zeros();
        t[0] = 1.0;
        t[1] = setup.t_seg;
        v[0] = 0.0;
        v[1] = 1.0;
        for j in 2..N {
            t[j] = (2.0 * setup.t_seg).mul_add(t[j - 1], -t[j - 2]);
            v[j] = 2.0f64.mul_add(t[j - 1], (2.0 * setup.t_seg).mul_add(v[j - 1], -v[j - 2]));
        }

        let mut pos = Vector3::zeros();
        let mut vel = Vector3::zeros();
        for ix in 0..3 {
            let mut psum = 0.0;
            let mut vsum = 0.0;
            for k in 0..N {
                let coeff = self.cheby[(setup.offset0 + N * ix + k, setup.int_num)];
                psum += coeff * t[k];
                vsum += coeff * v[k];
            }
            pos[ix] = psum;
            vel[ix] = vsum;
        }

        Ok((
            pos * 1.0e3,
            vel * 2.0e3 * setup.nsubint as f64 / self.jd_step / 86400.0,
        ))
    }

    /// Return the position of the given body in
    /// Geocentric coordinate system
    ///
    /// # Arguments
    ///  * body - the solar system body for which to return position
    ///  * tm - The time at which to return position
    ///
    /// # Return
    ///    3-vector of cartesian Geocentric position in meters
    ///
    fn geocentric_pos(&self, body: SolarSystem, tm: &Instant) -> Result<Vector3> {
        if body == SolarSystem::Moon {
            self.barycentric_pos(body, tm)
        } else {
            let emb: Vector3 = self.barycentric_pos(SolarSystem::EMB, tm)?;
            let moon: Vector3 = self.barycentric_pos(SolarSystem::Moon, tm)?;
            let b: Vector3 = self.barycentric_pos(body, tm)?;

            // Compute the position of the body relative to the Earth-moon
            // barycenter, then "correct" to Earth-center by accounting
            // for moon position and Earth/moon mass ratio
            Ok(b - emb + moon / (1.0 + self.emrat))
        }
    }

    /// Return the position and velocity of the given body in
    ///  Geocentric coordinate system
    ///
    /// # Arguments
    ///
    ///  * `body` - the solar system body for which to return position
    ///  * `tm` - The time at which to return position
    ///
    /// # Return
    ///   * 2-element tuple with following elements:
    ///     * 3-vector of cartesian Geocentric position in meters
    ///     * 3-vector of cartesian Geocentric velocity in meters / second
    ///       Note: velocity is relative to Earth
    ///
    fn geocentric_state(&self, body: SolarSystem, tm: &Instant) -> Result<(Vector3, Vector3)> {
        if body == SolarSystem::Moon {
            self.barycentric_state(body, tm)
        } else {
            let emb: (Vector3, Vector3) = self.barycentric_state(SolarSystem::EMB, tm)?;
            let moon: (Vector3, Vector3) = self.barycentric_state(SolarSystem::Moon, tm)?;
            let b: (Vector3, Vector3) = self.barycentric_state(body, tm)?;

            // Compute the position of the body relative to the Earth-moon
            // barycenter, then "correct" to Earth-center by accounting
            // for moon position and Earth/moon mass ratio
            Ok((
                b.0 - emb.0 + moon.0 / (1.0 + self.emrat),
                b.1 - emb.1 + moon.1 / (1.0 + self.emrat),
            ))
        }
    }
}

pub fn consts(s: &str) -> Option<&f64> {
    jplephem_singleton().as_ref().unwrap().consts(s)
}

/// Return the position of the given body in the Barycentric
/// coordinate system (origin is solarsystem barycenter)
///
/// # Arguments
///  * `body` - the solar system body for which to return position
///  * `tm` - The time at which to return position
///
/// # Returns
///    3-vector of cartesian Heliocentric position in meters
///
///
/// # Notes:
///  * Positions for all bodies are natively relative to solar system barycenter,
///    with exception of moon, which is computed in Geocentric system
///  * EMB (2) is the Earth-Moon barycenter
///  * The sun position is relative to the solar system barycenter
///    (it will be close to origin)
pub fn barycentric_pos<T: TimeLike>(body: SolarSystem, tm: &T) -> Result<Vector3> {
    let tm = tm.as_instant();
    jplephem_singleton()
        .as_ref()
        .unwrap()
        .barycentric_pos(body, &tm)
}

/// Return the position and velocity of the given body in
///  Geocentric coordinate system
///
/// # Arguments
///  * `body` - the solar system body for which to return position
///  * `tm` - The time at which to return position
///
/// # Returns
///   * two-element tuple with following elements:
///     * 3-vector of cartesian Geocentric position in meters
///     * 3-vector of cartesian Geocentric velocity in meters / second
///       Note: velocity is relative to Earth
///
pub fn geocentric_state<T: TimeLike>(body: SolarSystem, tm: &T) -> Result<(Vector3, Vector3)> {
    let tm = tm.as_instant();
    jplephem_singleton()
        .as_ref()
        .unwrap()
        .geocentric_state(body, &tm)
}

/// Return the position of the given body in
/// Geocentric coordinate system
///
/// # Arguments
///  * `body` - the solar system body for which to return position
///  * `tm` - The time at which to return position
///
/// # Returns
///    3-vector of Cartesian Geocentric position in meters
///
pub fn geocentric_pos<T: TimeLike>(body: SolarSystem, tm: &T) -> Result<Vector3> {
    let tm = tm.as_instant();
    jplephem_singleton()
        .as_ref()
        .unwrap()
        .geocentric_pos(body, &tm)
}

/// Return the position & velocity the given body in the barycentric coordinate system
/// (origin is solar system barycenter)
///
/// # Arguments
///  * `body` - the solar system body for which to return position
///  * `tm` - The time at which to return position
///
/// # Returns
///  * two-element tuple with following values:
///    * 3-vector of cartesian Heliocentric position in meters
///    * 3-vector of cartesian Heliocentric velocity in meters / second
///
///
/// # Notes:
///  * Positions for all bodies are natively relative to solar system barycenter,
///    with exception of moon, which is computed in Geocentric system
///  * EMB (2) is the Earth-Moon barycenter
///  * The sun position is relative to the solar system barycenter
///    (it will be close to origin)
pub fn barycentric_state<T: TimeLike>(body: SolarSystem, tm: &T) -> Result<(Vector3, Vector3)> {
    let tm = tm.as_instant();
    jplephem_singleton()
        .as_ref()
        .unwrap()
        .barycentric_state(body, &tm)
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::utils::test;
    use std::io::{self, BufRead};

    #[test]
    fn load_test() {
        //let tm = &Instant::from_date(2010, 3, 1);
        let jpl = jplephem_singleton().as_ref().unwrap();

        let tm = Instant::from_jd_with_scale(2451545.0, TimeScale::TT);
        //let tm = &Instant::from_jd(2451545.0, Scale::UTC);
        let (_, _): (Vector3, Vector3) = jpl.geocentric_state(SolarSystem::Moon, &tm).unwrap();
        println!("au = {:.20}", jpl._au);
    }

    /// Load the test vectors that come with the JPL ephemeris files
    /// and compare calculated positions to test vectors.
    #[test]
    fn testvecs() {
        // Read in test vectors from the NASA JPL web site

        let jpl = jplephem_singleton().as_ref().unwrap();

        let testvecfile = test::get_testvec_dir()
            .unwrap()
            .join("jplephem")
            .join("testpo.440");

        if !testvecfile.is_file() {
            println!(
                "Required JPL ephemeris test vectors file: \"{}\" does not exist
                clone test vectors repo at
                https://github.com/StevenSamirMichael/satkit-testvecs.git
                from root of repo or set \"SATKIT_TESTVEC_ROOT\"
                to point to directory",
                testvecfile.to_string_lossy()
            );
            return;
        }

        let file = std::fs::File::open(testvecfile).unwrap();
        let b = io::BufReader::new(file);

        for rline in b.lines().skip(14) {
            let line = match rline {
                Ok(v) => v,
                Err(_) => continue,
            };
            let s: Vec<&str> = line.split_whitespace().collect();
            assert!(s.len() >= 7);
            let jd: f64 = s[2].parse().unwrap();
            let tar: i32 = s[3].parse().unwrap();
            let src: i32 = s[4].parse().unwrap();
            let coord: usize = s[5].parse().unwrap();
            let truth: f64 = s[6].parse().unwrap();
            let tm = Instant::from_jd_with_scale(jd, TimeScale::TT);
            if tar <= 10 && src <= 10 && coord <= 6 {
                let (mut tpos, mut tvel) = jpl
                    .geocentric_state(SolarSystem::try_from(tar - 1).unwrap(), &tm)
                    .unwrap();
                let (mut spos, mut svel) = jpl
                    .geocentric_state(SolarSystem::try_from(src - 1).unwrap(), &tm)
                    .unwrap();

                // in test vectors, index 3 is not EMB, but rather Earth
                // (this took me a long time to figure out...)
                if tar == 3 {
                    tpos = Vector3::zeros();
                    let (_mpos, mvel): (Vector3, Vector3) = jplephem_singleton()
                        .as_ref()
                        .unwrap()
                        .geocentric_state(SolarSystem::Moon, &tm)
                        .unwrap();
                    // Earth velocity from EMB velocity minus scaled
                    // moon velocity
                    tvel -= mvel / (1.0 + jplephem_singleton().as_ref().unwrap().emrat);
                }
                if src == 3 {
                    spos = Vector3::zeros();
                    let (_mpos, mvel): (Vector3, Vector3) = jplephem_singleton()
                        .as_ref()
                        .unwrap()
                        .geocentric_state(SolarSystem::Moon, &tm)
                        .unwrap();
                    // Earth velocity from EMB velocity minus scaled
                    // moon velocity
                    svel -= mvel / (1.0 + jplephem_singleton().as_ref().unwrap().emrat);
                }


                // Comparing positions
                if coord <= 3 {
                    let calc = (tpos - spos)[coord - 1] / jpl._au / 1.0e3;
                    // These should be very exact
                    // Allow for errors of only ~ 1e-12
                    let maxerr = 1.0e-12;
                    let err = ((truth - calc) / truth).abs();
                    assert!(err < maxerr);
                }
                // Comparing velocities
                else {
                    let calc = (tvel - svel)[coord - 4] / jpl._au / 1.0e3 * 86400.0;
                    let maxerr: f64 = 1.0e-12;
                    let err: f64 = ((truth - calc) / truth).abs();
                    assert!(err < maxerr);
                }
            }
        }
    }
}