amlich_core/
sun.rs

1/**
2 * Sun Longitude Calculations
3 *
4 * Algorithm from "Astronomical Algorithms" by Jean Meeus, 1998
5 */
6use std::f64::consts::PI;
7
8/// Compute the longitude of the sun at any time
9///
10/// # Arguments
11/// * `jdn` - Julian day number (can be fractional)
12///
13/// # Returns
14/// Sun longitude in radians, normalized to [0, 2π)
15pub fn sun_longitude(jdn: f64) -> f64 {
16    // Time in Julian centuries from 2000-01-01 12:00:00 GMT
17    let t = (jdn - 2451545.0) / 36525.0;
18    let t2 = t * t;
19    let dr = PI / 180.0; // degree to radian
20
21    // Sun's mean anomaly, in degrees
22    let m = 357.52910 + 35999.05030 * t - 0.0001559 * t2 - 0.00000048 * t * t2;
23
24    // Sun's mean longitude, in degrees
25    let l0 = 280.46645 + 36000.76983 * t + 0.0003032 * t2;
26
27    // Sun's equation of center
28    let mut dl = (1.914600 - 0.004817 * t - 0.000014 * t2) * (dr * m).sin();
29    dl += (0.019993 - 0.000101 * t) * (dr * 2.0 * m).sin();
30    dl += 0.000290 * (dr * 3.0 * m).sin();
31
32    // Sun's true longitude, in degrees
33    let l = l0 + dl;
34
35    // Convert to radians
36    let mut l_rad = l * dr;
37
38    // Normalize to [0, 2π)
39    l_rad = l_rad - PI * 2.0 * (l_rad / (PI * 2.0)).floor();
40
41    l_rad
42}
43
44/// Compute sun position at midnight of the day with the given Julian day number
45///
46/// The function returns a number between 0 and 11.
47/// From the day after March equinox and the 1st major term after March equinox, 0 is returned.
48/// After that, return 1, 2, 3 ...
49///
50/// # Arguments
51/// * `day_number` - Julian day number
52/// * `time_zone` - Time zone offset (e.g., 7.0 for UTC+7:00)
53///
54/// # Returns
55/// Sun longitude index (0-11)
56pub fn get_sun_longitude(day_number: i32, time_zone: f64) -> i32 {
57    let jdn = day_number as f64 - 0.5 - time_zone / 24.0;
58    let longitude = sun_longitude(jdn);
59    (longitude / PI * 6.0).floor() as i32
60}
61
62#[cfg(test)]
63mod tests {
64    use super::*;
65
66    #[test]
67    fn test_sun_longitude_range() {
68        // Test that sun_longitude returns values in [0, 2π)
69        let test_jds = vec![
70            2451545.0, // 2000-01-01
71            2460349.0, // 2024-02-10
72            2460703.0, // 2025-01-29
73        ];
74
75        for jd in test_jds {
76            let sl = sun_longitude(jd);
77            assert!(
78                sl >= 0.0 && sl < 2.0 * PI,
79                "Sun longitude {} not in [0, 2π)",
80                sl
81            );
82        }
83    }
84
85    #[test]
86    fn test_get_sun_longitude_range() {
87        // Test that get_sun_longitude returns values in [0, 11]
88        let test_days = vec![
89            2451545, // 2000-01-01
90            2460349, // 2024-02-10
91            2460703, // 2025-01-29
92        ];
93
94        for day in test_days {
95            let sl = get_sun_longitude(day, 7.0);
96            assert!(
97                sl >= 0 && sl <= 11,
98                "Sun longitude index {} not in [0, 11]",
99                sl
100            );
101        }
102    }
103
104    #[test]
105    fn test_sun_longitude_consistency() {
106        // Test that sun longitude changes gradually over time
107        let jd_start = 2451545.0; // 2000-01-01
108
109        for i in 1..365 {
110            let jd = jd_start + i as f64;
111            let sl = sun_longitude(jd);
112
113            // Sun longitude should be between 0 and 2π
114            assert!(sl >= 0.0 && sl < 2.0 * PI);
115        }
116    }
117
118    #[test]
119    fn test_march_equinox_approximation() {
120        // Around March 20-21, sun longitude should be near 0 (spring equinox)
121        // Using 2024-03-20 as reference
122        let jd = 2460388; // 2024-03-20
123        let sl_index = get_sun_longitude(jd, 7.0);
124
125        // Should be 0 or close to transition
126        assert!(
127            sl_index == 0 || sl_index == 11,
128            "March equinox should have sun longitude index near 0, got {}",
129            sl_index
130        );
131    }
132}
amlich_core/sun.rs

amlich_core/
sun.rs
