#!/usr/bin/env python3
"""
Solar–Lunar Local Time Convergence Demo

This script illustrates the formalism in the pack:

- local solar phase:  theta_sun(t, lambda)  = (Theta(t) + lambda - alpha_sun(t)) mod 2pi
- local lunar phase:  theta_moon(t, lambda) = (Theta(t) + lambda - alpha_moon(t)) mod 2pi
- synodic phase:      phi(t)               = (alpha_moon(t) - alpha_sun(t)) mod 2pi
- matching identity:  theta_sun - theta_moon = phi(t) mod 2pi (independent of lambda)

It also shows equidistribution of phi_n under daily sampling in the ideal mean-motion model
when the daily increment is irrational.

Outputs:
- outputs/phi_hist.png
- outputs/phi_timeseries.png
"""

import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path

TAU = 2 * np.pi

def mod2pi(x):
    return np.mod(x, TAU)

def circ_dist(a, b):
    # minimal circular distance between angles a and b mod 2pi
    d = mod2pi(a - b)
    return np.minimum(d, TAU - d)

def main():
    outdir = Path("outputs")
    outdir.mkdir(exist_ok=True)

    # Sampling
    D = 86400.0  # mean solar day (s)
    N = 4000     # number of daily samples

    # Simplified mean-motion model (not a full ephemeris).
    omega_E = TAU / 86400.0           # Earth mean solar rotation phase rate
    omega_sun = TAU / (365.2422 * D)  # Sun inertial drift
    omega_moon = TAU / (27.32166 * D) # Moon inertial drift (sidereal month)

    t0 = 0.0
    n = np.arange(N)
    t = t0 + n * D

    Theta = mod2pi(omega_E * t)
    alpha_sun = mod2pi(omega_sun * t)
    alpha_moon = mod2pi(omega_moon * t)

    lam = 1.234  # longitude (radians)
    theta_sun = mod2pi(Theta + lam - alpha_sun)
    theta_moon = mod2pi(Theta + lam - alpha_moon)
    phi = mod2pi(alpha_moon - alpha_sun)

    diff = mod2pi(theta_sun - theta_moon)
    max_circ = float(np.max(circ_dist(diff, phi)))
    print("Max circular distance between (theta_sun - theta_moon) and phi:", max_circ)

    rho = float(mod2pi((omega_moon - omega_sun) * D))
    print("Daily increment rho (radians):", rho)
    print("rho / (2pi):", rho / float(TAU))

    # Histogram for equidistribution heuristic
    plt.figure()
    plt.hist(phi, bins=60, density=True)
    plt.xlabel("phi (radians)")
    plt.ylabel("density")
    plt.title("Histogram of sampled synodic phase phi_n")
    plt.tight_layout()
    plt.savefig(outdir / "phi_hist.png", dpi=200)

    # Timeseries segment
    plt.figure()
    plt.plot(n[:600], phi[:600])
    plt.xlabel("day index n")
    plt.ylabel("phi_n (radians)")
    plt.title("Synodic phase samples (first 600 days)")
    plt.tight_layout()
    plt.savefig(outdir / "phi_timeseries.png", dpi=200)

    print("Wrote:", outdir / "phi_hist.png")
    print("Wrote:", outdir / "phi_timeseries.png")

if __name__ == "__main__":
    main()
