#!/usr/bin/env python3
"""
Test whether solar p-mode oscillation data (BiSON) shows the tridiagonal
A_N matrix eigenvalue signature λ_k = 2cos(kπ/(N+1)).

Also tests for golden-ratio (φ) signatures in frequency ratios and spacings.

Data source: Davies et al. (2014), 22 years of BiSON observations.
"""

import math
import os

PHI = (1 + math.sqrt(5)) / 2  # 1.6180339887...
ALPHA_INV = 137.035999084

# ─── Load BiSON data ───────────────────────────────────────────────────────────

def load_bison(path):
    """Load BiSON frequency table: n, l, freq (μHz), error."""
    data = []
    with open(path) as f:
        for line in f:
            parts = line.split()
            if len(parts) >= 4:
                n, l, freq, err = int(parts[0]), int(parts[1]), float(parts[2]), float(parts[3])
                data.append((n, l, freq, err))
    return data

script_dir = os.path.dirname(os.path.abspath(__file__))
bison = load_bison(os.path.join(script_dir, "davies2014.txt"))

# ─── Organize by degree l ──────────────────────────────────────────────────────

modes = {}  # l -> sorted list of (n, freq)
for n, l, freq, err in bison:
    modes.setdefault(l, []).append((n, freq))
for l in modes:
    modes[l].sort()

print("=" * 72)
print("SOLAR P-MODE EIGENVALUE ANALYSIS")
print("Data: BiSON Davies et al. (2014), 22-year dataset")
print("=" * 72)

# ─── TEST 1: Large frequency separation (Δν) ─────────────────────────────────
# The large separation Δν_n,l = ν(n+1,l) - ν(n,l) is the fundamental
# spacing of the p-mode spectrum. For a homogeneous sphere, modes are
# equally spaced; deviations encode interior structure.

print("\n── TEST 1: Large Frequency Separation Δν ──")
print(f"{'l':>2}  {'n→n+1':>7}  {'Δν (μHz)':>10}  {'Δν/Δν_mean':>10}")
print("-" * 45)

for l in sorted(modes.keys()):
    freqs = modes[l]
    if len(freqs) < 2:
        continue
    separations = []
    for i in range(len(freqs) - 1):
        n1, f1 = freqs[i]
        n2, f2 = freqs[i + 1]
        if n2 == n1 + 1:
            separations.append((n1, n2, f2 - f1))
    if separations:
        mean_sep = sum(s[2] for s in separations) / len(separations)
        for n1, n2, dv in separations:
            print(f"{l:>2}  {n1:>2}→{n2:<2}   {dv:>10.3f}  {dv/mean_sep:>10.5f}")
        print(f"    mean Δν(l={l}) = {mean_sep:.3f} μHz")
        print()

# ─── TEST 2: Frequency ratios and φ ──────────────────────────────────────────
# Check whether consecutive frequency ratios approach φ or simple
# functions of φ.

print("\n── TEST 2: Frequency Ratios vs Golden Ratio φ ──")
print(f"φ = {PHI:.10f}")
print(f"1/φ = {1/PHI:.10f}")
print(f"φ² = {PHI**2:.10f}")
print()

for l in sorted(modes.keys()):
    freqs = modes[l]
    if len(freqs) < 2:
        continue
    print(f"  l = {l}:")
    print(f"  {'n₂/n₁':>8}  {'ν₂/ν₁':>10}  {'closest φ^p':>12}  {'p':>4}  {'error':>10}")
    print(f"  {'-'*50}")
    for i in range(len(freqs) - 1):
        n1, f1 = freqs[i]
        n2, f2 = freqs[i + 1]
        ratio = f2 / f1
        # Find closest φ^p for integer or half-integer p
        best_p, best_err = None, 1e10
        for p_num in range(-8, 9):
            p = p_num / 4.0
            val = PHI ** p
            if abs(ratio - val) < abs(best_err):
                best_p = p
                best_err = ratio - val
        print(f"  {n1:>2}/{n2:<2}     {ratio:>10.6f}  {PHI**best_p:>12.6f}  {best_p:>4.2f}  {best_err:>+10.6f}")
    print()

# ─── TEST 3: A_N eigenvalue fit ──────────────────────────────────────────────
# For a chain of N sites, eigenvalues are λ_k = 2cos(kπ/(N+1)).
# We normalize the observed frequencies and test whether they match
# the A_N pattern for some N.
#
# Strategy: use the l=0 radial modes. The "large separation" is roughly
# constant (~135 μHz), so the frequencies are approximately linear in n.
# But the DEVIATIONS from linearity (the "second differences") are what
# encode the A_N structure.

print("\n── TEST 3: Second Differences vs A_N Eigenvalue Pattern ──")
print("The second difference δ²ν = ν(n+1) - 2ν(n) + ν(n-1) probes")
print("the tridiagonal structure. For A_N: δ² is related to the")
print("discrete Laplacian eigenvalues.\n")

for l in sorted(modes.keys()):
    freqs_list = [f for _, f in modes[l]]
    ns = [n for n, _ in modes[l]]
    if len(freqs_list) < 3:
        continue
    print(f"  l = {l}: second differences δ²ν(n)")
    print(f"  {'n':>4}  {'ν(n) μHz':>12}  {'δ²ν (μHz)':>10}")
    print(f"  {'-'*32}")
    d2_vals = []
    for i in range(1, len(freqs_list) - 1):
        if ns[i+1] == ns[i] + 1 and ns[i] == ns[i-1] + 1:
            d2 = freqs_list[i + 1] - 2 * freqs_list[i] + freqs_list[i - 1]
            d2_vals.append((ns[i], freqs_list[i], d2))
            print(f"  {ns[i]:>4}  {freqs_list[i]:>12.3f}  {d2:>+10.3f}")
    print()

    if len(d2_vals) >= 2:
        # Test: do the second differences match 2cos(kπ/(N+1)) pattern?
        N = len(d2_vals)
        d2_only = [v[2] for v in d2_vals]
        d2_mean = sum(d2_only) / N
        d2_range = max(d2_only) - min(d2_only)

        # Normalize and compare to A_N eigenvalues
        if d2_range > 0:
            d2_normed = [(d - min(d2_only)) / d2_range for d in d2_only]
            an_eigenvalues = [2 * math.cos(k * math.pi / (N + 1)) for k in range(1, N + 1)]
            an_normed = [(e - min(an_eigenvalues)) / (max(an_eigenvalues) - min(an_eigenvalues))
                         for e in an_eigenvalues]
            rms = math.sqrt(sum((a - b) ** 2 for a, b in zip(d2_normed, an_normed)) / N)
            print(f"  A_{N} fit (N={N}): RMS residual = {rms:.4f}")
            print(f"  Normalized data:  {['%.3f' % v for v in d2_normed]}")
            print(f"  A_{N} eigenvals:  {['%.3f' % v for v in an_normed]}")
            print()

# ─── TEST 4: Small separation ratios and φ ────────────────────────────────────
# The small separation d_02 = ν(n,0) - ν(n-1,2) is sensitive to the
# core structure (i.e., the tachocline region).

print("\n── TEST 4: Small Separations d₀₂ = ν(n,0) − ν(n−1,2) ──")
print("These probe the deep interior / tachocline region.\n")

l0_dict = {n: f for n, f in modes.get(0, [])}
l2_dict = {n: f for n, f in modes.get(2, [])}

d02_vals = []
for n in sorted(l0_dict.keys()):
    if n - 1 in l2_dict:
        d02 = l0_dict[n] - l2_dict[n - 1]
        d02_vals.append((n, d02))
        print(f"  n={n:>2}: d₀₂ = {l0_dict[n]:.3f} - {l2_dict[n-1]:.3f} = {d02:>+8.3f} μHz")

if len(d02_vals) >= 2:
    print()
    d02_only = [v[1] for v in d02_vals]
    for i in range(len(d02_only) - 1):
        r = d02_only[i + 1] / d02_only[i] if d02_only[i] != 0 else float('inf')
        print(f"  d₀₂(n={d02_vals[i+1][0]})/d₀₂(n={d02_vals[i][0]}) = {r:.6f}  "
              f"(φ⁻¹ = {1/PHI:.6f}, diff = {r - 1/PHI:+.6f})")

# ─── TEST 5: Ratio of large separation to fundamental constants ──────────────

print("\n\n── TEST 5: Large Separation vs Fundamental Constants ──")

l0_freqs = sorted(modes.get(0, []))
if len(l0_freqs) >= 2:
    seps = [l0_freqs[i+1][1] - l0_freqs[i][1]
            for i in range(len(l0_freqs)-1)
            if l0_freqs[i+1][0] == l0_freqs[i][0] + 1]
    if seps:
        mean_dnu = sum(seps) / len(seps)
        print(f"  Mean Δν(l=0) = {mean_dnu:.3f} μHz")
        print(f"  Δν / α⁻¹     = {mean_dnu / ALPHA_INV:.6f}")
        print(f"  Δν / φ        = {mean_dnu / PHI:.6f}")
        print(f"  Δν × φ        = {mean_dnu * PHI:.6f}")
        print(f"  Δν / (φ × α⁻¹)= {mean_dnu / (PHI * ALPHA_INV):.6f}")
        print(f"  Δν / 135      = {mean_dnu / 135:.6f}  (canonical Δν ≈ 135 μHz)")

# ─── TEST 6: Explicit A_N eigenvalue overlay ─────────────────────────────────
# For a given N, compute λ_k = 2cos(kπ/(N+1)) and scale to match the
# frequency range. Try N = 7 (8 radial modes for l=0) and N = 8.

print("\n\n── TEST 6: Direct A_N Eigenvalue Overlay ──")
print("Scaling A_N eigenvalues to the observed frequency range.\n")

for l in [0, 1]:
    freqs_sorted = sorted(modes.get(l, []))
    if len(freqs_sorted) < 3:
        continue
    N = len(freqs_sorted)
    f_min = freqs_sorted[0][1]
    f_max = freqs_sorted[-1][1]

    an_raw = [2 * math.cos(k * math.pi / (N + 1)) for k in range(1, N + 1)]
    an_min, an_max = min(an_raw), max(an_raw)

    # Affine map: A_N eigenvalue range → observed frequency range
    scaled = [f_min + (f_max - f_min) * (e - an_min) / (an_max - an_min)
              for e in an_raw]

    print(f"  l = {l}, N = {N} modes")
    print(f"  {'k':>4}  {'A_N eig':>10}  {'scaled (μHz)':>12}  {'observed':>12}  {'Δ (μHz)':>10}")
    print(f"  {'-'*55}")
    rms = 0
    for i, (n, f_obs) in enumerate(freqs_sorted):
        diff = scaled[i] - f_obs
        rms += diff**2
        print(f"  {i+1:>4}  {an_raw[i]:>+10.5f}  {scaled[i]:>12.3f}  {f_obs:>12.3f}  {diff:>+10.3f}")
    rms = math.sqrt(rms / N)
    print(f"  RMS residual: {rms:.3f} μHz")
    print(f"  Frequency range: {f_max - f_min:.3f} μHz")
    print(f"  Relative RMS: {rms / (f_max - f_min) * 100:.2f}%")
    print()

# ─── TEST 7: Consecutive spacing ratios ──────────────────────────────────────
# For a pure A_N spectrum, consecutive eigenvalue spacings have a specific
# pattern determined by sin functions. Check if the observed spacing
# ratios match.

print("\n── TEST 7: Consecutive Spacing Ratios ──")
print("For A_N: Δλ_k / Δλ_{k-1} follows a specific sine pattern.\n")

for l in [0, 1, 2]:
    freqs_sorted = sorted(modes.get(l, []))
    if len(freqs_sorted) < 4:
        continue
    N = len(freqs_sorted)
    spacings = []
    for i in range(len(freqs_sorted) - 1):
        n1, f1 = freqs_sorted[i]
        n2, f2 = freqs_sorted[i + 1]
        spacings.append(f2 - f1)

    # A_N predicted spacing ratios
    an_raw = [2 * math.cos(k * math.pi / (N + 1)) for k in range(1, N + 1)]
    an_spacings = [an_raw[i + 1] - an_raw[i] for i in range(N - 1)]

    print(f"  l = {l}:")
    print(f"  {'k':>4}  {'Δν_obs':>10}  {'ratio_obs':>10}  {'ratio_A_N':>10}  {'diff':>10}")
    print(f"  {'-'*50}")
    for i in range(1, len(spacings)):
        r_obs = spacings[i] / spacings[i - 1] if spacings[i - 1] != 0 else float('inf')
        r_an = an_spacings[i] / an_spacings[i - 1] if i < len(an_spacings) and an_spacings[i - 1] != 0 else float('nan')
        diff = r_obs - r_an if not math.isnan(r_an) else float('nan')
        print(f"  {i+1:>4}  {spacings[i]:>10.3f}  {r_obs:>10.5f}  {r_an:>10.5f}  {diff:>+10.5f}")
    print()

print("\n" + "=" * 72)
print("ANALYSIS COMPLETE")
print("=" * 72)