"""
Neutrino Mass Prediction - Compact Version
For Supplementary Material Submission
"""

import numpy as np
from scipy import stats

# Constants
alpha_inv = 137.035999084  # α^{-1}
dm21_sq = 7.53e-5  # Δm21² [eV²]
dm31_sq = 2.453e-3  # |Δm31²| [eV²]

# Empirical relations from paper
ratio_21 = 2.003  # m2/m1
ratio_31 = alpha_inv / 13.5  # m3/m1 = α^{-1}/13.5

# Calculate masses
m1 = np.sqrt(dm21_sq / (ratio_21**2 - 1))  # [eV]
m2 = ratio_21 * m1
m3 = ratio_31 * m1

# Verify with oscillation data
dm21_sq_pred = m2**2 - m1**2
dm31_sq_pred = m3**2 - m1**2
diff_21 = abs(dm21_sq_pred - dm21_sq) / dm21_sq * 100
diff_31 = abs(dm31_sq_pred - dm31_sq) / dm31_sq * 100

# Experimental observables
theta12 = np.radians(33.44)
theta13 = np.radians(8.57)
U_e1 = np.cos(theta12) * np.cos(theta13)
U_e2 = np.sin(theta12) * np.cos(theta13)
U_e3 = np.sin(theta13)

m_ee = abs(m1 * U_e1**2 + m2 * U_e2**2 + m3 * U_e3**2)
m_beta = np.sqrt(m1**2 * U_e1**2 + m2**2 * U_e2**2 + m3**2 * U_e3**2)
sum_m_nu = m1 + m2 + m3

# Print results
print("NEUTRINO MASS PREDICTIONS")
print("=" * 50)
print(f"m1 = {m1*1000:.2f} meV")
print(f"m2 = {m2*1000:.2f} meV")
print(f"m3 = {m3*1000:.2f} meV")
print(f"Σm_ν = {sum_m_nu*1000:.1f} meV")
print()
print("VERIFICATION WITH OSCILLATION DATA")
print(f"Δm21²: predicted = {dm21_sq_pred:.3e}, experimental = {dm21_sq:.3e}")
print(f"Difference: {diff_21:.1f}%")
print(f"|Δm31²|: predicted = {dm31_sq_pred:.3e}, experimental = {dm31_sq:.3e}")
print(f"Difference: {diff_31:.1f}%")
print()
print("EXPERIMENTAL OBSERVABLES")
print(f"0νββ m_ee = {m_ee*1000:.2f} meV")
print(f"KATRIN m_β = {m_beta*1000:.2f} meV")
print()
print("CONSISTENCY CHECKS")
print(f"Planck bound Σm_ν < 120 meV: {'✓' if sum_m_nu*1000 < 120 else '✗'}")
print(f"KATRIN bound m_β < 800 meV: {'✓' if m_beta*1000 < 800 else '✗'}")
print(f"KamLAND-Zen range 36-156 meV for m_ee: {'✓' if 36 < m_ee*1000 < 156 else '✗'}")

# Save to file
with open("neutrino_results.txt", "w") as f:
    f.write("Parameter,Value,Units\n")
    f.write(f"m1,{m1*1000:.2f},meV\n")
    f.write(f"m2,{m2*1000:.2f},meV\n")
    f.write(f"m3,{m3*1000:.2f},meV\n")
    f.write(f"sum_m_nu,{sum_m_nu*1000:.1f},meV\n")
    f.write(f"m_ee,{m_ee*1000:.2f},meV\n")
    f.write(f"m_beta,{m_beta*1000:.2f},meV\n")
    f.write(f"dm21_diff,{diff_21:.1f},%\n")
    f.write(f"dm31_diff,{diff_31:.1f},%\n")

print("\nResults saved to neutrino_results.txt")