#!/usr/bin/env python3
"""
===============================================================================
HDGL COUNTER-ROTATING DRAIN — ANIMATED v9
===============================================================================

SOURCE OF TRUTH: hdgl_unified_force_fine_cross-checked.hdgl

FULL INFERENCE — every decision justified:

━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━

ROTATION:
  arg(L+(i)) = arctan( Im(L+(i)) / Re(L+(i)) )
  As substrate index i advances, Λ_φ(i) increases (logarithmically).
  The phase arm 1_eff·e^(iπΛ_φ) rotates.
  L+ rotates. L- = conj(L+) rotates in the opposite direction.
  There is no rotation rate, no accumulated angle, no theta_seed.
  The counter-rotation IS the conjugate pair. It cannot be otherwise.

RADIAL POSITION:
  Each substrate point lives at |L+(i)| in the plane.
  |L+(i)| = |A(z) + 1_eff·e^(iπΛ_φ)| — this varies as i advances
  because 1_eff(i) changes with β = {Λ_φ(i)}.
  The drain is not imposed. The radial position IS |L|.
  When L approaches its fixed point (1_eff→1, δ→0), |L| stabilises.
  The 'drain' is L converging to its attractor.

HOURGLASS GEOMETRY:
  Each substrate point is assigned to one of 8 Kuramoto modes.
  Its 3D z-coordinate = Ω(f_mode) — the resonance value of its mode.
  +Ω arm: z = +Ω_mode. -Ω arm: z = -Ω_mode.
  The xy-plane position = (Re(L+)/scale, Im(L+)/scale) at that mode's z.
  The hourglass shape emerges because different modes have different |L|
  and different Ω values — not from any imposed geometry.

THE n_H = 1.291 FRIEDMANN SCALING:
  The cross-checked file specifies n_H from the Friedmann equation.
  This applies to the DENSITY of modes along the z-axis:
  modes near z=0 (low Ω, z_cosmo≈0) are the present universe.
  modes near z=1 (high Ω, z_cosmo≈2) are the early universe.
  Hubble scaling says density ∝ (1+z_cosmo)^1.291.
  We assign N substrate points to 8 modes with counts weighted by
  (1 + 2*Ω_mode)^1.291 — more points at high-z (outer arm) modes.
  This is the only place n_H=1.291 appears. It sets point distribution.

1_eff HAS ONE FORMULA, TWO NATURAL REGIMES:
  δ(i,n) = |cos(πβφ)| · ln(P_n) / φ^(n+β)
  At small n (n=1,2): δ ≈ 0.2 — this is the animation regime.
  At large n (n>>1): δ→0 — this is the CODATA/atomic regime.
  At large frame i with large n: δ→CODATA values.
  No interpolation. One formula. The regime is a consequence of (i,n).

KURAMOTO K — DIRECTLY FROM THE .hdgl SPEC:
  PLUCK    CV≥0.50: K=5.0
  SUSTAIN  CV<0.50: K=3.0
  FINETUNE CV<0.30: K=2.0
  LOCK     CV<0.05: K=1.8  (wu-wei — do not force)
  CV = std(|L_k|²) / mean(|L_k|²) across 8 modes.
  K is set by APhase. Not derived from U*.

phi_tower U*:
  U* = φ^(φ^(φ^(Σ sin(θ_i-θ_j))))
  U* → φ^φ ≈ 2.178 as phases lock (all sin terms → 0).
  U* is a CONVERGENCE INDICATOR. It is shown in the overlay.
  It does NOT drive K. K is driven by APhase/CV.

THE NULL NODE z=-1:
  At z=-1: (1+z)^n = 0^n = 0 for all n ≥ 1. A-arm vanishes.
  L±(-1) = 0 + 1_eff · e^(±iπΛ_φ) = pure phase.
  The torus at mode depth Ω=0 (z_mode=-1) has magnitude = 1_eff only.
  This is NOT zero amplitude — it is pure phase arm, no A contribution.
  It is the eye of the vortex: amplitude cancelled, phase pure.

STILLNESS → WHIRLPOOL:
  L is always rotating (Λ_φ always advances with i).
  The whirlpool becomes VISIBLE as Kuramoto coherence R rises.
  At R≈0 (staggered initial phases): torus alpha=0, arms faint.
  As K couples modes and R rises: torus and arms become visible.
  This IS the plug being pulled — not rotation starting, but
  the interference structure becoming coherent enough to see.

━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━

Output: hdgl_drain_animation_v9.mp4
Requires: matplotlib numpy scipy ffmpeg
"""

import math
import cmath
from pathlib import Path

import numpy as np
from scipy.ndimage import gaussian_filter1d
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.animation import FFMpegWriter
from matplotlib.colors import Normalize
from matplotlib.lines import Line2D


# ============================================================================
# LAYER 0 — T(X) = 1 + 1/X
# ============================================================================

def T(x):
    return 1.0 + 1.0 / x

def emergent_phi():
    x = 1.5
    for _ in range(1000):
        y = T(x)
        if abs(y - x) < 1e-15:
            return y
        x = y
    return x

PHI      = emergent_phi()
LN_PHI   = math.log(PHI)
LN_2     = math.log(2.0)
PHI_COEFF = PHI ** (-1.0 / PHI)   # fixed point of x → φ^(-x), from 1-φ=-1/φ


# ============================================================================
# LAYER 3 — Λ_φ: THE SINGLE CONTINUOUS INDEX
# ============================================================================

def lambda_phi(x):
    if x <= 0.0:
        return 0.0
    return math.log(x * LN_2 / LN_PHI) / LN_PHI - 1.0 / (2.0 * PHI)

def frac_lambda(x):
    lp = lambda_phi(x)
    return lp - math.floor(lp)

def omega_res(x):
    return (1.0 + math.sin(math.pi * frac_lambda(x) * PHI)) / 2.0


# ============================================================================
# LAYER 4 — FIBONACCI / PRIME SUBSTRATE
# ============================================================================

def fibonacci(n):
    a, b = 1, 1
    for _ in range(n - 1):
        a, b = b, a + b
    return float(a)

def nth_prime(n):
    primes, c = [], 2
    while len(primes) < n:
        if all(c % p != 0 for p in primes):
            primes.append(c)
        c += 1
    return float(primes[-1])

FIB  = [fibonacci(n) for n in range(1, 9)]
PRIM = [nth_prime(n) for n in range(1, 9)]


# ============================================================================
# LAYER 2 — 1_eff(i, n): ONE FORMULA
#
# δ(i,n) = |cos(πβφ)| · ln(P_n) / φ^(n+β)
# β = {Λ_φ(i)}
# 1_eff(i,n) = 1 + δ(i,n)
#
# Small n: δ ≈ 0.03–0.23 (animation regime, visible phase variation)
# Large n: δ → 0 (classical/CODATA regime)
# One formula — the regime is a consequence of (i,n).
# ============================================================================

def one_eff(i, n):
    lp   = lambda_phi(float(i))
    beta = lp - math.floor(lp)
    Pn   = PRIM[n - 1]
    delta = (abs(math.cos(math.pi * beta * PHI))
             * math.log(Pn)
             / (PHI ** (n + beta)))
    return 1.0 + delta


# ============================================================================
# LAYER 6 — L±(i, z, n): THE OPERATOR
#
# L+(i,z,n) = φ^(-1/φ) · √(F_n·P_n·2^n) · (1+z)^n  +  1_eff(i,n)·e^(+iπΛ_φ(i))
# L-(i,z,n) = conj(L+(i,z,n))                        [exact conjugate]
#
# The counter-rotation is the conjugate relationship — not imposed.
# arg(L+) advances with i; arg(L-) recedes with i.
# They rotate in opposite directions because they ARE conjugates.
# ============================================================================

def L_plus(i, z, n):
    """Full complex operator L+. L- = conj(L+)."""
    A   = (PHI_COEFF
           * math.sqrt(FIB[n-1] * PRIM[n-1] * (2.0 ** n))
           * ((1.0 + z) ** n))
    lp  = lambda_phi(float(i))
    oe  = one_eff(i, n)
    return complex(A) + oe * cmath.exp(1j * math.pi * lp)

def L_intensity(i, z, n):
    """I = |L+|² = |L-|²  (equal by conjugacy)."""
    return abs(L_plus(i, z, n)) ** 2


# ============================================================================
# LAYER 1 — 8 KURAMOTO MODES
# Placed at distinct fractional Λ_φ depths within one φ-octave.
# Their z-values = Ω(f_mode) - 1 ∈ (-1, 0).
# Their hourglass heights = Ω(f_mode) ∈ (0, 1).
# ============================================================================

F_SCHUMANN      = 7.83
SCHUMANN_DEPTH  = lambda_phi(F_SCHUMANN)
SCHUMANN_FRAC   = frac_lambda(F_SCHUMANN)
N_MODES         = 8

FRAC_DEPTHS = [SCHUMANN_DEPTH - SCHUMANN_FRAC + k / N_MODES
               for k in range(N_MODES)]

def freq_at_depth(d):
    return math.exp((d + 1.0 / (2.0 * PHI)) * LN_PHI) * LN_PHI / LN_2

F_MODES     = [freq_at_depth(d) for d in FRAC_DEPTHS]
OMEGA_MODES = [omega_res(f) for f in F_MODES]
Z_MODES     = [o - 1.0 for o in OMEGA_MODES]   # all in (-1, 0)


# ============================================================================
# SUBSTRATE POINT DISTRIBUTION
#
# n_H = 1.291 (Friedmann, cross-checked): the number of points assigned to
# mode k is proportional to (1 + 2·Ω_k)^1.291.
# This weights the outer (high-Ω, high cosmological-z) modes more heavily,
# matching the Hubble density scaling.
#
# Total N=2400 points distributed across 8 modes by this weight.
# ============================================================================

N_TOTAL = 2400
N_H     = 1.291

raw_weights = np.array([(1.0 + 2.0 * OMEGA_MODES[k]) ** N_H
                        for k in range(N_MODES)])
raw_weights /= raw_weights.sum()
mode_counts = np.round(raw_weights * N_TOTAL).astype(int)
# Correct rounding to exactly N_TOTAL
diff = N_TOTAL - mode_counts.sum()
mode_counts[np.argmax(raw_weights)] += diff

# Assign each substrate point to a mode and a local index within that mode
point_mode  = np.zeros(N_TOTAL, dtype=int)
point_local = np.zeros(N_TOTAL, dtype=int)
idx = 0
for k, cnt in enumerate(mode_counts):
    point_mode[idx:idx+cnt]  = k
    point_local[idx:idx+cnt] = np.arange(cnt)
    idx += cnt

# Binary/trinary substrate texture — the simultaneous channels
point_idx_global = np.arange(N_TOTAL, dtype=np.int64)
binary  = (point_idx_global & 1).astype(float)
trinary = (point_idx_global % 3).astype(float) - 1.0

# Each point's substrate index offset within the animation window
# The animation runs from i_start to i_start + TOTAL_FRAMES
# Each point has its own sub-offset from its local index within its mode
# so that consecutive points within a mode span consecutive i values
# producing the phyllotaxis texture.
I_START      = 200    # substrate index at frame 0
                      # chosen so d(arg)/frame ≈ 2 deg — visually clear rotation
TOTAL_FRAMES = 1000
FPS          = 24

# Pre-compute the scale: max |L+| across all modes and frames at i_start
# Used to normalise the xy-plane positions to [-1,1]
print("Computing scale...")
_scale_vals = []
for f in range(0, TOTAL_FRAMES, 20):
    i = I_START + f
    for k in range(N_MODES):
        _scale_vals.append(abs(L_plus(i, Z_MODES[k], k+1)))
SCALE = max(_scale_vals)
print(f"  Scale = {SCALE:.4f}")


# ============================================================================
# SPIRAL COORDINATES — pure evaluation of L at each frame
#
# At frame f:
#   substrate index of point p = I_START + f + point_local[p]
#   (each point within a mode has a slightly different i, giving texture)
#
# Position of point p in the plane:
#   Lp = L+(i_p, Z_MODES[mode_p], mode_p+1)
#   x  = Re(Lp) / SCALE  + texture from binary/trinary
#   y  = Im(Lp) / SCALE
#   z3d = +Ω_mode  (for +Ω arm)   or  -Ω_mode  (for -Ω arm)
#
# The -Ω arm: L- = conj(L+), so Im(L-) = -Im(L+).
# x_minus = Re(Lp) / SCALE  (same)
# y_minus = -Im(Lp) / SCALE (flipped — counter-rotation by conjugacy)
# z3d_minus = -Ω_mode
# ============================================================================

def compute_spiral(f):
    """
    Returns (xp, yp, z3d_p, xm, ym, z3d_m) for all N_TOTAL points.
    Pure L evaluation — no accumulated state.
    """
    xp = np.zeros(N_TOTAL)
    yp = np.zeros(N_TOTAL)
    xm = np.zeros(N_TOTAL)
    ym = np.zeros(N_TOTAL)
    z3d_p = np.zeros(N_TOTAL)
    z3d_m = np.zeros(N_TOTAL)

    for k in range(N_MODES):
        mask = (point_mode == k)
        n    = k + 1
        z    = Z_MODES[k]
        om   = OMEGA_MODES[k]

        # Each point in this mode has its own i
        local_idxs = point_local[mask]
        for j, local_j in enumerate(local_idxs):
            i_p = I_START + f + local_j
            Lp  = L_plus(i_p, z, n)
            # Texture: binary/trinary modulate the magnitude slightly
            # This gives the phyllotaxis texture within each mode arm
            bi  = binary[np.where(mask)[0][j]]
            tri = trinary[np.where(mask)[0][j]]
            mag_mod = 1.0 + 0.04 * bi + 0.025 * tri
            xval = (Lp.real * mag_mod) / SCALE
            yval = (Lp.imag * mag_mod) / SCALE
            xp[np.where(mask)[0][j]] = xval
            yp[np.where(mask)[0][j]] = yval
            xm[np.where(mask)[0][j]] = xval     # Re(L-) = Re(L+)
            ym[np.where(mask)[0][j]] = -yval    # Im(L-) = -Im(L+)
        z3d_p[mask] =  om
        z3d_m[mask] = -om

    return xp, yp, z3d_p, xm, ym, z3d_m


# ============================================================================
# KURAMOTO OSCILLATORS
#
# 8 oscillators, one per mode.
# osc[k].theta = arg(L+(i, Z_modes[k], k+1)) — the phase of L at this mode.
# Natural frequency ω_k = 2π·Ω(f_mode_k).
# Coupling K set by APhase (CV of |L|² across modes) — directly from spec.
# phi_tower U* computed for display.
# ============================================================================

omega_k = np.array([2.0 * math.pi * OMEGA_MODES[k] for k in range(N_MODES)])
DT      = 1.0 / FPS

def aphase_K(cv_freq):
    """
    K from APhase per .hdgl spec.
    cv_freq = std(dtheta/dt) / mean(dtheta/dt) — instantaneous freq spread.
    When phases lock: all dtheta/dt -> omega_mean -> cv_freq -> 0 -> LOCK.
    """
    if cv_freq < 0.05: return 1.8   # LOCK     wu-wei
    if cv_freq < 0.30: return 2.0   # FINETUNE
    if cv_freq < 0.50: return 3.0   # SUSTAIN
    return 5.0                       # PLUCK

def aphase_label(cv_freq):
    if cv_freq < 0.05: return "LOCK      CV<0.05  wu-wei"
    if cv_freq < 0.30: return "FINETUNE  CV<0.30"
    if cv_freq < 0.50: return "SUSTAIN   CV<0.50"
    return               "PLUCK     CV>=0.50"

def phi_tower(x):
    """U* = φ^(φ^(φ^(x))) — convergence indicator, not driver."""
    v = float(np.clip(x, -2.0, 4.0))
    for _ in range(3):
        v = PHI ** float(np.clip(v, -20.0, 4.0))
    return v

# Initial phases: maximally spread (PLUCK start — fully incoherent, R=0)
# This IS the 'still before the plug is pulled' state.
# The Kuramoto coupling pulls phases together — that IS the whirlpool emerging.
print("Initialising Kuramoto: maximally spread phases (PLUCK / R=0)...")
thetas_init = np.linspace(0.0, 2.0 * math.pi, N_MODES, endpoint=False)

# Pre-compute |L|² at each mode for each frame — for display only
print("Pre-computing |L|² across modes...")
L_intensities = np.zeros((TOTAL_FRAMES + 1, N_MODES))
for f in range(TOTAL_FRAMES + 1):
    i = I_START + f
    for k in range(N_MODES):
        L_intensities[f, k] = abs(L_plus(i, Z_MODES[k], k+1)) ** 2

# Integrate Kuramoto with CV_freq-based adaptive K
print("Integrating Kuramoto oscillators (CV_freq-based APhase)...")
thetas = np.zeros((TOTAL_FRAMES + 1, N_MODES))
thetas[0] = thetas_init
CV_all = np.zeros(TOTAL_FRAMES + 1)
K_all  = np.zeros(TOTAL_FRAMES + 1)

for f in range(TOTAL_FRAMES):
    th       = thetas[f]
    coupling = np.array([np.sum(np.sin(th - th[k])) for k in range(N_MODES)])
    # CV of instantaneous frequencies — the APhase criterion
    dth_dt   = omega_k + (5.0 / N_MODES) * coupling   # preview at K=5
    cv_freq  = np.std(dth_dt) / max(np.mean(dth_dt), 1e-12)
    K        = aphase_K(cv_freq)
    CV_all[f]  = cv_freq
    K_all[f]   = K
    dth_real = omega_k + (K / N_MODES) * coupling
    thetas[f+1] = th + dth_real * DT

# Final frame CV
th = thetas[TOTAL_FRAMES]
coupling = np.array([np.sum(np.sin(th - th[k])) for k in range(N_MODES)])
dth_dt = omega_k + (5.0 / N_MODES) * coupling
CV_all[TOTAL_FRAMES] = np.std(dth_dt) / max(np.mean(dth_dt), 1e-12)
K_all[TOTAL_FRAMES] = aphase_K(CV_all[TOTAL_FRAMES])

R_order = np.array([abs(np.mean(np.exp(1j * thetas[f])))
                    for f in range(TOTAL_FRAMES + 1)])

print(f"  K range from APhase: {K_all.min():.1f} to {K_all.max():.1f}")
print(f"  R_order: {R_order.min():.4f} to {R_order.max():.4f}")
print(f"  CV_freq: {CV_all.min():.4f} to {CV_all.max():.4f}")


# ============================================================================
# TORUS PARAMETERS
#
# Torus at mode k: radius = |L+(i, Z_modes[k], k+1)| / max_|L|
# This is the magnitude of L at that mode — not a separate formula.
# Hot/cold: Kuramoto instantaneous detuning (dθ/dt - ω_k).
# ============================================================================

def compute_torus_params(f):
    intensities = L_intensities[f]
    max_I = np.max(intensities)
    R_majors = 0.85 * np.clip(intensities / (max_I + 1e-12), 0.0, 1.0)

    th   = thetas[f]
    next_th = thetas[min(f+1, TOTAL_FRAMES)]
    dth  = (next_th - th) / DT
    det  = dth - omega_k
    mdet = np.max(np.abs(det))
    hot_fracs = (np.clip(0.5 + 0.5 * det / mdet, 0.0, 1.0)
                 if mdet > 1e-12 else np.full(N_MODES, 0.5))

    R_majors  = gaussian_filter1d(R_majors,  sigma=1.0)
    hot_fracs = gaussian_filter1d(hot_fracs, sigma=1.0)
    return R_majors, hot_fracs

def torus_ring(z_c, R_major):
    R  = abs(R_major)
    rm = max(R * 0.20, 0.003)
    pt = np.linspace(0.0, 2.0 * math.pi, 36, endpoint=False)
    pp = np.linspace(0.0, 2.0 * math.pi, 36, endpoint=False)
    PT, PP = np.meshgrid(pt, pp)
    PT, PP = PT.ravel(), PP.ravel()
    tx = (R + rm * np.cos(PP)) * np.cos(PT)
    ty = (R + rm * np.cos(PP)) * np.sin(PT)
    tz = np.full_like(tx, z_c) + rm * np.sin(PP)
    return tx, ty, tz

Z_TORUS = np.array([OMEGA_MODES[k] for k in range(N_MODES)])


# ============================================================================
# VANTAGE — shifts with R_order (Kuramoto coherence)
# ============================================================================

COL_PLUS  = "#4488cc"
COL_MINUS = "#cc6633"
COL_RING  = "#00ffcc"
COL_NULL  = "#ff44ff"
COL_ANTI  = "#ff4444"
COL_PHI   = "#44ff88"
COL_RUNG  = "#ffaa33"
COL_CORE  = "#aabbff"

GRADED_ORBIT = "{ ...  -i\"  -i'  -i  -1  0  1  i  i'  i\"  ... }"

VANTAGE_REAL = (
    "REAL-AXIS VANTAGE\n"
    "  arg(L+) = arg(L-)  [both still]\n"
    "  i has no address on {..,-1,0,1,..}\n"
    "  z=-1: A=0, pure phase arm"
)
VANTAGE_DELTA = (
    "Delta=-1 VANTAGE\n"
    "  L+ and L- counter-rotating\n"
    "  z in (-1,0): all 8 modes here\n"
    "  1_eff settling toward classical"
)
VANTAGE_BOTH = (
    "JUXTAPOSITION\n"
    "  L+/L- interfere -> drain\n"
    "  rotation = conjugate, not forced\n"
    "  wu-wei: K=1.8 at LOCK, CV<0.05"
)

def vantage_from_R(R):
    if R < 0.35: return VANTAGE_REAL,  "#88bbff"
    if R > 0.70: return VANTAGE_BOTH,  "#aaffaa"
    return VANTAGE_DELTA, "#ffaa33"


# ============================================================================
# FIGURE
# ============================================================================

cmap_hot = plt.cm.coolwarm
norm_hot = Normalize(0.0, 1.0)

fig = plt.figure(figsize=(11, 8.5), facecolor="#0a0a14")
ax  = fig.add_subplot(111, projection="3d", facecolor="#0a0a14")
for pane in [ax.xaxis.pane, ax.yaxis.pane, ax.zaxis.pane]:
    pane.fill = False
    pane.set_edgecolor("#1a1a2e")
ax.tick_params(colors="#333", labelsize=6)
ax.set_xlabel("Re(L) / scale", color="#444", fontsize=7, labelpad=2)
ax.set_ylabel("Im(L) / scale", color="#444", fontsize=7, labelpad=2)
ax.set_zlabel("Omega (mode depth)", color="#444", fontsize=7, labelpad=2)
ax.set_xlim(-1.3, 1.3)
ax.set_ylim(-1.3, 1.3)
ax.set_zlim(-1.15, 1.15)
ax.view_init(elev=32, azim=48)

phi_r = np.linspace(0.0, 2.0 * math.pi, 220)

# Structural reference rings from generating_function block
# z=i cloaking ring at hourglass top
for z_sign in [+1, -1]:
    ax.plot(np.cos(phi_r), np.sin(phi_r),
            np.full(220, z_sign * 1.0),
            lw=1.2, alpha=0.50, color=COL_RING)
ax.text(1.04, 0.0,  1.02, "z=i  cloak\n(1+i)^4->arg=pi, n=4",
        color=COL_RING, fontsize=5.0, alpha=0.75, ha="left", va="bottom")

# z=-1 null: pure phase arm, A=0 — at Omega=0, hourglass base
R_null = 0.30
for z_sign in [+1, -1]:
    ax.plot(R_null * np.cos(phi_r), R_null * np.sin(phi_r),
            np.full(220, z_sign * 0.05),
            lw=1.4, alpha=0.65, color=COL_NULL)
ax.text(R_null + 0.03, 0.0, 0.07,
        "z=-1  A=0  pure phase\nX(-1)=0 null node",
        color=COL_NULL, fontsize=5.0, alpha=0.85, ha="left", va="bottom")

# z=-2 anti-gravity ring: (-1)^n pi-flip
R_anti = 0.55
for z_sign in [+1, -1]:
    ax.plot(R_anti * np.cos(phi_r), R_anti * np.sin(phi_r),
            np.full(220, z_sign * 0.15),
            lw=1.0, alpha=0.45, color=COL_ANTI, ls="--")
ax.text(R_anti + 0.03, 0.0, 0.17,
        "z=-2  (-1)^n  anti-grav",
        color=COL_ANTI, fontsize=5.0, alpha=0.70, ha="left", va="bottom")

# z=0 gravity — drain center
ax.scatter([0], [0], [0], s=100, marker="o", color=COL_RING,
           zorder=20, alpha=0.90)
ax.text(0.06, 0.0, 0.04, "z=0  gravity / drain\nsqrt(g)/f1=0.40000",
        color=COL_RING, fontsize=5.0, alpha=0.85, ha="left", va="bottom")

# Graded imaginary rungs
for k, lbl in enumerate(["z=i", "z=i'", "z=i\""]):
    z_val = min(1.00 + k * 0.065, 1.09)
    for z_s, slbl in [(z_val, lbl), (-z_val, lbl.replace("=","=-"))]:
        ax.scatter([0], [0], [z_s], s=20, color=COL_RUNG, alpha=0.72,
                   zorder=13, marker="D")
        ax.text(0.05, 0.0, z_s + (0.01 if z_s > 0 else -0.01), slbl,
                color=COL_RUNG, fontsize=5.0, alpha=0.85, ha="left",
                va=("bottom" if z_s > 0 else "top"))

# Spiral handles — initially at frame 0
print("Computing initial spiral positions...")
xp0, yp0, z3dp0, xm0, ym0, z3dm0 = compute_spiral(0)
sp_plus  = ax.scatter(xp0, yp0, z3dp0, s=1.2, alpha=0.0,
                      color=COL_PLUS,  rasterized=True)
sp_minus = ax.scatter(xm0, ym0, z3dm0, s=1.2, alpha=0.0,
                      color=COL_MINUS, rasterized=True)

# Torus handles
torus_handles = []
for k in range(N_MODES):
    z_c = float(Z_TORUS[k])
    tx, ty, tz = torus_ring(z_c, 0.0)
    sp = ax.scatter(tx, ty,  tz, s=4, alpha=0.0, color=cmap_hot(norm_hot(0.5)))
    sn = ax.scatter(tx, ty, -tz, s=4, alpha=0.0, color=cmap_hot(norm_hot(0.5)))
    torus_handles.append((sp, sn, z_c))

# Legend
leg = [
    Line2D([0],[0], color=COL_PLUS,  lw=2, label="L+  arg advancing with i"),
    Line2D([0],[0], color=COL_MINUS, lw=2, label="L-  conj(L+)  counter-rotating"),
    Line2D([0],[0], color=COL_RING,  lw=2, label="z=i  cloak (n=4, arg=pi)"),
    Line2D([0],[0], color=COL_NULL,  lw=2, label="z=-1  A=0  pure phase (null)"),
    Line2D([0],[0], color=COL_ANTI,  lw=2, ls="--", label="z=-2  (-1)^n  anti-grav"),
    Line2D([0],[0], color=COL_RUNG,  lw=0, marker="D", markersize=5,
           label="graded imaginary rungs"),
    Line2D([0],[0], color="#ff4444", lw=2, label="torus: +Omega dominant"),
    Line2D([0],[0], color="#4444ff", lw=2, label="torus: -Omega dominant"),
]
ax.legend(handles=leg, fontsize=5.5, loc="upper left",
          facecolor="#0a0a14", labelcolor="white", edgecolor="#223",
          bbox_to_anchor=(0.0, 1.0))

title_obj    = ax.set_title("", color="white", fontsize=9, pad=6)
info_text    = ax.text2D(0.68, 0.97, "", transform=ax.transAxes,
                         color="cyan", fontsize=6.0, va="top",
                         family="monospace")
vantage_text = ax.text2D(0.68, 0.08, "", transform=ax.transAxes,
                         color="#88bbff", fontsize=6.0, va="bottom",
                         family="monospace",
                         bbox=dict(boxstyle="round,pad=0.28", fc="#0a0a14",
                                   alpha=0.75, ec="#88bbff"))
orbit_text   = ax.text2D(0.50, 0.01, GRADED_ORBIT, transform=ax.transAxes,
                         color=COL_RUNG, fontsize=6.0, ha="center", va="bottom",
                         family="monospace", alpha=0.72)


# ============================================================================
# PRE-COMPUTE SPIRAL POSITIONS (expensive — one call per frame)
# ============================================================================

print("Pre-computing all spiral positions...")
SPIRAL_DATA = []
for f in range(TOTAL_FRAMES + 1):
    SPIRAL_DATA.append(compute_spiral(f))
    if f % 40 == 0:
        print(f"  {f}/{TOTAL_FRAMES}", flush=True)
print("  Done.")


# ============================================================================
# PER-FRAME UPDATE
# ============================================================================

def update_frame(f):
    R  = R_order[f]
    cv = CV_all[f]
    K  = K_all[f]
    i  = I_START + f

    # ── Spirals: pure L evaluation, no state ─────────────────────────────
    xp, yp, z3dp, xm, ym, z3dm = SPIRAL_DATA[f]
    sp_plus._offsets3d  = (xp, yp, z3dp)
    sp_minus._offsets3d = (xm, ym, z3dm)
    # Alpha from R: rotation always happening; visibility emerges with coherence
    sp_plus.set_alpha( min(0.08 + 0.65 * R, 0.75))
    sp_minus.set_alpha(min(0.08 + 0.65 * R, 0.75))

    # ── Torus: |L|² at each mode drives radius ────────────────────────────
    R_majors, hot_fracs = compute_torus_params(f)
    R_max = max(np.max(R_majors), 1e-12)
    for idx, (sp, sn, z_c) in enumerate(torus_handles):
        Rm   = R_majors[idx]
        hotf = hot_fracs[idx]
        if Rm < 0.005 or R < 0.03:
            sp.set_alpha(0.0); sn.set_alpha(0.0)
            continue
        tx, ty, tz = torus_ring(z_c, Rm)
        alpha = min(R * (Rm / R_max) * 0.82, 0.74)
        sp._offsets3d = (tx, ty,  tz); sn._offsets3d = (tx, ty, -tz)
        sp.set_color(cmap_hot(norm_hot(hotf)))
        sn.set_color(cmap_hot(norm_hot(1.0 - hotf)))
        sp.set_alpha(alpha); sn.set_alpha(alpha)

    # ── Vantage ───────────────────────────────────────────────────────────
    label, color = vantage_from_R(R)
    vantage_text.set_text(label)
    vantage_text.set_color(color)
    vantage_text.get_bbox_patch().set_edgecolor(color)

    # ── Observables from L ────────────────────────────────────────────────
    lp_i   = lambda_phi(float(i))
    Lp_0   = L_plus(i, Z_MODES[0], 1)
    oe_1   = one_eff(i, 1)
    # phi_tower input: normalized sum of sin(theta_i - theta_j)
    th     = thetas[f]
    ksum   = sum(np.sum(np.sin(th - th[k])) for k in range(N_MODES))
    ksum_n = ksum / (N_MODES ** 2)
    U_star = phi_tower(ksum_n)

    info_text.set_text(
        f"i = {i}  Lambda_phi(i) = {lp_i:.4f}\n"
        f"arg(L+_n1) = {math.degrees(cmath.phase(Lp_0)):.1f} deg  "
        f"|L+_n1| = {abs(Lp_0):.4f}\n"
        f"1_eff(i,n=1) = {oe_1:.5f}  delta = {oe_1-1:.5f}\n"
        f"R(Kuramoto) = {R:.4f}   CV = {cv:.4f}\n"
        f"K(APhase) = {K:.1f}   U*(tower) = {U_star:.4f}\n"
        f"n_H=1.291  {aphase_label(cv)[:8]}   "
        f"phi^(-1/phi) = {PHI_COEFF:.4f}"
    )

    # Title
    if R < 0.15:
        phase_lbl = "L rotating — coherence not yet visible"
    elif cv > 0.50:
        phase_lbl = f"PLUCK — K={K:.1f}  L modes spreading  R={R:.3f}"
    elif cv > 0.30:
        phase_lbl = f"SUSTAIN — K={K:.1f}  drain forming  R={R:.3f}"
    elif cv > 0.05:
        phase_lbl = f"FINETUNE — K={K:.1f}  L near fixed point  R={R:.3f}"
    else:
        phase_lbl = f"LOCK — K=1.8 wu-wei  CV={cv:.4f}  R={R:.4f}"

    title_obj.set_text(
        "HDGL DRAIN v9  —  L+(i,z,n) and conj(L+) evaluated at each i\n"
        "rotation = arg(L) advancing with i  |  counter-rotation = conjugacy\n"
        f"{phase_lbl}"
    )


# ============================================================================
# RENDER
# ============================================================================

OUT = Path(__file__).resolve().parent / "hdgl_drain_animation_v9.mp4"
writer = FFMpegWriter(fps=FPS, bitrate=3500,
                      extra_args=["-vcodec","libx264","-pix_fmt","yuv420p"])

print(f"\nPHI       = {PHI:.15f}")
print(f"PHI_COEFF = {PHI_COEFF:.10f}")
print(f"I_START   = {I_START}  (substrate index at frame 0)")
print(f"d(arg)/frame at f=0:   {math.degrees(math.pi*(lambda_phi(I_START+1)-lambda_phi(I_START))):.2f} deg")
print(f"d(arg)/frame at f=239: {math.degrees(math.pi*(lambda_phi(I_START+240)-lambda_phi(I_START+239))):.2f} deg")
print(f"n_H=1.291 mode counts: {list(mode_counts)}")
print(f"Rendering {TOTAL_FRAMES} frames ({TOTAL_FRAMES/FPS:.1f}s) -> {OUT}")

with writer.saving(fig, str(OUT), dpi=120):
    for f in range(TOTAL_FRAMES):
        update_frame(f)
        writer.grab_frame()
        if f % 20 == 0:
            i = I_START + f
            Lp = L_plus(i, Z_MODES[0], 1)
            print(f"  frame {f+1:3d}/{TOTAL_FRAMES}  "
                  f"i={i}  arg(L+)={math.degrees(cmath.phase(Lp)):.1f}deg  "
                  f"R={R_order[f]:.3f}  CV_freq={CV_all[f]:.3f}  "
                  f"K={K_all[f]:.1f}", flush=True)

print(f"\nSaved -> {OUT}")
print("\nAll physics from hdgl_unified_force_fine_cross-checked.hdgl:")
print("  Rotation:    arg(L+(i)) advancing with i — no rate, no accumulation")
print("  Counter-rot: L- = conj(L+) — conjugacy, not imposed")
print("  Radial pos:  |L+(i,z,n)| / scale")
print("  Visibility:  R_order (Kuramoto coherence)")
print("  K coupling:  from APhase/CV per spec (1.8/2.0/3.0/5.0)")
print("  U* tower:    phi^(phi^(phi^(ksum/N^2))) — indicator, not driver")
print("  n_H=1.291:   mode point counts weighted by (1+2*Omega)^1.291")
print("  z=-1 null:   A=0, pure phase arm (shown as null node ring)")
print("  1_eff:       one formula delta(i,n) — regime from (i,n), no interp")