helioxide/
equatorial.rs

1// SPDX-License-Identifier: MIT OR Apache-2.0
2// helioxide - Rust implementation of NREL Solar Position Algorithm (SPA)
3// Copyright (c) 2026 Juan Luis Leal Contreras (Kuenlun)
4
5//! Sun geocentric right ascension (`α`) and declination (`δ`). Sections 3.9 and 3.10.
6
7use crate::helper::limit_degrees;
8
9/// `α = atan2(sin λ · cos ε − tan β · sin ε, cos λ)` (degrees, wrapped into
10/// `[0°, 360°)`). Equation 30.
11#[inline]
12#[must_use]
13pub fn geocentric_right_ascension(
14    apparent_sun_longitude: f64,
15    geocentric_latitude: f64,
16    true_obliquity: f64,
17) -> f64 {
18    let (sin_lambda, cos_lambda) = apparent_sun_longitude.to_radians().sin_cos();
19    let (sin_epsilon, cos_epsilon) = true_obliquity.to_radians().sin_cos();
20    let tan_beta = geocentric_latitude.to_radians().tan();
21
22    let numerator = sin_lambda.mul_add(cos_epsilon, -(tan_beta * sin_epsilon));
23    limit_degrees(numerator.atan2(cos_lambda).to_degrees())
24}
25
26/// `δ = arcsin(sin β · cos ε + cos β · sin ε · sin λ)` (degrees, signed,
27/// in `[-90°, 90°]`). Equation 31.
28#[inline]
29#[must_use]
30pub fn geocentric_declination(
31    apparent_sun_longitude: f64,
32    geocentric_latitude: f64,
33    true_obliquity: f64,
34) -> f64 {
35    let (sin_beta, cos_beta) = geocentric_latitude.to_radians().sin_cos();
36    let (sin_epsilon, cos_epsilon) = true_obliquity.to_radians().sin_cos();
37    let sin_lambda = apparent_sun_longitude.to_radians().sin();
38
39    (cos_beta * sin_epsilon)
40        .mul_add(sin_lambda, sin_beta * cos_epsilon)
41        .asin()
42        .to_degrees()
43}
44
45#[cfg(test)]
46#[cfg_attr(coverage_nightly, coverage(off))]
47mod tests {
48    use super::{geocentric_declination, geocentric_right_ascension};
49    use crate::apparent::{aberration_correction, apparent_sun_longitude};
50    use crate::geocentric::{geocentric_latitude, geocentric_longitude};
51    use crate::heliocentric::{
52        earth_heliocentric_latitude, earth_heliocentric_longitude, earth_radius_vector,
53    };
54    use crate::nutation::nutation_in_longitude_and_obliquity;
55    use crate::obliquity::true_obliquity_of_ecliptic;
56    use crate::test_fixtures::{reference_jce, reference_jme};
57
58    fn reference_lambda_beta_epsilon() -> (f64, f64, f64) {
59        let jme = reference_jme();
60        let theta = geocentric_longitude(earth_heliocentric_longitude(jme));
61        let beta = geocentric_latitude(earth_heliocentric_latitude(jme));
62        let (delta_psi, delta_epsilon) = nutation_in_longitude_and_obliquity(reference_jce());
63        let epsilon = true_obliquity_of_ecliptic(jme, delta_epsilon);
64        let delta_tau = aberration_correction(earth_radius_vector(jme));
65        let lambda = apparent_sun_longitude(theta, delta_psi, delta_tau);
66        (lambda, beta, epsilon)
67    }
68
69    #[test]
70    fn geocentric_right_ascension_matches_table_a5_1() {
71        let (lambda, beta, epsilon) = reference_lambda_beta_epsilon();
72        let alpha = geocentric_right_ascension(lambda, beta, epsilon);
73        assert!((alpha - 202.227_41).abs() < 1e-4);
74    }
75
76    #[test]
77    fn geocentric_declination_matches_table_a5_1() {
78        let (lambda, beta, epsilon) = reference_lambda_beta_epsilon();
79        let delta = geocentric_declination(lambda, beta, epsilon);
80        assert!((delta - -9.314_34).abs() < 1e-4);
81    }
82
83    #[test]
84    fn right_ascension_collapses_to_longitude_when_beta_and_epsilon_zero() {
85        for &lambda in &[15.0_f64, 75.0, 165.0, 200.0, 269.5, 350.0, -45.0, 720.5] {
86            let alpha = geocentric_right_ascension(lambda, 0.0, 0.0);
87            assert!((alpha - lambda.rem_euclid(360.0)).abs() < 1e-10);
88        }
89    }
90
91    #[test]
92    fn right_ascension_resolves_atan2_quadrant() {
93        for &(lambda, expected) in &[(180.0_f64, 180.0_f64), (270.0, 270.0)] {
94            let alpha = geocentric_right_ascension(lambda, 0.0, 23.44);
95            assert!((alpha - expected).abs() < 1e-10);
96        }
97    }
98
99    #[test]
100    fn right_ascension_pins_tan_beta_correction() {
101        for &(beta, eps) in &[(0.001_f64, 23.0), (-0.001, 23.0), (0.5, 45.0), (-0.5, 45.0)] {
102            let alpha = geocentric_right_ascension(0.0, beta, eps);
103            let expected = (-beta.to_radians().tan() * eps.to_radians().sin())
104                .atan()
105                .to_degrees()
106                .rem_euclid(360.0);
107            assert!((alpha - expected).abs() < 1e-12);
108        }
109    }
110
111    #[test]
112    fn right_ascension_wraps_into_zero_360() {
113        for &lambda in &[-720.0_f64, -10.0, 250.0, 1000.0] {
114            let alpha = geocentric_right_ascension(lambda, 0.000_1, 23.44);
115            assert!((0.0..360.0).contains(&alpha));
116        }
117    }
118
119    #[test]
120    fn declination_collapses_to_beta_when_epsilon_zero() {
121        for &beta in &[-30.0_f64, -1e-3, 1e-3, 30.0, 89.999] {
122            for &lambda in &[0.0_f64, 90.0, 180.0, 270.0] {
123                let delta = geocentric_declination(lambda, beta, 0.0);
124                assert!((delta - beta).abs() < 1e-10);
125            }
126        }
127    }
128
129    #[test]
130    fn declination_pins_sin_eps_sin_lambda_when_beta_zero() {
131        for &lambda in &[15.0_f64, 90.0, 150.0, 210.0, 300.0] {
132            for &eps in &[5.0_f64, 23.44, 45.0, 89.0] {
133                let delta = geocentric_declination(lambda, 0.0, eps);
134                let expected = (eps.to_radians().sin() * lambda.to_radians().sin())
135                    .asin()
136                    .to_degrees();
137                assert!((delta - expected).abs() < 1e-12);
138            }
139        }
140    }
141}
helioxide/equatorial.rs

helioxide/
equatorial.rs
