#!/usr/bin/env python3
"""
A_N Tridiagonal Eigenvalue Analysis for Stellar Rotation

Computes the eigenvalues of the path-graph adjacency matrix A_N for small N,
demonstrates that φ = (1+√5)/2 emerges exactly at N=4, and lays out the
eigenvalue ratios √2, φ, φ/√2 that arise from the first three nontrivial
cases.

This structure is forced by the Recognition Composition Law (d'Alembert
functional equation) restricted to a finite chain of N sites.

References:
  - McQuillan et al. 2014 (ApJS 211, 24)  — 34,030 Kepler rotation periods
  - Dungee et al. 2022 (ApJ 938, 118)     — M67 K dwarfs from K2
  - Lu et al. 2024 (AJ 167, 159)          — 94,627 TESS rotation periods
  - Schou et al. 1998 (ApJ 505, 390)      — helioseismic tachocline rotation
"""

import math

PHI = (1 + math.sqrt(5)) / 2
SQRT2 = math.sqrt(2)
SQRT3 = math.sqrt(3)

# ═══════════════════════════════════════════════════════════════════════
# PART 1: A_N EIGENVALUES FROM FIRST PRINCIPLES
# ═══════════════════════════════════════════════════════════════════════

print("=" * 72)
print("PART 1: A_N EIGENVALUE STRUCTURE")
print("Eigenvalues of the N×N path-graph adjacency matrix")
print("λ_k = 2cos(kπ/(N+1)),  k = 1, ..., N")
print("=" * 72)

for N in range(1, 9):
    eigs = [2 * math.cos(k * math.pi / (N + 1)) for k in range(1, N + 1)]
    print(f"\n  N = {N}:")
    for k, e in enumerate(eigs, 1):
        # Identify algebraic form
        label = ""
        if abs(e - PHI) < 1e-10: label = "= φ"
        elif abs(e + PHI) < 1e-10: label = "= -φ"
        elif abs(e - 1/PHI) < 1e-10: label = "= 1/φ"
        elif abs(e + 1/PHI) < 1e-10: label = "= -1/φ"
        elif abs(e - SQRT2) < 1e-10: label = "= √2"
        elif abs(e + SQRT2) < 1e-10: label = "= -√2"
        elif abs(e - SQRT3) < 1e-10: label = "= √3"
        elif abs(e + SQRT3) < 1e-10: label = "= -√3"
        elif abs(e - 1) < 1e-10: label = "= 1"
        elif abs(e + 1) < 1e-10: label = "= -1"
        elif abs(e) < 1e-10: label = "= 0"
        print(f"    λ_{k} = {e:+.10f}  {label}")

# ═══════════════════════════════════════════════════════════════════════
# PART 2: THE GOLDEN RATIO EMERGES AT N=4
# ═══════════════════════════════════════════════════════════════════════

print("\n" + "=" * 72)
print("PART 2: WHY φ APPEARS AT N = 4")
print("=" * 72)

print("""
  The characteristic polynomial of A_4 is:
    det(A_4 - λI) = λ⁴ - 3λ² + 1

  Setting x = λ², this factors as:
    x² - 3x + 1 = 0
    x = (3 ± √5)/2

  So λ² = (3+√5)/2 = φ²  or  λ² = (3-√5)/2 = 1/φ²

  Therefore λ = ±φ  or  λ = ±1/φ

  The golden ratio is NOT assumed. It is the algebraic solution
  of the A_4 characteristic equation. The same second-order
  recursion is the discrete form of the d'Alembert/RCL functional
  equation restricted to a 4-site chain (a modeling step).
""")

# Verify
print("  Verification:")
print(f"    φ        = {PHI:.10f}")
print(f"    2cos(π/5)= {2*math.cos(math.pi/5):.10f}  ← λ₁ of A₄")
print(f"    Match: {abs(PHI - 2*math.cos(math.pi/5)) < 1e-14}")
print(f"    1/φ      = {1/PHI:.10f}")
print(f"    2cos(2π/5)={2*math.cos(2*math.pi/5):.10f}  ← λ₂ of A₄")
print(f"    Match: {abs(1/PHI - 2*math.cos(2*math.pi/5)) < 1e-14}")

# ═══════════════════════════════════════════════════════════════════════
# PART 3: THE THREE CRITICAL EIGENVALUE RATIOS
# ═══════════════════════════════════════════════════════════════════════

print("\n" + "=" * 72)
print("PART 3: EIGENVALUE RATIOS — √2, φ, φ/√2")
print("=" * 72)

lambda_max = {}
for N in range(2, 7):
    lambda_max[N] = 2 * math.cos(math.pi / (N + 1))

print("\n  Dominant eigenvalues λ₁(A_N):")
for N in range(2, 7):
    print(f"    N={N}: λ₁ = {lambda_max[N]:.10f}")

print("\n  Key ratios (zero free parameters):")

r_32 = lambda_max[3] / lambda_max[2]
r_42 = lambda_max[4] / lambda_max[2]
r_43 = lambda_max[4] / lambda_max[3]

print(f"    λ₁(A₃)/λ₁(A₂) = {r_32:.10f}  =  √2  = {SQRT2:.10f}  "
      f"(error: {abs(r_32-SQRT2):.1e})")
print(f"    λ₁(A₄)/λ₁(A₂) = {r_42:.10f}  =  φ   = {PHI:.10f}  "
      f"(error: {abs(r_42-PHI):.1e})")
print(f"    λ₁(A₄)/λ₁(A₃) = {r_43:.10f}  =  φ/√2= {PHI/SQRT2:.10f}  "
      f"(error: {abs(r_43-PHI/SQRT2):.1e})")

print(f"\n  Note: φ/√2 = {PHI/SQRT2:.10f} ≈ 1.1441")
print(f"        φ²   = {PHI**2:.10f} ≈ 2.6180")

# ═══════════════════════════════════════════════════════════════════════
# PART 4: PHYSICAL PREDICTION — STELLAR ROTATION PERIODS
# ═══════════════════════════════════════════════════════════════════════

print("\n" + "=" * 72)
print("PART 4: STELLAR ROTATION PERIOD PREDICTIONS")
print("=" * 72)

print("""
  Model: A star with N coupled, differentially rotating convective shells
  has its rotation dynamics governed by the A_N tridiagonal matrix.

  The dominant eigenmode (k=1) sets the long-period convergence edge.
  If the period scales as P ∝ 1/λ₁(A_N) (eigenvalue sets the coupling
  strength, stronger coupling → faster equilibration → shorter period):

  Assignment (from stellar structure):
    N=2: M dwarfs  (0.35-0.55 M☉) — 2 convective shells
    N=3: K dwarfs  (0.55-0.85 M☉) — 3 convective shells
    N=4: G dwarfs  (0.85-1.10 M☉) — 4 convective shells (Sun-like)
""")

# If we normalize so that N=4 (Sun-like) gives the solar rotation period
P_sun_sidereal = 25.38  # Carrington sidereal period, days
P_tachocline = 26.9     # interior/tachocline rotation, approximately

print("  Using solar equatorial sidereal period = 25.38 days as the N=4 anchor:")
print()

for N, label, mass_range in [
    (2, "M dwarfs", "0.35-0.55 M☉"),
    (3, "K dwarfs", "0.55-0.85 M☉"),
    (4, "G dwarfs (Sun)", "0.85-1.10 M☉"),
    (5, "F dwarfs", "1.10-1.40 M☉"),
    (6, "A dwarfs", "1.40-2.00 M☉"),
]:
    # P ∝ 1/λ₁  (stronger coupling = shorter period)
    P_predicted_inv = P_sun_sidereal * lambda_max[4] / lambda_max[N]
    # P ∝ λ₁  (dominant mode period scales with eigenvalue)
    P_predicted_dir = P_sun_sidereal * lambda_max[N] / lambda_max[4]
    print(f"  N={N} ({label}, {mass_range}):")
    print(f"    If P ∝ 1/λ₁: P = {P_predicted_inv:.2f} days  "
          f"(ratio to Sun = {lambda_max[4]/lambda_max[N]:.6f})")
    print(f"    If P ∝  λ₁:  P = {P_predicted_dir:.2f} days  "
          f"(ratio to Sun = {lambda_max[N]/lambda_max[4]:.6f})")
    print()

# ═══════════════════════════════════════════════════════════════════════
# PART 5: PERIOD RATIOS — THE ZERO-PARAMETER TEST
# ═══════════════════════════════════════════════════════════════════════

print("=" * 72)
print("PART 5: ZERO-PARAMETER PERIOD RATIO PREDICTIONS")
print("These ratios are exact algebraic numbers, independent of any")
print("anchor period or physical scale.")
print("=" * 72)

print("""
  If period P ∝ 1/λ₁(A_N):

    P(M dwarf, N=2) / P(K dwarf, N=3)  = λ₁(A₃)/λ₁(A₂) = √2    ≈ 1.4142
    P(M dwarf, N=2) / P(G dwarf, N=4)  = λ₁(A₄)/λ₁(A₂) = φ     ≈ 1.6180
    P(K dwarf, N=3) / P(G dwarf, N=4)  = λ₁(A₄)/λ₁(A₃) = φ/√2  ≈ 1.1441

  If period P ∝ λ₁(A_N):

    P(K dwarf, N=3) / P(M dwarf, N=2)  = λ₁(A₃)/λ₁(A₂) = √2    ≈ 1.4142
    P(G dwarf, N=4) / P(M dwarf, N=2)  = λ₁(A₄)/λ₁(A₂) = φ     ≈ 1.6180
    P(G dwarf, N=4) / P(K dwarf, N=3)  = λ₁(A₄)/λ₁(A₃) = φ/√2  ≈ 1.1441
""")

# ═══════════════════════════════════════════════════════════════════════
# PART 6: KNOWN OBSERVATIONAL BENCHMARKS
# ═══════════════════════════════════════════════════════════════════════

print("=" * 72)
print("PART 6: OBSERVATIONAL BENCHMARKS")
print("(Published values from the cited papers)")
print("=" * 72)

print("""
  Solar rotation:
    Equatorial sidereal:  25.38 days  (Carrington)
    Tachocline (interior): ~26.9 days  (Schou et al. 1998)
    Ratio surface/interior: {:.4f}

  McQuillan et al. 2014 (Kepler):
    34,030 rotation periods, 0.2-70 days
    Upper envelope (slow rotators) at ~4.5 Gyr:
      M dwarfs (0.4 M☉):  ~35-40 days (slow rotator convergence)
      K dwarfs (0.7 M☉):  ~28-32 days (slow rotator convergence)
      G dwarfs (1.0 M☉):  ~25-28 days (slow rotator convergence)

  Dungee et al. 2022 (M67 cluster, 4 Gyr):
    K dwarfs: ~26-30 days median period

  If the long-period convergence edge is:
    M dwarf (N=2): ~37 days
    K dwarf (N=3): ~28 days
    G dwarf (N=4, Sun): ~25.4 days

  Then observed ratios:
    P_M / P_K ≈ 37/28 ≈ 1.32   (predicted √2 ≈ 1.414)
    P_M / P_G ≈ 37/25.4 ≈ 1.46  (predicted φ  ≈ 1.618)
    P_K / P_G ≈ 28/25.4 ≈ 1.10  (predicted φ/√2 ≈ 1.144)
""".format(25.38/26.9))

# ═══════════════════════════════════════════════════════════════════════
# PART 7: SENSITIVITY ANALYSIS
# ═══════════════════════════════════════════════════════════════════════

print("=" * 72)
print("PART 7: WHAT PERIOD VALUES WOULD GIVE EXACT MATCHES?")
print("=" * 72)

# If we anchor at P_G = 25.38 days:
P_G = 25.38
P_K_pred = P_G * SQRT2 / (PHI / SQRT2)  # Wait, let me think about this more carefully

# If P ∝ 1/λ₁:
# P(N=2)/P(N=4) = λ₁(N=4)/λ₁(N=2) = φ
# So P(N=2) = φ × P(N=4)
P_M_pred = PHI * P_G
P_K_pred = (PHI / SQRT2) * P_G  # P(N=3)/P(N=4) = λ₁(A₄)/λ₁(A₃) = φ/√2

print(f"\n  Anchoring at P(Sun, N=4) = {P_G:.2f} days:")
print(f"    Predicted P(M dwarf, N=2) = φ × {P_G:.2f}     = {P_M_pred:.2f} days")
print(f"    Predicted P(K dwarf, N=3) = (φ/√2) × {P_G:.2f} = {P_K_pred:.2f} days")
print()
print(f"  Cross-check ratios:")
print(f"    P_M / P_K = {P_M_pred/P_K_pred:.10f}  should be √2 = {SQRT2:.10f}")
print(f"    P_M / P_G = {P_M_pred/P_G:.10f}  should be φ  = {PHI:.10f}")
print(f"    P_K / P_G = {P_K_pred/P_G:.10f}  should be φ/√2= {PHI/SQRT2:.10f}")

print(f"\n  Anchoring at P(tachocline) = 26.9 days:")
P_tach = 26.9
P_M_tach = PHI * P_tach
P_K_tach = (PHI / SQRT2) * P_tach
print(f"    Predicted P(M dwarf) = φ × {P_tach} = {P_M_tach:.2f} days")
print(f"    Predicted P(K dwarf) = (φ/√2) × {P_tach} = {P_K_tach:.2f} days")

# ═══════════════════════════════════════════════════════════════════════
# PART 8: THE RS FORCING CHAIN
# ═══════════════════════════════════════════════════════════════════════

print("\n\n" + "=" * 72)
print("PART 8: RECOGNITION SCIENCE FORCING CHAIN")
print("=" * 72)

print("""
  The complete derivation from RS axioms:

  1. A₁+A₂+A₃ (Recognition Composition Law)
     ↓
  2. J(xy) + J(x/y) = 2J(x)J(y) + 2J(x) + 2J(y)
     ↓  [log substitution x=eᵗ, y=eᵘ, H=G+1]
  3. H(t+u) + H(t-u) = 2H(t)H(u)   [d'Alembert]
     ↓  [restrict to integer steps on finite chain]
  4. v_{j-1} + v_{j+1} = λ v_j       [A_N eigenvalue equation]
     ↓  [solve characteristic polynomial]
  5. λ_k = 2cos(kπ/(N+1))            [A_N eigenvalues]
     ↓  [evaluate at N=4]
  6. λ₁ = 2cos(π/5) = φ              [golden ratio FORCED]
     ↓  [ratios between N=2,3,4]
  7. Period ratios = √2, φ, φ/√2      [zero free parameters]

  Steps 1-3 are proved in Lean 4 (0 sorry):
    - IndisputableMonolith/Cost/FunctionalEquation.lean
    - IndisputableMonolith/Cost/AczelTheorem.lean
    - IndisputableMonolith/Foundation/ContinuumLimit.lean

  Steps 4-7 are standard linear algebra (Toeplitz tridiagonal diagonalization).
""")

print("=" * 72)
print("ANALYSIS COMPLETE")
print("=" * 72)