solar_positioning/
grena3.rs

1//! Grena3 solar position algorithm implementation.
2//!
3//! This follows the no. 3 algorithm described by Roberto Grena of ENEA (2012).
4//!
5//! The algorithm is designed for the years 2010 to 2110, with a maximum error of 0.01 degrees.
6//! It's approximately 10x faster than the SPA algorithm but with reduced accuracy and time range.
7//!
8//! Reference: Grena, R. (2012). Five new algorithms for the computation of sun position from 2010 to 2110.
9//! Solar Energy, 86(5), 1323-1337. DOI: <http://dx.doi.org/10.1016/j.solener.2012.01.024>
10
11#![allow(clippy::unreadable_literal)]
12#![allow(clippy::many_single_char_names)]
13#![allow(clippy::suboptimal_flops)]
14
15use crate::error::check_coordinates;
16use crate::math::{
17    PI, asin, atan2, cos, degrees_to_radians, normalize_degrees_0_to_360, radians_to_degrees, sin,
18    sqrt, tan,
19};
20use crate::{RefractionCorrection, Result, SolarPosition};
21use chrono::{DateTime, Datelike, TimeZone, Timelike};
22
23/// Calculate solar position using the Grena3 algorithm.
24///
25/// This is a simplified algorithm designed for years 2010-2110 with maximum error of 0.01°.
26/// It's much faster than SPA but less accurate and has a limited time range.
27///
28/// # Arguments
29/// * `datetime` - Timezone-aware date and time
30/// * `latitude` - Observer latitude in degrees (-90 to +90)
31/// * `longitude` - Observer longitude in degrees (-180 to +180)
32/// * `delta_t` - ΔT in seconds (difference between TT and UT1)
33///
34/// # Returns
35/// Solar position or error
36///
37/// # Errors
38/// Returns error for invalid coordinates (latitude outside ±90°, longitude outside ±180°)
39///
40/// # Example
41/// ```rust
42/// use solar_positioning::{grena3, RefractionCorrection};
43/// use chrono::{DateTime, FixedOffset};
44///
45/// let datetime = "2023-06-21T12:00:00-07:00".parse::<DateTime<FixedOffset>>().unwrap();
46///
47/// // Without refraction correction
48/// let position = grena3::solar_position(
49///     datetime,
50///     37.7749,     // San Francisco latitude
51///     -122.4194,   // San Francisco longitude
52///     69.0,        // deltaT (seconds)
53///     None,        // no atmospheric refraction correction
54/// ).unwrap();
55///
56/// // With refraction correction
57/// let position_with_refraction = grena3::solar_position(
58///     datetime,
59///     37.7749,     // San Francisco latitude
60///     -122.4194,   // San Francisco longitude
61///     69.0,        // deltaT (seconds)
62///     Some(RefractionCorrection::standard()), // standard atmospheric conditions
63/// ).unwrap();
64///
65/// println!("Azimuth: {:.3}°", position.azimuth());
66/// println!("Elevation: {:.3}°", position.elevation_angle());
67/// ```
68#[allow(clippy::needless_pass_by_value)]
69pub fn solar_position<Tz: TimeZone>(
70    datetime: DateTime<Tz>,
71    latitude: f64,
72    longitude: f64,
73    delta_t: f64,
74    refraction: Option<RefractionCorrection>,
75) -> Result<SolarPosition> {
76    check_coordinates(latitude, longitude)?;
77
78    // Calculate t (days since 2000-01-01 12:00:00 TT)
79    let t = calc_t(&datetime);
80    let t_e = t + 1.1574e-5 * delta_t;
81    let omega_at_e = 0.0172019715 * t_e;
82
83    // Calculate apparent sun longitude (lambda)
84    let lambda = -1.388803
85        + 1.720279216e-2 * t_e
86        + 3.3366e-2 * sin(omega_at_e - 0.06172)
87        + 3.53e-4 * sin(2.0 * omega_at_e - 0.1163);
88
89    // Calculate obliquity of ecliptic (epsilon)
90    let epsilon = 4.089567e-1 - 6.19e-9 * t_e;
91
92    let s_lambda = sin(lambda);
93    let c_lambda = cos(lambda);
94    let s_epsilon = sin(epsilon);
95    let c_epsilon = sqrt(1.0 - s_epsilon * s_epsilon);
96
97    // Calculate right ascension (alpha)
98    let mut alpha = atan2(s_lambda * c_epsilon, c_lambda);
99    if alpha < 0.0 {
100        alpha += 2.0 * PI;
101    }
102
103    // Calculate declination (delta)
104    let delta = asin(s_lambda * s_epsilon);
105
106    // Calculate hour angle (H)
107    let mut h = 1.7528311 + 6.300388099 * t + degrees_to_radians(longitude) - alpha;
108    h = ((h + PI) % (2.0 * PI)) - PI;
109    if h < -PI {
110        h += 2.0 * PI;
111    }
112
113    // Calculate topocentric coordinates
114    let s_phi = sin(degrees_to_radians(latitude));
115    let c_phi = sqrt(1.0 - s_phi * s_phi);
116    let s_delta = sin(delta);
117    let c_delta = sqrt(1.0 - s_delta * s_delta);
118    let s_h = sin(h);
119    let c_h = cos(h);
120
121    let s_epsilon0 = s_phi * s_delta + c_phi * c_delta * c_h;
122    let e_p = asin(s_epsilon0) - 4.26e-5 * sqrt(1.0 - s_epsilon0 * s_epsilon0);
123    let gamma = atan2(s_h, c_h * s_phi - s_delta * c_phi / c_delta);
124
125    // Apply refraction correction if provided and sun is visible
126    let delta_re = refraction.map_or(0.0, |correction| {
127        if e_p > 0.0 {
128            let pressure = correction.pressure();
129            let temperature = correction.temperature();
130            (0.08422 * (pressure / 1000.0))
131                / ((273.0 + temperature) * tan(e_p + 0.003138 / (e_p + 0.08919)))
132        } else {
133            0.0
134        }
135    });
136
137    let z = PI / 2.0 - e_p - delta_re;
138
139    let azimuth = normalize_degrees_0_to_360(radians_to_degrees(gamma + PI));
140    let zenith = radians_to_degrees(z);
141
142    SolarPosition::new(azimuth, zenith)
143}
144
145/// Calculate t parameter (days since 2000-01-01 12:00:00 TT)
146fn calc_t<Tz: TimeZone>(datetime: &DateTime<Tz>) -> f64 {
147    // Convert to UTC for proper astronomical calculations
148    let utc_datetime = datetime.with_timezone(&chrono::Utc);
149    let mut m = i32::try_from(utc_datetime.month()).expect("month should fit in i32");
150    let mut y = utc_datetime.year();
151    let d = i32::try_from(utc_datetime.day()).expect("day should fit in i32");
152    let h = f64::from(utc_datetime.hour())
153        + f64::from(utc_datetime.minute()) / 60.0
154        + f64::from(utc_datetime.second()) / 3600.0;
155
156    if m <= 2 {
157        m += 12;
158        y -= 1;
159    }
160
161    f64::from((365.25 * f64::from(y - 2000)) as i32)
162        + f64::from((30.6001 * f64::from(m + 1)) as i32)
163        - f64::from((0.01 * f64::from(y)) as i32)
164        + f64::from(d)
165        + 0.0416667 * h
166        - 21958.0
167}
168
169#[cfg(test)]
170mod tests {
171    use super::*;
172    use chrono::{DateTime, FixedOffset};
173
174    #[test]
175    fn test_grena3_basic_functionality() {
176        let datetime = "2023-06-21T12:00:00-07:00"
177            .parse::<DateTime<FixedOffset>>()
178            .unwrap();
179
180        let result = solar_position(datetime, 37.7749, -122.4194, 69.0, None);
181
182        assert!(result.is_ok());
183        let position = result.unwrap();
184        assert!(position.azimuth() >= 0.0 && position.azimuth() <= 360.0);
185        assert!(position.zenith_angle() >= 0.0 && position.zenith_angle() <= 180.0);
186    }
187
188    #[test]
189    fn test_grena3_with_refraction() {
190        let datetime = "2023-06-21T12:00:00-07:00"
191            .parse::<DateTime<FixedOffset>>()
192            .unwrap();
193
194        let result = solar_position(
195            datetime,
196            37.7749,
197            -122.4194,
198            69.0,
199            Some(RefractionCorrection::new(1013.25, 15.0).unwrap()),
200        );
201
202        assert!(result.is_ok());
203        let position = result.unwrap();
204        assert!(position.azimuth() >= 0.0 && position.azimuth() <= 360.0);
205        assert!(position.zenith_angle() >= 0.0 && position.zenith_angle() <= 180.0);
206    }
207
208    #[test]
209    fn test_grena3_coordinate_validation() {
210        let datetime = "2023-06-21T12:00:00-07:00"
211            .parse::<DateTime<FixedOffset>>()
212            .unwrap();
213
214        // Invalid latitude
215        assert!(solar_position(datetime, 95.0, 0.0, 0.0, None).is_err());
216
217        // Invalid longitude
218        assert!(solar_position(datetime, 0.0, 185.0, 0.0, None).is_err());
219    }
220
221    #[test]
222    fn test_calc_t() {
223        // Test with a known date
224        let datetime = "2023-06-21T12:00:00-07:00"
225            .parse::<DateTime<FixedOffset>>()
226            .unwrap();
227        let t = calc_t(&datetime);
228
229        // The Grena3 algorithm uses a specific reference point that may result in negative values
230        // This is correct behavior - just ensure the calculation is consistent
231        assert!(t.is_finite(), "t should be finite");
232
233        // Test that the calculation is stable
234        let t2 = calc_t(&datetime);
235        assert!(
236            (t - t2).abs() < f64::EPSILON,
237            "calc_t should be deterministic"
238        );
239
240        // Test that different dates give different results
241        let datetime2 = "2023-06-22T12:00:00Z"
242            .parse::<DateTime<FixedOffset>>()
243            .unwrap();
244        let t3 = calc_t(&datetime2);
245        assert!(
246            (t - t3).abs() > 0.5,
247            "Different dates should give different t values"
248        );
249    }
250}
solar_positioning/grena3.rs

solar_positioning/
grena3.rs
