# -*- coding: utf-8 -*-
"""
Friedmann-Analyse des Trigger-Verstärker-Mechanismus (Begleitcode zur Hauptarbeit, Abschnitt "Friedmann-Integration und w(z)")

Modell (dimensionslos, H0 = 1, Dichten in Einheiten von rho_crit,0):
  Feldgleichung (x = ln a, ' = d/dx):
     E^2 psi'' + (E E' + 3 E^2) psi' + mt^2 psi = q * exp(-4x)      (Photonquelle)
  Portal:      xi = psi^2
  Lambda-Boost: rho_boost = K * zeta^-2 * xi^2,  zeta^2 = 1 + xi
  Friedmann:   E^2 = Om*e^-3x + Or*e^-4x + rho_DE(x)

Drei Regime:
  S1: mt >> 1  (Tracking)      -> analytisch: boost ~ a^-8, w = +5/3
  S2: mt <~ 1  (eingefroren, erwacht heute) -> numerisch, w(z)
  S3: mt ~ 50  (oszillierend)  -> mittelt zu Staub, w ~ 0
Zusätzlich: Early-Dark-Energy-Schranke für das Tracking-Regime.
"""
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

Om, Or = 0.315, 9.0e-5          # Planck-2018-artig
Ode_target = 1.0 - Om - Or       # ~0.685
zstar = 1090.0                   # Rekombination
xi0, x0 = np.log(1e-7), 0.0      # Integration von z=1e7 bis heute

# ----------------------------------------------------------------------
# Hintergrund-E(x) (Iteration 0: LCDM)
def E2_lcdm(x):
    return Om*np.exp(-3*x) + Or*np.exp(-4*x) + Ode_target

def dlnE_dx(E2fun, x, h=1e-4):
    return 0.5*(np.log(E2fun(x+h)) - np.log(E2fun(x-h)))/(2*h)*2  # d ln E/dx

# ----------------------------------------------------------------------
# Feldintegration: linear in psi -> Amplitude q frei skalierbar
def integrate_field(mt, q, E2fun, x_ini=xi0, n=4000):
    def rhs(x, y):
        psi, dpsi = y
        E2 = E2fun(x)
        dlnE = dlnE_dx(E2fun, x)
        ddpsi = (q*np.exp(-4*x) - mt**2*psi)/E2 - (dlnE + 3.0)*dpsi
        return [dpsi, ddpsi]
    xs = np.linspace(x_ini, x0, n)
    sol = solve_ivp(rhs, (x_ini, x0), [0.0, 0.0], t_eval=xs,
                    method="LSODA", rtol=1e-9, atol=1e-14)
    return sol.t, sol.y[0], sol.y[1]

def w_eff(xs, rho):
    lnr = np.log(np.clip(rho, 1e-300, None))
    dlnr = np.gradient(lnr, xs)
    return -1.0 - dlnr/3.0

# ======================================================================
# S1 Tracking-Regime (mt >> 1): analytisch
#   psi_eq = q e^{-4x}/mt^2  ->  boost = K psi^4 = A e^{-8x}
#   w = -1 - (1/3)(-8) = +5/3;  Early-DE-Schranke:
rho_tot_star = Om*(1+zstar)**3 + Or*(1+zstar)**4
A_max = 0.01*rho_tot_star/((1+zstar)**8)   # boost(z*) < 1% der Gesamtdichte
print("=== S1: Tracking (m >> H) ===")
print(f"w_eff = +5/3 (steigt wie a^-8, schneller als Strahlung)")
print(f"Early-DE-Schranke: boost heute < {A_max:.2e} rho_crit")
print(f"  -> Faktor {Ode_target/A_max:.1e} unter der DE-Dichte. Regime tot.\n")

# ======================================================================
# S2 Eingefrorenes Feld, erwacht heute (mt <= 1)
print("=== S2: eingefroren/erwachend ===")
results = {}
for mt in [0.3, 1.0, 3.0]:
    xs, psi, dpsi = integrate_field(mt, q=1.0, E2fun=E2_lcdm)
    # Fall (a): Lambda-Boost dominiert DE: rho_DE ~ K psi^4, K so dass heute 0.685
    K = Ode_target/max(psi[-1]**4, 1e-300)
    rho_a = K*psi**4
    # Fall (b): Feldenergie dominiert: rho_chi = s(0.5 E^2 psi'^2 + 0.5 mt^2 psi^2)
    E2s = E2_lcdm(xs)
    rho_chi_raw = 0.5*E2s*dpsi**2 + 0.5*mt**2*psi**2
    s = Ode_target/max(rho_chi_raw[-1], 1e-300)
    rho_b = s*rho_chi_raw
    wa, wb = w_eff(xs, rho_a), w_eff(xs, rho_b)
    z = np.exp(-xs) - 1.0
    results[mt] = (z, rho_a, rho_b, wa, wb, psi)
    # Selbstkonsistenz-Iteration (eine Runde) für Fall (a)
    from scipy.interpolate import interp1d
    rhoDE_i = interp1d(xs, rho_a, bounds_error=False,
                       fill_value=(rho_a[0], rho_a[-1]))
    def E2_it(x): return Om*np.exp(-3*x) + Or*np.exp(-4*x) + float(rhoDE_i(x))
    xs2, psi2, dpsi2 = integrate_field(mt, 1.0, E2_it)
    K2 = Ode_target/max(psi2[-1]**4, 1e-300)
    wa2 = w_eff(xs2, K2*psi2**4)
    i0 = np.argmin(np.abs(np.exp(-xs)-1.0))       # z=0
    i05 = np.argmin(np.abs(np.exp(-xs)-1.5))      # z=0.5
    print(f" mt={mt:4.1f}:  w0(a)={wa[-1]:+.2f} (iter: {wa2[-1]:+.2f})   "
          f"w(z=0.5)(a)={wa[i05]:+.2f}   w0(b)={wb[-1]:+.2f}")

# ======================================================================
# S3 Oszillierend (mt = 50): Mittelung -> Staub
print("\n=== S3: oszillierend (m ~ 50 H0) ===")
xs3, psi3, dpsi3 = integrate_field(50.0, q=1.0, E2fun=E2_lcdm, n=20000)
E2s3 = E2_lcdm(xs3)
rho3 = 0.5*E2s3*dpsi3**2 + 0.5*50.0**2*psi3**2
# Mittleres w über letzte e-Faltung:
mask = xs3 > -1.0
w3 = w_eff(xs3, rho3)
print(f" <w> (letzte e-Faltung) = {np.mean(w3[mask]):+.2f}  (Staub-artig, w~0)")

# ======================================================================
# Plots
fig, axes = plt.subplots(1, 3, figsize=(15, 4.4))

# (1) Feldentwicklung xi(z) im erwachenden Regime
ax = axes[0]
for mt, c in zip([0.3, 1.0, 3.0], ["C0", "C1", "C2"]):
    z, ra, rb, wa, wb, psi = results[mt]
    ax.loglog(1+z, psi**2/max(psi[-1]**2, 1e-300), c, label=fr"$m/H_0={mt}$")
ax.set_xlabel(r"$1+z$"); ax.set_ylabel(r"$\xi(z)/\xi(0)$")
ax.set_title("S2: Anregung erwacht erst heute")
ax.set_xlim(1, 1e4); ax.legend(); ax.grid(alpha=0.3)

# (2) w(z) aller Regime vs. Beobachtung
ax = axes[1]
zz = np.linspace(0, 3, 100)
ax.fill_between(zz, -1.1, -0.9, color="0.85",
                label="beobachtungsvertr. Band")
ax.axhline(5/3, color="C3", ls="--", label=r"S1 Tracking: $w=+5/3$")
ax.axhline(0.0, color="C4", ls=":", label=r"S3 oszill.: $w\approx 0$")
for mt, c in zip([0.3, 1.0, 3.0], ["C0", "C1", "C2"]):
    z, ra, rb, wa, wb, psi = results[mt]
    m = (z >= 0) & (z <= 3)
    ax.plot(z[m], wa[m], c, label=fr"S2 Boost, $m/H_0={mt}$")
ax.set_xlabel(r"$z$"); ax.set_ylabel(r"$w_{\rm eff}(z)$")
ax.set_ylim(-6, 2.2); ax.set_title(r"Zustandsgleichung $w(z)$")
ax.legend(fontsize=7, loc="lower right"); ax.grid(alpha=0.3)

# (3) Early-DE-Katastrophe im Tracking-Regime
ax = axes[2]
zz2 = np.logspace(0, 3.2, 200)
for frac, c in zip([1.0, 1e-6, 1e-12, A_max/Ode_target], ["C3", "C1", "C0", "C2"]):
    boost = frac*Ode_target*zz2**8
    lbl = (r"$\rho_{\rm boost,0}=\rho_{\rm DE}$" if frac == 1.0 else
           fr"$10^{{{int(np.log10(frac))}}}\,\rho_{{\rm DE}}$")
    ax.loglog(zz2, boost, c, label=lbl)
rho_tot = Om*zz2**3 + Or*zz2**4 + Ode_target
ax.loglog(zz2, rho_tot, "k", lw=2, label=r"$\rho_{\rm total}$ (Std.)")
ax.axvline(1+zstar, color="0.5", ls=":")
ax.text(1+zstar, 1e3, " CMB", rotation=90, va="bottom", fontsize=8)
ax.set_xlabel(r"$1+z$"); ax.set_ylabel(r"$\rho/\rho_{c,0}$")
ax.set_title(r"S1: Boost $\propto (1+z)^8$ sprengt frühe Kosmologie")
ax.set_ylim(1e-8, 1e30); ax.legend(fontsize=7, loc="upper left")
ax.grid(alpha=0.3)

plt.tight_layout()
plt.savefig("/home/claude/friedmann_trigger.png", dpi=150)
print("\nPlot gespeichert.")

# Kernzahlen für den Bericht
z05 = {}
for mt in [0.3, 1.0, 3.0]:
    z, ra, rb, wa, wb, psi = results[mt]
    i05 = np.argmin(np.abs(z-0.5))
    z05[mt] = (wa[-1], wa[i05], wb[-1])
print("\nKernzahlen S2 (Boost-Fall):")
for mt, (w0, w05, w0b) in z05.items():
    print(f"  m/H0={mt}: w0={w0:+.2f}, w(0.5)={w05:+.2f}; Feldenergie-Fall w0={w0b:+.2f}")
print(f"\nS1 Early-DE: max. erlaubter Boost heute = {A_max:.1e} rho_c "
      f"= {A_max/Ode_target:.1e} * rho_DE")