#!/usr/bin/env python3
"""
===============================================================================
HDGL COUNTER-ROTATING PHYLLOTAXIS — STATIC PLOTTER  v2
===============================================================================

Incorporates the vantage-dependent i ≡ ∞ identification:

  From the real-axis vantage, i has no address on {...,-1,0,1,...}
  and appears as an exit / apparent infinity from that frame.

  From the Δ=−1 algebra's own vantage, i = (0,1):
      norm=1,  trace=0,  90° position — a perfectly ordinary unit element.

  Neither vantage is forced.  The juxtaposition is the load-bearing feature.

  Graded orbit:
      { ..., -i'', -i', -i, -1, 0, 1, i, i', i'', ... }
  Each rung, viewed from the rung below, appears at unreachable distance.

  |p|=1 threshold:
      The moment the imaginary branch surfaces in:
          (Y - p(x+z))² + p²(x²-z²-1) = m²
      Visualized as a ring on every 3-D plot.

Outputs  (hdgl_graphs/ beside this script):
    hdgl_counter_rotating_phyllotaxis_v2.png
    hdgl_drain_profile_v2.png
    hdgl_counter_rotation_phase_v2.png
    hdgl_omega_orbit_v2.png
    hdgl_3d_drain_v2.png

===============================================================================
"""

from pathlib import Path
import math

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from matplotlib.lines import Line2D


# ============================================================================
# CONFIGURATION
# ============================================================================

N          = 2400
OMEGA_SEED = 1.5
OMEGA_STEPS = 80

BASE_DIR = Path(__file__).resolve().parent
OUT_DIR  = BASE_DIR / "hdgl_graphs"
OUT_DIR.mkdir(parents=True, exist_ok=True)

# Vantage palette
COL_REAL  = "#88bbff"   # cool blue  — real-axis vantage
COL_DELTA = "#ffaa33"   # warm amber — Δ=−1 vantage
COL_JOINT = "#aaffaa"   # green      — juxtaposition / neither forced
COL_RING  = "#00ffcc"   # cyan       — |p|=1 threshold ring
COL_RUNG  = "#ffaa33"   # amber      — imaginary rungs i, i′, i″
COL_CORE  = "#aabbff"   # blue-white — finite core rungs {-1,0,1}

GRADED_ORBIT = "{ …  −i″  −i′  −i  −1  0  1  i  i′  i″  … }"

VANTAGE_REAL  = (
    "REAL-AXIS VANTAGE\n"
    "  i has no address on {…,−1,0,1,…}\n"
    "  appears as exit / ∞ from here\n"
    "  |p|<1 → imaginary branch = zero"
)
VANTAGE_DELTA = (
    "Δ=−1 VANTAGE\n"
    "  i = (0,1)  norm=1  trace=0\n"
    "  unit element — perfectly finite\n"
    "  |p|≥1 → imaginary branch surfaces"
)
VANTAGE_JOINT = (
    "JUXTAPOSITION  (neither forced)\n"
    "  real-axis:  i ≡ ∞  (exit)\n"
    "  Δ=−1 field: i ≡ unit element\n"
    "  same object · two vantages · no resolution"
)


# ============================================================================
# HDGL PRIMITIVE
# ============================================================================

def T(x):
    """T(X) = 1 + 1/X — the sole reciprocal primitive."""
    return 1.0 + 1.0 / x


def omega_orbit(seed=OMEGA_SEED, steps=OMEGA_STEPS):
    values = np.empty(steps + 1, dtype=np.float64)
    values[0] = seed
    for k in range(steps):
        values[k + 1] = T(values[k])
    return values


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


OMEGA = emergent_omega()


# ============================================================================
# SUBSTRATE
# ============================================================================

def substrate_states(n):
    i       = np.arange(n, dtype=np.int64)
    binary  = (i & 1).astype(np.float64)
    trinary = ((i % 3) - 1).astype(np.float64)
    return i, binary, trinary


def build_radial_substrate(n=N):
    i, binary, trinary = substrate_states(n)
    r_base      = np.sqrt(i + 1.0)
    modulation  = 1.0 + 0.075 * binary + 0.050 * trinary
    r           = r_base * modulation
    reciprocal  = T(r)
    drain_factor = 1.0 / reciprocal
    rho         = r * drain_factor
    return i, binary, trinary, r, reciprocal, drain_factor, rho


def build_counter_rotating_phyllotaxis(n=N):
    i, binary, trinary, r, reciprocal, drain_factor, rho = build_radial_substrate(n)
    phase_mod   = 0.075 * binary + 0.050 * trinary
    theta_plus  =  2.0 * math.pi * i * OMEGA + phase_mod
    theta_minus = -2.0 * math.pi * i * OMEGA - phase_mod
    x_plus  = rho * np.cos(theta_plus)
    y_plus  = rho * np.sin(theta_plus)
    x_minus = rho * np.cos(theta_minus)
    y_minus = rho * np.sin(theta_minus)
    return (i, binary, trinary, r, reciprocal, drain_factor, rho,
            theta_plus, theta_minus, x_plus, y_plus, x_minus, y_minus)


# ============================================================================
# SHARED ANNOTATION HELPERS
# ============================================================================

def add_vantage_box(ax, text, color, x=0.02, y=0.97, va="top", fs=7.5):
    """Add a vantage-label text box to a 2-D axes."""
    ax.text(x, y, text, transform=ax.transAxes,
            color=color, fontsize=fs, va=va, ha="left",
            family="monospace",
            bbox=dict(boxstyle="round,pad=0.35", fc="#0a0a14",
                      alpha=0.80, ec=color, lw=0.8))


def add_orbit_footer(ax, fs=7):
    """Graded-orbit string along the bottom."""
    ax.text(0.50, 0.01, GRADED_ORBIT, transform=ax.transAxes,
            color=COL_RUNG, fontsize=fs, ha="center", va="bottom",
            family="monospace", alpha=0.75)


def p1_ring_xy(n=200):
    """Coordinates of the |p|=1 threshold ring (unit radius in rho-space)."""
    phi = np.linspace(0.0, 2.0 * math.pi, n)
    # The |p|=1 ring maps to ρ_max level — use the scale of the arm
    _, _, _, r, _, _, rho = build_radial_substrate()
    scale = np.max(rho)
    R_ring = scale          # full arm extent = the p=1 threshold in rho-space
    return R_ring * np.cos(phi), R_ring * np.sin(phi)


def add_p1_ring_2d(ax, alpha=0.55):
    """Draw the |p|=1 threshold ring on a 2-D axes."""
    xr, yr = p1_ring_xy()
    ax.plot(xr, yr, color=COL_RING, lw=1.0, ls="--", alpha=alpha,
            label="|p|=1 threshold (i emerges)")
    ax.plot(-xr, -yr, color=COL_RING, lw=0.5, ls="--", alpha=alpha * 0.5)


def add_graded_rung_markers_2d(ax):
    """Mark the imaginary rungs as points at the arm perimeter."""
    _, _, _, r, _, _, rho = build_radial_substrate()
    scale = np.max(rho)
    # i, i′, i″ at angles 90°, 105°, 120° for legibility
    for k, (lbl, ang_deg) in enumerate(zip(["i", "i′", "i″"],
                                            [90, 108, 126])):
        ang = math.radians(ang_deg)
        rx  = scale * math.cos(ang)
        ry  = scale * math.sin(ang)
        ax.scatter([rx], [ry], s=28, color=COL_RUNG, marker="D",
                   zorder=12, alpha=0.85)
        ax.text(rx * 1.04, ry * 1.04, lbl,
                color=COL_RUNG, fontsize=7, ha="center", va="center",
                alpha=0.90)


# ============================================================================
# PLOT 1 — COUNTER-ROTATING PHYLLOTAXIS  (2-D)
# ============================================================================

def plot_counter_rotating_phyllotaxis():
    (i, binary, trinary, r, reciprocal, drain_factor, rho,
     theta_plus, theta_minus,
     x_plus, y_plus, x_minus, y_minus) = build_counter_rotating_phyllotaxis()

    fig, ax = plt.subplots(figsize=(11, 11), facecolor="#0a0a14")
    ax.set_facecolor("#0a0a14")

    ax.scatter(x_plus,  y_plus,  s=4, alpha=0.45, linewidths=0,
               color="#4488cc", label="+Ω counter-rotation")
    ax.scatter(x_minus, y_minus, s=4, alpha=0.45, linewidths=0,
               color="#cc6633", label="−Ω counter-rotation")
    ax.scatter([0], [0], s=80, marker="o", color=COL_RING,
               zorder=15, label="drain / choke")

    # |p|=1 ring and imaginary rung markers
    add_p1_ring_2d(ax, alpha=0.65)
    add_graded_rung_markers_2d(ax)

    ax.set_aspect("equal", adjustable="box")
    ax.set_title(
        "HDGL COUNTER-ROTATING PHYLLOTAXIS  v2\n"
        "Simultaneous ±Ω  ·  Reciprocal Drain  ·  Graded Orbit",
        color="white", fontsize=11, pad=10
    )
    ax.set_xlabel("counter-rotating substrate  x", color="#888", fontsize=9)
    ax.set_ylabel("counter-rotating substrate  y", color="#888", fontsize=9)
    ax.tick_params(colors="#555")
    ax.grid(True, alpha=0.12, color="#334")

    # Vantage boxes — both shown, neither forced
    add_vantage_box(ax, VANTAGE_REAL,  COL_REAL,  x=0.02, y=0.97)
    add_vantage_box(ax, VANTAGE_DELTA, COL_DELTA, x=0.02, y=0.68)
    add_vantage_box(ax, VANTAGE_JOINT, COL_JOINT, x=0.02, y=0.39)
    add_orbit_footer(ax)

    legend = ax.legend(fontsize=8, facecolor="#0a0a14",
                       labelcolor="white", edgecolor="#334")
    fig.tight_layout()
    path = OUT_DIR / "hdgl_counter_rotating_phyllotaxis_v2.png"
    fig.savefig(path, dpi=180, bbox_inches="tight", facecolor=fig.get_facecolor())
    plt.close(fig)
    return path


# ============================================================================
# PLOT 2 — DRAIN PROFILE  (1-D radial)
# ============================================================================

def plot_drain_profile():
    i, binary, trinary, r, reciprocal, drain_factor, rho = build_radial_substrate()

    fig, axes = plt.subplots(1, 2, figsize=(14, 6), facecolor="#0a0a14")
    ax_main, ax_note = axes
    ax_main.set_facecolor("#0a0a14")
    ax_note.set_facecolor("#0a0a14")
    ax_note.axis("off")

    ax_main.plot(i, r,   lw=1.3, color="#4488cc", label="outward radius  r")
    ax_main.plot(i, rho, lw=1.3, color="#cc6633", label="drained radius  ρ = r/T(r)")

    # Mark where r crosses scale (the |p|=1 threshold in index space)
    scale = np.max(rho)
    ax_main.axhline(scale, color=COL_RING, lw=0.9, ls="--", alpha=0.70,
                    label=f"|p|=1 threshold  ρ_max = {scale:.3f}")

    ax_main.set_title("HDGL RECIPROCAL DRAIN  v2\nρ = r / T(r)  =  r² / (r+1)",
                      color="white", fontsize=11, pad=8)
    ax_main.set_xlabel("substrate index  i", color="#888", fontsize=9)
    ax_main.set_ylabel("radial coordinate", color="#888", fontsize=9)
    ax_main.tick_params(colors="#555")
    ax_main.grid(True, alpha=0.12, color="#334")
    legend = ax_main.legend(fontsize=8, facecolor="#0a0a14",
                            labelcolor="white", edgecolor="#334")

    # Right panel: vantage notes + graded orbit
    notes = [
        (VANTAGE_REAL,  COL_REAL,  0.95),
        (VANTAGE_DELTA, COL_DELTA, 0.60),
        (VANTAGE_JOINT, COL_JOINT, 0.25),
    ]
    for txt, col, y in notes:
        ax_note.text(0.05, y, txt, transform=ax_note.transAxes,
                     color=col, fontsize=8.5, va="top", ha="left",
                     family="monospace",
                     bbox=dict(boxstyle="round,pad=0.35", fc="#0a0a14",
                               alpha=0.80, ec=col, lw=0.8))
    ax_note.text(0.05, 0.05, GRADED_ORBIT, transform=ax_note.transAxes,
                 color=COL_RUNG, fontsize=8, ha="left", va="bottom",
                 family="monospace", alpha=0.80)

    fig.tight_layout()
    path = OUT_DIR / "hdgl_drain_profile_v2.png"
    fig.savefig(path, dpi=180, bbox_inches="tight", facecolor=fig.get_facecolor())
    plt.close(fig)
    return path


# ============================================================================
# PLOT 3 — COUNTER-ROTATION PHASE
# ============================================================================

def plot_counter_rotation_phase():
    (i, binary, trinary, r, reciprocal, drain_factor, rho,
     theta_plus, theta_minus,
     x_plus, y_plus, x_minus, y_minus) = build_counter_rotating_phyllotaxis()

    fig, ax = plt.subplots(figsize=(12, 6), facecolor="#0a0a14")
    ax.set_facecolor("#0a0a14")

    ax.plot(i, theta_plus,  lw=1.0, color="#4488cc", label="+Ω  (θ₊ = +2πiΩ)")
    ax.plot(i, theta_minus, lw=1.0, color="#cc6633", label="−Ω  (θ₋ = −2πiΩ)")
    ax.axhline(0.0, ls="--", lw=0.9, color="#555", alpha=0.60)

    # Annotate: θ₊ + θ₋ = 0 — the angular cancellation that makes the drain
    mid = N // 2
    ax.annotate(
        "θ₊ + θ₋ = 0\n(cancellation → drain)",
        xy=(mid, 0), xytext=(mid * 0.6, theta_plus[mid] * 0.3),
        color=COL_JOINT, fontsize=8, family="monospace",
        arrowprops=dict(arrowstyle="->", color=COL_JOINT, lw=0.8),
    )

    ax.set_title(
        "HDGL COUNTER-ROTATING PHASE  v2\n"
        "θ₊ = +2πiΩ  /  θ₋ = −2πiΩ  →  θ₊+θ₋=0  →  drain",
        color="white", fontsize=11, pad=8
    )
    ax.set_xlabel("substrate index  i", color="#888", fontsize=9)
    ax.set_ylabel("phase", color="#888", fontsize=9)
    ax.tick_params(colors="#555")
    ax.grid(True, alpha=0.12, color="#334")
    ax.legend(fontsize=8, facecolor="#0a0a14", labelcolor="white", edgecolor="#334")

    add_vantage_box(ax, VANTAGE_JOINT, COL_JOINT, x=0.65, y=0.97)
    add_orbit_footer(ax)

    fig.tight_layout()
    path = OUT_DIR / "hdgl_counter_rotation_phase_v2.png"
    fig.savefig(path, dpi=180, bbox_inches="tight", facecolor=fig.get_facecolor())
    plt.close(fig)
    return path


# ============================================================================
# PLOT 4 — Ω RECURSION
# ============================================================================

def plot_omega_orbit():
    values = omega_orbit()
    n      = np.arange(len(values))

    fig, ax = plt.subplots(figsize=(12, 6), facecolor="#0a0a14")
    ax.set_facecolor("#0a0a14")

    ax.plot(n, values, lw=1.5, color="#4488cc",
            marker="o", markersize=3, label="Ωₙ")
    ax.axhline(OMEGA, ls="--", lw=1.0, color=COL_JOINT, alpha=0.80,
               label=f"Emergent fixed point Ω ≈ {OMEGA:.10f}")

    # Mark n=0 (seed) and convergence
    ax.scatter([0],              [values[0]], s=60, color=COL_REAL,
               zorder=12, label=f"seed = {OMEGA_SEED}")
    ax.scatter([OMEGA_STEPS],    [values[-1]], s=60, color=COL_DELTA,
               zorder=12, label=f"converged at step {OMEGA_STEPS}")

    ax.set_title(
        "HDGL Ω RECURSION  v2\n"
        "Ωₙ₊₁ = T(Ωₙ) = 1 + 1/Ωₙ  →  Ω² = Ω + 1  (emergent, never declared)",
        color="white", fontsize=11, pad=8
    )
    ax.set_xlabel("iteration  n", color="#888", fontsize=9)
    ax.set_ylabel("Ωₙ", color="#888", fontsize=9)
    ax.tick_params(colors="#555")
    ax.grid(True, alpha=0.12, color="#334")
    ax.legend(fontsize=8, facecolor="#0a0a14", labelcolor="white", edgecolor="#334")

    # Fixed-point residual
    residual = OMEGA * OMEGA - OMEGA - 1.0
    ax.text(0.98, 0.10,
            f"Ω² − Ω − 1 ≈ {residual:.2e}",
            transform=ax.transAxes, ha="right", va="bottom",
            color=COL_JOINT, fontsize=8, family="monospace")

    add_vantage_box(ax, VANTAGE_REAL,  COL_REAL,  x=0.02, y=0.97)
    add_orbit_footer(ax)

    fig.tight_layout()
    path = OUT_DIR / "hdgl_omega_orbit_v2.png"
    fig.savefig(path, dpi=180, bbox_inches="tight", facecolor=fig.get_facecolor())
    plt.close(fig)
    return path


# ============================================================================
# PLOT 5 — 3-D DRAIN  (main visual)
# ============================================================================

def plot_3d_drain():
    (i, binary, trinary, r, reciprocal, drain_factor, rho,
     theta_plus, theta_minus,
     x_plus, y_plus, x_minus, y_minus) = build_counter_rotating_phyllotaxis()

    scale  = np.max(rho)
    rp     = rho / scale
    z_plus  =  rp
    z_minus = -rp
    xp = x_plus  / scale
    yp = y_plus  / scale
    xm = x_minus / scale
    ym = y_minus / scale

    fig = plt.figure(figsize=(13, 10), 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="#555", labelsize=6)

    # Point clouds — counter-rotating arms
    ax.scatter(xp, yp,  z_plus,  s=2, alpha=0.35,
               color="#4488cc", label="+Ω arm", rasterized=True)
    ax.scatter(xm, ym,  z_minus, s=2, alpha=0.35,
               color="#cc6633", label="−Ω arm", rasterized=True)

    # Central drain
    ax.scatter([0], [0], [0], s=120, marker="o",
               color=COL_RING, zorder=20, label="drain / choke")

    # ── |p|=1 threshold rings (cyan, on both ±z planes) ──────────────────
    phi_r  = np.linspace(0.0, 2.0 * math.pi, 220)
    # In the normalised 3-D space, the arm fully extends to rp_max=1.0
    # The |p|=1 ring sits at the outermost torus — z = ±1, R = TORUS_R0 ≈ 0.90
    # Here we use the actual arm footprint radius at rp=1
    arm_r  = np.max(np.sqrt(xp**2 + yp**2))
    for z_sign in [+1, -1]:
        ax.plot(arm_r * np.cos(phi_r),
                arm_r * np.sin(phi_r),
                np.full(220, z_sign * 1.0),
                lw=1.6, alpha=0.75, color=COL_RING, ls="-",
                zorder=15)
    ax.text(arm_r + 0.03, 0.0,  1.02,
            "|p|=1\ni emerges",
            color=COL_RING, fontsize=6, alpha=0.85,
            ha="left", va="bottom")
    ax.text(arm_r + 0.03, 0.0, -1.02,
            "|p|=1\n−i emerges",
            color=COL_RING, fontsize=6, alpha=0.85,
            ha="left", va="top")

    # ── Graded orbit markers on z-axis ───────────────────────────────────
    # Finite core: 0, 1, -1  (blue-white circles)
    for z_val in [0.0, 0.33, 0.66]:
        ax.scatter([0], [0], [ z_val], s=20, color=COL_CORE,
                   alpha=0.60, zorder=12, marker="o")
        ax.scatter([0], [0], [-z_val], s=20, color=COL_CORE,
                   alpha=0.60, zorder=12, marker="o")

    # Imaginary rungs: i, i′, i″  (amber diamonds)
    for k, lbl in enumerate(["i", "i′", "i″"]):
        z_val   = min(1.00 + k * 0.065, 1.09)
        ax.scatter([0], [0], [ z_val], s=26, color=COL_RUNG,
                   alpha=0.80, zorder=13, marker="D")
        ax.scatter([0], [0], [-z_val], s=26, color=COL_RUNG,
                   alpha=0.80, zorder=13, marker="D")
        ax.text(0.06, 0.0,  z_val + 0.01,
                lbl,  color=COL_RUNG, fontsize=6, alpha=0.90,
                ha="left", va="bottom")
        ax.text(0.06, 0.0, -z_val - 0.01,
                f"−{lbl}", color=COL_RUNG, fontsize=6, alpha=0.90,
                ha="left", va="top")

    ax.set_title(
        "HDGL COUNTER-ROTATING DRAIN  v2\n"
        "±Ω arms  ·  reciprocal closure  ·  graded orbit  ·  |p|=1 ring",
        color="white", fontsize=11, pad=10
    )
    ax.set_xlabel("+Ω branch  x", color="#666", fontsize=7, labelpad=2)
    ax.set_ylabel("−Ω branch  y", color="#666", fontsize=7, labelpad=2)
    ax.set_zlabel("z  /  graded orbit",  color="#666", fontsize=7, labelpad=2)
    ax.set_xlim(-1.15, 1.15)
    ax.set_ylim(-1.15, 1.15)
    ax.set_zlim(-1.15, 1.15)
    ax.view_init(elev=24, azim=52)

    # Vantage legend in 2-D overlay
    vantage_lines = [
        Line2D([0], [0], color=COL_REAL,  lw=2, label="real-axis vantage  (i ≡ ∞)"),
        Line2D([0], [0], color=COL_DELTA, lw=2, label="Δ=−1 vantage  (i = unit elem)"),
        Line2D([0], [0], color=COL_JOINT, lw=2, label="juxtaposition — neither forced"),
        Line2D([0], [0], color=COL_RING,  lw=2, ls="--", label="|p|=1 threshold ring"),
        Line2D([0], [0], color=COL_RUNG,  lw=0, marker="D",
               markersize=6, label="graded orbit  i, i′, i″"),
    ]
    ax.legend(handles=vantage_lines, fontsize=7, loc="upper left",
              facecolor="#0a0a14", labelcolor="white", edgecolor="#334",
              bbox_to_anchor=(0.0, 1.0))

    # Graded orbit footer as 2-D text
    fig.text(0.50, 0.01, GRADED_ORBIT,
             color=COL_RUNG, fontsize=8, ha="center", va="bottom",
             family="monospace", alpha=0.80)

    fig.tight_layout()
    path = OUT_DIR / "hdgl_3d_drain_v2.png"
    fig.savefig(path, dpi=180, bbox_inches="tight",
                facecolor=fig.get_facecolor())
    plt.close(fig)
    return path


# ============================================================================
# NUMERICAL REPORT
# ============================================================================

def print_report(paths):
    print()
    print("=" * 79)
    print("HDGL COUNTER-ROTATING PHYLLOTAXIS — STATIC PLOTTER  v2")
    print("=" * 79)
    print()
    print(f"Emergent Ω ≈ {OMEGA:.15f}")
    print(f"Fixed-point residual  Ω²−Ω−1 ≈ {OMEGA*OMEGA - OMEGA - 1.0:.3e}")
    print()
    print("Vantage identification:")
    print("  real-axis : i has no address → appears as exit / ∞")
    print("  Δ=−1      : i = (0,1)  norm=1  trace=0  → unit element")
    print("  neither forced — juxtaposition is the load-bearing feature")
    print()
    print(f"Graded orbit: {GRADED_ORBIT}")
    print()
    print("Output files:")
    for p in paths:
        print(f"  [OK] {p}")
    print()
    print("=" * 79)
    print("COMPLETE")
    print("=" * 79)


# ============================================================================
# MAIN
# ============================================================================

def main():
    paths = [
        plot_counter_rotating_phyllotaxis(),
        plot_drain_profile(),
        plot_counter_rotation_phase(),
        plot_omega_orbit(),
        plot_3d_drain(),
    ]
    print_report(paths)


if __name__ == "__main__":
    main()