import numpy as np
# ============ CONSTANTS (no fudge: standard values) ============
M_Pl   = 1.220910e19      # full Planck mass, GeV
Mbar   = 2.435323e18      # reduced Planck mass, GeV
alpha  = 1/137.035999
s158   = np.sqrt(158.0)   # = 3N/2 - 1, N=2(k_max-dim M7)=2(60-7)=106 ; S^3 spectral, FIXED
s0     = 2.2255e-38       # entropy density today, GeV^3  (2891 cm^-3)
rhoc_h2= 8.0992e-47       # rho_crit,0 / h^2, GeV^4
gstar  = 106.75           # relativistic dof at T_osc (T_osc >> EW for these masses)

# ============ 1. MASS  m_chi = sqrt158 * Mbar * alpha^(16-n)  (see-saw: 16-n=2*8-n) ============
def m_chi(n): return s158 * Mbar * alpha**(16-n)            # GeV
def f_chi(n): return Mbar * alpha**n                         # GeV  (the open exponent)
print("MASS (forced once n is fixed; see-saw 16-n, sqrt158 from S^3):")
for n in [5,3]:
    m=m_chi(n); f=f_chi(n)
    print(f"  n={n}: m_chi = {m*1e9:.3g} eV = {m*1e6:.3g} keV    f_chi = {f:.3g} GeV")

# ============ 2. ABUNDANCE  (rigorous misalignment, constant-mass ALP, RD era) ============
def Omega_h2(n, theta=1.0):
    m=m_chi(n); f=f_chi(n)
    Tosc = np.sqrt(m*Mbar/( (3*np.sqrt(np.pi**2*gstar/90.0)) /1.0 * np.sqrt(1) ))  # placeholder
    # 3H=m, H=sqrt(pi^2 g/90) T^2/Mbar  -> Tosc^2 = m*Mbar/(3*sqrt(pi^2 g/90))
    Tosc = np.sqrt( m*Mbar/(3*np.sqrt(np.pi**2*gstar/90.0)) )
    n_over_s = (0.5*m*f**2*theta**2)/((2*np.pi**2/45.0)*gstar*Tosc**3)
    rho0 = m * n_over_s * s0
    return rho0/rhoc_h2, m, Tosc
print("\nABUNDANCE (rigorous misalignment, theta_i=1):")
for n in [5,3]:
    O,m,T=Omega_h2(n,1.0)
    print(f"  n={n}: Omega_chi h^2 = {O:.3g}  (T_osc={T:.2g} GeV)   {'<< 0.12 FAILS' if O<0.01 else ''}")
# theta needed for 0.12 at n=3
O3,_,_=Omega_h2(3,1.0); th=np.sqrt(0.12/O3)
print(f"\n  n=3 hits the OBSERVED 0.12 at theta_i = {th:.3f}  (natural O(1) misalignment angle)")
print(f"  (n=5 underproduces by ~{0.12/Omega_h2(5,1.0)[0]:.1e}x even at theta=pi -> chi NOT dark matter at n=5)")

# ============ 3. DETECTION (forced, given the branch) ============
print("\nDETECTION signatures (derived, no freedom):")
for n,label in [(3,'DM branch'),(5,'X-ray branch')]:
    m=m_chi(n); f=f_chi(n)
    g = alpha/(2*np.pi*f)                       # axion-photon coupling ~ alpha/(2pi f), GeV^-1
    Egamma = m/2*1e9                            # chi->gamma gamma photon energy, eV
    print(f"  n={n} ({label}): m={m*1e9:.2g} eV, photon E={Egamma:.2g} eV, g_chi-gamma ~ {g:.2g} GeV^-1")
print("  (CAST bound g<6.6e-11; ADMX/IAXO future reach ~1e-12..1e-15)")

print("\n============ HONEST VERDICT ============")
print("FORCED: mass (sqrt158*Mbar*alpha^(16-n)), detection signature (g=alpha/2pi f).")
print("FORCED *given n*: abundance Omega(n=3) ~ 0.12 for theta_i~0.27 (O(1)).")
print("OPEN (the ONE lemma): is the decay-constant exponent n=3 (=dim S^3) or n=5 (chiral)?")
print("  v4.0 integer catalog HAS 3=dim S^3, has NO standalone 5. But Toeplitz proof unbuilt.")
