#!/usr/bin/env python3
"""
chi_matter_verify.py  --  canonical numerical check of every quantity in
Appendix AV (chi-matter) of the DFD v4.0 corpus.

Run:  python3 chi_matter_verify.py

Reproduces, from first principles + measured (alpha, Planck-mass) constants only:
  f_chi = Mbar_P * alpha^3                       (decay constant, d_chi = 3)
  m_chi = sqrt(158) * Mbar_P * alpha^13           (see-saw, Lambda = Mbar_P alpha^8)
  Omega_chi h^2 = 1.62 * theta_i^2                (misalignment relic)
  theta_i for Omega_chi h^2 = 0.12                (the one O(1) free input)
  g_chi_gamma = alpha/(2 pi f_chi)                (two-photon coupling)
  lambda_res  = hc / m_chi                         (resonant conversion line)

Expected output (the numbers that appear in Appendix AV):
  f_chi ~ 9.46e11 GeV ; m_chi ~ 5.09 eV ; Omega(theta=1) ~ 1.62 ;
  theta_i ~ 0.27 ; g ~ 1.23e-15 GeV^-1 ; lambda ~ 243 nm.
"""
import math

# --- measured / fixed constants -------------------------------------------------
alpha = 1/137.035999084          # fine-structure constant (CODATA)
Mbar  = 2.435e18                 # reduced Planck mass, GeV  = (8 pi G)^-1/2
MPl   = 1.22091e19               # full Planck mass, GeV
T0    = 2.35e-13                 # CMB temperature today, GeV (2.725 K)
gs0   = 3.91                     # g_*s today
rho_crit_over_h2 = 8.0992e-47    # GeV^4

# --- Step 3: decay constant (d_chi = 3 = dim sl(2,C)) ---------------------------
f_chi = Mbar * alpha**3
print(f"alpha^3                  = {alpha**3:.4e}")
print(f"f_chi = Mbar*alpha^3     = {f_chi:.4e} GeV     (App AV: ~9.46e11)")

# --- Step 4: see-saw mass  m = Lambda^2 / f_chi,  Lambda = Mbar*alpha^8 ----------
#     => m_chi = sqrt(158) * Mbar * alpha^(16-3) = sqrt(158)*Mbar*alpha^13
m_chi = math.sqrt(158) * Mbar * alpha**13     # GeV
print(f"alpha^13                 = {alpha**13:.4e}")
print(f"m_chi = sqrt158*Mbar*a13 = {m_chi*1e9:.3f} eV       (App AV: ~5.1 eV)")

# --- Step 5: misalignment relic abundance --------------------------------------
def gstar(T):
    return 106.75 if T > 0.3 else 10.75
# oscillation onset: 3H = m_chi, radiation era H = 1.66 sqrt(g*) T^2 / MPl
T = 1e4
for _ in range(80):
    T = math.sqrt((m_chi/3) * MPl / (1.66 * math.sqrt(gstar(T))))
gosc = gstar(T)
aratio3 = (gs0 * T0**3) / (gosc * T**3)            # (a_osc/a0)^3, entropy-conserved
rho_per_th2 = 0.5 * m_chi**2 * f_chi**2 * aratio3  # GeV^4 per theta_i^2
Om_h2_per_th2 = rho_per_th2 / rho_crit_over_h2
theta_i = math.sqrt(0.12 / Om_h2_per_th2)
print(f"T_osc                    = {T:.3e} GeV   (g*={gosc})")
print(f"Omega_chi h^2 (theta=1)  = {Om_h2_per_th2:.3f}        (App AV: ~1.62)")
print(f"theta_i for Omega=0.12   = {theta_i:.3f}         (App AV: ~0.27, natural O(1))")

# --- Step 7: detection ---------------------------------------------------------
g = alpha / (2*math.pi*f_chi)
lam = 1239.84 / (m_chi*1e9)        # nm, with E in eV
tau = 64*math.pi / (g**2 * (m_chi)**3)   # GeV^-1; convert below
print(f"g_chi_gamma              = {g:.3e} GeV^-1   (App AV: ~1.2e-15)")
print(f"resonant wavelength      = {lam:.1f} nm        (App AV: ~244)")
# lifetime in seconds: hbar = 6.582e-25 GeV*s
tau_s = tau * 6.582e-25
print(f"lifetime tau_chi         = {tau_s:.1e} s         (App AV: ~1e32, stable)")
