ako 0.0.3

Ako is a Rust crate that offers a practical and human-friendly approach to creating, manipulating, formatting and converting dates, times and timestamps.
Documentation
use core::f64::consts::PI;
use std::f64::consts::TAU;

use crate::astronomy::{eval_polynomial, solar};

/// Returns the Julian ephemeris day (JDE) of the March equinox for the ISO `year`.
pub fn march(year: i32) -> f64 {
    longitude(year, 0.0)
}

/// Returns the Julian ephemeris day (JDE) of the June solstice for the ISO `year`.
pub fn june(year: i32) -> f64 {
    longitude(year, PI / 2.0)
}

/// Returns the Julian ephemeris day (JDE) of the September equinox for the ISO `year`.
pub fn september(year: i32) -> f64 {
    longitude(year, PI)
}

/// Returns the Julian ephemeris day (JDE) of the December solstice for the ISO `year`.
pub fn december(year: i32) -> f64 {
    longitude(year, PI * 3.0 / 2.0)
}

/// Returns the Julian ephemeris day (JDE) for a given ISO `year`
/// at a given geocentric solar longitude `lon`.
fn longitude(mut year: i32, mut long: f64) -> f64 {
    let terms = if year < 1000 {
        &TABLE_1000
    } else {
        year -= 2000;
        &TABLE_2000
    };

    long %= TAU;

    let term_index = if long < PI / 2.0 {
        0
    } else if long < PI {
        1
    } else if long < PI * 3.0 / 2.0 {
        2
    } else {
        3
    };

    let term = &terms[term_index];
    let mut j0 = eval_polynomial(f64::from(year) * 0.001, term);

    loop {
        let a = solar::position_apparent(j0);
        let c = 58.0 * (long - a.longitude).sin();

        j0 += c;

        if c.abs() < 0.000005 {
            break;
        }
    }

    j0
}

// Tables from Meeus, page 178
// Chapter 27: Equinoxes and Solstices

// Table 27.A, for the years -1000 to 1000
const TABLE_1000: [[f64; 5]; 4] = [
    [1721139.29189, 365242.13740, 0.06134, 0.00111, -0.00071],
    [1721233.25401, 365241.72562, -0.05323, 0.00907, 0.00025],
    [1721325.70455, 365242.49558, -0.11677, -0.00297, 0.00074],
    [1721414.39987, 365242.88257, -0.00769, -0.00933, -0.00006],
];

// Table 27.B, for the years 1000 to 3000
const TABLE_2000: [[f64; 5]; 4] = [
    [2451623.80984, 365242.37404, 0.05169, -0.00411, -0.00057],
    [2451716.56767, 365241.62603, 0.00325, 0.00888, -0.00030],
    [2451810.21715, 365242.01767, -0.11575, 0.00337, 0.00078],
    [2451900.05952, 365242.74049, -0.06223, -0.00823, 0.00032],
];

#[cfg(test)]
mod tests {
    use float_cmp::assert_approx_eq;
    use test_case::test_case;

    use super::{december, june, march, september};

    #[test_case(450, 1885498.2660674439)]
    #[test_case(1875, 2405968.5147716254)]
    #[test_case(2020, 2458928.6602445147)]
    #[test_case(2024, 2460389.6302345656)]
    #[test_case(2150, 2506410.170935882)]
    fn expect_march(year: i32, date: f64) {
        assert_approx_eq!(f64, march(year), date);
    }

    #[test_case(450, 1885592.0227529237)]
    #[test_case(1875, 2406061.365985682)]
    #[test_case(2020, 2459021.4061180665)]
    #[test_case(2024, 2460482.3695361647)]
    #[test_case(2150, 2506502.8153652167)]
    fn expect_june(year: i32, date: f64) {
        assert_approx_eq!(f64, june(year), date);
    }

    #[test_case(450, 1885684.8014704678)]
    #[test_case(1875, 2406154.968756302)]
    #[test_case(2020, 2459115.0637382525)]
    #[test_case(2024, 2460576.031104249)]
    #[test_case(2150, 2506596.523182565)]
    fn expect_september(year: i32, date: f64) {
        assert_approx_eq!(f64, september(year), date);
    }

    #[test_case(450, 1885773.6904944384)]
    #[test_case(1875, 2406244.719369422)]
    #[test_case(2020, 2459204.919078868)]
    #[test_case(2024, 2460665.890071845)]
    #[test_case(2150, 2506686.4764768207)]
    fn expect_december(year: i32, date: f64) {
        assert_approx_eq!(f64, december(year), date);
    }
}