Script: Navigational Star Catalog and Charts Generator

#!/usr/bin/env python3
"""
Generate navigational star catalog and seasonal star charts for the Recovery Library.

Produces:
  ../tables-stars.md        — Star catalog, monthly visibility, rise/transit/set tables
  ../site/images/stars-january.png   — Evening sky chart, Jan 15 21:00 NZDT
  ../site/images/stars-april.png     — Evening sky chart, Apr 15 21:00 NZST
  ../site/images/stars-july.png      — Evening sky chart, Jul 15 21:00 NZST
  ../site/images/stars-october.png   — Evening sky chart, Oct 15 21:00 NZDT

Observer location: Wellington, New Zealand (41.3°S, 174.8°E)
Epoch: 2026.0 for SHA/Dec values

Usage:
    scripts/.venv/bin/python generate_star_catalog.py

Requirements:
    skyfield, matplotlib (both available in project venv)
    Hipparcos star catalog will be downloaded automatically on first run (~8 MB)
"""

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

# ---------------------------------------------------------------------------
# The 57 Nautical Almanac navigational stars
# Format: (common_name, bayer_designation, constellation, hip_id_or_None)
# HIP IDs from the Hipparcos catalog — used to look up precise positions.
# ---------------------------------------------------------------------------
NAV_STARS = [
    ("Acamar",      "θ Eri",  "Eridanus",    13847),
    ("Achernar",    "α Eri",  "Eridanus",     7588),
    ("Acrux",       "α Cru",  "Crux",        60718),
    ("Adhara",      "ε CMa",  "Canis Major", 33579),
    ("Aldebaran",   "α Tau",  "Taurus",      21421),
    ("Alioth",      "ε UMa",  "Ursa Major",  62956),
    ("Alkaid",      "η UMa",  "Ursa Major",  67301),
    ("Al Na'ir",    "α Gru",  "Grus",       109268),
    ("Alnilam",     "ε Ori",  "Orion",       26311),
    ("Alphard",     "α Hya",  "Hydra",       46390),
    ("Alphecca",    "α CrB",  "Corona Borealis", 76267),
    ("Alpheratz",   "α And",  "Andromeda",     677),
    ("Altair",      "α Aql",  "Aquila",      97649),
    ("Ankaa",       "α Phe",  "Phoenix",      2081),
    ("Antares",     "α Sco",  "Scorpius",    80763),
    ("Arcturus",    "α Boo",  "Bootes",      69673),
    ("Atria",       "α TrA",  "Triangulum Australe", 82273),
    ("Avior",       "ε Car",  "Carina",      41037),
    ("Bellatrix",   "γ Ori",  "Orion",       25336),
    ("Betelgeuse",  "α Ori",  "Orion",       27989),
    ("Canopus",     "α Car",  "Carina",      30438),
    ("Capella",     "α Aur",  "Auriga",      24608),
    ("Deneb",       "α Cyg",  "Cygnus",     102098),
    ("Denebola",    "β Leo",  "Leo",         57632),
    ("Diphda",      "β Cet",  "Cetus",        3419),
    ("Dubhe",       "α UMa",  "Ursa Major",  54061),
    ("Elnath",      "β Tau",  "Taurus",      25428),
    ("Eltanin",     "γ Dra",  "Draco",       87833),
    ("Enif",        "ε Peg",  "Pegasus",    107315),
    ("Fomalhaut",   "α PsA",  "Piscis Austrinus", 113368),
    ("Gacrux",      "γ Cru",  "Crux",        61084),
    ("Gienah",      "γ Crv",  "Corvus",      59803),
    ("Hadar",       "β Cen",  "Centaurus",   68702),
    ("Hamal",       "α Ari",  "Aries",        9884),
    ("Kaus Australis", "ε Sgr", "Sagittarius", 90185),
    ("Kochab",      "β UMi",  "Ursa Minor",  72607),
    ("Markab",      "α Peg",  "Pegasus",    113963),
    ("Menkar",      "α Cet",  "Cetus",       14135),
    ("Menkent",     "θ Cen",  "Centaurus",   68933),
    ("Miaplacidus", "β Car",  "Carina",      45238),
    ("Mirfak",      "α Per",  "Perseus",     15863),
    ("Nunki",       "σ Sgr",  "Sagittarius", 92855),
    ("Peacock",     "α Pav",  "Pavo",       100751),
    ("Pollux",      "β Gem",  "Gemini",      37826),
    ("Procyon",     "α CMi",  "Canis Minor", 37279),
    ("Rasalhague",  "α Oph",  "Ophiuchus",   86032),
    ("Regulus",     "α Leo",  "Leo",         49669),
    ("Rigel",       "β Ori",  "Orion",       24436),
    ("Rigil Kentaurus", "α Cen", "Centaurus", 71683),
    ("Sabik",       "η Oph",  "Ophiuchus",   84012),
    ("Schedar",     "α Cas",  "Cassiopeia",   3179),
    ("Shaula",      "λ Sco",  "Scorpius",    85927),
    ("Sirius",      "α CMa",  "Canis Major", 32349),
    ("Spica",       "α Vir",  "Virgo",       65474),
    ("Suhail",      "λ Vel",  "Vela",        44816),
    ("Vega",        "α Lyr",  "Lyra",        91262),
    ("Zubenelgenubi", "α Lib", "Libra",      72622),
]

# Wellington, NZ coordinates
WELLINGTON_LAT = -41.3      # degrees (negative = south)
WELLINGTON_LON = 174.8      # degrees east
WELLINGTON_ELEV = 25        # metres

# New Zealand timezone offsets
NZST_OFFSET = +12           # hours, standard time (April–September)
NZDT_OFFSET = +13           # hours, daylight time (October–March)


def get_nz_offset(month):
    """Return NZ UTC offset in hours for a given month."""
    # NZDT: last Sun in September to first Sun in April
    # Approximate by month: Oct–Mar = NZDT (+13), Apr–Sep = NZST (+12)
    if month in (10, 11, 12, 1, 2, 3):
        return NZDT_OFFSET
    else:
        return NZST_OFFSET


def load_skyfield():
    """Import skyfield modules and load ephemeris + Hipparcos catalog."""
    try:
        from skyfield.api import load, Star, wgs84
        from skyfield.data import hipparcos
    except ImportError:
        print("ERROR: skyfield not installed.")
        sys.exit(1)

    ts = load.timescale()
    eph = load('de421.bsp')
    earth = eph['earth']

    # Load Hipparcos catalog
    print("Loading Hipparcos catalog...")
    with load.open(hipparcos.URL) as f:
        hip_df = hipparcos.load_dataframe(f)

    return ts, earth, hip_df, Star, wgs84


def build_star_objects(hip_df, Star):
    """Build a list of (name, bayer, constellation, magnitude, Star object) for all 57 stars."""
    stars_out = []
    for (name, bayer, const, hip_id) in NAV_STARS:
        if hip_id not in hip_df.index:
            print(f"  WARNING: HIP {hip_id} ({name}) not found in catalog")
            continue
        row = hip_df.loc[hip_id]
        mag = row['magnitude']
        star_obj = Star.from_dataframe(hip_df.loc[hip_id:hip_id])
        stars_out.append((name, bayer, const, mag, star_obj))
    return stars_out


def compute_sha_dec(ts, earth, star_obj, epoch_year=2026):
    """
    Compute SHA and Declination for a star at a given epoch.
    SHA = 360 - RA (in degrees).
    Returns (sha_deg, dec_deg).
    """
    t = ts.tt_jd(ts.utc(epoch_year, 1, 1, 12).tt)
    astrometric = earth.at(t).observe(star_obj)
    ra, dec, _ = astrometric.radec(epoch='date')
    ra_h = ra.hours
    dec_d = dec.degrees
    # Skyfield may return numpy arrays when Star is from a dataframe
    if hasattr(ra_h, '__len__'):
        ra_h = ra_h[0] if len(ra_h) == 1 else ra_h
    if hasattr(dec_d, '__len__'):
        dec_d = dec_d[0] if len(dec_d) == 1 else dec_d
    ra_deg = float(ra_h) * 15.0
    sha_deg = 360.0 - ra_deg
    if sha_deg >= 360.0:
        sha_deg -= 360.0
    return sha_deg, float(dec_d)


def _scalar(val):
    """Convert a Skyfield angle attribute (possibly a numpy array) to a Python float."""
    if hasattr(val, '__len__'):
        return float(val[0]) if len(val) == 1 else float(val)
    return float(val)


def fmt_angle(deg, signed=True):
    """Format a decimal degree angle as degrees and arcminutes."""
    if signed:
        sign = "N" if deg >= 0 else "S"
        d = int(abs(deg))
        m = (abs(deg) - d) * 60.0
        return f"{sign}{d:02d}° {m:04.1f}'"
    else:
        d = int(deg)
        m = (deg - d) * 60.0
        return f"{d:3d}° {m:04.1f}'"


def build_catalog_table(ts, earth, star_list):
    """Build the 57-star catalog table as a markdown string."""
    rows = []
    for (name, bayer, const, mag, star_obj) in star_list:
        sha, dec = compute_sha_dec(ts, earth, star_obj)
        sha_str = fmt_angle(sha, signed=False)
        dec_str = fmt_angle(dec, signed=True)
        rows.append((sha, name, bayer, sha_str, dec_str, mag, const))

    # Sort by SHA
    rows.sort(key=lambda r: r[0])

    lines = []
    lines.append("| # | Name | Bayer | SHA | Dec | Mag | Constellation |")
    lines.append("|---|------|-------|-----|-----|-----|---------------|")
    for i, (sha_raw, name, bayer, sha_str, dec_str, mag, const) in enumerate(rows, 1):
        lines.append(f"| {i:2d} | {name} | {bayer} | {sha_str} | {dec_str} | {mag:.1f} | {const} |")

    return "\n".join(lines), rows


def compute_monthly_visibility(ts, earth, star_list, lat=WELLINGTON_LAT, lon=WELLINGTON_LON):
    """
    For each month, find which of the 57 stars are above 10° altitude
    at 21:00 local time on the 15th, as seen from Wellington.
    Returns dict: month_num -> list of (name, alt_deg, az_deg)
    """
    from skyfield.api import wgs84
    observer = wgs84.latlon(lat, lon, elevation_m=WELLINGTON_ELEV)

    monthly = {}
    for month in range(1, 13):
        utc_offset = get_nz_offset(month)
        # 21:00 local = 21:00 - offset in UTC
        utc_hour = 21 - utc_offset
        utc_day = 15
        utc_month = month
        if utc_hour < 0:
            utc_hour += 24
            utc_day = 14  # previous day in UTC

        t = ts.utc(2026, utc_month, utc_day, utc_hour, 0, 0)
        visible = []
        for (name, bayer, const, mag, star_obj) in star_list:
            astrometric = (earth + observer).at(t).observe(star_obj)
            alt, az, _ = astrometric.apparent().altaz()
            if _scalar(alt.degrees) > 10.0:
                visible.append((name, _scalar(alt.degrees), _scalar(az.degrees), mag))

        # Sort by altitude descending
        visible.sort(key=lambda x: -x[1])
        monthly[month] = visible

    return monthly


def build_visibility_table(monthly_visibility):
    """Build markdown for monthly visibility section."""
    month_names = [
        "January", "February", "March", "April", "May", "June",
        "July", "August", "September", "October", "November", "December"
    ]
    lines = []
    for month in range(1, 13):
        tz_name = "NZDT" if get_nz_offset(month) == NZDT_OFFSET else "NZST"
        visible = monthly_visibility[month]
        lines.append(f"### {month_names[month-1]} (21:00 {tz_name}, 15th)")
        lines.append(f"*{len(visible)} stars above 10° altitude*\n")
        lines.append("| Star | Altitude | Azimuth | Mag |")
        lines.append("|------|----------|---------|-----|")
        for (name, alt, az, mag) in visible:
            lines.append(f"| {name} | {alt:.1f}° | {az:.1f}° | {mag:.1f} |")
        lines.append("")

    return "\n".join(lines)


def compute_rise_transit_set(ts, earth, star_list, lat=WELLINGTON_LAT, lon=WELLINGTON_LON):
    """
    Compute rise, transit, set times for the 10 brightest navigational stars
    for the 15th of each month from Wellington.
    Returns dict: star_name -> list of (month, rise_str, transit_str, set_str)
    """
    from skyfield.api import wgs84
    from skyfield.almanac import find_discrete, risings_and_settings

    observer_loc = wgs84.latlon(lat, lon, elevation_m=WELLINGTON_ELEV)
    observer = earth + observer_loc

    # Pick 10 brightest
    by_mag = sorted(star_list, key=lambda x: x[3])[:10]
    by_mag_names = [s[0] for s in by_mag]

    month_names = ["Jan", "Feb", "Mar", "Apr", "May", "Jun",
                   "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"]

    results = {}
    for (name, bayer, const, mag, star_obj) in by_mag:
        results[name] = []
        for month in range(1, 13):
            utc_offset = get_nz_offset(month)
            # Search window: 00:00 to 24:00 local time on the 15th
            t0 = ts.utc(2026, month, 15, -utc_offset, 0, 0)
            t1 = ts.utc(2026, month, 16, -utc_offset, 0, 0)

            try:
                f = risings_and_settings(earth, star_obj, observer_loc,
                                         horizon_degrees=-0.5667)
                times, events = find_discrete(t0, t1, f)

                rise_str = "—"
                set_str = "—"
                for ti, ev in zip(times, events):
                    local_dt = ti.utc_datetime().replace(tzinfo=timezone.utc)
                    local_hour = local_dt.hour + utc_offset
                    local_min = local_dt.minute
                    if local_hour >= 24:
                        local_hour -= 24
                    time_str = f"{local_hour:02d}:{local_min:02d}"
                    if ev == 1:   # rising
                        rise_str = time_str
                    else:         # setting
                        set_str = time_str

            except Exception:
                rise_str = "circ"
                set_str = "circ"

            # Compute transit: when star reaches highest altitude
            # Sample every 10 minutes through the day to find max
            try:
                best_alt = -90.0
                best_t = None
                for h in range(0, 144):  # 0 to 24h in 10-min steps
                    frac_h = h * 10.0 / 60.0
                    t_sample = ts.utc(2026, month, 15,
                                      int(frac_h) - utc_offset,
                                      int((frac_h % 1) * 60), 0)
                    alt, _, _ = observer.at(t_sample).observe(star_obj).apparent().altaz()
                    if _scalar(alt.degrees) > best_alt:
                        best_alt = _scalar(alt.degrees)
                        best_t = t_sample

                if best_t is not None and best_alt > 0:
                    local_dt = best_t.utc_datetime().replace(tzinfo=timezone.utc)
                    local_hour = local_dt.hour + utc_offset
                    local_min = local_dt.minute
                    if local_hour >= 24:
                        local_hour -= 24
                    transit_str = f"{local_hour:02d}:{local_min:02d}"
                else:
                    transit_str = "—"
            except Exception:
                transit_str = "—"

            results[name].append((month_names[month-1], rise_str, transit_str, set_str))

    return results, by_mag_names


def build_rise_transit_set_table(rts_data, star_names):
    """Build markdown for rise/transit/set section."""
    lines = []
    for name in star_names:
        if name not in rts_data:
            continue
        lines.append(f"### {name}")
        lines.append("| Month | Rise | Transit | Set |")
        lines.append("|-------|------|---------|-----|")
        for (month_str, rise, transit, setting) in rts_data[name]:
            lines.append(f"| {month_str} | {rise} | {transit} | {setting} |")
        lines.append("")
    return "\n".join(lines)


# ---------------------------------------------------------------------------
# Star chart generation
# ---------------------------------------------------------------------------

CHART_DATES = [
    ("january",  2026, 1,  15, "January 2026",  "NZDT"),
    ("april",    2026, 4,  15, "April 2026",    "NZST"),
    ("july",     2026, 7,  15, "July 2026",     "NZST"),
    ("october",  2026, 10, 15, "October 2026",  "NZDT"),
]


def make_star_chart(ts, earth, star_list, month_tag, year, month, day,
                    title_month, tz_name, output_path, hip_df=None, StarClass=None):
    """
    Generate a circular alt/az star chart for the evening sky from Wellington.
    Center = zenith, outer ring = horizon. Uses altitude/azimuth projection.
    """
    try:
        import matplotlib
        matplotlib.use('Agg')
        import matplotlib.pyplot as plt
        import matplotlib.patches as mpatches
        import numpy as np
    except ImportError:
        print("ERROR: matplotlib not installed.")
        return

    from skyfield.api import wgs84

    utc_offset = get_nz_offset(month)
    # 21:00 local time
    utc_hour = 21 - utc_offset
    utc_day_actual = day
    if utc_hour < 0:
        utc_hour += 24
        utc_day_actual = day - 1

    t = ts.utc(year, month, utc_day_actual, utc_hour, 0, 0)
    observer = wgs84.latlon(WELLINGTON_LAT, WELLINGTON_LON, elevation_m=WELLINGTON_ELEV)

    # Compute alt/az for each star
    chart_stars = []
    for (name, bayer, const, mag, star_obj) in star_list:
        astrometric = (earth + observer).at(t).observe(star_obj)
        alt, az, _ = astrometric.apparent().altaz()
        if _scalar(alt.degrees) > -5.0:  # include stars just below horizon for labeling
            chart_stars.append((name, _scalar(alt.degrees), _scalar(az.degrees), mag))

    # ---------- matplotlib figure ----------
    fig, ax = plt.subplots(figsize=(9, 9), facecolor='#090920')
    ax.set_facecolor('#090920')
    ax.set_xlim(-1.15, 1.15)
    ax.set_ylim(-1.15, 1.15)
    ax.set_aspect('equal')
    ax.axis('off')

    import numpy as np

    # Draw horizon circle
    theta_ring = np.linspace(0, 2 * np.pi, 360)
    ax.plot(np.cos(theta_ring), np.sin(theta_ring),
            color='#4488cc', linewidth=1.5, zorder=2)

    # Draw altitude grid circles (30°, 60°)
    for alt_line in (30, 60):
        r = 1.0 - alt_line / 90.0
        ax.plot(r * np.cos(theta_ring), r * np.sin(theta_ring),
                color='#223355', linewidth=0.6, linestyle='--', zorder=1)
        ax.text(0, r + 0.02, f"{alt_line}°",
                color='#334466', fontsize=6, ha='center', va='bottom')

    # Draw azimuth lines every 45°
    for az_line in range(0, 360, 45):
        az_rad = math.radians(az_line)
        # In our projection: N is up (az=0), E is left (az=90)
        # We map: x = sin(az), y = cos(az) — standard compass-to-Cartesian
        # But to preserve astronomical convention (E on left when facing south)
        # we flip x: x = -sin(az), y = cos(az)
        x_end = -math.sin(az_rad)
        y_end = math.cos(az_rad)
        ax.plot([0, x_end], [0, y_end],
                color='#1a2a3a', linewidth=0.5, zorder=1)

    # Cardinal direction labels
    # N up, E left (southern hemisphere sky view facing up = south at top is
    # conventional for SH star charts; but we use N-up for simplicity)
    cardinals = [
        (0,   "N",  0,    1.10),
        (90,  "E", -1.10, 0),
        (180, "S",  0,   -1.10),
        (270, "W",  1.10, 0),
    ]
    for (az_c, label, x_c, y_c) in cardinals:
        ax.text(x_c, y_c, label, color='#88aacc', fontsize=12,
                fontweight='bold', ha='center', va='center', zorder=5)

    # Plot stars
    for (name, alt_deg, az_deg, mag) in chart_stars:
        if alt_deg < 0:
            continue
        # r = 0 at zenith (alt=90°), r = 1 at horizon (alt=0°)
        r = 1.0 - alt_deg / 90.0
        az_rad = math.radians(az_deg)
        # Flip x for astronomical convention
        x = -r * math.sin(az_rad)
        y = r * math.cos(az_rad)

        # Star size: brighter = larger dot (reduced to ~50% of original)
        # Magnitude scale: mag 1 -> size ~27, mag 5 -> size ~4
        size = max(2, 45 - 9 * mag)

        # Color by magnitude
        if mag < 1.0:
            color = '#ffffff'
        elif mag < 2.0:
            color = '#ffffee'
        elif mag < 3.0:
            color = '#eeeedd'
        else:
            color = '#ccccbb'

        ax.scatter([x], [y], s=size, c=color, zorder=4, edgecolors='none')

        # Labels: only above horizon, avoid crowding
        if alt_deg > 5.0:
            offset_x = 0.015 if x >= 0 else -0.015
            ha = 'left' if x >= 0 else 'right'
            ax.text(x + offset_x, y + 0.015, name,
                    color='#ffee88', fontsize=5.5,
                    ha=ha, va='bottom', zorder=5,
                    fontfamily='monospace')

    # --- Constellation lines ---
    # Build a lookup from star name to (x, y) chart position for visible stars
    star_xy = {}
    for (name, alt_deg, az_deg, mag) in chart_stars:
        if alt_deg < 0:
            continue
        r = 1.0 - alt_deg / 90.0
        az_rad = math.radians(az_deg)
        sx = -r * math.sin(az_rad)
        sy = r * math.cos(az_rad)
        star_xy[name] = (sx, sy)

    # Extra (non-navigational) stars needed for constellation lines
    # These are not in the 57-star NAV_STARS list
    extra_stars_hip = {
        "Mimosa":       62434,   # Beta Crucis (Crux short arm)
        "Delta Crucis": 59747,   # Delta Crucis (Crux short arm)
        "Alnitak":      26727,   # Zeta Orionis (belt)
        "Mintaka":      25930,   # Delta Orionis (belt)
        "Saiph":        27366,   # Kappa Orionis (foot)
        "Sargas":       86228,   # Theta Scorpii (curve)
    }

    # Compute alt/az for extra stars and add them to star_xy if above horizon
    if hip_df is not None and StarClass is not None:
        for extra_name, extra_hip in extra_stars_hip.items():
            if extra_name in star_xy:
                continue  # Already present from nav star list
            if extra_hip not in hip_df.index:
                continue
            extra_star_obj = StarClass.from_dataframe(hip_df.loc[extra_hip:extra_hip])
            astrometric = (earth + observer).at(t).observe(extra_star_obj)
            alt_e, az_e, _ = astrometric.apparent().altaz()
            alt_val = _scalar(alt_e.degrees)
            az_val = _scalar(az_e.degrees)
            if alt_val > 0:
                r = 1.0 - alt_val / 90.0
                az_rad = math.radians(az_val)
                star_xy[extra_name] = (-r * math.sin(az_rad), r * math.cos(az_rad))

    # Define constellation line segments: (star_a, star_b)
    constellation_lines = [
        # Southern Cross (Crux)
        ("Acrux", "Gacrux"),            # long arm
        ("Mimosa", "Delta Crucis"),     # short arm
        # Orion
        ("Betelgeuse", "Bellatrix"),    # shoulders
        ("Rigel", "Saiph"),            # feet
        ("Alnilam", "Alnitak"),        # belt segment 1
        ("Alnilam", "Mintaka"),        # belt segment 2
        # Scorpius
        ("Antares", "Sargas"),         # head through curve
        ("Sargas", "Shaula"),          # curve to stinger
        # Centaurus pointer line
        ("Rigil Kentaurus", "Hadar"),  # Alpha Cen to Beta Cen
    ]

    # Draw lines only when both endpoint stars are visible (above horizon)
    for (star_a, star_b) in constellation_lines:
        if star_a in star_xy and star_b in star_xy:
            xa, ya = star_xy[star_a]
            xb, yb = star_xy[star_b]
            ax.plot([xa, xb], [ya, yb],
                    color='white', alpha=0.4, linewidth=0.8, zorder=3)

    # Zenith marker
    ax.scatter([0], [0], s=20, c='#8888ff', zorder=6, marker='+')
    ax.text(0, 0.04, "Zenith", color='#8888ff', fontsize=6, ha='center', va='bottom')

    # Title
    ax.set_title(
        f"Evening Sky from Wellington — {title_month}\n"
        f"21:00 {tz_name} · Observer: 41.3°S, 174.8°E",
        color='white', fontsize=11, pad=14,
        fontweight='bold'
    )

    # Legend for altitude rings
    legend_elements = [
        mpatches.Patch(facecolor='#090920', edgecolor='#4488cc',
                       label='Horizon (0°)'),
        mpatches.Patch(facecolor='#090920', edgecolor='#223355',
                       linestyle='--', label='30° / 60° altitude'),
    ]
    ax.legend(handles=legend_elements, loc='lower right',
              facecolor='#0d1530', edgecolor='#334466',
              labelcolor='#aabbcc', fontsize=7, framealpha=0.8)

    # Attribution footer
    fig.text(0.5, 0.01,
             "Star positions: Hipparcos catalogue via Skyfield. "
             "Ephemeris: JPL DE421. Computed for epoch 2026.0.",
             ha="center", fontsize=6, color="#7799bb", style="italic")

    os.makedirs(os.path.dirname(output_path), exist_ok=True)
    fig.savefig(output_path, dpi=200, bbox_inches='tight',
                facecolor='#090920', edgecolor='none')
    plt.close(fig)
    print(f"  Saved: {output_path}")


# ---------------------------------------------------------------------------
# Main
# ---------------------------------------------------------------------------

def main():
    script_dir = os.path.dirname(os.path.abspath(__file__))
    md_out = os.path.join(script_dir, "..", "tables-stars.md")
    img_dir = os.path.join(script_dir, "..", "site", "images")

    print("Loading Skyfield libraries and ephemeris...")
    ts, earth, hip_df, Star, wgs84 = load_skyfield()

    print("Building star objects from Hipparcos catalog...")
    star_list = build_star_objects(hip_df, Star)
    print(f"  Loaded {len(star_list)} of {len(NAV_STARS)} stars")

    # --- Section 1: SHA/Dec catalog ---
    print("Computing SHA and Declination for epoch 2026.0...")
    catalog_table, sorted_rows = build_catalog_table(ts, earth, star_list)

    # --- Section 2: Monthly visibility ---
    print("Computing monthly visibility (21:00 local, Wellington)...")
    monthly_vis = compute_monthly_visibility(ts, earth, star_list)
    visibility_md = build_visibility_table(monthly_vis)

    # --- Section 3: Rise/transit/set for 10 brightest ---
    print("Computing rise/transit/set times for 10 brightest stars...")
    rts_data, bright_names = compute_rise_transit_set(ts, earth, star_list)
    rts_md = build_rise_transit_set_table(rts_data, bright_names)

    # --- Assemble markdown ---
    now_str = datetime.now().strftime("%Y-%m-%d")

    content = f"""---
title: "Navigational Star Catalog and Charts"
subtitle: "Recovery Library — Computed Reference Data"
date: "{now_str}"
---

# Navigational Star Catalog and Charts

*Recovery Library — Doc #015: Star Atlas*

**Computation method:** All values computed using the [Skyfield](https://rhodesmill.org/skyfield/) astronomy library with the JPL DE421 planetary ephemeris and the Hipparcos star catalog (ESA, 1997). SHA and Declination are computed for epoch 2026.0. Rise/transit/set times and visibility tables are computed for Wellington, New Zealand (41.3°S, 174.8°E).

**SHA definition:** Sidereal Hour Angle (SHA) = 360° − Right Ascension (in degrees). This is the standard convention used in the Nautical Almanac. SHA increases westward from the vernal equinox.

**Accuracy:** Star positions are accurate to better than 1 arcminute for 2026, well within sextant observation precision. Proper motion corrections are applied via the Hipparcos catalog data.

**Verification:** Canopus (α Carinae) should show approximately S52° 42' declination. Sirius (α CMa) should show approximately S16° 43' declination.

**Observer location for visibility and rise/set tables:** Wellington, New Zealand — latitude 41.3°S, longitude 174.8°E, elevation 25 m. Times use NZDT (UTC+13) for October–March and NZST (UTC+12) for April–September.

---

## Seasonal Sky Charts

These charts show the evening sky (21:00 local time) from Wellington for the 15th of a representative month in each season. N is up, E is left (the standard orientation when facing south in the Southern Hemisphere, or equivalently the orientation of the sky above you when lying on your back with head pointing north).

Altitude rings are shown at 30° and 60°. The outer circle is the horizon.

| Season | Chart |
|--------|-------|
| Summer (January) | ![January Evening Sky](images/stars-january.png) |
| Autumn (April) | ![April Evening Sky](images/stars-april.png) |
| Winter (July) | ![July Evening Sky](images/stars-july.png) |
| Spring (October) | ![October Evening Sky](images/stars-october.png) |

---

## Part 1: The 57 Navigational Stars — SHA and Declination, Epoch 2026.0

Stars are sorted by SHA (increasing westward from the vernal equinox). This matches the ordering used in the Nautical Almanac.

{catalog_table}

---

## Part 2: Monthly Evening Visibility from Wellington (41.3°S)

Stars visible above 10° altitude at 21:00 local time on the 15th of each month.
Stars below 10° altitude are subject to significant atmospheric refraction and difficult to
observe reliably with a sextant. Stars below the horizon are omitted.

{visibility_md}

---

## Part 3: Rise, Transit, and Set Times — 10 Brightest Navigational Stars

Times are local Wellington time (NZDT UTC+13 for Oct–Mar; NZST UTC+12 for Apr–Sep),
computed for the 15th of each month 2026. "—" indicates the star does not rise or set
on that date (circumpolar or below horizon all day). "circ" indicates a circumpolar star.

Transit time is the moment of maximum altitude (meridian passage), computed by sampling
altitude at 10-minute intervals and finding the peak. Times are accurate to approximately
±10 minutes by this method.

{rts_md}

---

## Notes on Use for Celestial Navigation

**Selecting stars for a fix:** Ideally choose 3 stars with azimuths separated by approximately
120°, all at altitudes between 15° and 70°. Stars below 15° suffer from uncertain atmospheric
refraction; stars above 70° are difficult to observe accurately with a marine sextant.

**SHA and GHA relationship:** GHA of a star = GHA of Aries + SHA of star (mod 360°).
GHA of Aries can be read from the nautical almanac or computed from sidereal time.

**Magnitude and naked-eye visibility:** Stars of magnitude 1.0 or brighter are easily visible
under urban skies. Stars to magnitude 3.0 are visible under suburban conditions. All 57
navigational stars are visible to the naked eye under reasonably dark skies.

**Southern Hemisphere considerations:** From Wellington at 41.3°S, stars with declination
south of about −48.7° are circumpolar (never set). Stars with declination north of about
+48.7° never rise. The north celestial pole (near Polaris) is permanently below the horizon;
there is no bright southern pole star — sigma Octantis (mag 5.4) is too faint for reliable
naked-eye observation.

---

*Generated by `scripts/generate_star_catalog.py` using Skyfield with JPL DE421 ephemeris
and the Hipparcos star catalog. Run date: {now_str}.*
"""

    os.makedirs(os.path.dirname(md_out), exist_ok=True)
    with open(md_out, 'w') as f:
        f.write(content)
    print(f"\nMarkdown written to: {md_out}")
    print(f"  Lines: {content.count(chr(10))}")

    # --- Generate star chart PNGs ---
    print("\nGenerating seasonal star charts...")
    os.makedirs(img_dir, exist_ok=True)
    for (tag, yr, mo, dy, title_mo, tz_name) in CHART_DATES:
        print(f"  Rendering {tag} chart...")
        out_png = os.path.join(img_dir, f"stars-{tag}.png")
        make_star_chart(ts, earth, star_list,
                        tag, yr, mo, dy, title_mo, tz_name, out_png,
                        hip_df=hip_df, StarClass=Star)

    print("\nDone.")
    print(f"Markdown: {os.path.abspath(md_out)}")
    print(f"Images:   {os.path.abspath(img_dir)}/stars-{{january,april,july,october}}.png")


if __name__ == "__main__":
    main()