A living almanac of seasons, soil, and the quiet knowledge that used to be common. Rust axum with minijinja and Vite.
agriculturealmanacaxumfolk-knowledgegardeningminijinjarustseasonalvite
1//! Astronomical math for the almanac. Self-contained, no external API calls.
2//!
3//! All formulas trace to Jean Meeus, *Astronomical Algorithms* (2nd ed.):
4//! ch. 7 julian day
5//! ch. 22 mean obliquity of the ecliptic
6//! ch. 25 sun apparent position
7//! ch. 28 equation of time
8//! ch. 47 moon position and phase (table 47.A perturbation series, truncated)
9//!
10//! Conventions:
11//! * Time variable `t` is always Julian centuries since J2000.0
12//! (JD 2451545.0 = 2000-01-01 12:00 TT). One Julian century = 36525 days.
13//! We treat UTC as TT; the offset is ~70 sec, well below render precision.
14//! * Angle <-> time conversions: 1 deg of arc = 4 min of time, 15 deg = 1 hr
15//! (earth rotates 360 deg per 24*60 min).
16//! * Latitude is north-positive, longitude east-positive (Raleigh's lon is
17//! therefore negative).
18
19use chrono::{DateTime, Datelike, Duration, TimeZone, Timelike, Utc};
20use chrono_tz::Tz;
21
22/// Latitude of Raleigh, NC (zone 7a target).
23const LAT: f64 = 35.78;
24/// Longitude of Raleigh, NC. East-positive, so this is negative.
25const LON: f64 = -78.64;
26/// Reference epoch for every "centuries since" calculation.
27const J2000_JD: f64 = 2451545.0;
28/// Days in one Julian century (the unit of `t` everywhere below).
29const JULIAN_CENTURY: f64 = 36525.0;
30/// Mean length of the lunar synodic month (new moon to new moon), in days.
31const SYNODIC_MONTH: f64 = 29.53058867;
32
33/// Julian Day for a UTC instant. Meeus formula 7.1, Gregorian calendar only
34/// (fine for any modern date). Constants explained inline below.
35fn julian_day(dt_utc: DateTime<Utc>) -> f64 {
36 let mut y = dt_utc.year();
37 let mut m = dt_utc.month() as i32;
38 let d = dt_utc.day() as f64
39 + (dt_utc.hour() as f64
40 + (dt_utc.minute() as f64 + dt_utc.second() as f64 / 60.0) / 60.0)
41 / 24.0;
42 // Meeus's calendar trick: Jan/Feb are treated as months 13/14 of the
43 // previous year so a single polynomial covers the whole shifted year.
44 if m <= 2 {
45 y -= 1;
46 m += 12;
47 }
48 // `b` is the Gregorian leap-day correction (number of "missing" leap days
49 // in the Gregorian rules vs. the pure Julian 365.25 average).
50 let a = (y as f64 / 100.0).floor() as i32;
51 let b = 2 - a + (a as f64 / 4.0).floor() as i32;
52 // 365.25 = avg days/year (Julian)
53 // 30.6001 = avg days/month over the shifted Mar->Feb calendar
54 // 4716 = year offset that keeps the result non-negative for all dates
55 // -1524.5 = anchors JD 0 to noon on -4712 Jan 1 (Julian proleptic), the
56 // astronomical zero point.
57 (365.25 * (y as f64 + 4716.0)).floor()
58 + (30.6001 * (m as f64 + 1.0)).floor()
59 + d
60 + b as f64
61 - 1524.5
62}
63
64/// (age in days through the synodic month, illumination 0..1).
65/// Truncated form of Meeus chapter 47 with the table 47.A perturbation series.
66fn moon_state(date: DateTime<Tz>) -> (f64, f64) {
67 let dt_utc = date.with_timezone(&Utc);
68 let t = (julian_day(dt_utc) - J2000_JD) / JULIAN_CENTURY;
69 // Meeus 47.2 - 47.5: the four fundamental angles of the lunar problem
70 // (all degrees). Each polynomial is "value at J2000.0 + linear drift per
71 // Julian century"; higher-order terms are below our render precision.
72 // d = mean elongation of moon from sun
73 // ms = sun's mean anomaly
74 // mm = moon's mean anomaly
75 // f = argument of latitude (moon's distance from ascending node)
76 let d = (297.8501921_f64 + 445267.1114034 * t).rem_euclid(360.0);
77 let ms = (357.5291092_f64 + 35999.0502909 * t).rem_euclid(360.0);
78 let mm = (134.9633964_f64 + 477198.8675055 * t).rem_euclid(360.0);
79 let f = (93.2720950_f64 + 483202.0175233 * t).rem_euclid(360.0);
80 let dr = d.to_radians();
81 let msr = ms.to_radians();
82 let mmr = mm.to_radians();
83 let fr = f.to_radians();
84 // Meeus table 47.A: longitude perturbations of the moon (degrees). First
85 // 13 terms are above the arcminute level; the rest of the table is too
86 // small to matter for percent-illumination output. Coefficients come from
87 // the table as-published.
88 let dl_moon = 6.288774 * mmr.sin()
89 + 1.274027 * (2.0 * dr - mmr).sin()
90 + 0.658314 * (2.0 * dr).sin()
91 + 0.213618 * (2.0 * mmr).sin()
92 - 0.185116 * msr.sin()
93 - 0.114332 * (2.0 * fr).sin()
94 + 0.058793 * (2.0 * dr - 2.0 * mmr).sin()
95 + 0.057066 * (2.0 * dr - msr - mmr).sin()
96 + 0.053322 * (2.0 * dr + mmr).sin()
97 + 0.045758 * (2.0 * dr - msr).sin()
98 - 0.040923 * (msr - mmr).sin()
99 - 0.034720 * dr.sin()
100 - 0.030383 * (msr + mmr).sin();
101 // Sun's longitude correction (the equation of centre, kept to 3 terms).
102 let dl_sun = 1.914602 * msr.sin()
103 + 0.019993 * (2.0 * msr).sin()
104 + 0.000289 * (3.0 * msr).sin();
105 // Elongation: how far around the cycle the moon is from the sun (degrees).
106 let elong = (d + dl_moon - dl_sun).rem_euclid(360.0);
107 // Map elongation linearly onto days through the mean synodic month.
108 let age = elong / 360.0 * SYNODIC_MONTH;
109 // Half-angle identity: fraction of disk illuminated, exact for a sphere.
110 let illum = (1.0 - elong.to_radians().cos()) / 2.0;
111 (age, illum)
112}
113
114pub fn moon_phase(date: DateTime<Tz>) -> f64 {
115 moon_state(date).0
116}
117
118pub fn moon_illumination(date: DateTime<Tz>) -> f64 {
119 moon_state(date).1
120}
121
122/// Name for a given lunar age in days. The four "moment" phases (new, first
123/// quarter, full, last quarter) get a thin ~1.85-day window around the exact
124/// moment (0.92 days each side); the four crescent/gibbous phases fill the
125/// rest. Boundaries are placed so the cycle divides into eight pieces with
126/// the moments centred on 0, 7.38 (cycle/4), 14.77 (cycle/2), 22.15 (3*cycle/4).
127pub fn moon_name(phase: f64) -> &'static str {
128 if phase < 1.85 {
129 "new moon"
130 } else if phase < 7.38 {
131 "waxing crescent"
132 } else if phase < 9.23 {
133 "first quarter"
134 } else if phase < 14.77 {
135 "waxing gibbous"
136 } else if phase < 16.61 {
137 "full moon"
138 } else if phase < 22.15 {
139 "waning gibbous"
140 } else if phase < 23.99 {
141 "last quarter"
142 } else if phase < 27.68 {
143 "waning crescent"
144 } else {
145 "new moon"
146 }
147}
148
149// Standard sun-event altitudes, in degrees of the sun's centre below the
150// horizon. Sunrise/sunset uses -0.8333 deg = -50' total, which is the sum of
151// average atmospheric refraction at the horizon (~34') and the sun's apparent
152// semi-diameter (~16'); the upper limb appears tangent to the horizon at this
153// geometric altitude. Twilight thresholds are the usual conventions.
154const SUN_GEOMETRIC_ALT: f64 = -0.8333;
155const CIVIL_DUSK_ALT: f64 = -6.0;
156const NAUTICAL_DUSK_ALT: f64 = -12.0;
157const ASTRO_DUSK_ALT: f64 = -18.0;
158
159#[derive(Clone, Copy, Debug)]
160pub struct SunTimes {
161 pub sunrise: f64,
162 pub sunset: f64,
163 pub civil_dusk: f64,
164 pub nautical_dusk: f64,
165 pub astronomical_dusk: f64,
166 pub daylight_hours: f64,
167}
168
169/// Apparent geocentric position of the sun, plus equation of time.
170/// Returns (right ascension deg, declination rad, equation of time minutes).
171/// Implements Meeus chapter 25 with low-accuracy nutation (chapter 22 / 25.8).
172fn sun_apparent(jde: f64) -> (f64, f64, f64) {
173 // `t` = Julian centuries since J2000.0, the unit every polynomial below
174 // expects. All polynomial coefficients come straight from Meeus.
175 let t = (jde - J2000_JD) / JULIAN_CENTURY;
176 // 25.2: sun's mean longitude L0 (deg). Linear term ~36000 deg/century =
177 // one full cycle/year. Quadratic term is tiny secular drift.
178 let l0 = (280.46646 + 36000.76983 * t + 0.0003032 * t * t).rem_euclid(360.0);
179 // 25.3: sun's mean anomaly M (deg). Same order of magnitude as L0; differs
180 // because mean longitude tracks "where the sun would be on a circle" while
181 // anomaly tracks position relative to perihelion.
182 let m = (357.52911 + 35999.05029 * t - 0.0001537 * t * t).rem_euclid(360.0);
183 let mr = m.to_radians();
184 // Meeus 25, "equation of the centre" series. Three-term truncation gives
185 // sub-arcsecond accuracy. Result is a correction to L0 to get the sun's
186 // true (geometric) longitude.
187 let c = (1.914602 - 0.004817 * t - 0.000014 * t * t) * mr.sin()
188 + (0.019993 - 0.000101 * t) * (2.0 * mr).sin()
189 + 0.000289 * (3.0 * mr).sin();
190 let true_long = l0 + c;
191 // Omega = longitude of ascending node of the moon's mean orbit (deg).
192 // Required by the low-accuracy nutation corrections in 25.8 and 22.
193 let omega_r = (125.04 - 1934.136 * t).to_radians();
194 // 25.8: apparent longitude (corrects for aberration -0.00569 and the dominant
195 // nutation term -0.00478*sin(Omega), both in degrees).
196 let lambda = true_long - 0.00569 - 0.00478 * omega_r.sin();
197 // 22.2: mean obliquity of the ecliptic (deg). Constant term =
198 // 23 deg 26' 21.448" = 23.439291 deg; secular terms ~0.013 deg/century.
199 let eps0 = 23.439291 - 0.0130042 * t - 1.64e-7 * t * t + 5.04e-7 * t * t * t;
200 // Apparent obliquity = mean + dominant nutation term in obliquity.
201 let eps = eps0 + 0.00256 * omega_r.cos();
202 let lr = lambda.to_radians();
203 let er = eps.to_radians();
204 // 25.6 / 25.7: rotate ecliptic (lambda, beta=0) into equatorial coordinates.
205 let alpha = (er.cos() * lr.sin())
206 .atan2(lr.cos())
207 .to_degrees()
208 .rem_euclid(360.0);
209 let decl = (er.sin() * lr.sin()).asin();
210 // 28.1 (simplified): equation of time = L0 - aberration - alpha (deg),
211 // wrapped into (-180, 180], then 4 min per degree gives minutes of time.
212 // The 0.0057183 deg = 20.4965" is the standard aberration constant for
213 // light's annual deflection.
214 let mut diff = (l0 - 0.0057183 - alpha) % 360.0;
215 if diff > 180.0 {
216 diff -= 360.0;
217 }
218 if diff < -180.0 {
219 diff += 360.0;
220 }
221 (alpha, decl, 4.0 * diff)
222}
223
224/// Half-day length, in hours: how long before/after solar noon the sun is at
225/// altitude `alt_rad` for declination `decl`. Returns None if the sun never
226/// reaches that altitude on this day (polar day or polar night). Standard
227/// "sunrise equation": cos(H) = (sin(h) - sin(phi)*sin(delta)) / (cos(phi)*cos(delta)).
228/// Convert from degrees to hours by dividing by 15 (earth turns 15 deg/hr).
229fn hour_angle_at_alt(decl: f64, alt_rad: f64) -> Option<f64> {
230 let lat_r = LAT.to_radians();
231 let cos_h = (alt_rad.sin() - lat_r.sin() * decl.sin()) / (lat_r.cos() * decl.cos());
232 if !(-1.0..=1.0).contains(&cos_h) {
233 return None;
234 }
235 Some(cos_h.acos().to_degrees() / 15.0)
236}
237
238/// Local-tz fractional hour of the rise (`is_rise=true`) or set event when the
239/// sun is at altitude `alt_deg`. Returns `None` if the sun never crosses that
240/// altitude on the given date (polar day/night). Two-pass algorithm: first
241/// solve using the sun's position at noon UTC, then re-solve at the predicted
242/// event time so the answer accounts for declination drift across the day.
243/// One refinement is enough for sub-second accuracy at temperate latitudes.
244fn sun_event_hours(
245 jd_noon: f64,
246 utc_midnight: DateTime<Utc>,
247 tz: Tz,
248 alt_deg: f64,
249 is_rise: bool,
250) -> Option<f64> {
251 let alt_r = alt_deg.to_radians();
252
253 // Pass 1: estimate event using sun position at local noon.
254 let (_, decl0, eot0) = sun_apparent(jd_noon);
255 let h0 = hour_angle_at_alt(decl0, alt_r)?;
256 // Apparent solar noon (UTC, fractional hours):
257 // 12.0 = noon UTC at the prime meridian, ignoring everything else
258 // -LON/15 = longitude correction (15 deg = 1 hr; west is negative LON,
259 // so this term shifts solar noon later in UTC)
260 // -eot/60 = equation of time, converted from minutes to hours; this is
261 // how far apparent (sundial) noon is from mean (clock) noon
262 let solar_noon_0 = 12.0 - LON / 15.0 - eot0 / 60.0;
263 let evt0 = if is_rise { solar_noon_0 - h0 } else { solar_noon_0 + h0 };
264
265 // Pass 2: recompute sun position at the predicted event time. `jd_noon`
266 // is at 12:00 UTC, so `jd_noon - 0.5` is 00:00 UTC the same date; adding
267 // `evt0 / 24` of a day lands at the event itself.
268 let (_, decl1, eot1) = sun_apparent(jd_noon - 0.5 + evt0 / 24.0);
269 let h1 = hour_angle_at_alt(decl1, alt_r)?;
270 let solar_noon_1 = 12.0 - LON / 15.0 - eot1 / 60.0;
271 let evt1 = if is_rise { solar_noon_1 - h1 } else { solar_noon_1 + h1 };
272
273 // Convert UTC fractional hours back to a local-tz fractional hour.
274 let event_utc = utc_midnight + Duration::nanoseconds((evt1 * 3600.0 * 1e9) as i64);
275 let local = event_utc.with_timezone(&tz);
276 Some(local.hour() as f64 + local.minute() as f64 / 60.0 + local.second() as f64 / 3600.0)
277}
278
279/// Sunrise, sunset, and three twilight ends (civil/nautical/astronomical dusk)
280/// for the given local date. Uses Meeus chapter 25 sun position with one
281/// refinement pass; matches USNO to a few seconds. Hours are local fractional
282/// hours within the calendar date of `local_date`.
283pub fn sun_times(local_date: DateTime<Tz>) -> SunTimes {
284 let tz = local_date.timezone();
285 let local_noon = tz
286 .with_ymd_and_hms(
287 local_date.year(),
288 local_date.month(),
289 local_date.day(),
290 12,
291 0,
292 0,
293 )
294 .unwrap();
295 let noon_utc = local_noon.with_timezone(&Utc);
296 let utc_midnight = Utc
297 .with_ymd_and_hms(noon_utc.year(), noon_utc.month(), noon_utc.day(), 0, 0, 0)
298 .unwrap();
299 let jd_noon = julian_day(noon_utc);
300 let evt = |alt, is_rise| {
301 sun_event_hours(jd_noon, utc_midnight, tz, alt, is_rise).unwrap_or(0.0)
302 };
303 let sunrise = evt(SUN_GEOMETRIC_ALT, true);
304 let sunset = evt(SUN_GEOMETRIC_ALT, false);
305 SunTimes {
306 sunrise,
307 sunset,
308 civil_dusk: evt(CIVIL_DUSK_ALT, false),
309 nautical_dusk: evt(NAUTICAL_DUSK_ALT, false),
310 astronomical_dusk: evt(ASTRO_DUSK_ALT, false),
311 daylight_hours: sunset - sunrise,
312 }
313}
314
315pub fn format_hm(hours: f64) -> String {
316 let mut h = hours.trunc() as i64;
317 let mut m = ((hours - h as f64) * 60.0).round() as i64;
318 if m == 60 {
319 h += 1;
320 m = 0;
321 }
322 format!("{h}h {m}m")
323}
324
325pub fn format_clock(hours: f64) -> String {
326 let mut h = hours.trunc() as i64;
327 let mut m = ((hours - h as f64) * 60.0).round() as i64;
328 if m == 60 {
329 h += 1;
330 m = 0;
331 }
332 let suffix = if h >= 12 { "pm" } else { "am" };
333 let display = if h > 12 {
334 h - 12
335 } else if h == 0 {
336 12
337 } else {
338 h
339 };
340 format!("{display}:{m:02} {suffix}")
341}
342
343pub fn sky_data_lines(now: DateTime<Tz>) -> Vec<String> {
344 let phase = moon_phase(now);
345 let name = moon_name(phase);
346 let illum = (moon_illumination(now) * 100.0).round() as i64;
347 let st = sun_times(now);
348 let yesterday = now - Duration::days(1);
349 let gained = (st.daylight_hours - sun_times(yesterday).daylight_hours) * 60.0;
350 let sign = if gained > 0.0 { "+" } else { "" };
351 vec![
352 format!("<strong>{name}</strong>, {illum}% lit"),
353 format!(
354 "sunrise <strong>{}</strong> \u{00b7} sunset <strong>{}</strong>",
355 format_clock(st.sunrise),
356 format_clock(st.sunset)
357 ),
358 format!(
359 "<strong>{}</strong> of daylight ({sign}{:.1} minutes from yesterday)",
360 format_hm(st.daylight_hours),
361 gained
362 ),
363 format!(
364 "civil dusk <strong>{}</strong> \u{00b7} sailor's dark <strong>{}</strong> \u{00b7} true dark <strong>{}</strong>",
365 format_clock(st.civil_dusk),
366 format_clock(st.nautical_dusk),
367 format_clock(st.astronomical_dusk)
368 ),
369 ]
370}
371
372#[cfg(test)]
373mod tests {
374 use super::*;
375 use chrono_tz::America::New_York;
376
377 fn fixed_now() -> DateTime<Tz> {
378 New_York.with_ymd_and_hms(2026, 5, 6, 19, 0, 0).unwrap()
379 }
380
381 #[test]
382 fn moon_matches_python() {
383 let now = fixed_now();
384 let phase = moon_phase(now);
385 let illum = moon_illumination(now);
386 assert!((phase - 19.474334776875626).abs() < 1e-9, "phase={phase}");
387 assert!((illum - 0.7693359060289093).abs() < 1e-9, "illum={illum}");
388 assert_eq!(moon_name(phase), "waning gibbous");
389 }
390
391 #[test]
392 fn sun_matches_usno() {
393 // 2026-05-06 USNO/Naval Observatory values for lat 35.78, lon -78.64 EDT:
394 // sunrise 06:17, sunset 20:06, civil dusk 20:34, nautical 21:07, astro 21:43.
395 // Tolerance is 30 sec since both USNO and Meeus round display values.
396 let now = fixed_now();
397 let st = sun_times(now);
398 let approx = |actual: f64, h: i32, m: i32, label: &str| {
399 let expected = h as f64 + m as f64 / 60.0;
400 let diff = (actual - expected).abs();
401 assert!(diff < 1.0 / 60.0, "{label}: got {actual}, expected ~{h}:{m:02}");
402 };
403 approx(st.sunrise, 6, 17, "sunrise");
404 approx(st.sunset, 20, 6, "sunset");
405 approx(st.civil_dusk, 20, 34, "civil dusk");
406 approx(st.nautical_dusk, 21, 7, "nautical dusk");
407 approx(st.astronomical_dusk, 21, 43, "astro dusk");
408 assert!(
409 (st.daylight_hours - 13.81).abs() < 0.02,
410 "daylight={}",
411 st.daylight_hours
412 );
413 }
414}