sky-cam/season_info.py

#!/usr/bin/env python3
"""
season_info.py — output shell variables describing the current (or given) date's
astronomical season and movement within that season.

Season boundaries are the actual equinoxes/solstices computed via the Jean Meeus
approximation (Astronomical Algorithms, Ch. 27 Table 27.a) and converted to the
LOCAL timezone so the calendar date matches what the camera sees on the ground.
DST transitions and year-to-year variation in the exact equinox moment are both
handled correctly.

Usage:
    eval "$(python3 season_info.py [YYYY-MM-DD] [Timezone])"

Arguments (both optional):
    YYYY-MM-DD   Date to evaluate; defaults to today.
    Timezone     IANA timezone name, e.g. America/New_York.
                 Can also be set via the TIMEZONE environment variable.
                 Defaults to America/New_York (matches the sunrise scripts).

Outputs (shell-sourceable):
    SEASON         Spring | Summer | Autumn | Winter
    MVT_NUM        1 | 2 | 3
    DAY_OF_MVT     day number within the movement (1-based)
    DAYS_IN_MVT    total days in this movement
    MVT_START      YYYY-MM-DD of movement start (local date)
    MVT_END        YYYY-MM-DD of movement end   (local date)
    SEASON_START   YYYY-MM-DD of season start   (local date)
    SEASON_END     YYYY-MM-DD of season end     (local date)
    IS_LAST_DAY    true | false
"""

import sys
import os
import datetime

DEFAULT_TZ = "America/New_York"

# ---------------------------------------------------------------------------
# Timezone helpers — prefer pytz (already installed for sunrise.py),
# fall back to zoneinfo (Python ≥ 3.9), then to UTC with a warning.
# ---------------------------------------------------------------------------
def _load_tz(name: str):
    try:
        import pytz
        return pytz.timezone(name)
    except ImportError:
        pass
    try:
        from zoneinfo import ZoneInfo
        return ZoneInfo(name)
    except ImportError:
        pass
    print(f"WARNING: neither pytz nor zoneinfo available; using UTC instead of {name}", file=sys.stderr)
    return datetime.timezone.utc


def _utc_datetime_from_jde(jde: float) -> datetime.datetime:
    """Convert Julian Day Number to a UTC datetime (including time of day)."""
    jde_s = jde + 0.5          # shift so integer part = calendar day
    Z = int(jde_s)
    F = jde_s - Z              # fractional day = fraction of 24 h after midnight

    if Z >= 2299161:
        alpha = int((Z - 1867216.25) / 36524.25)
        A = Z + 1 + alpha - alpha // 4
    else:
        A = Z
    B = A + 1524
    C = int((B - 122.1) / 365.25)
    D = int(365.25 * C)
    E = int((B - D) / 30.6001)

    day_frac   = B - D - int(30.6001 * E) + F
    day        = int(day_frac)
    frac_hours = (day_frac - day) * 24
    hour       = int(frac_hours)
    frac_mins  = (frac_hours - hour) * 60
    minute     = int(frac_mins)
    second     = int((frac_mins - minute) * 60)

    month = E - 1 if E < 14 else E - 13
    year  = C - 4716 if month > 2 else C - 4715

    return datetime.datetime(year, month, day, hour, minute, second,
                             tzinfo=datetime.timezone.utc)


def _season_starts_local(year: int, tz) -> dict:
    """
    Return the LOCAL calendar date (not UTC date) of each equinox/solstice
    for the given year.  'tz' is a pytz/zoneinfo timezone object.
    """
    y = (year - 2000) / 1000.0

    jdes = {
        "Spring": 2451623.80984 + 365242.37404*y + 0.05169*y**2 - 0.00411*y**3 - 0.00057*y**4,
        "Summer": 2451716.56767 + 365241.62603*y + 0.00325*y**2 + 0.00888*y**3 - 0.00030*y**4,
        "Autumn": 2451810.21715 + 365242.01767*y - 0.11575*y**2 + 0.00337*y**3 + 0.00078*y**4,
        "Winter": 2451900.05952 + 365242.74049*y - 0.06223*y**2 - 0.00823*y**3 + 0.00032*y**4,
    }

    result = {}
    for name, jde in jdes.items():
        utc_dt   = _utc_datetime_from_jde(jde)
        local_dt = utc_dt.astimezone(tz)
        result[name] = local_dt.date()
    return result


def get_info(date: datetime.date, tz_name: str = DEFAULT_TZ) -> dict:
    tz   = _load_tz(tz_name)
    year = date.year

    cur  = _season_starts_local(year,     tz)
    prev = _season_starts_local(year - 1, tz)
    nxt  = _season_starts_local(year + 1, tz)

    # Determine season and its inclusive start/end dates.
    # Winter straddles the year boundary (Dec solstice → next Mar equinox).
    if date < cur["Spring"]:
        season, s_start, s_end = "Winter", prev["Winter"], cur["Spring"] - datetime.timedelta(days=1)
    elif date < cur["Summer"]:
        season, s_start, s_end = "Spring", cur["Spring"], cur["Summer"] - datetime.timedelta(days=1)
    elif date < cur["Autumn"]:
        season, s_start, s_end = "Summer", cur["Summer"], cur["Autumn"] - datetime.timedelta(days=1)
    elif date < cur["Winter"]:
        season, s_start, s_end = "Autumn", cur["Autumn"], cur["Winter"] - datetime.timedelta(days=1)
    else:
        season, s_start, s_end = "Winter", cur["Winter"], nxt["Spring"] - datetime.timedelta(days=1)

    days_in_season = (s_end - s_start).days + 1

    # Split season into 3 movements.  Mvt 1 and 2 get floor(n/3) days;
    # Mvt 3 takes the remainder so no days are lost to integer rounding.
    m1 = days_in_season // 3
    m2 = days_in_season // 3
    m3 = days_in_season - m1 - m2

    m1_end = s_start   + datetime.timedelta(days=m1 - 1)
    m2_end = m1_end    + datetime.timedelta(days=m2)
    m3_end = s_end

    if date <= m1_end:
        mvt, mvt_start, mvt_end, days_in_mvt = 1, s_start,                        m1_end, m1
    elif date <= m2_end:
        mvt, mvt_start, mvt_end, days_in_mvt = 2, m1_end + datetime.timedelta(1), m2_end, m2
    else:
        mvt, mvt_start, mvt_end, days_in_mvt = 3, m2_end + datetime.timedelta(1), m3_end, m3

    return {
        "SEASON":       season,
        "MVT_NUM":      mvt,
        "DAY_OF_MVT":   (date - mvt_start).days + 1,
        "DAYS_IN_MVT":  days_in_mvt,
        "MVT_START":    mvt_start.isoformat(),
        "MVT_END":      mvt_end.isoformat(),
        "SEASON_START": s_start.isoformat(),
        "SEASON_END":   s_end.isoformat(),
        "IS_LAST_DAY":  str(date == mvt_end).lower(),
    }


if __name__ == "__main__":
    args    = [a for a in sys.argv[1:] if not a.startswith("-")]
    date_s  = args[0] if len(args) > 0 else None
    tz_name = args[1] if len(args) > 1 else os.environ.get("TIMEZONE", DEFAULT_TZ)

    date = datetime.date.fromisoformat(date_s) if date_s else datetime.date.today()

    for k, v in get_info(date, tz_name).items():
        print(f'{k}="{v}"')