Moon phases

use moon_phase.qweruse moon_phase, moon_phases from moon_phase      // Quillo 0.6: only these

chrono = { version = "0.4", default-features = false, features = ["std", "clock"] }

// A fishing and night-sky calendar: tonight's moon, the month's phases, and when the moon is up.////   quillo run moon_phase.qwe      then open http://localhost:8080/moonuse "moon_phase.qwer"route GET /moon ? { at: str = "2026-08-28T04:00:00Z" } = moon_phase(at)?// August 2026: the total solar eclipse of 12 August is at a new moon, the partial lunar eclipse of 28 August at a full one.route GET /month = moon_phases("2026-08-01", "2026-08-31")?// Moonrise and moonset in Amsterdam (UTC+2 in summer).route GET /tonight ? { date: str = "2026-08-28" } = moon_times(52.37, 4.89, date, 120)?route GET / {    let m = moon_phase("now")?    <html><body>        <h1>{m.emoji} {m.phase}</h1>        <p>{round(m.illumination * 100)}% lit, {m.age_days} days old. Next full moon: {m.next_full}.</p>    </body></html>}

//! The moon for a date and a place: its phase and how much of it is lit, its age, the next new, first quarter, full and last quarter moons, the phases in a range, and moonrise and moonset. For tide, fishing, gardening and night-sky apps.//! version: 1.0.0//!//! [dependencies]//! chrono = { version = "0.4", default-features = false, features = ["std", "clock"] }use chrono::{DateTime, NaiveDate, NaiveDateTime, Utc};use serde_json::{json, Value};const J2000: f64 = 2451545.0;const SYNODIC: f64 = 29.530588861;const DELTA_T_DAYS: f64 = 69.0 / 86400.0;fn sin(d: f64) -> f64 { d.to_radians().sin() }fn cos(d: f64) -> f64 { d.to_radians().cos() }fn norm(d: f64) -> f64 { d.rem_euclid(360.0) }fn jd_of(t: DateTime<Utc>) -> f64 {    t.timestamp_millis() as f64 / 86_400_000.0 + 2440587.5}fn time_of(jd: f64) -> DateTime<Utc> {    DateTime::from_timestamp_millis(((jd - 2440587.5) * 86_400_000.0).round() as i64).unwrap_or_default()}fn parse_time(s: &str) -> Result<DateTime<Utc>, String> {    let t = s.trim();    if t.is_empty() || t == "now" { return Ok(Utc::now()) }    if let Ok(d) = DateTime::parse_from_rfc3339(t) { return Ok(d.with_timezone(&Utc)) }    for f in ["%Y-%m-%dT%H:%M:%S", "%Y-%m-%dT%H:%M", "%Y-%m-%d %H:%M"] {        if let Ok(d) = NaiveDateTime::parse_from_str(t, f) { return Ok(d.and_utc()) }    }    if let Ok(d) = NaiveDate::parse_from_str(t, "%Y-%m-%d") { return Ok(d.and_hms_opt(12, 0, 0).unwrap().and_utc()) }    Err(format!("`{}` is not a date or time (2026-10-29, 2026-10-29T21:00:00Z)", s))}fn iso(t: DateTime<Utc>, offset_min: i64) -> String {    let off = chrono::FixedOffset::east_opt((offset_min * 60) as i32).unwrap_or(chrono::FixedOffset::east_opt(0).unwrap());    t.with_timezone(&off).format("%Y-%m-%dT%H:%M%:z").to_string()}// Low-precision positions (Astronomical Almanac): the moon to about 0.3°, the sun to 0.01°.fn moon(jd: f64) -> (f64, f64, f64) {    let t = (jd - J2000) / 36525.0;    let l = 218.32 + 481267.881 * t + 6.29 * sin(135.0 + 477198.87 * t) - 1.27 * sin(259.3 - 413335.36 * t)        + 0.66 * sin(235.7 + 890534.22 * t) + 0.21 * sin(269.9 + 954397.74 * t) - 0.19 * sin(357.5 + 35999.05 * t)        - 0.11 * sin(186.5 + 966404.03 * t);    let b = 5.13 * sin(93.3 + 483202.02 * t) + 0.28 * sin(228.2 + 960400.89 * t) - 0.28 * sin(318.3 + 6003.15 * t) - 0.17 * sin(217.6 - 407332.21 * t);    let p = 0.9508 + 0.0518 * cos(135.0 + 477198.87 * t) + 0.0095 * cos(259.3 - 413335.36 * t) + 0.0078 * cos(235.7 + 890534.22 * t) + 0.0028 * cos(269.9 + 954397.74 * t);    (norm(l), b, p)}fn sun(jd: f64) -> f64 {    let n = jd - J2000;    let g = 357.528 + 0.9856003 * n;    norm(280.460 + 0.9856474 * n + 1.915 * sin(g) + 0.020 * sin(2.0 * g))}fn elongation(jd: f64) -> f64 {    norm(moon(jd).0 - sun(jd))}// Meeus, Astronomical Algorithms ch. 49: the instant of a phase (0 new, 1 first quarter, 2 full, 3 last quarter) for lunation k.fn phase_jde(k: f64, q: u8) -> f64 {    let k = k + q as f64 * 0.25;    let t = k / 1236.85;    let mut jde = 2451550.09766 + SYNODIC * k + 0.00015437 * t * t - 0.00000015 * t * t * t + 0.00000000073 * t.powi(4);    let e = 1.0 - 0.002516 * t - 0.0000074 * t * t;    let m = 2.5534 + 29.10535670 * k - 0.0000014 * t * t;    let mp = 201.5643 + 385.81693528 * k + 0.0107582 * t * t;    let f = 160.7108 + 390.67050284 * k - 0.0016118 * t * t;    let om = 124.7746 - 1.56375588 * k + 0.0020672 * t * t;    let c = match q {        0 | 2 => {            let (a, b) = if q == 0 { (-0.40720, 0.17241) } else { (-0.40614, 0.17302) };            a * sin(mp) + b * e * sin(m) + 0.01608 * sin(2.0 * mp) + 0.01039 * sin(2.0 * f) + 0.00739 * e * sin(mp - m)                - 0.00514 * e * sin(mp + m) + 0.00208 * e * e * sin(2.0 * m) - 0.00111 * sin(mp - 2.0 * f) - 0.00057 * sin(mp + 2.0 * f)                + 0.00056 * e * sin(2.0 * mp + m) - 0.00042 * sin(3.0 * mp) + 0.00042 * e * sin(m + 2.0 * f) + 0.00038 * e * sin(m - 2.0 * f)                - 0.00024 * e * sin(2.0 * mp - m) - 0.00017 * sin(om)        }        _ => {            let w = 0.00306 - 0.00038 * e * cos(m) + 0.00026 * cos(mp) - 0.00002 * cos(mp - m) + 0.00002 * cos(mp + m) + 0.00002 * cos(2.0 * f);            -0.62801 * sin(mp) + 0.17172 * e * sin(m) - 0.01183 * e * sin(mp + m) + 0.00862 * sin(2.0 * mp) + 0.00804 * sin(2.0 * f)                + 0.00454 * e * sin(mp - m) + 0.00204 * e * e * sin(2.0 * m) - 0.00180 * sin(mp - 2.0 * f) - 0.00070 * sin(mp + 2.0 * f)                - 0.00040 * sin(3.0 * mp) - 0.00034 * e * sin(2.0 * mp - m) + 0.00032 * e * sin(m + 2.0 * f) + 0.00032 * e * sin(m - 2.0 * f)                - 0.00028 * e * e * sin(mp + 2.0 * m) + 0.00027 * e * sin(2.0 * mp + m) - 0.00017 * sin(om)                + if q == 1 { w } else { -w }        }    };    jde += c;    jde - DELTA_T_DAYS}const NAMES: [&str; 4] = ["new moon", "first quarter", "full moon", "last quarter"];// Every principal phase from jd0 to jd1, in order.fn phases_between(jd0: f64, jd1: f64) -> Vec<(f64, u8)> {    let mut k = ((jd0 - 2451550.09766) / SYNODIC).floor() - 1.0;    let mut out = Vec::new();    loop {        for q in 0..4u8 {            let j = phase_jde(k, q);            if j > jd1 { return out }            if j >= jd0 { out.push((j, q)) }        }        k += 1.0;    }}fn phase_name(age: f64) -> (&'static str, &'static str) {    let q = SYNODIC / 4.0;    let near = 1.0;    if age < near || age > SYNODIC - near { ("new moon", "🌑") }    else if age < q - near { ("waxing crescent", "🌒") }    else if age < q + near { ("first quarter", "🌓") }    else if age < 2.0 * q - near { ("waxing gibbous", "🌔") }    else if age < 2.0 * q + near { ("full moon", "🌕") }    else if age < 3.0 * q - near { ("waning gibbous", "🌖") }    else if age < 3.0 * q + near { ("last quarter", "🌗") }    else { ("waning crescent", "🌘") }}/// The moon at a moment (`"2026-10-29T21:00:00Z"`, a date for its noon UTC, or `"now"`): `{ phase, emoji, illumination, age_days, waxing, next_new, next_first_quarter, next_full, next_last_quarter }`./// `illumination` is the lit fraction from 0 to 1; `age_days` the days since the last new moon; the next phases in UTC, to about a minute.pub fn moon_phase(at: &str) -> Result<Value, String> {    let t = parse_time(at)?;    let jd = jd_of(t);    let (l, b, _) = moon(jd);    let psi = (cos(b) * cos(l - sun(jd))).clamp(-1.0, 1.0).acos().to_degrees();    let illum = (1.0 - cos(psi)) / 2.0;    let last_new = phases_between(jd - 31.0, jd).into_iter().filter(|p| p.1 == 0).last().map(|p| p.0).unwrap_or(jd);    let age = jd - last_new;    let (name, emoji) = phase_name(age);    let next: Vec<(f64, u8)> = phases_between(jd, jd + 31.0);    let next_of = |q: u8| next.iter().find(|p| p.1 == q).map(|p| time_of(p.0).format("%Y-%m-%dT%H:%MZ").to_string());    Ok(json!({        "phase": name,        "emoji": emoji,        "illumination": (illum * 1000.0).round() / 1000.0,        "age_days": (age * 100.0).round() / 100.0,        "waxing": elongation(jd) < 180.0,        "next_new": next_of(0),        "next_first_quarter": next_of(1),        "next_full": next_of(2),        "next_last_quarter": next_of(3),    }))}/// The principal phases (new moon, first quarter, full moon, last quarter) from `from` to `to` (dates, at most 3 years apart): a list of `{ phase, time }`, times in UTC.pub fn moon_phases(from: &str, to: &str) -> Result<Vec<Value>, String> {    let (a, b) = (jd_of(parse_time(from)?) - 0.5, jd_of(parse_time(to)?) + 0.5);    if b < a { return Err("`to` is before `from`".into()) }    if b - a > 3.0 * 366.0 { return Err("at most 3 years at a time".into()) }    Ok(phases_between(a, b).into_iter().map(|(j, q)| json!({ "phase": NAMES[q as usize], "time": time_of(j).format("%Y-%m-%dT%H:%MZ").to_string() })).collect())}fn altitude(jd: f64, lat: f64, lon: f64) -> f64 {    let (l, b, p) = moon(jd);    let t = (jd - J2000) / 36525.0;    let eps = 23.439291 - 0.0130042 * t;    let ra = (sin(l) * cos(eps) - (b.to_radians().tan()) * sin(eps)).atan2(cos(l)).to_degrees();    let dec = (sin(b) * cos(eps) + cos(b) * sin(eps) * sin(l)).clamp(-1.0, 1.0).asin().to_degrees();    let gmst = norm(280.46061837 + 360.98564736629 * (jd - J2000));    let h = gmst + lon - ra;    let alt = (sin(lat) * sin(dec) + cos(lat) * cos(dec) * cos(h)).clamp(-1.0, 1.0).asin().to_degrees();    // Risen when the top of the disc clears the horizon, refraction and parallax counted.    alt - (0.7275 * p - 0.5667)}/// Moonrise and moonset at a place on `date` (`2026-10-29`), as local times at `utc_offset_minutes` (60 for UTC+1): `{ rise, set, always_up, always_down }`. A day can have no rise or no set (`null`); near the poles the moon can stay up or down all day.pub fn moon_times(lat: f64, lon: f64, date: &str, utc_offset_minutes: i64) -> Result<Value, String> {    if !(-90.0..=90.0).contains(&lat) || !(-180.0..=180.0).contains(&lon) {        return Err(format!("({}, {}) is not a latitude and longitude", lat, lon));    }    let d = NaiveDate::parse_from_str(date.trim(), "%Y-%m-%d").map_err(|_| format!("`{}` is not a date (2026-10-29)", date))?;    let start = jd_of(d.and_hms_opt(0, 0, 0).unwrap().and_utc()) - utc_offset_minutes as f64 / 1440.0;    let step = 10.0 / 1440.0;    let (mut rise, mut set) = (None, None);    let mut prev = altitude(start, lat, lon);    let (mut up, mut down) = (prev > 0.0, prev <= 0.0);    for i in 1..=144 {        let j1 = start + i as f64 * step;        let a1 = altitude(j1, lat, lon);        if a1 > 0.0 { down = false } else { up = false }        if (prev <= 0.0) != (a1 <= 0.0) {            let (mut lo, mut hi) = (j1 - step, j1);            for _ in 0..20 {                let mid = (lo + hi) / 2.0;                if (altitude(mid, lat, lon) <= 0.0) == (prev <= 0.0) { lo = mid } else { hi = mid }            }            let when = iso(time_of((lo + hi) / 2.0), utc_offset_minutes);            if prev <= 0.0 { if rise.is_none() { rise = Some(when) } } else if set.is_none() { set = Some(when) }        }        prev = a1;    }    Ok(json!({ "rise": rise, "set": set, "always_up": up, "always_down": down }))}