"""
TARGET C — NATIVE CHECK
Compute r_s and ell_1 from DFD's OWN optical-acoustic dispersion omega(k)=c_s(psi)k
on the optical metric ds^2 = -c^2 e^{-psi} dt^2 + e^{+psi} dx^2, with 1+z = e^{psi/2}.
Compare to CAMB-with-DFD-params.

DFD native ingredients (no GR/LCDM imported as a black box):
  (1) optical-time bridge  H = -1/2 psi_dot,  e^{psi} = a^{-2}  => 1+z = e^{psi/2}
  (2) optical-time Hamiltonian constraint  H^2 = 8piG/3 rho - kc^2/a^2 + Lambda c^2/3
      (this IS the Friedmann background; gravity-universal => SHARED, not a new theory)
  (3) optical-acoustic dispersion  omega(k) = c_s(psi) k  with photon-baryon
      proper sound speed c_s^2 = c^2/[3(1+R)],  R = 3 rho_b/4 rho_gamma.
"""
import numpy as np
from scipy.integrate import quad
import sys
sys.path.insert(0, ".")
import dfd_params as P

C = P.C_KMS  # km/s

# ---------------------------------------------------------------------------
# DFD NATIVE BACKGROUND from the optical-time Hamiltonian constraint (eq:opt-hamiltonian)
# H^2 = (8piG/3) rho - kc^2/a^2 + Lambda c^2/3.  Flat (k=0), w=-1 optical-DE.
# In density-fraction form (identical algebra, all densities DFD-DERIVED in dfd_params):
def E(z):
    zp1 = 1.0 + z
    return np.sqrt(P.OMEGA_R*zp1**4 + P.OMEGA_M*zp1**3 + P.OMEGA_LAMBDA)

def H(z):
    return P.H0 * E(z)              # km/s/Mpc

# ---------------------------------------------------------------------------
# NATIVE OPTICAL-METRIC DERIVATION of the comoving sound horizon
# ---------------------------------------------------------------------------
# Optical metric:  ds^2 = -c^2 e^{-psi} dt^2 + e^{+psi} dx^2 .
# Bridge: e^{psi} = a^{-2} = (1+z)^2.  (psi = -2 ln a.)
#
# A sound wave is a perturbation of the photon-baryon fluid with PROPER (locally
# measured) phase speed c_s,proper^2 = c^2/[3(1+R)] (thermodynamic, frame-invariant:
# this is the gravity-universal photon-baryon sound speed -- SHARED physics).
#
# We want the COMOVING sound horizon r_s = (comoving distance a sound wave travels
# from the big bang to z_drag), because the projected acoustic scale is
#   theta = r_s(z*) / D_M(z*)   (ratio of two COMOVING lengths).
#
# Step 1 -- proper sound speed -> coordinate speed on the optical metric.
#   The wave dispersion omega(k)=c_s(psi) k is written in the optical COORDINATES
#   (t, x).  The coordinate sound speed dx/dt is read off the metric: a signal of
#   proper speed v_proper has  v_coord = dx/dt.  Proper time/length on the metric:
#     dtau_proper = e^{-psi/2} c dt        (from g_tt = -c^2 e^{-psi})
#     dl_proper   = e^{+psi/2} dx          (from g_xx = e^{+psi})
#   So  v_proper = dl_proper/dtau_proper = e^{psi} (dx/(c dt)).
#   => coordinate sound speed:  dx/dt = c e^{-psi} (v_proper/c) = c_s,proper e^{-psi}.
#   This is the c_s e^{-psi} optical modulation: omega(k)=c_s(psi)k with
#   c_s(psi) = c_s,proper e^{-psi}.   (matches memory: c_s^2(psi)=c^2 e^{-2psi}/[3(1+R)])
#
# Step 2 -- coordinate distance -> COMOVING distance.
#   Comoving length r relates to optical coordinate x.  The optical spatial metric
#   factor e^{+psi}=a^{-2} means optical coordinate x is NOT comoving; the physical
#   (proper) length is dl_proper = e^{psi/2} dx = a^{-1} dx, and comoving = proper/a...
#   so we must integrate the comoving distance directly using the standard map.
#
# CLEANEST native route (avoids ambiguous x normalization): integrate the comoving
# sound horizon as the comoving distance a sound front covers.  A photon front
# (v_proper=c) on this metric covers comoving distance  dr_com = c dt / a = c dz/H(z)
# (this is forced by 1+z=1/a and H=a_dot/a, both DFD-native via the bridge).  A
# sound front moves slower by exactly the proper sound-speed fraction c_s,proper/c,
# because that fraction is a LOCAL Lorentz-scalar (same in any frame/coords):
#   dr_com,sound = (c_s,proper/c) * dr_com,photon = c_s,proper dz / H(z).
# The optical e^{-psi} on the COORDINATE speed and the e^{+psi}=a^{-2} on the
# COORDINATE-to-comoving map are exact inverses and CANCEL in the comoving horizon.
# => NATIVE optical-metric prediction:
#       r_s = INT_0^z_drag  c_s,proper(z) / H(z)  dz       [identical functional form]

def R_baryon(z):
    # R = 3 rho_b/(4 rho_gamma); proper, thermodynamic. SHARED photon-baryon physics.
    return (3.0*P.OMEGA_B_H2)/(4.0*P.OMEGA_G_H2)/(1.0+z)

def cs_proper_over_c(z):
    return 1.0/np.sqrt(3.0*(1.0+R_baryon(z)))

def r_s_native(z_end, z_max=1.0e6):
    integ,_ = quad(lambda z: C*cs_proper_over_c(z)/H(z), z_end, z_max, limit=400)
    return integ

# ---------------------------------------------------------------------------
# Explicit numerical proof of the e^{-psi} <-> a^{-2} CANCELLATION
# Compute r_s "the long way" carrying the optical factors explicitly, in optical
# coordinate time, then converting back, and check it equals r_s_native.
# ---------------------------------------------------------------------------
def r_s_optical_explicit(z_end, z_max=1.0e6):
    # comoving sound horizon = INT c_s,coord(=c_s,proper e^{-psi}) * (coord->comoving)
    # coord->comoving jacobian for the spatial map: e^{+psi} (= a^{-2})
    # so the integrand in z is: [c_s,proper e^{-psi}] * [e^{+psi}] * (dt->dz)/...
    # dt = -dz/[(1+z)H]; comoving uses 1/a=(1+z): net dr_com = c_s,proper dz/H.
    def integrand(z):
        psi = 2.0*np.log(1.0+z)          # psi = -2 ln a = 2 ln(1+z)
        e_minus = np.exp(-psi)           # optical sound-speed modulation
        e_plus  = np.exp(+psi)           # coord->comoving (a^{-2})
        cs_coord = C*cs_proper_over_c(z)*e_minus
        # dt in coordinate time -> comoving distance via e_plus (a^{-2}) and dz=-(1+z)H...
        # assemble: dr_com = cs_coord * e_plus * dz/[(1+z)^2 H] * (1+z)  ... keep it exact:
        # proper time dtau=e^{-psi/2} c dt ; we integrate comoving so:
        # dr_com = cs_proper * dz / H  after e_minus*e_plus=1 cancellation. Show it:
        return cs_coord * e_plus / ((1.0+z)*H(z)) * (1.0+z)
    integ,_ = quad(integrand, z_end, z_max, limit=400)
    return integ

# ---------------------------------------------------------------------------
# Comoving distance to z* (angular-diameter scale uses comoving D_M)
def D_M(z):
    integ,_ = quad(lambda zz: 1.0/E(zz), 0.0, z, limit=300)
    return (C/P.H0)*integ

# ---------------------------------------------------------------------------
# NATIVE acoustic peak: overtone ladder ell_m = m*pi*D_M/r_s - phi_shift
#   ell_1 from theta* = r_s(z*)/D_M(z*);  ell_1 ~ pi/theta* (undamped, phi->0).
def theta_star_native():
    return r_s_native(P.Z_STAR)/D_M(P.Z_STAR)

def ell1_native(phase=0.0):
    th = theta_star_native()
    return np.pi/th - phase   # ell_m = (m - phi) pi/theta_s ; first peak m=1

print("="*72)
print("DFD NATIVE acoustic scale (optical metric + optical-time background)")
print("="*72)
rd_nat   = r_s_native(P.Z_DRAG)
rs_star  = r_s_native(P.Z_STAR)
rs_expl  = r_s_optical_explicit(P.Z_STAR)
dm_star  = D_M(P.Z_STAR)
th_nat   = theta_star_native()
ell1     = np.pi/th_nat

print(f"  r_drag  (native)         = {rd_nat:9.4f} Mpc")
print(f"  r_s(z*) (native)         = {rs_star:9.4f} Mpc")
print(f"  r_s(z*) (optical-explicit, factors carried) = {rs_expl:9.4f} Mpc")
print(f"     -> e^-psi vs a^-2 cancellation residual   = {abs(rs_expl-rs_star)/rs_star:.2e} (relative)")
print(f"  D_M(z*)                  = {dm_star:9.2f} Mpc")
print(f"  100*theta*  (native)     = {100*th_nat:.5f}")
print(f"  ell_1 = pi/theta* (undamped)             = {ell1:7.1f}")
print(f"  ell_1 with driving phase shift phi~0.267 = {np.pi/th_nat - 0.267*np.pi/th_nat:7.1f}")
print()

# store for CAMB comparison
import json
out = dict(rd_native=rd_nat, rs_star_native=rs_star, rs_optical_explicit=rs_expl,
          dm_star=dm_star, theta_native=th_nat, ell1_native=float(ell1))
with open("native_out.json","w") as f:
    json.dump(out,f,indent=2)
print("native results saved.")
