"""
Prediction #38 — Mass as a function of Compton wavelength across Z = 1–92

Spectrum hypothesis: when plotting atomic number Z against log(Compton wavelength),
the residuals from a smooth trend cluster at the nuclear magic numbers
Z ∈ {2, 8, 20, 28, 50, 82}.

Falsifier: chi-square / Welch t-test on |residuals| at magic Z vs other Z.
Threshold: p < 0.01 two-sided.
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from mendeleev import element

# Physical constants (CODATA 2018)
h = 6.62607015e-34         # Planck constant [J·s]
c = 299_792_458            # speed of light [m/s]
u_to_kg = 1.66053906660e-27  # atomic mass unit [kg]

MAGIC_Z = {2, 8, 20, 28, 50, 82}

def atomic_mass_kg(Z: int) -> float:
    """Standard atomic weight in kg (natural-abundance-weighted)."""
    el = element(Z)
    return el.atomic_weight * u_to_kg

def compton_wavelength(mass_kg: float) -> float:
    return h / (mass_kg * c)

# Build dataset for Z = 1 .. 92
Z_arr = np.arange(1, 93)
mass_arr = np.array([atomic_mass_kg(int(z)) for z in Z_arr])
lambda_C = compton_wavelength(mass_arr)
log_lambda = np.log10(lambda_C)

# --- Step 1: smooth trend fit ---
# Mass roughly scales as A ~ 2-2.5 × Z, so log(λ_C) is monotone but curved.
# We fit a polynomial of degree 3 to capture the bulk trend without
# overfitting local shell-structure variations.
poly_deg = 3
coeffs = np.polyfit(Z_arr, log_lambda, deg=poly_deg)
trend = np.polyval(coeffs, Z_arr)
residuals = log_lambda - trend

# --- Step 2: separate residuals at magic Z vs other Z ---
mask_magic = np.array([z in MAGIC_Z for z in Z_arr])
res_magic = residuals[mask_magic]
res_other = residuals[~mask_magic]

abs_magic = np.abs(res_magic)
abs_other = np.abs(res_other)

# --- Step 3: Welch's t-test on |residuals| ---
t_stat, p_value = stats.ttest_ind(abs_magic, abs_other, equal_var=False)

# Sign test: are magic-Z residuals systematically positive?
# (More binding → less mass → larger λ_C → positive residual)
sign_test = stats.binomtest(int(np.sum(res_magic > 0)), n=len(res_magic), p=0.5)

# --- Print numerical summary ---
print("=" * 60)
print("PREDICTION #38 — Compton wavelength magic-number signature")
print("=" * 60)
print(f"Polynomial trend degree: {poly_deg}")
print(f"Magic Z tested: {sorted(MAGIC_Z)}")
print()
print(f"Residuals at magic Z:")
for z, r in zip(Z_arr[mask_magic], res_magic):
    print(f"  Z = {z:3d}  residual = {r:+.4e}")
print()
print(f"|residual| mean (magic):  {abs_magic.mean():.4e}")
print(f"|residual| mean (other):  {abs_other.mean():.4e}")
print(f"Ratio (magic / other):    {abs_magic.mean()/abs_other.mean():.3f}×")
print()
print(f"Welch's t-test on |residuals|:")
print(f"  t-statistic = {t_stat:+.3f}")
print(f"  p-value     = {p_value:.4f}")
print(f"  Significant at p<0.01? {'YES' if p_value < 0.01 else 'NO'}")
print()
print(f"Sign test on magic-Z residuals (positive means binding-energy effect):")
print(f"  positives: {int(np.sum(res_magic > 0))} / {len(res_magic)}")
print(f"  p-value:   {sign_test.pvalue:.4f}")
print("=" * 60)

# --- Step 4: plot ---
fig, axes = plt.subplots(2, 1, figsize=(10, 8), sharex=True)

# Top: log(λ_C) vs Z with trend
ax = axes[0]
ax.plot(Z_arr, log_lambda, 'o-', color='#0a0a1a', markersize=4,
        linewidth=0.8, label=r'log$_{10}$(λ$_C$) [m]', zorder=2)
ax.plot(Z_arr, trend, '--', color='#c9a84c', linewidth=1.5,
        label=f'polynomial trend (deg {poly_deg})', zorder=3)
for z in sorted(MAGIC_Z):
    ax.axvline(z, color='#00d4ff', alpha=0.25, linewidth=1, zorder=1)
ax.set_ylabel(r'log$_{10}$(Compton wavelength) [m]')
ax.set_title('Prediction #38 — Periodic Table as Frequency Spectrum')
ax.legend(loc='upper right')
ax.grid(alpha=0.2)

# Bottom: residuals
ax = axes[1]
colors = ['#00d4ff' if z in MAGIC_Z else '#888' for z in Z_arr]
ax.bar(Z_arr, residuals, color=colors, edgecolor='none', width=0.8)
ax.axhline(0, color='#0a0a1a', linewidth=0.5)
for z in sorted(MAGIC_Z):
    ax.axvline(z, color='#00d4ff', alpha=0.15, linewidth=1, zorder=0)
ax.set_xlabel('Atomic number Z')
ax.set_ylabel('Residual from trend')
ax.set_title(f'Residuals · Welch t={t_stat:+.2f}, p={p_value:.3f} · '
             f'magic |res|/other |res| = {abs_magic.mean()/abs_other.mean():.2f}×')
ax.grid(alpha=0.2)

plt.tight_layout()
out_path = '/Users/maraldbes/Projects/spectrum-video-production/Theorie/experimenten/voorspelling_38_residuals.png'
plt.savefig(out_path, dpi=150, bbox_inches='tight', facecolor='white')
print(f"\nPlot saved → {out_path}")