sky-cam/moon_phase.py

#!/usr/bin/env python3
"""moon_phase.py — moon ephemeris lookups for sky-cam.

Provides:
  - full_moons_in_range(start, end)           list of UTC datetimes of full moons
  - nearest_full_moon(when)                   UTC datetime of full moon nearest 'when'
  - illumination(when)                        0.0..1.0 fraction lit
  - phase_angle(when)                         degrees, 0=full, 180=new
  - waxing(when)                              True if moon is waxing
  - altaz(when, lat, lon)                     (altitude_deg, azimuth_deg)
  - parallactic_angle(when, lat, lon)         degrees, rotation to put celestial north up

Used as a library by other scripts and as a CLI:

  python3 moon_phase.py next-full              # next full moon UTC + local
  python3 moon_phase.py altaz <ISO-UTC>        # alt/az from configured location
  python3 moon_phase.py info <ISO-UTC>         # everything for one timestamp

All times assume UTC unless tagged otherwise.  Local timezone comes from
sky-cam.conf TIMEZONE for display only.
"""
from __future__ import annotations

import math
import os
import pathlib
import re
import sys
from datetime import datetime, timedelta, timezone

_here = pathlib.Path(__file__).resolve().parent


def _read_conf(path):
    conf = {}
    try:
        with open(path) as f:
            for line in f:
                line = line.strip()
                if not line or line.startswith('#') or '=' not in line:
                    continue
                k, v = line.split('=', 1)
                k = k.strip()
                v = v.strip()
                v = re.sub(r'\s+#.*$', '', v)
                v = v.strip('"').strip("'")
                m = re.match(r'^\$\{[^}]+:-([^}]*)\}$', v)
                if m:
                    v = m.group(1).strip('"').strip("'")
                conf[k] = v
    except FileNotFoundError:
        pass
    return conf


_conf = _read_conf(_here / 'sky-cam.conf')
_conf.update(_read_conf(_here / '.env'))

LATITUDE = float(_conf.get('LATITUDE') or 0.0)
LONGITUDE = float(_conf.get('LONGITUDE') or 0.0)
TIMEZONE = _conf.get('TIMEZONE', 'UTC')

# Cached ephemeris — bootstrap.sh pre-downloads de421.bsp into SCRIPT_DIR so
# this never needs to hit the network at run time.
_EPH_PATH = _here / 'de421.bsp'

_ts = None
_eph = None
_observer = None
_topos = None  # bare geographic position; sunrise_sunset() adds earth itself


def _lazy():
    """Defer skyfield import + ephemeris load until first use.

    Lets the module be imported by tests and by --help paths even when
    skyfield is missing or the ephemeris hasn't been downloaded yet.
    """
    global _ts, _eph, _observer, _topos
    if _ts is not None:
        return
    from skyfield.api import Loader, wgs84
    loader = Loader(str(_here), verbose=False)
    _ts = loader.timescale()
    if _EPH_PATH.exists():
        _eph = loader('de421.bsp')
    else:
        _eph = loader('de421.bsp')  # downloads on first run
    _topos = wgs84.latlon(LATITUDE, LONGITUDE)
    _observer = _eph['earth'] + _topos


def _to_utc(when: datetime) -> datetime:
    if when.tzinfo is None:
        return when.replace(tzinfo=timezone.utc)
    return when.astimezone(timezone.utc)


def _t(when: datetime):
    _lazy()
    w = _to_utc(when)
    return _ts.from_datetime(w)


def full_moons_in_range(start: datetime, end: datetime) -> list[datetime]:
    """Return UTC datetimes of every full moon between start and end (inclusive)."""
    return phase_events_in_range(start, end, 2)


def crescent_times_in_range(
    start: datetime,
    end: datetime,
    target_illum: float = 0.28,
    waxing_side: bool = True,
) -> list[datetime]:
    """Return UTC datetimes when illumination crosses target_illum on the given side.

    Scans in 2-hour steps and records each moment illumination passes through
    target_illum while waxing (waxing_side=True) or waning (waxing_side=False).
    Yields one event per lunar cycle, used for scheduling just like
    phase_events_in_range() is used for quarter/full events.
    """
    _lazy()
    step = timedelta(hours=2)
    results: list[datetime] = []
    t = _to_utc(start)
    end_dt = _to_utc(end)
    prev_illum: float | None = None
    prev_t: datetime | None = None

    while t <= end_dt:
        illum_val = illumination(t)
        is_wax = waxing(t)
        if is_wax == waxing_side:
            if prev_illum is not None:
                if (prev_illum - target_illum) * (illum_val - target_illum) < 0:
                    # Linear interpolation to the crossing moment
                    frac = (target_illum - prev_illum) / (illum_val - prev_illum)
                    results.append(prev_t + frac * (t - prev_t))
            prev_illum = illum_val
            prev_t = t
        else:
            prev_illum = None
            prev_t = None
        t += step

    return results


def phase_events_in_range(start: datetime, end: datetime, phase_index: int) -> list[datetime]:
    """Return UTC datetimes of every occurrence of `phase_index` between start and end.

    skyfield phase indices: 0 = new, 1 = first quarter, 2 = full, 3 = last quarter.
    """
    _lazy()
    from skyfield.almanac import find_discrete, moon_phases
    t0 = _t(start)
    t1 = _t(end)
    times, phases = find_discrete(t0, t1, moon_phases(_eph))
    return [t.utc_datetime() for t, p in zip(times, phases) if p == phase_index]


def nearest_full_moon(when: datetime) -> datetime:
    """Full moon UTC nearest 'when' — searches a 45-day window centred on it."""
    w = _to_utc(when)
    candidates = full_moons_in_range(w - timedelta(days=45), w + timedelta(days=45))
    return min(candidates, key=lambda d: abs(d - w))


def illumination(when: datetime) -> float:
    """Fraction of the moon's disk that is illuminated (0..1)."""
    _lazy()
    from skyfield.almanac import fraction_illuminated
    return float(fraction_illuminated(_eph, 'moon', _t(when)))


def phase_angle(when: datetime) -> float:
    """Sun-Moon-Earth phase angle in degrees: 0=full, 90=quarter, 180=new."""
    _lazy()
    earth = _eph['earth']
    sun = _eph['sun']
    moon = _eph['moon']
    t = _t(when)
    e = earth.at(t)
    s_from_moon = (sun - moon).at(t)
    e_from_moon = (e.position.au - moon.at(t).position.au)
    # angle between sun→moon→earth
    a = s_from_moon.position.au
    b = e_from_moon
    cosang = (a[0]*b[0] + a[1]*b[1] + a[2]*b[2]) / (
        math.sqrt(a[0]**2 + a[1]**2 + a[2]**2)
        * math.sqrt(b[0]**2 + b[1]**2 + b[2]**2)
    )
    cosang = max(-1.0, min(1.0, cosang))
    return math.degrees(math.acos(cosang))


def waxing(when: datetime) -> bool:
    """True if moon is waxing (illumination growing)."""
    now = illumination(when)
    later = illumination(when + timedelta(hours=6))
    return later > now


def altaz(when: datetime, lat: float | None = None, lon: float | None = None):
    """Apparent altitude/azimuth in degrees as seen from (lat, lon).

    Falls back to configured LATITUDE/LONGITUDE if not specified.
    """
    _lazy()
    from skyfield.api import wgs84
    if lat is None and lon is None:
        obs = _observer
    else:
        obs = _eph['earth'] + wgs84.latlon(
            LATITUDE if lat is None else lat,
            LONGITUDE if lon is None else lon,
        )
    t = _t(when)
    alt, az, _ = obs.at(t).observe(_eph['moon']).apparent().altaz()
    return float(alt.degrees), float(az.degrees)


def sun_events_in_range(
    search_start: datetime,
    search_end: datetime,
) -> list[tuple[datetime, bool]]:
    """Return every sunrise/sunset transition between search_start and search_end.

    Each item is (utc_datetime, is_rise): is_rise=True for sunrise, False for sunset.
    Uses skyfield's almanac — accurate to within a minute at the configured location.
    """
    _lazy()
    from skyfield.almanac import find_discrete, sunrise_sunset
    t0 = _t(search_start)
    t1 = _t(search_end)
    times, values = find_discrete(t0, t1, sunrise_sunset(_eph, _topos))
    return [(t.utc_datetime(), bool(v)) for t, v in zip(times, values)]


def sun_altaz(when: datetime, lat: float | None = None, lon: float | None = None):
    """Apparent altitude/azimuth of the sun in degrees as seen from (lat, lon)."""
    _lazy()
    from skyfield.api import wgs84
    if lat is None and lon is None:
        obs = _observer
    else:
        obs = _eph['earth'] + wgs84.latlon(
            LATITUDE if lat is None else lat,
            LONGITUDE if lon is None else lon,
        )
    t = _t(when)
    alt, az, _ = obs.at(t).observe(_eph['sun']).apparent().altaz()
    return float(alt.degrees), float(az.degrees)


def parallactic_angle(when: datetime, lat: float | None = None, lon: float | None = None) -> float:
    """Parallactic angle in degrees — rotates a moon image so celestial north is up.

    sin(q) = sin(H) * cos(lat) / cos(alt)
    where H is the hour angle and alt is the altitude.

    Returned value is the angle to rotate the lunar texture clockwise (in the
    image sense, y-down) so celestial north points up in the camera frame.
    """
    _lazy()
    from skyfield.api import wgs84
    if lat is None:
        lat = LATITUDE
    if lon is None:
        lon = LONGITUDE
    obs = _eph['earth'] + wgs84.latlon(lat, lon)
    t = _t(when)
    apparent = obs.at(t).observe(_eph['moon']).apparent()
    alt, _, _ = apparent.altaz()
    ra, dec, _ = apparent.radec(epoch='date')
    # local sidereal time at observer's longitude → hour angle
    lst = t.gast * 15.0 + lon  # gast in hours → degrees, plus longitude
    H = math.radians(lst - ra._degrees)
    phi = math.radians(lat)
    delta = math.radians(dec.degrees)
    # Standard parallactic angle formula
    q = math.atan2(math.sin(H), math.tan(phi) * math.cos(delta) - math.sin(delta) * math.cos(H))
    return math.degrees(q)


def _format_local(dt_utc: datetime) -> str:
    try:
        import pytz
        tz = pytz.timezone(TIMEZONE)
        return dt_utc.astimezone(tz).strftime('%Y-%m-%d %H:%M:%S %Z')
    except Exception:
        return dt_utc.strftime('%Y-%m-%d %H:%M:%S UTC')


def _cli():
    if len(sys.argv) < 2:
        print(__doc__.strip())
        return 1
    cmd = sys.argv[1]
    if cmd == 'next-full':
        now = datetime.now(timezone.utc)
        fm = full_moons_in_range(now, now + timedelta(days=45))
        if not fm:
            print('No full moon found in the next 45 days?', file=sys.stderr)
            return 1
        next_fm = fm[0]
        print(f'next_full_moon_utc={next_fm.strftime("%Y-%m-%dT%H:%M:%SZ")}')
        print(f'next_full_moon_local={_format_local(next_fm)}')
        return 0
    if cmd == 'nearest-full':
        when = datetime.fromisoformat(sys.argv[2].replace('Z', '+00:00'))
        fm = nearest_full_moon(when)
        print(f'nearest_full_moon_utc={fm.strftime("%Y-%m-%dT%H:%M:%SZ")}')
        return 0
    if cmd == 'altaz':
        when = datetime.fromisoformat(sys.argv[2].replace('Z', '+00:00'))
        alt, az = altaz(when)
        print(f'altitude_deg={alt:.3f}')
        print(f'azimuth_deg={az:.3f}')
        return 0
    if cmd == 'info':
        when = datetime.fromisoformat(sys.argv[2].replace('Z', '+00:00'))
        alt, az = altaz(when)
        print(f'utc={when.strftime("%Y-%m-%dT%H:%M:%SZ")}')
        print(f'local={_format_local(when)}')
        print(f'illumination={illumination(when):.4f}')
        print(f'phase_angle_deg={phase_angle(when):.2f}')
        print(f'waxing={waxing(when)}')
        print(f'altitude_deg={alt:.3f}')
        print(f'azimuth_deg={az:.3f}')
        print(f'parallactic_angle_deg={parallactic_angle(when):.3f}')
        return 0
    print(f'unknown command: {cmd}', file=sys.stderr)
    return 1


if __name__ == '__main__':
    sys.exit(_cli())