//! Just enough astronomy for lunar calendars: where the Sun and Moon are, and
//! when the Moon reaches a phase. Formulas are the shorter series from Jean
//! Meeus' "Astronomical Algorithms", good to a few minutes of time.
use chrono::{DateTime, Duration, NaiveDate, NaiveDateTime, Utc};
/// Julian day of a UTC time
pub fn jd(t: NaiveDateTime) -> f64 {
let epoch = NaiveDate::from_ymd_opt(2000, 1, 1)
.and_then(|d| d.and_hms_opt(12, 0, 0))
.unwrap_or_default();
2451545.0 + (t - epoch).num_seconds() as f64 / 86400.0
}
/// UTC time of a Julian day
pub fn utc(jd: f64) -> DateTime<Utc> {
let secs = ((jd - 2451545.0) * 86400.0).round() as i64;
DateTime::<Utc>::from_timestamp(946728000 + secs, 0).unwrap_or_default()
}
fn norm(deg: f64) -> f64 {
deg.rem_euclid(360.0)
}
fn sin(deg: f64) -> f64 {
deg.to_radians().sin()
}
fn centuries(jd: f64) -> f64 {
(jd - 2451545.0) / 36525.0
}
/// The Sun's apparent ecliptic longitude, in degrees
pub fn sun_longitude(jd: f64) -> f64 {
let t = centuries(jd);
let l0 = 280.46646 + 36000.76983 * t + 0.0003032 * t * t;
let m = 357.52911 + 35999.05029 * t - 0.0001537 * t * t;
let c = (1.914602 - 0.004817 * t - 0.000014 * t * t) * sin(m)
+ (0.019993 - 0.000101 * t) * sin(2.0 * m)
+ 0.000289 * sin(3.0 * m);
let omega = 125.04 - 1934.136 * t;
norm(l0 + c - 0.00569 - 0.00478 * sin(omega))
}
/// The Moon's apparent ecliptic longitude, in degrees
pub fn moon_longitude(jd: f64) -> f64 {
let t = centuries(jd);
let t2 = t * t;
let t3 = t2 * t;
let t4 = t3 * t;
let lp = 218.3164477 + 481267.88123421 * t - 0.0015786 * t2 + t3 / 538841.0 - t4 / 65194000.0;
let d = 297.8501921 + 445267.1114034 * t - 0.0018819 * t2 + t3 / 545868.0 - t4 / 113065000.0;
let m = 357.5291092 + 35999.0502909 * t - 0.0001536 * t2 + t3 / 24490000.0;
let mp = 134.9633964 + 477198.8675055 * t + 0.0087414 * t2 + t3 / 69699.0 - t4 / 14712000.0;
let f = 93.2720950 + 483202.0175233 * t - 0.0036539 * t2 - t3 / 3526000.0 + t4 / 863310000.0;
let e = 1.0 - 0.002516 * t - 0.0000074 * t2;
// The largest terms of Meeus' table 47.A: multiples of D, M, M', F and the coefficient in 1e-6 degrees
#[rustfmt::skip]
const TERMS: [(f64, f64, f64, f64, f64); 34] = [
(0.0, 0.0, 1.0, 0.0, 6288774.0), (2.0, 0.0, -1.0, 0.0, 1274027.0),
(2.0, 0.0, 0.0, 0.0, 658314.0), (0.0, 0.0, 2.0, 0.0, 213618.0),
(0.0, 1.0, 0.0, 0.0, -185116.0), (0.0, 0.0, 0.0, 2.0, -114332.0),
(2.0, 0.0, -2.0, 0.0, 58793.0), (2.0, -1.0, -1.0, 0.0, 57066.0),
(2.0, 0.0, 1.0, 0.0, 53322.0), (2.0, -1.0, 0.0, 0.0, 45758.0),
(0.0, 1.0, -1.0, 0.0, -40923.0), (1.0, 0.0, 0.0, 0.0, -34720.0),
(0.0, 1.0, 1.0, 0.0, -30383.0), (2.0, 0.0, 0.0, -2.0, 15327.0),
(0.0, 0.0, 1.0, 2.0, -12528.0), (0.0, 0.0, 1.0, -2.0, 10980.0),
(4.0, 0.0, -1.0, 0.0, 10675.0), (0.0, 0.0, 3.0, 0.0, 10034.0),
(4.0, 0.0, -2.0, 0.0, 8548.0), (2.0, 1.0, -1.0, 0.0, -7888.0),
(2.0, 1.0, 0.0, 0.0, -6766.0), (1.0, 0.0, -1.0, 0.0, -5163.0),
(1.0, 1.0, 0.0, 0.0, 4987.0), (2.0, -1.0, 1.0, 0.0, 4036.0),
(2.0, 0.0, 2.0, 0.0, 3994.0), (4.0, 0.0, 0.0, 0.0, 3861.0),
(2.0, 0.0, -3.0, 0.0, 3665.0), (0.0, 1.0, -2.0, 0.0, -2689.0),
(2.0, 0.0, -1.0, 2.0, -2602.0), (2.0, -1.0, -2.0, 0.0, 2390.0),
(1.0, 0.0, 1.0, 0.0, -2348.0), (2.0, -2.0, 0.0, 0.0, 2236.0),
(0.0, 1.0, 2.0, 0.0, -2120.0), (0.0, 2.0, 0.0, 0.0, -2069.0),
];
let mut sum = 0.0;
for (cd, cm, cmp, cf, coeff) in TERMS {
// Terms with the Sun's anomaly shrink as Earth's orbit becomes less eccentric
let ecc = match cm.abs() as i32 {
1 => e,
2 => e * e,
_ => 1.0,
};
sum += coeff * ecc * sin(cd * d + cm * m + cmp * mp + cf * f);
}
let a1 = 119.75 + 131.849 * t;
let a2 = 53.09 + 479264.290 * t;
sum += 3958.0 * sin(a1) + 1962.0 * sin(lp - f) + 318.0 * sin(a2);
let omega = 125.04 - 1934.136 * t;
norm(lp + sum / 1e6 - 0.00478 * sin(omega))
}
/// How far the Moon is ahead of the Sun: 0 at new moon, 180 at full
pub fn elongation(jd: f64) -> f64 {
norm(moon_longitude(jd) - sun_longitude(jd))
}
/// The Lahiri ayanamsa: how far the sidereal zodiac lags the tropical one
fn ayanamsa(jd: f64) -> f64 {
23.853 + (jd - 2451545.0) / 365.25 * 0.013969
}
pub fn sidereal_sun_longitude(jd: f64) -> f64 {
norm(sun_longitude(jd) - ayanamsa(jd))
}
/// Every time between `from` and `to` that `f` passes `target` degrees, `f` always increasing
fn crossings(from: f64, to: f64, target: f64, step: f64, f: impl Fn(f64) -> f64) -> Vec<f64> {
let off = |jd: f64| (f(jd) - target + 540.0).rem_euclid(360.0) - 180.0;
let mut out = vec![];
let mut a = from;
while a < to {
let b = a + step;
if off(a) < 0.0 && off(b) >= 0.0 {
let (mut lo, mut hi) = (a, b);
for _ in 0..40 {
let mid = (lo + hi) / 2.0;
if off(mid) < 0.0 {
lo = mid
} else {
hi = mid
}
}
out.push((lo + hi) / 2.0);
}
a = b;
}
out
}
/// When the Moon reaches a phase angle (0 new, 90 first quarter, 180 full, 270 last quarter)
pub fn moon_phases(from: f64, to: f64, angle: f64) -> Vec<f64> {
crossings(from, to, angle, 1.0, elongation)
}
/// When the Sun enters a sidereal longitude, e.g. 270 for Makara (Capricorn)
pub fn sidereal_ingress(from: f64, to: f64, longitude: f64) -> Vec<f64> {
crossings(from, to, longitude, 5.0, sidereal_sun_longitude)
}
/// When the Sun reaches a tropical longitude: 0 and 180 are the equinoxes
pub fn sun_ingress(from: f64, to: f64, longitude: f64) -> Vec<f64> {
crossings(from, to, longitude, 5.0, sun_longitude)
}
/// Julian day of midnight UTC starting a date, shifted by hours
pub fn jd_of(date: NaiveDate, hours: f64) -> f64 {
jd(date.and_hms_opt(0, 0, 0).unwrap_or_default()) + hours / 24.0
}
/// Sunrise on a date at a latitude and longitude (east positive), as a Julian day.
/// The NOAA approximation, good to a minute or two away from the poles.
pub fn sunrise(date: NaiveDate, lat: f64, lon: f64) -> f64 {
let noon = jd_of(date, 12.0 - lon / 15.0);
let t = centuries(noon);
let l0 = norm(280.46646 + 36000.76983 * t);
let m = 357.52911 + 35999.05029 * t;
let e = 0.016708634 - 0.000042037 * t;
let eps = 23.439291 - 0.0130042 * t;
let decl = (sin(eps) * sin(sun_longitude(noon))).asin();
let y = (eps / 2.0).to_radians().tan().powi(2);
let (l0r, mr) = (l0.to_radians(), m.to_radians());
// Equation of time, in minutes
let eot = 4.0
* (y * (2.0 * l0r).sin() - 2.0 * e * mr.sin() + 4.0 * e * y * mr.sin() * (2.0 * l0r).cos()
- 0.5 * y * y * (4.0 * l0r).sin()
- 1.25 * e * e * (2.0 * mr).sin())
.to_degrees();
let phi = lat.to_radians();
let cos_h =
((-0.833f64).to_radians().sin() - phi.sin() * decl.sin()) / (phi.cos() * decl.cos());
let h = cos_h.clamp(-1.0, 1.0).acos().to_degrees();
let minutes_utc = 720.0 - 4.0 * lon - eot - 4.0 * h;
jd_of(date, minutes_utc / 60.0)
}
/// A UTC time as a date and time in a fixed offset from UTC
pub fn local(jd: f64, offset_hours: f64) -> NaiveDateTime {
utc(jd).naive_utc() + Duration::seconds((offset_hours * 3600.0) as i64)
}
#[cfg(test)]
mod tests {
use super::*;
fn at(y: i32, m: u32, d: u32, h: u32, min: u32) -> f64 {
jd(NaiveDate::from_ymd_opt(y, m, d)
.unwrap()
.and_hms_opt(h, min, 0)
.unwrap())
}
#[test]
fn phases_match_eclipses() {
// The 2024-04-08 total solar eclipse new moon was at 18:21 UTC
let nm = moon_phases(at(2024, 4, 1, 0, 0), at(2024, 4, 15, 0, 0), 0.0);
assert_eq!(nm.len(), 1);
assert!(
(nm[0] - at(2024, 4, 8, 18, 21)).abs() < 10.0 / 1440.0,
"{}",
utc(nm[0])
);
// The 2024-09-18 lunar eclipse full moon was at 02:34 UTC
let fm = moon_phases(at(2024, 9, 10, 0, 0), at(2024, 9, 25, 0, 0), 180.0);
assert!(
(fm[0] - at(2024, 9, 18, 2, 34)).abs() < 10.0 / 1440.0,
"{}",
utc(fm[0])
);
}
#[test]
fn sunrise_in_delhi() {
// 2026-06-21 05:24 and 2026-12-21 07:10 IST
let day = |m, d| NaiveDate::from_ymd_opt(2026, m, d).unwrap();
let near = |jd: f64, d: NaiveDate, h: u32, m: u32| {
let want = super::jd(d.and_hms_opt(h, m, 0).unwrap()) - 5.5 / 24.0;
(jd - want).abs() < 2.0 / 1440.0
};
assert!(near(sunrise(day(6, 21), 28.61, 77.21), day(6, 21), 5, 24));
assert!(near(sunrise(day(12, 21), 28.61, 77.21), day(12, 21), 7, 10));
}
}
⭐
Log in to mark this repository.
8.7 KiB · text
History
d9460a8