#!/usr/bin/env python3
"""
===============================================================================
HDGL COUNTER-ROTATING PHYLLOTAXIS — DRAIN EFFECT
===============================================================================

Two simultaneous phyllotactic fields:

    θ+ = +2π i Ω
    θ- = -2π i Ω

They counter-rotate while sharing one reciprocal radial closure:

    T(X) = 1 + 1/X

The drain is produced by reciprocal compression toward the common origin:

    d(r) = 1 / T(r)

    ρ(r) = r * d(r)

The binary/trinary substrate remains simultaneous:

    B_i ∈ {0,1}
    τ_i ∈ {-1,0,+1}

No stored φ is used.
Ω emerges from:

    Ω_(n+1) = T(Ω_n)
            = 1 + 1/Ω_n

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

from pathlib import Path
import math

import numpy as np

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt


# ============================================================================
# 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)


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

def T(x):
    """
    Reciprocal HDGL transformation:

        T(X) = 1 + 1/X
    """
    return 1.0 + 1.0 / x


def omega_orbit(seed=OMEGA_SEED, steps=OMEGA_STEPS):
    """
    Generate the emergent Ω orbit.
    """

    values = np.empty(steps + 1, dtype=np.float64)
    values[0] = seed

    for i in range(steps):
        values[i + 1] = T(values[i])

    return values


def emergent_omega():
    """
    Generate Ω dynamically from the reciprocal closure.
    """

    x = OMEGA_SEED

    for _ in range(1000):

        y = T(x)

        if abs(y - x) < 1e-15:
            break

        x = y

    return x


OMEGA = emergent_omega()


# ============================================================================
# SIMULTANEOUS BINARY / TRINARY 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


# ============================================================================
# COMMON RADIAL SUBSTRATE
# ============================================================================

def build_radial_substrate(n=N):

    i, binary, trinary = substrate_states(n)

    # ------------------------------------------------------------------------
    # Base outward organization.
    # ------------------------------------------------------------------------

    r_base = np.sqrt(i + 1.0)

    # ------------------------------------------------------------------------
    # Simultaneous binary/trinary modulation.
    #
    # Both channels act on the SAME substrate.
    # ------------------------------------------------------------------------

    modulation = (
        1.0
        + 0.075 * binary
        + 0.050 * trinary
    )

    r = r_base * modulation

    # ------------------------------------------------------------------------
    # Reciprocal drain.
    #
    # T(r) = 1 + 1/r
    #
    # Reciprocal factor:
    #
    # d(r) = 1/T(r)
    #      = r/(r+1)
    #
    # Therefore:
    #
    # rho = r*d
    #     = r²/(r+1)
    #
    # This remains finite and preserves the radial ordering while introducing
    # reciprocal closure.
    # ------------------------------------------------------------------------

    reciprocal = T(r)

    drain_factor = 1.0 / reciprocal

    rho = r * drain_factor

    return (
        i,
        binary,
        trinary,
        r,
        reciprocal,
        drain_factor,
        rho,
    )


# ============================================================================
# COUNTER-ROTATING PHYLLOTAXIS
# ============================================================================

def build_counter_rotating_phyllotaxis(n=N):

    (
        i,
        binary,
        trinary,
        r,
        reciprocal,
        drain_factor,
        rho,
    ) = build_radial_substrate(n)

    # ------------------------------------------------------------------------
    # Two simultaneous angular fields.
    # ------------------------------------------------------------------------

    theta_plus = 2.0 * math.pi * i * OMEGA
    theta_minus = -2.0 * math.pi * i * OMEGA

    # ------------------------------------------------------------------------
    # Shared binary/trinary phase perturbation.
    #
    # The perturbation is applied symmetrically so that the two branches
    # remain counter-rotating.
    # ------------------------------------------------------------------------

    phase_modulation = (
        0.075 * binary
        + 0.050 * trinary
    )

    theta_plus += phase_modulation
    theta_minus -= phase_modulation

    # ------------------------------------------------------------------------
    # Counter-rotating coordinates.
    # ------------------------------------------------------------------------

    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,
    )


# ============================================================================
# DRAIN STRENGTH
# ============================================================================

def drain_profile(n=N):

    (
        i,
        binary,
        trinary,
        r,
        reciprocal,
        drain_factor,
        rho,
    ) = build_radial_substrate(n)

    return i, r, reciprocal, drain_factor, rho


# ============================================================================
# PLOT 1 — COUNTER-ROTATING PHYLLOTAXIS
# ============================================================================

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))

    # ------------------------------------------------------------------------
    # Positive rotation.
    # ------------------------------------------------------------------------

    ax.scatter(
        x_plus,
        y_plus,
        s=4,
        alpha=0.45,
        linewidths=0,
        label="counter-rotation +Ω",
    )

    # ------------------------------------------------------------------------
    # Negative rotation.
    # ------------------------------------------------------------------------

    ax.scatter(
        x_minus,
        y_minus,
        s=4,
        alpha=0.45,
        linewidths=0,
        label="counter-rotation −Ω",
    )

    # ------------------------------------------------------------------------
    # Drain origin.
    # ------------------------------------------------------------------------

    ax.scatter(
        [0],
        [0],
        s=70,
        marker="o",
        label="drain",
    )

    ax.set_aspect("equal", adjustable="box")

    ax.set_title(
        "HDGL COUNTER-ROTATING PHYLLOTAXIS\n"
        "Simultaneous ±Ω with Reciprocal Drain"
    )

    ax.set_xlabel("counter-rotating substrate")
    ax.set_ylabel("counter-rotating substrate")

    ax.grid(True, alpha=0.20)
    ax.legend()

    fig.tight_layout()

    path = OUT_DIR / "hdgl_counter_rotating_phyllotaxis.png"

    fig.savefig(
        path,
        dpi=180,
        bbox_inches="tight",
    )

    plt.close(fig)

    return path


# ============================================================================
# PLOT 2 — DRAIN PROFILE
# ============================================================================

def plot_drain_profile():

    i, r, reciprocal, drain_factor, rho = drain_profile()

    fig, ax = plt.subplots(figsize=(11, 7))

    ax.plot(
        i,
        r,
        linewidth=1.2,
        label="outward radius r",
    )

    ax.plot(
        i,
        rho,
        linewidth=1.2,
        label="drained radius ρ",
    )

    ax.set_title(
        "HDGL RECIPROCAL DRAIN\n"
        "ρ = r / T(r)"
    )

    ax.set_xlabel("substrate index i")
    ax.set_ylabel("radial coordinate")

    ax.grid(True, alpha=0.20)
    ax.legend()

    fig.tight_layout()

    path = OUT_DIR / "hdgl_drain_profile.png"

    fig.savefig(
        path,
        dpi=180,
        bbox_inches="tight",
    )

    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=(11, 7))

    ax.plot(
        i,
        theta_plus,
        linewidth=1.0,
        label="+Ω",
    )

    ax.plot(
        i,
        theta_minus,
        linewidth=1.0,
        label="−Ω",
    )

    ax.axhline(
        0.0,
        linestyle="--",
        linewidth=1.0,
    )

    ax.set_title(
        "HDGL COUNTER-ROTATING PHASE\n"
        "θ+ = +2πiΩ   /   θ− = −2πiΩ"
    )

    ax.set_xlabel("substrate index i")
    ax.set_ylabel("phase")

    ax.grid(True, alpha=0.20)
    ax.legend()

    fig.tight_layout()

    path = OUT_DIR / "hdgl_counter_rotation_phase.png"

    fig.savefig(
        path,
        dpi=180,
        bbox_inches="tight",
    )

    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=(11, 7))

    ax.plot(
        n,
        values,
        linewidth=1.5,
        marker="o",
        markersize=3,
    )

    ax.axhline(
        OMEGA,
        linestyle="--",
        linewidth=1.0,
    )

    ax.set_title(
        "HDGL Ω RECURSION\n"
        "Ωₙ₊₁ = T(Ωₙ) = 1 + 1/Ωₙ"
    )

    ax.set_xlabel("iteration n")
    ax.set_ylabel("Ωₙ")

    ax.grid(True, alpha=0.20)

    ax.text(
        0.98,
        0.05,
        f"Emergent Ω ≈ {OMEGA:.15f}",
        transform=ax.transAxes,
        ha="right",
        va="bottom",
    )

    fig.tight_layout()

    path = OUT_DIR / "hdgl_omega_orbit.png"

    fig.savefig(
        path,
        dpi=180,
        bbox_inches="tight",
    )

    plt.close(fig)

    return path


# ============================================================================
# PLOT 5 — 3D DRAIN
# ============================================================================

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()

    # ------------------------------------------------------------------------
    # Use normalized radial coordinates to make the drain visible.
    # ------------------------------------------------------------------------

    scale = np.max(rho)

    rp = rho / scale

    # The two planar phyllotaxes are folded into opposing z-gradients.
    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=(12, 10))

    ax = fig.add_subplot(
        111,
        projection="3d",
    )

    # Positive branch.

    ax.scatter(
        xp,
        yp,
        z_plus,
        s=2,
        alpha=0.35,
        label="+Ω",
    )

    # Negative branch.

    ax.scatter(
        xm,
        ym,
        z_minus,
        s=2,
        alpha=0.35,
        label="−Ω",
    )

    # Central drain.

    ax.scatter(
        [0],
        [0],
        [0],
        s=80,
        marker="o",
        label="DRAIN",
    )

    ax.set_title(
        "HDGL COUNTER-ROTATING DRAIN\n"
        "±Ω → reciprocal closure"
    )

    ax.set_xlabel("+Ω branch")
    ax.set_ylabel("−Ω branch")
    ax.set_zlabel("radial closure")

    ax.legend()

    fig.tight_layout()

    path = OUT_DIR / "hdgl_3d_drain.png"

    fig.savefig(
        path,
        dpi=180,
        bbox_inches="tight",
    )

    plt.close(fig)

    return path


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

def print_report(paths):

    print()
    print("=" * 79)
    print("HDGL COUNTER-ROTATING PHYLLOTAXIS — DRAIN EFFECT")
    print("=" * 79)
    print()

    print("Primitive:")
    print("    Ω_(n+1) = T(Ω_n)")
    print("    T(X)    = 1 + 1/X")
    print()

    print("Emergent Ω:")
    print(f"    Ω ≈ {OMEGA:.15f}")
    print()

    print("Fixed-point residual:")
    print(
        f"    Ω² - Ω - 1 ≈ "
        f"{OMEGA * OMEGA - OMEGA - 1.0:.6e}"
    )
    print()

    print("Simultaneous substrate:")
    print("    B_i ∈ {0,1}")
    print("    τ_i ∈ {-1,0,+1}")
    print()

    print("Counter-rotation:")
    print("    θ+ = +2π i Ω")
    print("    θ- = -2π i Ω")
    print()

    print("Angular cancellation:")
    print("    Ω + (-Ω) = 0")
    print()

    print("Reciprocal drain:")
    print("    T(r) = 1 + 1/r")
    print("    d(r) = 1/T(r)")
    print("    ρ(r) = r/T(r)")
    print()

    print("Geometry:")
    print("    PHYLLOTAXIS+")
    print("           ↘")
    print("            DRAIN")
    print("           ↗")
    print("    PHYLLOTAXIS−")
    print()

    print("Output:")
    print(f"    {OUT_DIR}")
    print()

    for path in paths:
        print(f"    [OK] {path}")

    print()
    print("=" * 79)
    print("COMPLETE")
    print("=" * 79)
    print()


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

def main():

    paths = []

    paths.append(
        plot_counter_rotating_phyllotaxis()
    )

    paths.append(
        plot_drain_profile()
    )

    paths.append(
        plot_counter_rotation_phase()
    )

    paths.append(
        plot_omega_orbit()
    )

    paths.append(
        plot_3d_drain()
    )

    print_report(paths)


if __name__ == "__main__":
    main()