Script: NZ Tide Table Generator

#!/usr/bin/env python3
"""
Generate NZ tide tables for the Recovery Library.

Produces tables-tides.md with high/low water predictions for 2026–2030
at three NZ standard ports: Auckland, Wellington, and Nelson.

Uses the standard harmonic prediction formula:
    h(t) = Z0 + sum(An * fn * cos(wn*t + Vn + un - gn))

Harmonic constants are from LINZ published data as reproduced in
Doc #12 (Tide Tables), Appendix A. Only the principal constituents
(M2, S2, N2, K1, O1) are used — sufficient for predictions accurate
to ~15–20 min in time and 0.1–0.2 m in height.

Includes nodal modulation factors for the 18.61-year lunar node cycle.
Uses only Python standard library (no external dependencies).

Usage:
    python generate_tide_tables.py [start_year] [end_year] [output_file]
    Default: 2026 2030 ../tables-tides.md
"""

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

# ---------------------------------------------------------------------------
# Tidal constituent angular speeds (degrees per solar hour)
# From Schureman, Manual of Harmonic Analysis and Prediction of Tides (1958)
# ---------------------------------------------------------------------------
CONSTITUENT_SPEEDS = {
    'M2': 28.984104,   # Principal lunar semidiurnal
    'S2': 30.000000,   # Principal solar semidiurnal
    'N2': 28.439730,   # Larger lunar elliptic semidiurnal
    'K1': 15.041069,   # Luni-solar diurnal
    'O1': 13.943036,   # Principal lunar diurnal
}

# Periods in hours (for reference)
CONSTITUENT_PERIODS = {k: 360.0 / v for k, v in CONSTITUENT_SPEEDS.items()}

# ---------------------------------------------------------------------------
# Port harmonic constants
# From Doc #12 Appendix A (representative LINZ values)
# Z0 = mean water level above chart datum (metres)
# For each constituent: (amplitude_m, phase_g_degrees)
# ---------------------------------------------------------------------------
PORTS = {
    'Auckland (Waitematā)': {
        'Z0': 1.74,
        'constituents': {
            'M2': (1.07, 255),
            'S2': (0.21, 289),
            'N2': (0.21, 237),
            'K1': (0.07, 175),
            'O1': (0.04, 151),
        },
        'region': 'Auckland',
    },
    'Wellington': {
        'Z0': 0.75,
        'constituents': {
            'M2': (0.44, 174),
            'S2': (0.07, 219),
            'N2': (0.09, 154),
            'K1': (0.06,  21),
            'O1': (0.04, 350),
        },
        'region': 'Wellington',
    },
    'Nelson': {
        'Z0': 2.15,
        'constituents': {
            'M2': (1.44, 179),
            'S2': (0.29, 215),
            'N2': (0.28, 160),
            'K1': (0.07, 352),
            'O1': (0.04, 322),
        },
        'region': 'Tasman Bay',
    },
}

# ---------------------------------------------------------------------------
# Astronomical functions
# ---------------------------------------------------------------------------

def julian_date(year, month, day, hour=0.0):
    """Compute Julian Date from calendar date (UTC)."""
    if month <= 2:
        year -= 1
        month += 12
    A = int(year / 100)
    B = 2 - A + int(A / 4)
    JD = (int(365.25 * (year + 4716)) + int(30.6001 * (month + 1))
          + day + hour / 24.0 + B - 1524.5)
    return JD


def julian_centuries(jd):
    """Julian centuries from J2000.0 epoch."""
    return (jd - 2451545.0) / 36525.0


def lunar_node_longitude(T):
    """
    Longitude of the mean ascending node of the lunar orbit (degrees).
    T = Julian centuries from J2000.0.
    From Meeus, Astronomical Algorithms.
    """
    N = 125.0445479 - 1934.1362891 * T + 0.0020754 * T**2
    return N % 360.0


def mean_lunar_longitude(T):
    """Mean longitude of the Moon (s), degrees. From Schureman/Meeus."""
    s = 218.3164477 + 481267.88123421 * T - 0.0015786 * T**2
    return s % 360.0


def mean_solar_longitude(T):
    """Mean longitude of the Sun (h), degrees."""
    h = 280.46646 + 36000.76983 * T + 0.0003032 * T**2
    return h % 360.0


def mean_lunar_perigee(T):
    """Mean longitude of the lunar perigee (p), degrees."""
    p = 83.3532465 + 4069.0137287 * T - 0.0103200 * T**2
    return p % 360.0


def mean_solar_perigee(T):
    """Mean longitude of the solar perigee (p1), degrees."""
    p1 = 282.93768 + 1.71946 * T + 0.00045688 * T**2
    return p1 % 360.0


def nodal_factors(T):
    """
    Compute nodal modulation factors f and u for each constituent.
    T = Julian centuries from J2000.0.

    f = amplitude modulation factor (dimensionless, ~1.0)
    u = phase modulation (degrees)

    Simplified formulas from Schureman (1958) and Pugh & Woodworth (2014).
    """
    N = math.radians(lunar_node_longitude(T))
    cos_N = math.cos(N)
    sin_N = math.sin(N)

    # M2: f = 1.0004 - 0.0373*cos(N) + 0.0002*cos(2N)
    #     u = -2.14*sin(N) degrees
    f_M2 = 1.0004 - 0.0373 * cos_N + 0.0002 * math.cos(2 * N)
    u_M2 = -2.14 * math.sin(N)

    # S2: no nodal modulation (purely solar)
    f_S2 = 1.0
    u_S2 = 0.0

    # N2: same nodal factors as M2
    f_N2 = f_M2
    u_N2 = u_M2

    # K1: f = 1.0060 + 0.1150*cos(N) - 0.0088*cos(2N)
    #     u = -8.86*sin(N) degrees
    f_K1 = 1.0060 + 0.1150 * cos_N - 0.0088 * math.cos(2 * N)
    u_K1 = -8.86 * math.sin(N)

    # O1: f = 1.0089 + 0.1871*cos(N) - 0.0147*cos(2N)
    #     u = 10.80*sin(N) degrees
    f_O1 = 1.0089 + 0.1871 * cos_N - 0.0147 * math.cos(2 * N)
    u_O1 = 10.80 * math.sin(N)

    return {
        'M2': (f_M2, u_M2),
        'S2': (f_S2, u_S2),
        'N2': (f_N2, u_N2),
        'K1': (f_K1, u_K1),
        'O1': (f_O1, u_O1),
    }


def astronomical_arguments(T):
    """
    Compute equilibrium arguments V0 for each constituent at T.
    T = Julian centuries from J2000.0.

    V0 is the phase of the equilibrium tide at Greenwich at the given time.
    Formulas from Schureman (1958), Table 2.

    Uses astronomical variables:
        s = mean lunar longitude
        h = mean solar longitude
        p = mean lunar perigee longitude
    """
    s = mean_lunar_longitude(T)
    h = mean_solar_longitude(T)
    p = mean_lunar_perigee(T)

    # Equilibrium arguments (V0) for each constituent
    # These follow from the Doodson numbers
    V0 = {
        'M2': 2 * (h - s) % 360,           # 2h - 2s
        'S2': 0.0,                           # S2 V0 = 0 by definition
        'N2': (2 * h - 3 * s + p) % 360,    # 2h - 3s + p
        'K1': (h + 90) % 360,               # h + 90°
        'O1': (h - 2 * s - 90) % 360,       # h - 2s - 90°
    }
    return V0


def predict_height(dt, port_data, T_ref=None, nf=None, V0_ref=None):
    """
    Predict tide height at a given datetime for a port.

    dt: datetime object (UTC)
    port_data: dict with Z0 and constituents
    T_ref: Julian centuries at start of year (cached)
    nf: nodal factors (cached)
    V0_ref: astronomical arguments at start of year (cached)

    Returns height in metres above chart datum.
    """
    jd = julian_date(dt.year, dt.month, dt.day,
                     dt.hour + dt.minute / 60.0 + dt.second / 3600.0)

    if T_ref is None:
        T_ref = julian_centuries(julian_date(dt.year, 1, 1))
    if nf is None:
        nf = nodal_factors(T_ref)
    if V0_ref is None:
        V0_ref = astronomical_arguments(T_ref)

    # Hours since the start of the year
    jd_start = julian_date(dt.year, 1, 1)
    t_hours = (jd - jd_start) * 24.0

    h = port_data['Z0']
    for name, (amp, phase_g) in port_data['constituents'].items():
        f, u = nf[name]
        V0 = V0_ref[name]
        speed = CONSTITUENT_SPEEDS[name]  # degrees per hour
        angle = speed * t_hours + V0 + u - phase_g
        h += amp * f * math.cos(math.radians(angle))

    return h


def find_extrema(dt_start, dt_end, port_data, step_minutes=6):
    """
    Find all high and low water events between dt_start and dt_end.

    Uses a coarse scan at step_minutes intervals, then refines each
    turning point to 1-second precision using bisection.

    Returns list of (datetime, height, type) where type is 'HW' or 'LW'.
    """
    # Cache astronomical parameters for the year
    T_ref = julian_centuries(julian_date(dt_start.year, 1, 1))
    nf = nodal_factors(T_ref)
    V0_ref = astronomical_arguments(T_ref)

    step = timedelta(minutes=step_minutes)
    extrema = []

    t = dt_start
    h_prev = predict_height(t, port_data, T_ref, nf, V0_ref)
    t += step
    h_curr = predict_height(t, port_data, T_ref, nf, V0_ref)
    t += step
    direction_prev = 1 if h_curr > h_prev else -1

    while t <= dt_end:
        h_next = predict_height(t, port_data, T_ref, nf, V0_ref)
        direction_curr = 1 if h_next > h_curr else -1

        if direction_curr != direction_prev and direction_prev != 0:
            # Turning point between t-step and t
            t_lo = t - 2 * step
            t_hi = t
            ext_type = 'HW' if direction_prev > 0 else 'LW'

            # Refine to ~1 second
            for _ in range(20):
                t_mid = t_lo + (t_hi - t_lo) / 2
                h_lo = predict_height(
                    t_lo + (t_mid - t_lo) / 2, port_data, T_ref, nf, V0_ref)
                h_hi = predict_height(
                    t_mid + (t_hi - t_mid) / 2, port_data, T_ref, nf, V0_ref)

                if ext_type == 'HW':
                    if h_lo > h_hi:
                        t_hi = t_mid
                    else:
                        t_lo = t_mid
                else:
                    if h_lo < h_hi:
                        t_hi = t_mid
                    else:
                        t_lo = t_mid

            t_ext = t_lo + (t_hi - t_lo) / 2
            h_ext = predict_height(t_ext, port_data, T_ref, nf, V0_ref)
            extrema.append((t_ext, h_ext, ext_type))

        h_prev = h_curr
        h_curr = h_next
        direction_prev = direction_curr
        t += step

    return extrema


def format_time_nzst(dt_utc):
    """Convert UTC datetime to NZST (UT+12) and format as HH:MM."""
    dt_nzst = dt_utc + timedelta(hours=12)
    return dt_nzst.strftime('%H:%M'), dt_nzst


def day_of_week(dt):
    """Return 2-letter day abbreviation."""
    days = ['Mo', 'Tu', 'We', 'Th', 'Fr', 'Sa', 'Su']
    return days[dt.weekday()]


def moon_phase_approx(dt):
    """
    Approximate moon phase. Returns symbol or empty string.
    Based on the synodic month (29.53059 days) from a known new moon.
    """
    known_new_moon = datetime(2026, 1, 29, 12, 36)
    synodic = 29.53059
    days_since = (dt - known_new_moon).total_seconds() / 86400.0
    phase = (days_since % synodic) / synodic

    if abs(phase) < 0.02 or abs(phase - 1.0) < 0.02:
        return '●'  # New moon
    elif abs(phase - 0.5) < 0.02:
        return '○'  # Full moon
    return ''


def generate_month_table(year, month, port_name, port_data):
    """Generate tide table for one port for one month."""
    import calendar

    days_in_month = calendar.monthrange(year, month)[1]
    month_name = calendar.month_name[month].upper()

    dt_start = datetime(year, month, 1) - timedelta(hours=12)
    if month == 12:
        dt_end = datetime(year + 1, 1, 1) + timedelta(hours=12)
    else:
        dt_end = datetime(year, month + 1, 1) + timedelta(hours=12)

    extrema = find_extrema(dt_start, dt_end, port_data)

    # Group extrema by NZST date
    daily = {}
    for t_utc, h, etype in extrema:
        _, t_nzst = format_time_nzst(t_utc)
        day = t_nzst.day
        m = t_nzst.month
        y = t_nzst.year
        if y == year and m == month:
            if day not in daily:
                daily[day] = []
            time_str = t_nzst.strftime('%H:%M')
            daily[day].append((time_str, h, etype))

    lines = []
    lines.append(f'**{port_name.upper()} — {month_name} {year}**')
    lines.append(
        'Times in NZST (UT +12:00). Heights in metres above Chart Datum (LAT).')
    lines.append('')
    lines.append('```')
    lines.append(
        'Day   HW Time  Ht(m)   LW Time  Ht(m)   '
        'HW Time  Ht(m)   LW Time  Ht(m)')

    for day in range(1, days_in_month + 1):
        dt_day = datetime(year, month, day)
        dow = day_of_week(dt_day)
        moon = moon_phase_approx(dt_day)

        events = daily.get(day, [])
        events.sort(key=lambda x: x[0])

        line_parts = f'{day:2d} {dow}'
        for t, h, e in events[:4]:
            line_parts += f'  {t}  {h:5.1f}'

        if moon:
            line_parts += f'  {moon}'

        lines.append(line_parts)

    lines.append('```')
    lines.append('')
    lines.append('● = New Moon    ○ = Full Moon')
    lines.append('')
    return '\n'.join(lines)


def generate_tables(start_year=2026, end_year=2030, output_file=None):
    """Generate the complete tide tables document."""
    if output_file is None:
        output_file = os.path.join(
            os.path.dirname(os.path.abspath(__file__)),
            '..', 'tables-tides.md')

    print(f"Generating NZ tide tables for {start_year}–{end_year}...")

    lines = []
    lines.append('---')
    lines.append('title: "NZ Tide Tables: 2026–2030"')
    lines.append('subtitle: "Recovery Library — Computed Reference Data"')
    lines.append(f'date: "{datetime.now().strftime("%Y-%m-%d")}"')
    lines.append('---')
    lines.append('')
    lines.append('# NZ Tide Tables: 2026–2030')
    lines.append('')
    lines.append('*Recovery Library — Computed Reference Data*')
    lines.append('')
    lines.append('**Computation method:** Harmonic tidal prediction ...')
    # [Full header text omitted for brevity in this display;
    #  see the actual script for complete output formatting]
    lines.append('---')

    for port_name, port_data in PORTS.items():
        lines.append('')
        lines.append(f'## {port_name}')
        lines.append('')
        for year in range(start_year, end_year + 1):
            lines.append(f'### {year}')
            lines.append('')
            for month in range(1, 13):
                print(f"  {port_name} — {year}/{month:02d}")
                table = generate_month_table(
                    year, month, port_name, port_data)
                lines.append(table)

    content = '\n'.join(lines)

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

    print(f"\nTide tables written to {output_file}")


if __name__ == '__main__':
    start = int(sys.argv[1]) if len(sys.argv) > 1 else 2026
    end = int(sys.argv[2]) if len(sys.argv) > 2 else 2030
    out = sys.argv[3] if len(sys.argv) > 3 else None
    generate_tables(start, end, out)