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 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125

use crate::astronomy::eval_polynomial;
use crate::astronomy::julian::j2000_century;

/// Returns nutation in longitude for the given Julian ephemeris day (JDE).
pub fn nutation(date: f64) -> f64 {
    // https://ui.adsabs.harvard.edu/abs/1982CeMec..27...79S/abstract
    // https://neoprogrammics.com/nutations/Nutation_In_Obliquity.php
    // https://observablehq.com/@igor-a1/meeus-aa-v1#cell-405

    let t = j2000_century(date);

    // Mean elongation of the Moon from the Sun
    let d = eval_polynomial(t, &[297.85036, 445267.11148, -0.0019142, 1.0 / 189474.0]).to_radians();

    // Mean anomaly of the Sun (Earth)
    let m = eval_polynomial(t, &[357.52772, 35999.05034, -0.0001603, -1.0 / 300000.0]).to_radians();

    // Mean anomaly of the Moon
    let m1 = eval_polynomial(t, &[134.96298, 477198.867398, 0.0086972, 1.0 / 56250.0]).to_radians();

    // Moon's argument of latitude
    let f = eval_polynomial(t, &[93.27191, 483202.017538, -0.0036825, 1.0 / 327270.0]).to_radians();

    // Longitude of the ascending node of the Moon's mean orbit on the ecliptic, measured from the mean equinox of the date
    let w = eval_polynomial(t, &[125.04452, -1934.136261, 0.0020708, 1.0 / 450000.0]).to_radians();

    // Compute the nutation in longitude
    let mut long = 0.0;

    for ([a0, a1, a2, a3, a4], sin, _) in TABLE.iter().rev() {
        let a0 = *a0 as f64;
        let a1 = *a1 as f64;
        let a2 = *a2 as f64;
        let a3 = *a3 as f64;
        let a4 = *a4 as f64;

        let arg = a4.mul_add(w, a3.mul_add(f, a2.mul_add(m1, a0.mul_add(d, a1 * m))));

        long += arg.sin() * sin[1].mul_add(t, sin[0]);
    }

    (long * 0.0001 / 3600.0_f64).to_radians()
}

// https://observablehq.com/@igor-a1/meeus-aa-v1#cell-409
// From “Astronomical Algorithms” (2nd ed.), chapter 22
#[rustfmt::skip]
const TABLE: [([i8; 5], [f64; 2], [f64; 2]); 63] = [
    ([0, 0, 0, 0, 1], [-171996.0, -174.2], [92025.0, 8.9]),
    ([-2, 0, 0, 2, 2], [-13187.0, -1.6], [5736.0, -3.1]),
    ([0, 0, 0, 2, 2], [-2274.0, -0.2], [977.0, -0.5]),
    ([0, 0, 0, 0, 2], [2062.0, 0.2], [-895.0, 0.5]),
    ([0, 1, 0, 0, 0], [1426.0, -3.4], [54.0, -0.1]),
    ([0, 0, 1, 0, 0], [712.0, 0.1], [-7.0, 0.0]),
    ([-2, 1, 0, 2, 2], [-517.0, 1.2], [224.0, -0.6]),
    ([0, 0, 0, 2, 1], [-386.0, -0.4], [200.0, 0.0]),
    ([0, 0, 1, 2, 2], [-301.0, 0.0], [129.0, -0.1]),
    ([-2, -1, 0, 2, 2], [217.0, -0.5], [-95.0, 0.3]),
    ([-2, 0, 1, 0, 0], [-158.0, 0.0], [0.0, 0.0]),
    ([-2, 0, 0, 2, 1], [129.0, 0.1], [-70.0, 0.0]),
    ([0, 0, -1, 2, 2], [123.0, 0.0], [-53.0, 0.0]),
    ([2, 0, 0, 0, 0], [63.0, 0.0], [0.0, 0.0]),
    ([0, 0, 1, 0, 1], [63.0, 0.1], [-33.0, 0.0]),
    ([2, 0, -1, 2, 2], [-59.0, 0.0], [26.0, 0.0]),
    ([0, 0, -1, 0, 1], [-58.0, -0.1], [32.0, 0.0]),
    ([0, 0, 1, 2, 1], [-51.0, 0.0], [27.0, 0.0]),
    ([-2, 0, 2, 0, 0], [48.0, 0.0], [0.0, 0.0]),
    ([0, 0, -2, 2, 1], [46.0, 0.0], [-24.0, 0.0]),
    ([2, 0, 0, 2, 2], [-38.0, 0.0], [16.0, 0.0]),
    ([0, 0, 2, 2, 2], [-31.0, 0.0], [13.0, 0.0]),
    ([0, 0, 2, 0, 0], [29.0, 0.0], [0.0, 0.0]),
    ([-2, 0, 1, 2, 2], [29.0, 0.0], [-12.0, 0.0]),
    ([0, 0, 0, 2, 0], [26.0, 0.0], [0.0, 0.0]),
    ([-2, 0, 0, 2, 0], [-22.0, 0.0], [0.0, 0.0]),
    ([0, 0, -1, 2, 1], [21.0, 0.0], [-10.0, 0.0]),
    ([0, 2, 0, 0, 0], [17.0, -0.1], [0.0, 0.0]),
    ([2, 0, -1, 0, 1], [16.0, 0.0], [-8.0, 0.0]),
    ([-2, 2, 0, 2, 2], [-16.0, 0.1], [7.0, 0.0]),
    ([0, 1, 0, 0, 1], [-15.0, 0.0], [9.0, 0.0]),
    ([-2, 0, 1, 0, 1], [-13.0, 0.0], [7.0, 0.0]),
    ([0, -1, 0, 0, 1], [-12.0, 0.0], [6.0, 0.0]),
    ([0, 0, 2, -2, 0], [11.0, 0.0], [0.0, 0.0]),
    ([2, 0, -1, 2, 1], [-10.0, 0.0], [5.0, 0.0]),
    ([2, 0, 1, 2, 2], [-8.0, 0.0], [3.0, 0.0]),
    ([0, 1, 0, 2, 2], [7.0, 0.0], [-3.0, 0.0]),
    ([-2, 1, 1, 0, 0], [-7.0, 0.0], [0.0, 0.0]),
    ([0, -1, 0, 2, 2], [-7.0, 0.0], [3.0, 0.0]),
    ([2, 0, 0, 2, 1], [-7.0, 0.0], [3.0, 0.0]),
    ([2, 0, 1, 0, 0], [6.0, 0.0], [0.0, 0.0]),
    ([-2, 0, 2, 2, 2], [6.0, 0.0], [-3.0, 0.0]),
    ([-2, 0, 1, 2, 1], [6.0, 0.0], [-3.0, 0.0]),
    ([2, 0, -2, 0, 1], [-6.0, 0.0], [3.0, 0.0]),
    ([2, 0, 0, 0, 1], [-6.0, 0.0], [3.0, 0.0]),
    ([0, -1, 1, 0, 0], [5.0, 0.0], [0.0, 0.0]),
    ([-2, -1, 0, 2, 1], [-5.0, 0.0], [3.0, 0.0]),
    ([-2, 0, 0, 0, 1], [-5.0, 0.0], [3.0, 0.0]),
    ([0, 0, 2, 2, 1], [-5.0, 0.0], [3.0, 0.0]),
    ([-2, 0, 2, 0, 1], [4.0, 0.0], [0.0, 0.0]),
    ([-2, 1, 0, 2, 1], [4.0, 0.0], [0.0, 0.0]),
    ([0, 0, 1, -2, 0], [4.0, 0.0], [0.0, 0.0]),
    ([-1, 0, 1, 0, 0], [-4.0, 0.0], [0.0, 0.0]),
    ([-2, 1, 0, 0, 0], [-4.0, 0.0], [0.0, 0.0]),
    ([1, 0, 0, 0, 0], [-4.0, 0.0], [0.0, 0.0]),
    ([0, 0, 1, 2, 0], [3.0, 0.0], [0.0, 0.0]),
    ([0, 0, -2, 2, 2], [-3.0, 0.0], [0.0, 0.0]),
    ([-1, -1, 1, 0, 0], [-3.0, 0.0], [0.0, 0.0]),
    ([0, 1, 1, 0, 0], [-3.0, 0.0], [0.0, 0.0]),
    ([0, -1, 1, 2, 2], [-3.0, 0.0], [0.0, 0.0]),
    ([2, -1, -1, 2, 2], [-3.0, 0.0], [0.0, 0.0]),
    ([0, 0, 3, 2, 2], [-3.0, 0.0], [0.0, 0.0]),
    ([2, -1, 0, 2, 2], [-3.0, 0.0], [0.0, 0.0]),
];

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

    use super::nutation;

    #[test_case(2460665.885255802, -0.0000020569363459503144)]
    fn expect_nutation_longitude(date: f64, expected: f64) {
        assert_approx_eq!(f64, nutation(date), expected);
    }
}