heartwood every commit a ring

the sun keeps truer hours and dusk announces itself three times

0cbc0f7f by Isaac Bythewood · 2 days ago

the sun keeps truer hours and dusk announces itself three times

replace spencer fourier with meeus chapter 25 sun position; sunset
no longer drifts ~3 minutes from usno. add civil dusk, sailor's dark,
and true dark to the sky lines.
modified CLAUDE.md
@@ -62,9 +62,11 @@ strips a single wrapping `<p>...</p>` for inline contexts (matches Mistune'soutput for the original flask version); `render_block` keeps block tags.**Astronomical math:** `src/astro.rs`. Moon phase + illumination via Meeus'slunar series (table 47.A perturbations). Sunrise/sunset/daylight via NOAA'sSpencer fourier series. Locked to zone 7a (lat 35.78, lon -78.64). Local-tzhandling uses chrono-tz with America/New_York. Accurate to ~1 minute.lunar series (table 47.A perturbations). Sun position via Meeus chapter 25with one refinement pass at predicted event time; produces sunrise, sunset,civil/nautical/astronomical dusk (sun at -6/-12/-18°). Matches USNO to a fewseconds. Locked to zone 7a (lat 35.78, lon -78.64). Local-tz handling useschrono-tz with America/New_York.**Daily-stable seeded RNG:** `src/rng.rs` is a mulberry32 PRNG withJS-Math.imul-compatible semantics. Keyed by day-of-year so picks shift day
modified src/astro.rs
@@ -89,54 +89,123 @@ pub fn moon_name(phase: f64) -> &'static str {    }}fn is_leap(year: i32) -> bool {    (year % 4 == 0 && year % 100 != 0) || year % 400 == 0}// standard sun-event altitudes (degrees below horizon).// -0.8333 = atmospheric refraction (-34') + sun semi-diameter (-16').const SUN_GEOMETRIC_ALT: f64 = -0.8333;const CIVIL_DUSK_ALT: f64 = -6.0;const NAUTICAL_DUSK_ALT: f64 = -12.0;const ASTRO_DUSK_ALT: f64 = -18.0;/// (sunrise local hours, sunset local hours, day length hours)./// uses NOAA Spencer fourier series; accounts for longitude, equation of/// time, atmospheric refraction; converts to local-tz of `local_date`.pub fn sun_times(local_date: DateTime<Tz>) -> (f64, f64, f64) {    let year = local_date.year();    let doy = local_date.ordinal() as i32;    let year_len = if is_leap(year) { 366.0 } else { 365.0 };    let gamma = (2.0 * std::f64::consts::PI / year_len) * (doy as f64 - 1.0);    let eot = 229.18        * (0.000075 + 0.001868 * gamma.cos()            - 0.032077 * gamma.sin()            - 0.014615 * (2.0 * gamma).cos()            - 0.040849 * (2.0 * gamma).sin());    let decl = 0.006918 - 0.399912 * gamma.cos()        + 0.070257 * gamma.sin()        - 0.006758 * (2.0 * gamma).cos()        + 0.000907 * (2.0 * gamma).sin()        - 0.002697 * (3.0 * gamma).cos()        + 0.001480 * (3.0 * gamma).sin();    let lat_rad = LAT.to_radians();    let mut cos_ha = (90.833_f64.to_radians().cos() - lat_rad.sin() * decl.sin())        / (lat_rad.cos() * decl.cos());    cos_ha = cos_ha.clamp(-1.0, 1.0);    let ha = cos_ha.acos().to_degrees();    let solar_noon_min = 720.0 - 4.0 * LON - eot;    let sr_min = solar_noon_min - 4.0 * ha;    let ss_min = solar_noon_min + 4.0 * ha;#[derive(Clone, Copy, Debug)]pub struct SunTimes {    pub sunrise: f64,    pub sunset: f64,    pub civil_dusk: f64,    pub nautical_dusk: f64,    pub astronomical_dusk: f64,    pub daylight_hours: f64,}    let tz = local_date.timezone();    let base = Utc        .with_ymd_and_hms(year, local_date.month(), local_date.day(), 0, 0, 0)        .unwrap();/// (right ascension deg, declination rad, equation of time min). Meeus ch 25.fn sun_apparent(jde: f64) -> (f64, f64, f64) {    let t = (jde - 2451545.0) / 36525.0;    let l0 = (280.46646 + 36000.76983 * t + 0.0003032 * t * t).rem_euclid(360.0);    let m = (357.52911 + 35999.05029 * t - 0.0001537 * t * t).rem_euclid(360.0);    let mr = m.to_radians();    let c = (1.914602 - 0.004817 * t - 0.000014 * t * t) * mr.sin()        + (0.019993 - 0.000101 * t) * (2.0 * mr).sin()        + 0.000289 * (3.0 * mr).sin();    let true_long = l0 + c;    let omega_r = (125.04 - 1934.136 * t).to_radians();    let lambda = true_long - 0.00569 - 0.00478 * omega_r.sin();    let eps0 = 23.439291 - 0.0130042 * t - 1.64e-7 * t * t + 5.04e-7 * t * t * t;    let eps = eps0 + 0.00256 * omega_r.cos();    let lr = lambda.to_radians();    let er = eps.to_radians();    let alpha = (er.cos() * lr.sin())        .atan2(lr.cos())        .to_degrees()        .rem_euclid(360.0);    let decl = (er.sin() * lr.sin()).asin();    let mut diff = (l0 - 0.0057183 - alpha) % 360.0;    if diff > 180.0 {        diff -= 360.0;    }    if diff < -180.0 {        diff += 360.0;    }    (alpha, decl, 4.0 * diff)}    let sr = (base + Duration::nanoseconds((sr_min * 60.0 * 1e9) as i64)).with_timezone(&tz);    let ss = (base + Duration::nanoseconds((ss_min * 60.0 * 1e9) as i64)).with_timezone(&tz);fn hour_angle_at_alt(decl: f64, alt_rad: f64) -> Option<f64> {    let lat_r = LAT.to_radians();    let cos_h = (alt_rad.sin() - lat_r.sin() * decl.sin()) / (lat_r.cos() * decl.cos());    if !(-1.0..=1.0).contains(&cos_h) {        return None;    }    Some(cos_h.acos().to_degrees() / 15.0)}    let sr_h = sr.hour() as f64 + sr.minute() as f64 / 60.0 + sr.second() as f64 / 3600.0;    let ss_h = ss.hour() as f64 + ss.minute() as f64 / 60.0 + ss.second() as f64 / 3600.0;/// Local-tz fractional hour of the rise (`is_rise=true`) or set event when the/// sun is at altitude `alt_deg`. Returns `None` if the sun never crosses that/// altitude on the given date (no event, e.g. polar day/night).fn sun_event_hours(    jd_noon: f64,    utc_midnight: DateTime<Utc>,    tz: Tz,    alt_deg: f64,    is_rise: bool,) -> Option<f64> {    let alt_r = alt_deg.to_radians();    // first pass uses noon-UTC sun position    let (_, decl0, eot0) = sun_apparent(jd_noon);    let h0 = hour_angle_at_alt(decl0, alt_r)?;    let solar_noon_0 = 12.0 - LON / 15.0 - eot0 / 60.0;    let evt0 = if is_rise { solar_noon_0 - h0 } else { solar_noon_0 + h0 };    // refine once at the predicted event time    let (_, decl1, eot1) = sun_apparent(jd_noon - 0.5 + evt0 / 24.0);    let h1 = hour_angle_at_alt(decl1, alt_r)?;    let solar_noon_1 = 12.0 - LON / 15.0 - eot1 / 60.0;    let evt1 = if is_rise { solar_noon_1 - h1 } else { solar_noon_1 + h1 };    let event_utc = utc_midnight + Duration::nanoseconds((evt1 * 3600.0 * 1e9) as i64);    let local = event_utc.with_timezone(&tz);    Some(local.hour() as f64 + local.minute() as f64 / 60.0 + local.second() as f64 / 3600.0)}    (sr_h, ss_h, (ss_min - sr_min) / 60.0)/// Sunrise, sunset, and three twilight ends (civil/nautical/astronomical dusk)/// for the given local date. Uses Meeus chapter 25 sun position with one/// refinement pass; matches USNO to a few seconds. Hours are local fractional/// hours within the calendar date of `local_date`.pub fn sun_times(local_date: DateTime<Tz>) -> SunTimes {    let tz = local_date.timezone();    let local_noon = tz        .with_ymd_and_hms(            local_date.year(),            local_date.month(),            local_date.day(),            12,            0,            0,        )        .unwrap();    let noon_utc = local_noon.with_timezone(&Utc);    let utc_midnight = Utc        .with_ymd_and_hms(noon_utc.year(), noon_utc.month(), noon_utc.day(), 0, 0, 0)        .unwrap();    let jd_noon = julian_day(noon_utc);    let evt = |alt, is_rise| {        sun_event_hours(jd_noon, utc_midnight, tz, alt, is_rise).unwrap_or(0.0)    };    let sunrise = evt(SUN_GEOMETRIC_ALT, true);    let sunset = evt(SUN_GEOMETRIC_ALT, false);    SunTimes {        sunrise,        sunset,        civil_dusk: evt(CIVIL_DUSK_ALT, false),        nautical_dusk: evt(NAUTICAL_DUSK_ALT, false),        astronomical_dusk: evt(ASTRO_DUSK_ALT, false),        daylight_hours: sunset - sunrise,    }}pub fn format_hm(hours: f64) -> String {
@@ -167,22 +236,28 @@ pub fn sky_data_lines(now: DateTime<Tz>) -> Vec<String> {    let phase = moon_phase(now);    let name = moon_name(phase);    let illum = (moon_illumination(now) * 100.0).round() as i64;    let (sunrise, sunset, hours) = sun_times(now);    let st = sun_times(now);    let yesterday = now - Duration::days(1);    let gained = (hours - sun_times(yesterday).2) * 60.0;    let gained = (st.daylight_hours - sun_times(yesterday).daylight_hours) * 60.0;    let sign = if gained > 0.0 { "+" } else { "" };    vec![        format!("<strong>{name}</strong>, {illum}% lit"),        format!(            "sunrise <strong>{}</strong> \u{00b7} sunset <strong>{}</strong>",            format_clock(sunrise),            format_clock(sunset)            format_clock(st.sunrise),            format_clock(st.sunset)        ),        format!(            "<strong>{}</strong> of daylight ({sign}{:.1} minutes from yesterday)",            format_hm(hours),            format_hm(st.daylight_hours),            gained        ),        format!(            "civil dusk <strong>{}</strong> \u{00b7} sailor's dark <strong>{}</strong> \u{00b7} true dark <strong>{}</strong>",            format_clock(st.civil_dusk),            format_clock(st.nautical_dusk),            format_clock(st.astronomical_dusk)        ),    ]}
@@ -206,11 +281,26 @@ mod tests {    }    #[test]    fn sun_matches_python() {    fn sun_matches_usno() {        // 2026-05-06 USNO/Naval Observatory values for lat 35.78, lon -78.64 EDT:        //   sunrise 06:17, sunset 20:06, civil dusk 20:34, nautical 21:07, astro 21:43.        // Tolerance is 30 sec since both USNO and Meeus round display values.        let now = fixed_now();        let (sr, ss, hours) = sun_times(now);        assert!((sr - 6.298611111111111).abs() < 1e-6, "sr={sr}");        assert!((ss - 20.067777777777778).abs() < 1e-6, "ss={ss}");        assert!((hours - 13.768966853824192).abs() < 1e-9, "hours={hours}");        let st = sun_times(now);        let approx = |actual: f64, h: i32, m: i32, label: &str| {            let expected = h as f64 + m as f64 / 60.0;            let diff = (actual - expected).abs();            assert!(diff < 1.0 / 60.0, "{label}: got {actual}, expected ~{h}:{m:02}");        };        approx(st.sunrise, 6, 17, "sunrise");        approx(st.sunset, 20, 6, "sunset");        approx(st.civil_dusk, 20, 34, "civil dusk");        approx(st.nautical_dusk, 21, 7, "nautical dusk");        approx(st.astronomical_dusk, 21, 43, "astro dusk");        assert!(            (st.daylight_hours - 13.81).abs() < 0.02,            "daylight={}",            st.daylight_hours        );    }}