sun_times/
lib.rs

1use chrono::{DateTime, NaiveDate, TimeZone, Timelike, Utc};
2
3const UNIX_EPOCH: JulianDate = JulianDate(2440587.5);
4const SECONDS_PER_DAY: u64 = 24 * 60 * 60;
5const JAN_2000: JulianDate = JulianDate(2451545.0);
6const LEAP_SECONDS: JulianDate = JulianDate(0.0008);
7const OBLIQUITY_OF_THE_ECLIPTIC: f64 = 23.44;
8
9#[derive(Debug, Clone, Copy)]
10struct JulianDate(f64);
11
12impl JulianDate {
13    fn ceil_days(&self) -> f64 {
14        self.0.ceil()
15    }
16
17    fn to_datetime(self) -> Option<DateTime<Utc>> {
18        Utc.timestamp_opt(
19            ((self - UNIX_EPOCH).0 * SECONDS_PER_DAY as f64).round() as i64,
20            0,
21        )
22        .single()
23    }
24}
25
26impl From<DateTime<Utc>> for JulianDate {
27    fn from(date: DateTime<Utc>) -> Self {
28        Self((date.timestamp() as f64 / SECONDS_PER_DAY as f64) + UNIX_EPOCH.0)
29    }
30}
31
32impl std::ops::Sub<JulianDate> for JulianDate {
33    type Output = Self;
34
35    fn sub(self, rhs: JulianDate) -> Self::Output {
36        Self(self.0 - rhs.0)
37    }
38}
39
40impl std::ops::Add<JulianDate> for JulianDate {
41    type Output = Self;
42
43    fn add(self, rhs: JulianDate) -> Self::Output {
44        Self(self.0 + rhs.0)
45    }
46}
47
48/// Calculates the approximate sunset and sunrise times at a given latitude, longitude, and altitude
49///
50/// Note that elevation is used to correct for atmospheric refraction, so negative elevations are treated as being at
51/// sea level due to having minimal difference in refraction to being at sea level
52///
53/// # Arguments
54///
55/// * `date` - The date on which to calculate the sunset and sunrise, in UTC
56/// * `latitude` - The latitude at which to calculate the times. Expressed as degrees
57/// * `longitude` - The longitude at which to calculate the times. Expressed as degrees
58/// * `elevation` - The elevation at which to calculate the times. Expressed as meters above sea level. Negative values will be ignored
59///
60/// # Return value
61///
62/// Returns
63///  - `None` if the date is not representable in chrono (~5M years from now), or sunsets/rises cannot be calculated due to long arctic/antarctic day/night (outside ~±67° of latitude)
64///  - `Some((sunrise,sunset))` otherwise
65///
66/// # Examples
67///
68/// ```
69/// //Calculate the sunset and sunrise times today at Sheffield university's new computer science building
70/// let times = sun_times(Utc::today(),53.38,-1.48,100.0);
71/// println!("Sunrise: {}, Sunset: {}",times.0,times.1);
72/// ```
73pub fn sun_times(
74    date: NaiveDate,
75    latitude: f64,
76    longitude: f64,
77    elevation: f64,
78) -> Option<(DateTime<Utc>, DateTime<Utc>)> {
79    //see https://en.wikipedia.org/wiki/Sunrise_equation
80
81    const ARGUMENT_OF_PERIHELION: f64 = 102.9372;
82
83    let julian_date = JulianDate::from(
84        date.and_hms_opt(0, 0, 0)?
85            .and_local_timezone(Utc)
86            .single()?,
87    );
88
89    //elevations below sea level will have minimal atmospheric refraction + the
90    //calculation is broken below sea level, so treat negative elevations as being at sea level
91    let elevation = elevation.max(0.0);
92    let elevation_correction = -2.076 * (elevation.sqrt()) / 60.0;
93
94    let days_since_2000 = (julian_date - JAN_2000 + LEAP_SECONDS).ceil_days();
95
96    let mean_solar_time = days_since_2000 - (longitude / 360.0);
97    let solar_mean_anomaly = (357.5291 + 0.98560028 * mean_solar_time).rem_euclid(360.0);
98    let center = 1.9148 * solar_mean_anomaly.to_radians().sin()
99        + 0.0200 * (2.0 * solar_mean_anomaly).to_radians().sin()
100        + 0.0003 * (3.0 * solar_mean_anomaly).to_radians().sin();
101    let ecliptic_longitude =
102        (solar_mean_anomaly + center + 180.0 + ARGUMENT_OF_PERIHELION).rem_euclid(360.0);
103
104    let declination = (ecliptic_longitude.to_radians().sin()
105        * OBLIQUITY_OF_THE_ECLIPTIC.to_radians().sin())
106    .asin();
107    let event_hour_angle = (((-0.83 + elevation_correction).to_radians().sin()
108        - (latitude.to_radians().sin() * declination.sin()))
109        / (latitude.to_radians().cos() * declination.cos()))
110    .acos()
111    .to_degrees();
112
113    if event_hour_angle.is_nan() {
114        return None;
115    }
116
117    let solar_transit =
118        JAN_2000.0 + mean_solar_time + 0.0053 * solar_mean_anomaly.to_radians().sin()
119            - 0.0069 * (2.0 * ecliptic_longitude).to_radians().sin();
120    let solar_transit_julian = JulianDate(solar_transit);
121
122    let julian_rise = JulianDate(solar_transit_julian.0 - event_hour_angle / 360.0);
123    let julian_set = JulianDate(solar_transit_julian.0 + event_hour_angle / 360.0);
124    let rise = julian_rise.to_datetime();
125    let set = julian_set.to_datetime();
126    if let (Some(rise), Some(set)) = (rise, set) {
127        Some((rise, set))
128    } else {
129        None
130    }
131}
132
133/// Calculates the altitude (angle from the horizon) of the sun at a given place and moment
134/// # Arguments
135///
136/// * `date_time` - The date and time on which to calculate the sunset and sunrise
137/// * `latitude` - The latitude at which to calculate the times. Expressed as degrees
138/// * `longitude` - The longitude at which to calculate the times. Expressed as degrees
139///
140/// # Return value
141///
142/// Returns the altitude in degrees
143///
144/// # Notes
145///
146/// This currently does not include an altitude correction like the [sun_times] function does, but *is* usable within arctic regions
147///
148/// # Examples
149///
150/// ```
151/// //Calculate the sunset and sunrise times today at Sheffield university's new computer science building
152/// let altitude = altitude(Utc::now(),53.38,-1.48);
153/// println!("Altitude: {}",altitude);
154/// ```
155pub fn altitude(date_time: DateTime<Utc>, latitude: f64, longitude: f64) -> f64 {
156    //see https://en.wikipedia.org/wiki/Sunrise_equation
157    //see https://en.wikipedia.org/wiki/Astronomical_coordinate_systems
158    //see http://www.stargazing.net/kepler/altaz.html
159
160    const ARGUMENT_OF_PERIHELION: f64 = 102.9372;
161
162    let julian_date = JulianDate::from(date_time);
163
164    let days_since_2000 = (julian_date - JAN_2000 + LEAP_SECONDS).ceil_days();
165
166    let mean_solar_time = days_since_2000 - (longitude / 360.0);
167    let solar_mean_anomaly = (357.5291 + 0.98560028 * mean_solar_time).rem_euclid(360.0);
168    let center = 1.9148 * solar_mean_anomaly.to_radians().sin()
169        + 0.0200 * (2.0 * solar_mean_anomaly).to_radians().sin()
170        + 0.0003 * (3.0 * solar_mean_anomaly).to_radians().sin();
171    let ecliptic_longitude =
172        (solar_mean_anomaly + center + 180.0 + ARGUMENT_OF_PERIHELION).rem_euclid(360.0);
173
174    let sin_declination =
175        ecliptic_longitude.to_radians().sin() * OBLIQUITY_OF_THE_ECLIPTIC.to_radians().sin();
176    let declination = sin_declination.asin();
177
178    let right_ascension = (ecliptic_longitude.to_radians().sin()
179        * OBLIQUITY_OF_THE_ECLIPTIC.to_radians().cos())
180    .atan2(ecliptic_longitude.to_radians().cos())
181    .to_degrees();
182
183    let greenwich_sidereal_time = mean_solar_time + 0.0;
184    let local_sideral_time = greenwich_sidereal_time
185        + (date_time.time().hour() as f64
186            + (date_time.time().minute() as f64 / 60.0)
187            + (date_time.time().second() as f64 / 60.0 * 60.0))
188            * 15.0
189        + longitude.to_degrees();
190    let local_hour_angle = local_sideral_time - right_ascension;
191
192    let sin_altitude = (latitude.to_radians().sin() * declination.sin())
193        + (latitude.to_radians().cos() * declination.cos() * local_hour_angle.to_radians().cos());
194    sin_altitude.asin().to_degrees()
195}
196
197#[cfg(test)]
198mod tests {
199    use chrono::{Duration, NaiveDate};
200
201    #[test]
202    ///Test for https://github.com/Eroc33/sun-times/issues/1
203    fn sunrise_and_sunset_land_on_requested_day() {
204        let date_range =
205            std::iter::successors(Some(NaiveDate::from_ymd_opt(2022, 1, 1).unwrap()), |date| {
206                let next = *date + Duration::days(1);
207                if next > NaiveDate::from_ymd_opt(2022, 12, 31).unwrap() {
208                    None
209                } else {
210                    Some(next)
211                }
212            });
213        for date in date_range {
214            let times = super::sun_times(date, 53.38, -1.48, 0.0);
215            assert!(times.is_some());
216            let times = times.unwrap();
217            assert_eq!(date, times.0.naive_utc().date());
218            assert_eq!(date, times.1.naive_utc().date());
219        }
220    }
221}
sun_times/lib.rs

sun_times/
lib.rs
