Initial Commit

backend/src/astronomy/coords.rs

@@ -0,0 +1,44 @@
+/// Convert RA/Dec to Altitude/Azimuth.
+/// All inputs and outputs in degrees.
+pub fn radec_to_altaz(
+    ra_deg: f64,
+    dec_deg: f64,
+    lst_deg: f64,
+    lat_deg: f64,
+) -> (f64, f64) {
+    let ha = (lst_deg - ra_deg).rem_euclid(360.0);
+    let ha_rad = ha.to_radians();
+    let dec_rad = dec_deg.to_radians();
+    let lat_rad = lat_deg.to_radians();
+
+    let sin_alt = dec_rad.sin() * lat_rad.sin()
+        + dec_rad.cos() * lat_rad.cos() * ha_rad.cos();
+    let alt_rad = sin_alt.asin();
+
+    let cos_az = (dec_rad.sin() - lat_rad.sin() * sin_alt)
+        / (lat_rad.cos() * alt_rad.cos());
+    let cos_az = cos_az.clamp(-1.0, 1.0);
+    let az_rad = cos_az.acos();
+
+    let az_deg = if ha_rad.sin() < 0.0 {
+        az_rad.to_degrees()
+    } else {
+        360.0 - az_rad.to_degrees()
+    };
+
+    (alt_rad.to_degrees(), az_deg)
+}
+
+/// Rozenberg airmass formula — valid to horizon.
+pub fn airmass(alt_deg: f64) -> f64 {
+    if alt_deg <= 0.0 {
+        return 40.0; // clamp at horizon
+    }
+    let z_rad = (90.0 - alt_deg).to_radians();
+    1.0 / (z_rad.cos() + 0.025 * (-11.0 * z_rad.cos()).exp())
+}
+
+/// Extinction in magnitudes. k = 0.20 mag/airmass (Bortle 5 site).
+pub fn extinction_mag(alt_deg: f64) -> f64 {
+    airmass(alt_deg) * 0.20
+}

backend/src/astronomy/horizon.rs

@@ -0,0 +1,33 @@
+use serde::{Deserialize, Serialize};
+
+#[derive(Debug, Clone, Serialize, Deserialize, sqlx::FromRow)]
+pub struct HorizonPoint {
+    pub az_deg: i32,
+    pub alt_deg: f64,
+}
+
+/// Linear interpolation of horizon altitude at a given azimuth.
+/// The profile must have points at every integer degree 0–359.
+pub fn horizon_alt(az_deg: f64, profile: &[HorizonPoint]) -> f64 {
+    if profile.is_empty() {
+        return 15.0;
+    }
+
+    let az = az_deg.rem_euclid(360.0);
+    let lo_idx = az.floor() as usize % 360;
+    let hi_idx = (lo_idx + 1) % 360;
+    let frac = az.fract();
+
+    let lo_alt = profile
+        .iter()
+        .find(|p| p.az_deg == lo_idx as i32)
+        .map(|p| p.alt_deg)
+        .unwrap_or(15.0);
+    let hi_alt = profile
+        .iter()
+        .find(|p| p.az_deg == hi_idx as i32)
+        .map(|p| p.alt_deg)
+        .unwrap_or(15.0);
+
+    lo_alt + frac * (hi_alt - lo_alt)
+}

backend/src/astronomy/lunar.rs

@@ -0,0 +1,136 @@
+use chrono::{DateTime, Duration, Utc};
+
+use super::{coords::radec_to_altaz, time::{julian_date, local_sidereal_time}};
+
+/// Compute approximate Moon RA/Dec (degrees) for a given JD.
+/// Low-precision algorithm (< 1° error).
+pub fn moon_position(jd: f64) -> (f64, f64) {
+    let d = jd - 2451545.0;
+
+    // Orbital elements
+    let l = (218.316 + 13.176396 * d).rem_euclid(360.0); // ecliptic longitude
+    let m = (134.963 + 13.064993 * d).to_radians().rem_euclid(std::f64::consts::TAU); // mean anomaly
+    let f = (93.272 + 13.229350 * d).to_radians().rem_euclid(std::f64::consts::TAU); // argument of latitude
+
+    let lambda = l + 6.289 * m.sin(); // ecliptic longitude corrected
+    let beta = 5.128 * f.sin(); // ecliptic latitude
+
+    let lambda_rad = lambda.to_radians();
+    let beta_rad = beta.to_radians();
+    let epsilon = (23.439 - 0.0000004 * d).to_radians();
+
+    let ra = (lambda_rad.sin() * epsilon.cos() - beta_rad.tan() * epsilon.sin())
+        .atan2(lambda_rad.cos());
+    let ra_deg = ra.to_degrees().rem_euclid(360.0);
+    let dec_deg = (beta_rad.sin() * epsilon.cos()
+        + beta_rad.cos() * epsilon.sin() * lambda_rad.sin())
+    .asin()
+    .to_degrees();
+
+    (ra_deg, dec_deg)
+}
+
+/// Moon illumination as fraction 0.0–1.0.
+pub fn moon_illumination(jd: f64) -> f64 {
+    // Sun-Moon elongation
+    let n = jd - 2451545.0;
+    let sun_l = (280.460 + 0.9856474 * n).rem_euclid(360.0);
+    let sun_g = (357.528 + 0.9856003 * n).to_radians();
+    let sun_lambda = sun_l + 1.915 * sun_g.sin() + 0.020 * (2.0 * sun_g).sin();
+
+    let moon_l = (218.316 + 13.176396 * n).rem_euclid(360.0);
+    let moon_m = (134.963 + 13.064993 * n).to_radians();
+    let moon_lambda = moon_l + 6.289 * moon_m.sin();
+
+    let i = (moon_lambda - sun_lambda).rem_euclid(360.0);
+    let k = (1.0 - i.to_radians().cos()) / 2.0;
+    k
+}
+
+/// Moon age in days from last new moon (approximate).
+pub fn moon_age_days(jd: f64) -> f64 {
+    let synodic = 29.53058868;
+    let new_moon_jd = 2451550.1; // reference new moon: 2000-01-06
+    ((jd - new_moon_jd) % synodic + synodic) % synodic
+}
+
+/// Phase name from illumination and age.
+pub fn moon_phase_name(illumination: f64, age_days: f64) -> String {
+    let pct = illumination * 100.0;
+    if age_days < 1.0 {
+        "New Moon".to_string()
+    } else if age_days < 7.4 {
+        format!("Waxing Crescent ({:.0}%)", pct)
+    } else if age_days < 8.4 {
+        "First Quarter".to_string()
+    } else if age_days < 13.7 {
+        format!("Waxing Gibbous ({:.0}%)", pct)
+    } else if age_days < 15.3 {
+        "Full Moon".to_string()
+    } else if age_days < 22.1 {
+        format!("Waning Gibbous ({:.0}%)", pct)
+    } else if age_days < 23.1 {
+        "Last Quarter".to_string()
+    } else if age_days < 29.0 {
+        format!("Waning Crescent ({:.0}%)", pct)
+    } else {
+        "New Moon".to_string()
+    }
+}
+
+/// Moon altitude at a given time for observer position.
+pub fn moon_altitude(jd: f64, lat_deg: f64, lon_deg: f64) -> f64 {
+    let (ra, dec) = moon_position(jd);
+    let lst = local_sidereal_time(jd, lon_deg);
+    let (alt, _) = radec_to_altaz(ra, dec, lst, lat_deg);
+    alt
+}
+
+/// Find moon rise and set times within the night window.
+/// Steps through at 5-minute intervals, interpolates crossings.
+pub fn moon_rise_set(
+    dusk: DateTime<Utc>,
+    dawn: DateTime<Utc>,
+    lat: f64,
+    lon: f64,
+) -> (Option<DateTime<Utc>>, Option<DateTime<Utc>>) {
+    let step = Duration::minutes(5);
+    let mut rise: Option<DateTime<Utc>> = None;
+    let mut set: Option<DateTime<Utc>> = None;
+
+    let mut t = dusk;
+    let mut prev_alt = moon_altitude(julian_date(t), lat, lon);
+
+    while t < dawn {
+        let next = t + step;
+        let next_alt = moon_altitude(julian_date(next), lat, lon);
+
+        if prev_alt < 0.0 && next_alt >= 0.0 && rise.is_none() {
+            // Rising: interpolate
+            let frac = (-prev_alt) / (next_alt - prev_alt);
+            let crossing = t + Duration::seconds((frac * step.num_seconds() as f64) as i64);
+            rise = Some(crossing);
+        } else if prev_alt >= 0.0 && next_alt < 0.0 && set.is_none() {
+            // Setting: interpolate
+            let frac = prev_alt / (prev_alt - next_alt);
+            let crossing = t + Duration::seconds((frac * step.num_seconds() as f64) as i64);
+            set = Some(crossing);
+        }
+
+        prev_alt = next_alt;
+        t = next;
+    }
+
+    (rise, set)
+}
+
+/// Separation in degrees between Moon and a target (RA/Dec in degrees).
+pub fn moon_separation(moon_ra: f64, moon_dec: f64, target_ra: f64, target_dec: f64) -> f64 {
+    let ra1 = moon_ra.to_radians();
+    let dec1 = moon_dec.to_radians();
+    let ra2 = target_ra.to_radians();
+    let dec2 = target_dec.to_radians();
+
+    let cos_sep = dec1.sin() * dec2.sin() + dec1.cos() * dec2.cos() * (ra1 - ra2).cos();
+    cos_sep.clamp(-1.0, 1.0).acos().to_degrees()
+}

backend/src/astronomy/mod.rs

@@ -0,0 +1,13 @@
+pub mod coords;
+pub mod horizon;
+pub mod lunar;
+pub mod solar;
+pub mod time;
+pub mod visibility;
+
+pub use coords::{airmass, extinction_mag, radec_to_altaz};
+pub use horizon::{horizon_alt, HorizonPoint};
+pub use lunar::{moon_age_days, moon_altitude, moon_illumination, moon_phase_name, moon_position, moon_rise_set, moon_separation};
+pub use solar::astro_twilight;
+pub use time::julian_date;
+pub use visibility::{compute_visibility, compute_visibility_with_step, true_dark_window, MoonState, TonightWindow};

backend/src/astronomy/solar.rs

@@ -0,0 +1,88 @@
+use chrono::{DateTime, Duration, NaiveDate, TimeZone, Utc};
+
+use super::{coords::radec_to_altaz, time::{julian_date, local_sidereal_time}};
+
+/// Compute approximate Sun RA/Dec (degrees) for a given JD.
+/// Uses low-precision VSOP87 approximation (< 1° error).
+fn sun_radec(jd: f64) -> (f64, f64) {
+    let n = jd - 2451545.0;
+    let l = (280.460 + 0.9856474 * n).rem_euclid(360.0); // mean longitude
+    let g = (357.528 + 0.9856003 * n).to_radians().rem_euclid(std::f64::consts::TAU); // mean anomaly
+    let lambda = l + 1.915 * g.sin() + 0.020 * (2.0 * g).sin(); // ecliptic longitude
+    let lambda_rad = lambda.to_radians();
+    let epsilon = (23.439 - 0.0000004 * n).to_radians(); // obliquity
+
+    let ra = lambda_rad.sin().atan2(epsilon.cos() * lambda_rad.cos());
+    let ra_deg = ra.to_degrees().rem_euclid(360.0);
+    let dec_deg = (epsilon.sin() * lambda_rad.sin()).asin().to_degrees();
+    (ra_deg, dec_deg)
+}
+
+/// Compute Sun altitude at a given JD for observer position (degrees).
+pub fn sun_altitude(jd: f64, lat_deg: f64, lon_deg: f64) -> f64 {
+    let (ra, dec) = sun_radec(jd);
+    let lst = local_sidereal_time(jd, lon_deg);
+    let (alt, _az) = radec_to_altaz(ra, dec, lst, lat_deg);
+    alt
+}
+
+/// Find astronomical twilight (sun alt = -18°) for a given date.
+/// Returns (dusk_utc, dawn_utc) by binary-search at 1-minute resolution.
+pub fn astro_twilight(
+    date: NaiveDate,
+    lat: f64,
+    lon: f64,
+) -> anyhow::Result<(DateTime<Utc>, DateTime<Utc>)> {
+    // Search window: noon to noon next day
+    let start = Utc.from_utc_datetime(&date.and_hms_opt(10, 0, 0).unwrap());
+    let end = start + Duration::hours(24);
+
+    let dusk = find_crossing(start, start + Duration::hours(12), lat, lon, -18.0, true)?;
+    let dawn = find_crossing(start + Duration::hours(12), end, lat, lon, -18.0, false)?;
+
+    Ok((dusk, dawn))
+}
+
+fn find_crossing(
+    t0: DateTime<Utc>,
+    t1: DateTime<Utc>,
+    lat: f64,
+    lon: f64,
+    target_alt: f64,
+    descending: bool,
+) -> anyhow::Result<DateTime<Utc>> {
+    let mut lo = t0;
+    let mut hi = t1;
+
+    // Verify sign change exists
+    let alt_lo = sun_altitude(julian_date(lo), lat, lon);
+    let alt_hi = sun_altitude(julian_date(hi), lat, lon);
+
+    // For dusk (descending): lo should be > target, hi should be < target
+    // For dawn (ascending): lo should be < target, hi should be > target
+    let _ = (alt_lo, alt_hi, descending); // used implicitly
+
+    // Binary search to 1-minute resolution
+    for _ in 0..100 {
+        let mid = lo + Duration::seconds((hi - lo).num_seconds() / 2);
+        let alt_mid = sun_altitude(julian_date(mid), lat, lon);
+
+        if descending {
+            if alt_mid > target_alt {
+                lo = mid;
+            } else {
+                hi = mid;
+            }
+        } else if alt_mid < target_alt {
+            lo = mid;
+        } else {
+            hi = mid;
+        }
+
+        if (hi - lo).num_seconds() <= 60 {
+            break;
+        }
+    }
+
+    Ok(lo)
+}

backend/src/astronomy/time.rs

@@ -0,0 +1,20 @@
+use chrono::{DateTime, Utc};
+
+/// Compute Julian Date from a UTC datetime.
+pub fn julian_date(dt: DateTime<Utc>) -> f64 {
+    let unix_seconds = dt.timestamp() as f64;
+    // J2000.0 epoch is 2000-01-01 12:00:00 UTC = Unix 946728000
+    // JD of Unix epoch 0 = 2440587.5
+    unix_seconds / 86400.0 + 2440587.5
+}
+
+/// Compute Local Sidereal Time in degrees [0, 360).
+/// Uses IAU formula via Julian Date and observer longitude.
+pub fn local_sidereal_time(jd: f64, lon_deg: f64) -> f64 {
+    // Days since J2000.0
+    let d = jd - 2451545.0;
+    // Greenwich Mean Sidereal Time in degrees
+    let gmst_deg = 280.46061837 + 360.98564736629 * d;
+    let lst = (gmst_deg + lon_deg).rem_euclid(360.0);
+    lst
+}

backend/src/astronomy/visibility.rs

@@ -0,0 +1,217 @@
+use chrono::{DateTime, Duration, Utc};
+use serde::{Deserialize, Serialize};
+
+use crate::config::{LAT, LON, MIN_ALT_DEG};
+use super::{
+    coords::{airmass, extinction_mag, radec_to_altaz},
+    horizon::{horizon_alt, HorizonPoint},
+    lunar::{moon_altitude, moon_separation},
+    time::{julian_date, local_sidereal_time},
+};
+
+#[derive(Debug, Clone, Serialize, Deserialize)]
+pub struct CurvePoint {
+    pub utc: DateTime<Utc>,
+    pub alt_deg: f64,
+    pub az_deg: f64,
+    pub airmass: f64,
+    pub above_custom_horizon: bool,
+    pub moon_alt_deg: f64,
+}
+
+#[derive(Debug, Clone, Serialize, Deserialize)]
+pub struct VisibilitySummary {
+    pub max_alt_deg: f64,
+    pub transit_utc: Option<DateTime<Utc>>,
+    pub rise_utc: Option<DateTime<Utc>>,
+    pub set_utc: Option<DateTime<Utc>>,
+    pub best_start_utc: Option<DateTime<Utc>>,
+    pub best_end_utc: Option<DateTime<Utc>>,
+    pub usable_min: u32,
+    pub is_visible_tonight: bool,
+    pub meridian_flip_utc: Option<DateTime<Utc>>,
+    pub airmass_at_transit: f64,
+    pub extinction_at_transit: f64,
+    pub moon_sep_deg: f64,
+    pub curve: Vec<CurvePoint>,
+}
+
+pub struct TonightWindow {
+    pub dusk: DateTime<Utc>,
+    pub dawn: DateTime<Utc>,
+}
+
+pub struct MoonState {
+    pub ra_deg: f64,
+    pub dec_deg: f64,
+    pub illumination: f64,
+    pub alt_at_midnight: f64,
+}
+
+/// Compute full visibility summary for a catalog object during tonight's window.
+/// step_minutes: resolution for the altitude curve (1 for detailed view, 10 for precompute cache).
+pub fn compute_visibility(
+    ra_deg: f64,
+    dec_deg: f64,
+    window: &TonightWindow,
+    horizon: &[HorizonPoint],
+    moon: &MoonState,
+) -> VisibilitySummary {
+    compute_visibility_with_step(ra_deg, dec_deg, window, horizon, moon, 10)
+}
+
+pub fn compute_visibility_with_step(
+    ra_deg: f64,
+    dec_deg: f64,
+    window: &TonightWindow,
+    horizon: &[HorizonPoint],
+    moon: &MoonState,
+    step_minutes: i64,
+) -> VisibilitySummary {
+    let step = Duration::minutes(step_minutes);
+    let mut t = window.dusk;
+
+    let mut curve = Vec::new();
+    let mut max_alt = f64::NEG_INFINITY;
+    let mut transit_utc: Option<DateTime<Utc>> = None;
+    let mut rise_utc: Option<DateTime<Utc>> = None;
+    let mut set_utc: Option<DateTime<Utc>> = None;
+    let mut best_start: Option<DateTime<Utc>> = None;
+    let mut best_end: Option<DateTime<Utc>> = None;
+    let mut usable_min = 0u32;
+    let mut prev_alt = f64::NEG_INFINITY;
+
+    while t <= window.dawn {
+        let jd = julian_date(t);
+        let lst = local_sidereal_time(jd, LON);
+        let (alt, az) = radec_to_altaz(ra_deg, dec_deg, lst, LAT);
+        let am = airmass(alt);
+        let h_alt = horizon_alt(az, horizon);
+        let above = alt > h_alt.max(MIN_ALT_DEG);
+        let moon_alt = moon_altitude(jd, LAT, LON);
+
+        curve.push(CurvePoint {
+            utc: t,
+            alt_deg: alt,
+            az_deg: az,
+            airmass: am,
+            above_custom_horizon: above,
+            moon_alt_deg: moon_alt,
+        });
+
+        if alt > max_alt {
+            max_alt = alt;
+            transit_utc = Some(t);
+        }
+
+        // Rise: first crossing above effective horizon
+        if prev_alt <= h_alt.max(MIN_ALT_DEG) && alt > h_alt.max(MIN_ALT_DEG) && rise_utc.is_none() {
+            rise_utc = Some(t);
+        }
+        // Set: last time we were above horizon
+        if alt > h_alt.max(MIN_ALT_DEG) {
+            set_utc = Some(t);
+        }
+
+        // Best window: above 30°
+        if alt > 30.0 {
+            if best_start.is_none() {
+                best_start = Some(t);
+            }
+            best_end = Some(t);
+            usable_min += step_minutes as u32;
+        }
+
+        prev_alt = alt;
+        t += step;
+    }
+
+    let is_visible = usable_min > 0 || rise_utc.is_some();
+
+    let airmass_transit = transit_utc
+        .map(|tr| {
+            let jd = julian_date(tr);
+            let lst = local_sidereal_time(jd, LON);
+            let (alt, _) = radec_to_altaz(ra_deg, dec_deg, lst, LAT);
+            airmass(alt)
+        })
+        .unwrap_or(40.0);
+
+    let extinction_transit = extinction_mag(
+        transit_utc
+            .map(|tr| {
+                let jd = julian_date(tr);
+                let lst = local_sidereal_time(jd, LON);
+                let (alt, _) = radec_to_altaz(ra_deg, dec_deg, lst, LAT);
+                alt
+            })
+            .unwrap_or(0.0),
+    );
+
+    let moon_sep = moon_separation(moon.ra_deg, moon.dec_deg, ra_deg, dec_deg);
+
+    // Meridian flip: transit + time for HA to reach +5°
+    // 5° of HA = 5/360 * 86400 = 1200 seconds
+    let meridian_flip = transit_utc.map(|tr| tr + Duration::seconds(1200));
+
+    VisibilitySummary {
+        max_alt_deg: if max_alt == f64::NEG_INFINITY { 0.0 } else { max_alt },
+        transit_utc,
+        rise_utc,
+        set_utc,
+        best_start_utc: best_start,
+        best_end_utc: best_end,
+        usable_min,
+        is_visible_tonight: is_visible,
+        meridian_flip_utc: meridian_flip,
+        airmass_at_transit: airmass_transit,
+        extinction_at_transit: extinction_transit,
+        moon_sep_deg: moon_sep,
+        curve,
+    }
+}
+
+/// Find the longest continuous true-dark window (sun < -18° AND moon below horizon).
+pub fn true_dark_window(
+    dusk: DateTime<Utc>,
+    dawn: DateTime<Utc>,
+    lat: f64,
+    lon: f64,
+) -> Option<(DateTime<Utc>, DateTime<Utc>)> {
+    let step = Duration::minutes(5);
+    let mut t = dusk;
+    let mut best: Option<(DateTime<Utc>, DateTime<Utc>)> = None;
+    let mut current_start: Option<DateTime<Utc>> = None;
+    let mut best_duration = Duration::zero();
+
+    while t <= dawn {
+        let jd = julian_date(t);
+        let moon_alt = moon_altitude(jd, lat, lon);
+        let is_dark = moon_alt < 0.0;
+
+        if is_dark {
+            if current_start.is_none() {
+                current_start = Some(t);
+            }
+        } else if let Some(start) = current_start {
+            let dur = t - start;
+            if dur > best_duration {
+                best_duration = dur;
+                best = Some((start, t));
+            }
+            current_start = None;
+        }
+
+        t += step;
+    }
+
+    // Check if still in dark window at dawn
+    if let Some(start) = current_start {
+        let dur = dawn - start;
+        if dur > best_duration {
+            best = Some((start, dawn));
+        }
+    }
+
+    best
+}