Script: Solar Almanac Generator

#!/usr/bin/env python3
"""
Generate solar almanac data for the Recovery Library using Skyfield.

Produces tables-almanac.md with daily sun position data:
- GHA (Greenwich Hour Angle) at 00:00 UTC
- Declination at 00:00 UTC
- Equation of Time
- Semi-diameter

Uses the Skyfield library with JPL DE421 ephemeris for validated
astronomical accuracy (sub-arcsecond precision).

Usage:
    python generate_almanac.py [start_year] [end_year] [output_file]
    Default: 2026 2030 ../tables-almanac.md

Requirements:
    pip install skyfield
    (DE421 ephemeris will be downloaded automatically on first run, ~17MB)
"""

import sys
import os
import math
from datetime import datetime


def generate_almanac(start_year=2026, end_year=2030, output_file=None):
    # Import skyfield here so the script gives a clear error if not installed
    try:
        from skyfield.api import load
        from skyfield.almanac import find_discrete, seasons
    except ImportError:
        print("ERROR: skyfield not installed. Run: pip install skyfield")
        sys.exit(1)

    if output_file is None:
        output_file = os.path.join(
            os.path.dirname(os.path.abspath(__file__)), "..", "tables-almanac.md"
        )

    # Load ephemeris
    print("Loading JPL DE421 ephemeris...")
    eph_loader = load
    ts = load.timescale()
    eph = load('de421.bsp')
    sun = eph['sun']
    earth = eph['earth']

    num_years = end_year - start_year + 1
    total_days = num_years * 366  # approximate

    content = f"""---
title: "Solar Almanac Data: {start_year}–{end_year}"
subtitle: "Recovery Library — Computed Reference Data"
date: "{datetime.now().strftime('%Y-%m-%d')}"
---

# Solar Almanac Data: {start_year}–{end_year}

*Recovery Library — Computed Reference Data*

**Computation method:** All values computed using the Skyfield astronomy
library with the JPL DE421 planetary ephemeris.

**Reproduction:** Any computer with Python and the Skyfield library can
reproduce and extend these tables by running `scripts/generate_almanac.py`.

---

"""

    for year in range(start_year, end_year + 1):
        print(f"Computing {year}...")

        content += f"## {year}\n\n"
        content += "| Date | GHA Sun 00h | Dec | EoT | SD |\n"
        content += "|------|-------------|-----|-----|----|\n"

        import calendar
        days_in_year = 366 if calendar.isleap(year) else 365

        for doy in range(1, days_in_year + 1):
            from datetime import date as dt_date
            d = dt_date.fromordinal(dt_date(year, 1, 1).toordinal() + doy - 1)

            t = ts.utc(d.year, d.month, d.day, 0, 0, 0)
            astrometric = earth.at(t).observe(sun)
            ra, dec, distance = astrometric.apparent().radec()

            gast = t.gast
            gha_hours = gast - ra.hours
            if gha_hours < 0:
                gha_hours += 24
            gha_deg = gha_hours * 15
            if gha_deg >= 360:
                gha_deg -= 360

            gha_d = int(gha_deg)
            gha_m = (gha_deg - gha_d) * 60

            dec_deg = dec.degrees
            dec_sign = "N" if dec_deg >= 0 else "S"
            dec_abs = abs(dec_deg)
            dec_d = int(dec_abs)
            dec_m = (dec_abs - dec_d) * 60

            eot_deg = gha_deg - 180.0
            if eot_deg > 180:
                eot_deg -= 360
            if eot_deg < -180:
                eot_deg += 360
            eot_minutes = -eot_deg * 4
            eot_sign = "+" if eot_minutes >= 0 else "−"
            eot_abs = abs(eot_minutes)
            eot_m = int(eot_abs)
            eot_s = (eot_abs - eot_m) * 60

            sun_radius_km = 695700.0
            distance_km = distance.km
            sd_rad = math.atan(sun_radius_km / distance_km)
            sd_arcmin = math.degrees(sd_rad) * 60

            date_str = d.strftime("%b %d")
            content += (
                f"| {date_str} | {gha_d:3d}° {gha_m:04.1f}' | "
                f"{dec_sign} {dec_d:2d}° {dec_m:04.1f}' | "
                f"{eot_sign}{eot_m:02d}m {eot_s:04.1f}s | "
                f"{sd_arcmin:.1f}' |\n"
            )

        content += "\n---\n\n"

    # Annual summary with solstice/equinox dates
    content += "## Solstice, Equinox, and Perihelion Dates\n\n"
    content += "| Year | Perihelion | March Equinox | June Solstice "
    content += "| September Equinox | December Solstice |\n"
    content += "|------|------------|---------------|---------------"
    content += "|-------------------|-------------------|\n"

    for year in range(start_year, end_year + 1):
        t0 = ts.utc(year, 1, 1)
        t1 = ts.utc(year + 1, 1, 1)
        season_times, season_types = find_discrete(t0, t1, seasons(eph))

        season_names = {
            0: "mar_eq", 1: "jun_sol", 2: "sep_eq", 3: "dec_sol"
        }
        events = {}
        for st, stype in zip(season_times, season_types):
            name = season_names.get(stype, "unknown")
            events[name] = st.utc_strftime('%b %d %H:%M UTC')

        min_dist = float('inf')
        peri_date = ""
        for day in range(1, 15):
            t = ts.utc(year, 1, day, 12)
            _, _, dist = earth.at(t).observe(sun).apparent().radec()
            if dist.km < min_dist:
                min_dist = dist.km
                peri_date = f"Jan {day:02d} ~12h"

        content += (
            f"| {year} | {peri_date} | "
            f"{events.get('mar_eq', '—')} | "
            f"{events.get('jun_sol', '—')} | "
            f"{events.get('sep_eq', '—')} | "
            f"{events.get('dec_sol', '—')} |\n"
        )

    content += "\n---\n\n"
    content += (
        f"*Generated by `scripts/generate_almanac.py` using Skyfield "
        f"{_get_skyfield_version()} with JPL DE421 ephemeris.*\n"
    )

    with open(output_file, 'w') as f:
        f.write(content)

    print(f"Written to {output_file}")
    print(f"Total lines: {content.count(chr(10))}")


def _get_skyfield_version():
    try:
        import skyfield
        return skyfield.__version__
    except Exception:
        return "unknown"


if __name__ == "__main__":
    args = sys.argv[1:]
    start = int(args[0]) if len(args) > 0 else 2026
    end = int(args[1]) if len(args) > 1 else 2030
    out = args[2] if len(args) > 2 else None
    generate_almanac(start, end, out)