#!/usr/bin/env python3
"""tide_stats.py — 全地点の潮汐特性を一括計算する

static/tide/<年>/*.txt(気象庁の年間毎時潮位)と static/tide/stations.json から、
地点ごとの潮汐特性を計算して static/tide/stats.json に書き出す。

計算する項目:
  age_of_tide  潮齢[日]      新月・満月から、実際に潮差が最大になるまでの遅れ
  mhwi         平均高潮間隔[時] 月がその地点の南中を過ぎてから満潮になるまでの時間
  max_rate     年間最大の潮位変化速度[cm/h]
  mean_range   年間平均の日潮差[cm]  (その日の最高潮位 - 最低潮位 の年平均)
  max_range    年間最大の日潮差[cm]
  mean_level   平均水面[cm]

注意: mean_range / max_range は当サイト独自の定義(日ごとの実測差の集計)であり、
      気象庁が公表している「大潮差」(調和定数から定義される値)とは別物。

使い方:
  python tide_stats.py --year 2026
"""

import argparse
import json
import math
import re
from pathlib import Path

D2R = math.pi / 180
SYNODIC = 29.530588853
DAYMS = 86400.0            # 秒
J1970 = 2440588
J2000 = 2451545
OBLIQ = D2R * 23.4397


# ---------- 気象庁テキストのパース ----------
def parse_jma(text: str) -> list[tuple[float, int]]:
    """(エポック秒, 潮位cm) のリストを返す。時刻はJSTの各正時。"""
    import datetime as dt

    rows = []
    for line in text.splitlines():
        if len(line) < 80:
            continue
        try:
            yy = int(line[72:74])
            mm = int(line[74:76])
            dd = int(line[76:78])
        except ValueError:
            continue
        for h in range(24):
            chunk = line[h * 3:h * 3 + 3].strip()
            if not chunk:
                continue
            try:
                level = int(chunk)
            except ValueError:
                continue
            t = dt.datetime(2000 + yy, mm, dd, h,
                            tzinfo=dt.timezone(dt.timedelta(hours=9)))
            rows.append((t.timestamp(), level))
    rows.sort(key=lambda r: r[0])
    return rows


def gradient(y: list[int]) -> list[float]:
    """中央差分(端は前進/後退差分)。numpy.gradient と同じ考え方。"""
    n = len(y)
    if n < 2:
        return [0.0] * n
    g = [0.0] * n
    g[0] = float(y[1] - y[0])
    g[-1] = float(y[-1] - y[-2])
    for i in range(1, n - 1):
        g[i] = (y[i + 1] - y[i - 1]) / 2.0
    return g


# ---------- 天文計算 ----------
def moon_age(ts: float) -> float:
    """月齢[日]。月と太陽の黄経差から求める(主要な摂動項まで)。"""
    jd = ts / DAYMS + 2440587.5
    T = (jd - 2451545.0) / 36525
    Lm = 218.3164477 + 481267.88123421 * T
    Ls = 280.46646 + 36000.76983 * T
    Mm = 134.9633964 + 477198.8675055 * T
    Ms = 357.5291092 + 35999.0502909 * T
    Dm = 297.8501921 + 445267.1114034 * T
    lm = (Lm + 6.289 * math.sin(Mm * D2R) - 1.274 * math.sin((2 * Dm - Mm) * D2R)
          - 0.658 * math.sin(2 * Dm * D2R) - 0.214 * math.sin(2 * Mm * D2R)
          - 0.186 * math.sin(Ms * D2R) - 0.114 * math.sin(2 * (Lm - Ls) * D2R))
    ls = Ls + 1.915 * math.sin(Ms * D2R) + 0.020 * math.sin(2 * Ms * D2R)
    phase = ((lm - ls) % 360 + 360) % 360
    return phase / 360 * SYNODIC


def syzygies(t0: float, t1: float) -> list[float]:
    """期間内の朔(新月)と望(満月)の時刻を二分法で求める。"""
    out = []
    prev = moon_age(t0)
    t = t0 + 3600
    while t < t1:
        cur = moon_age(t)
        if cur < prev:                       # 月齢が29.5→0に戻る = 朔
            out.append(_refine(t - 3600, t, 0.0))
        if prev < SYNODIC / 2 <= cur:        # 望
            out.append(_refine(t - 3600, t, SYNODIC / 2))
        prev = cur
        t += 3600
    return out


def _refine(a: float, b: float, target: float) -> float:
    for _ in range(26):
        m = (a + b) / 2
        vm = moon_age(m)
        if target == 0.0:
            if vm > SYNODIC / 2:
                a = m
            else:
                b = m
        else:
            if vm < target:
                a = m
            else:
                b = m
    return (a + b) / 2


def moon_hour_angle(ts: float, lon: float) -> float:
    """月の時角(南中で0、-πからπに正規化)。"""
    d = ts / DAYMS - 0.5 + J1970 - J2000
    L = D2R * (218.316 + 13.176396 * d)
    M = D2R * (134.963 + 13.064993 * d)
    F = D2R * (93.272 + 13.229350 * d)
    l = L + D2R * 6.289 * math.sin(M)
    b = D2R * 5.128 * math.sin(F)
    ra = math.atan2(math.sin(l) * math.cos(OBLIQ) - math.tan(b) * math.sin(OBLIQ),
                    math.cos(l))
    H = D2R * (280.16 + 360.9856235 * d) - D2R * (-lon) - ra
    return ((H + math.pi) % (2 * math.pi) + 2 * math.pi) % (2 * math.pi) - math.pi


# ---------- 特性の計算 ----------
def daily_ranges(rows) -> list[tuple[float, int]]:
    """(その日の0時のエポック秒, 日潮差cm) のリスト。"""
    import datetime as dt

    jst = dt.timezone(dt.timedelta(hours=9))
    byday = {}
    for ts, lv in rows:
        d = dt.datetime.fromtimestamp(ts, jst).date()
        lo, hi = byday.get(d, (10 ** 9, -10 ** 9))
        byday[d] = (min(lo, lv), max(hi, lv))
    out = []
    for d in sorted(byday):
        lo, hi = byday[d]
        t0 = dt.datetime(d.year, d.month, d.day, tzinfo=jst).timestamp()
        out.append((t0, hi - lo))
    return out


def age_of_tide(rows) -> float | None:
    """潮齢: 朔望から日潮差が最大になる日までの遅れ[日]。放物線補間で小数日まで。"""
    days = daily_ranges(rows)
    if len(days) < 30:
        return None
    lags = []
    for s in syzygies(rows[0][0], rows[-1][0]):
        best_i, best_r = None, -1
        for i, (t0, rng) in enumerate(days):
            diff = (t0 - s) / DAYMS
            if -1 <= diff <= 6 and rng > best_r:
                best_r, best_i = rng, i
        if best_i is None or best_i == 0 or best_i == len(days) - 1:
            continue
        y0, y1, y2 = days[best_i - 1][1], days[best_i][1], days[best_i + 1][1]
        den = y0 - 2 * y1 + y2
        dx = 0.5 * (y0 - y2) / den if den != 0 else 0.0
        peak = days[best_i][0] + (0.5 + dx) * DAYMS   # その日の正午を基準
        lags.append((peak - s) / DAYMS)
    return sum(lags) / len(lags) if lags else None


def mean_high_water_interval(rows, lon: float, sample_days: int = 90) -> float | None:
    """平均高潮間隔: 月の南中から次の満潮までの時間[時]。"""
    levels = [lv for _, lv in rows]
    # 満潮時刻(毎時潮位の極大を放物線補間)
    highs = []
    for i in range(1, len(rows) - 1):
        if levels[i] >= levels[i - 1] and levels[i] > levels[i + 1]:
            y0, y1, y2 = levels[i - 1], levels[i], levels[i + 1]
            den = y0 - 2 * y1 + y2
            dx = 0.5 * (y0 - y2) / den if den != 0 else 0.0
            highs.append(rows[i][0] + dx * 3600)
    if not highs:
        return None

    lags = []
    start = rows[0][0]
    for k in range(sample_days):
        day = start + k * DAYMS
        prev = moon_hour_angle(day, lon)
        for j in range(1, 145):
            t = day + j * 600
            cur = moon_hour_angle(t, lon)
            if prev < 0 <= cur:                       # 南中を通過
                a, b = t - 600, t
                for _ in range(22):
                    m = (a + b) / 2
                    if moon_hour_angle(m, lon) < 0:
                        a = m
                    else:
                        b = m
                tr = (a + b) / 2
                for hw in highs:
                    if hw > tr:
                        lag = (hw - tr) / 3600
                        if lag < 12.5:
                            lags.append(lag)
                        break
            prev = cur
    return sum(lags) / len(lags) if lags else None


def stats_for(rows, lon: float | None) -> dict:
    levels = [lv for _, lv in rows]
    g = gradient(levels)
    ranges = [r for _, r in daily_ranges(rows)]
    out = {
        "max_rate": round(max(abs(x) for x in g), 1),
        "mean_range": round(sum(ranges) / len(ranges), 1),
        "max_range": max(ranges),
        "mean_level": round(sum(levels) / len(levels), 1),
    }
    a = age_of_tide(rows)
    out["age_of_tide"] = round(a, 2) if a is not None else None
    if lon is not None:
        m = mean_high_water_interval(rows, lon)
        out["mhwi"] = round(m, 2) if m is not None else None
    else:
        out["mhwi"] = None
    return out


def main() -> None:
    p = argparse.ArgumentParser()
    p.add_argument("--year", type=int, required=True)
    p.add_argument("--tide-dir", default="static/tide")
    p.add_argument("--only", help="地点記号をカンマ区切りで指定(動作確認用)")
    args = p.parse_args()

    tide_dir = Path(args.tide_dir)
    stations = json.loads((tide_dir / "stations.json").read_text(encoding="utf-8"))
    targets = set(args.only.split(",")) if args.only else None

    result = {}
    codes = sorted(stations)
    for i, code in enumerate(codes, 1):
        if targets and code not in targets:
            continue
        f = tide_dir / str(args.year) / f"{code}.txt"
        if not f.exists():
            print(f"  [{i}/{len(codes)}] {code} データなし、スキップ")
            continue
        info = stations[code]
        if isinstance(info, str):          # 旧形式(名前だけ)
            info = {"name": info, "lat": None, "lon": None}
        rows = parse_jma(f.read_text(encoding="ascii", errors="replace"))
        if len(rows) < 24 * 300:
            print(f"  [{i}/{len(codes)}] {code} {info['name']} データ不足、スキップ")
            continue
        s = stats_for(rows, info.get("lon"))
        s["name"] = info["name"]
        s["lat"] = info.get("lat")
        s["lon"] = info.get("lon")
        result[code] = s
        print(f"  [{i}/{len(codes)}] {code} {info['name']}: "
              f"潮齢 {s['age_of_tide']}日 / 高潮間隔 {s['mhwi']}時 / "
              f"最大 {s['max_rate']}cm/h / 平均潮差 {s['mean_range']}cm")

    outfile = tide_dir / "stats.json"
    outfile.write_text(json.dumps(result, ensure_ascii=False, indent=1),
                       encoding="utf-8")
    print(f"\n完了: {len(result)} 地点 → {outfile}")


if __name__ == "__main__":
    main()