"""
Test 7: Does the "cooling-by-expansion" brane mechanism produce a cosmology
that resembles observations?

Setup:
- The 3-brane has internal energy density rho_int(t) and scale factor a(t).
- EM radiation is confined to the brane, so the brane cannot radiatively cool
  into the 4D bulk.
- Hypothesis: the brane sheds internal energy adiabatically by expanding,
  i.e., dU = -p dV.

Two scenarios we test:
  (A) Ordinary adiabatic expansion: dU = -p dV with p = w*rho
       where w is an equation-of-state parameter.
  (B) "Blocked radiation" modification: if radiative channel is blocked,
       all energy loss must be mechanical. Expansion rate is determined by
       the energy disposal rate.

We compare against:
  - LambdaCDM expansion history (Omega_m ~ 0.3, Omega_Lambda ~ 0.7)
  - Observed acceleration onset at z ~ 0.7
  - Effective equation of state at late times (observed w ~ -1)

Key diagnostic: does the mechanism naturally produce
    (i) accelerating expansion at late times,
    (ii) an equation of state close to w = -1,
    (iii) an onset redshift consistent with observation?
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# ==============================================================
# SCENARIO A: Standard adiabatic cooling with various EOS
# ==============================================================
# First law: d(rho * V)/dt = -p * dV/dt
# With V = a^3, p = w*rho:
#   d(rho * a^3)/dt = -w*rho * d(a^3)/dt
#   => rho_dot * a^3 + 3*rho*a^2*a_dot = -3*w*rho*a^2*a_dot
#   => rho_dot = -3*(1+w)*rho*H   where H = a_dot/a
#
# This is the standard continuity equation. For it to drive expansion
# (Friedmann), we need a field equation. In LCDM:
#   H^2 = (8*pi*G/3) * rho_total
#
# For Jay's brane picture, we need to ASK: what sets H?
# His proposal: H is set by the rate at which the brane NEEDS to shed energy.
# If the brane can't radiate, it must expand fast enough to dissipate its
# excess energy via pdV work.

# Let me formalize this. Suppose the brane has a "thermalization" timescale
# tau such that dU/dt = -U/tau (exponential cooling, if unblocked).
# If EM cooling is blocked, this must be accomplished via expansion:
#   dU/dt = -p dV/dt = -w*rho * 3*a^2*a_dot = -3*w*rho*a^3 * H
# Setting these equal (the rate that "would have been" radiative becomes mechanical):
#   -U/tau = -3*w*rho*V*H
#   rho*V/tau = 3*w*rho*V*H
#   H = 1/(3*w*tau)
# 
# This gives a CONSTANT H if w and tau are constant. Constant H is de Sitter!
# That's exactly what dark energy looks like observationally.

# Let me test this more carefully.

print("=" * 70)
print("SCENARIO A: Forced mechanical cooling gives constant H")
print("=" * 70)

# If radiative cooling rate (per unit volume) is L, and it's blocked,
# then H must satisfy: 3*w*rho*H = L
# If L is proportional to rho (blackbody-like: L ~ sigma*T^4, T ~ rho^(1/4)),
# then actually L ~ rho, so rho cancels and H is constant.
#
# But if L ~ rho^n for some other power, H would depend on rho.

# Test: simulate the cosmology for various assumptions about L(rho)

def cooling_expansion_cosmology(L_exponent, w, tau0, t_max=20.0, rho0=1.0, a0=1.0):
    """
    Simulate expansion driven by blocked radiative cooling.
    
    L(rho) = (rho/rho0)^L_exponent / tau0  [energy per unit volume per unit time]
    Forced: 3*w*rho*H = L  =>  H = L / (3*w*rho)
    Continuity: rho_dot = -3*(1+w)*rho*H
    
    Returns t, a, rho, H arrays.
    """
    def dydt(t, y):
        a, rho = y
        if rho <= 0:
            return [0, 0]
        L = (rho/rho0)**L_exponent / tau0
        H = L / (3.0 * w * rho)
        a_dot = H * a
        rho_dot = -3.0 * (1.0 + w) * rho * H
        return [a_dot, rho_dot]
    
    sol = solve_ivp(dydt, [0, t_max], [a0, rho0], dense_output=True,
                    rtol=1e-8, atol=1e-10, max_step=0.01)
    t = sol.t
    a = sol.y[0]
    rho = sol.y[1]
    H = np.gradient(np.log(a), t)
    return t, a, rho, H

# Case 1: L ~ rho (blackbody-like with T ~ rho^(1/4), i.e., L ~ T^4 ~ rho)
# Actually, blackbody: L ~ sigma*T^4, and for radiation rho ~ T^4, so L ~ rho.
# For non-relativistic matter, L depends on microphysics.

fig, axes = plt.subplots(2, 2, figsize=(12, 9))

colors = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728']
exponents = [1.0, 1.25, 0.75, 0.5]
labels = [r'$L \propto \rho$ (blackbody)', r'$L \propto \rho^{1.25}$',
          r'$L \propto \rho^{0.75}$', r'$L \propto \rho^{0.5}$']

w_eos = 1/3  # radiation-like EOS for the internal energy being shed
tau0 = 1.0

for i, (n, lbl, c) in enumerate(zip(exponents, labels, colors)):
    t, a, rho, H = cooling_expansion_cosmology(n, w_eos, tau0, t_max=10)
    axes[0,0].plot(t, a, color=c, label=lbl)
    axes[0,1].plot(t, rho, color=c, label=lbl)
    axes[1,0].plot(t, H, color=c, label=lbl)
    # Acceleration: a_ddot = d(H*a)/dt = H_dot*a + H^2*a = a*(H_dot + H^2)
    # q = -a_ddot*a/a_dot^2 = -(H_dot + H^2)/H^2 = -1 - H_dot/H^2
    H_dot = np.gradient(H, t)
    q = -1.0 - H_dot / (H**2 + 1e-30)
    axes[1,1].plot(t, q, color=c, label=lbl)

axes[0,0].set_xlabel('time')
axes[0,0].set_ylabel('scale factor a(t)')
axes[0,0].set_title('Scale factor evolution')
axes[0,0].set_yscale('log')
axes[0,0].legend(fontsize=8)
axes[0,0].grid(alpha=0.3)

axes[0,1].set_xlabel('time')
axes[0,1].set_ylabel(r'energy density $\rho$')
axes[0,1].set_title('Energy density')
axes[0,1].set_yscale('log')
axes[0,1].legend(fontsize=8)
axes[0,1].grid(alpha=0.3)

axes[1,0].set_xlabel('time')
axes[1,0].set_ylabel('H(t)')
axes[1,0].set_title('Hubble rate')
axes[1,0].legend(fontsize=8)
axes[1,0].grid(alpha=0.3)

axes[1,1].set_xlabel('time')
axes[1,1].set_ylabel('deceleration parameter q')
axes[1,1].set_title('q = -1: de Sitter / pure DE | q = 0: coasting | q > 0: decelerating')
axes[1,1].axhline(0, color='k', linestyle=':', alpha=0.5)
axes[1,1].axhline(-1, color='r', linestyle=':', alpha=0.5, label='q=-1 (de Sitter)')
axes[1,1].legend(fontsize=8)
axes[1,1].grid(alpha=0.3)

plt.tight_layout()
plt.savefig('/home/claude/cooling_expansion_test.png', dpi=120)
print("Saved figure: cooling_expansion_test.png")

# Analysis: for each case, what's the asymptotic behavior?
print("\nAsymptotic behavior for each case:")
for n, lbl in zip(exponents, labels):
    t, a, rho, H = cooling_expansion_cosmology(n, w_eos, tau0, t_max=10)
    # Late-time H
    H_late = np.mean(H[-100:])
    # Deceleration parameter at late time
    H_dot = np.gradient(H, t)
    q_late = -1.0 - np.mean(H_dot[-100:]) / (np.mean(H[-100:])**2 + 1e-30)
    # Effective equation of state: w_eff = -1 - (2/3)*H_dot/H^2
    w_eff_late = -1.0 - (2.0/3.0) * np.mean(H_dot[-100:]) / (np.mean(H[-100:])**2 + 1e-30)
    print(f"  {lbl}: H_late = {H_late:.3f}, q_late = {q_late:.3f}, w_eff = {w_eff_late:.3f}")

# ==============================================================
# KEY RESULT ANALYSIS
# ==============================================================
print("\n" + "=" * 70)
print("KEY ANALYTICAL RESULT")
print("=" * 70)
print("""
For L(rho) = (rho/rho0)^n / tau0 and forced mechanical cooling:
    H = L / (3*w*rho) = (rho/rho0)^(n-1) / (3*w*tau0*rho0)

CASE n = 1 (blackbody-like, L ~ rho):
    H = 1 / (3*w*tau0*rho0)  <-- CONSTANT, independent of rho
    => pure de Sitter expansion
    => w_eff = -1 exactly
    => q = -1 exactly
    
This is the closest thing to observed dark energy.

CASE n != 1:
    H depends on rho, and since rho decreases as a grows, H evolves.
    - n > 1: H increases as rho decreases => phantom-like (w < -1)
    - n < 1: H decreases => quintessence-like (w > -1)
    
Observed cosmology: w = -1.03 +/- 0.03 (Planck 2018 + SN + BAO)
    => Current data strongly consistent with n = 1 (exact de Sitter)
    => Allowing small deviations from n = 1 gives room for
       quintessence or phantom fits within observational error
""")

# ==============================================================
# COMPARISON WITH STANDARD LCDM
# ==============================================================
print("=" * 70)
print("COMPARISON: LambdaCDM vs cooling-expansion with n=1")
print("=" * 70)

# LCDM: H^2 = H0^2 * (Omega_m/a^3 + Omega_L)
# Cooling-expansion: H = const (for n=1), no matter contribution

# The cooling-expansion gives PURE de Sitter from the start.
# LCDM has matter domination early, then transitions to DE around z ~ 0.7.

# For the cooling-expansion model to match LCDM, we'd need:
# - Some mechanism that gives matter-like behavior early
# - Then transitions to the blocked-radiation regime

# One possibility: the "L ~ rho" relation only holds for the NON-radiative
# component of the energy density. If there's an ordinary matter component
# that dilutes as a^(-3), it would dominate early and de-Sitter component 
# dominates late.

# Let's test a TWO-component model:

def two_component_model(t_max=30.0, rho_m0=1.0, rho_de0=0.3, H_de=0.1):
    """
    rho_m: ordinary matter, dilutes as a^-3
    rho_de: 'blocked thermal' component, drives constant-H expansion
    
    Total expansion: H^2 = (8*pi*G/3)*(rho_m + rho_de) in standard picture
    In cooling-expansion: H is forced by the de component only.
    
    Let's use the standard Friedmann picture for comparison.
    """
    G_eff = 1.0
    def dydt(t, y):
        a, rho_m = y
        if a <= 0:
            return [0, 0]
        H = np.sqrt((8*np.pi*G_eff/3) * (rho_m + rho_de0))
        a_dot = H * a
        rho_m_dot = -3 * H * rho_m
        return [a_dot, rho_m_dot]
    
    sol = solve_ivp(dydt, [0, t_max], [1.0, rho_m0], dense_output=True,
                    rtol=1e-8, atol=1e-10, max_step=0.01)
    t = sol.t
    a = sol.y[0]
    rho_m = sol.y[1]
    H = np.sqrt((8*np.pi*G_eff/3) * (rho_m + rho_de0))
    H_dot = np.gradient(H, t)
    q = -1.0 - H_dot / (H**2 + 1e-30)
    return t, a, rho_m, H, q

# Pure blocked-thermal (single component, n=1)
def blocked_thermal_single():
    # H = constant = 1/(3*w*tau0*rho0_norm)
    # rho decays exponentially
    w = 1/3
    tau0 = 1.0
    rho0 = 1.0
    H_const = 1.0 / (3*w*tau0*rho0)  # For n=1
    # Actually with n=1: H = (rho/rho0)^0 / (3*w*tau0*rho0) = 1/(3*w*tau0*rho0)
    # rho_dot = -3*(1+w)*rho*H
    t = np.linspace(0, 20, 2000)
    rho = rho0 * np.exp(-3*(1+w)*H_const*t)
    a = np.exp(H_const * t)
    H = np.full_like(t, H_const)
    q = np.full_like(t, -1.0)
    return t, a, rho, H, q

fig, axes = plt.subplots(1, 2, figsize=(13, 5))

# Comparison plot: a(t) for different models
t1, a1, rho1, H1, q1 = blocked_thermal_single()
t2, a2, rho_m2, H2, q2 = two_component_model()

axes[0].plot(t1, a1, 'b-', label='Pure blocked-thermal (n=1)', linewidth=2)
axes[0].plot(t2, a2, 'r-', label='Matter + blocked-thermal (LCDM-like)', linewidth=2)
axes[0].set_xlabel('cosmic time')
axes[0].set_ylabel('scale factor a(t)')
axes[0].set_title('Scale factor: pure de Sitter vs mixed model')
axes[0].set_yscale('log')
axes[0].legend()
axes[0].grid(alpha=0.3)

axes[1].plot(t1, q1, 'b-', label='Pure blocked-thermal', linewidth=2)
axes[1].plot(t2, q2, 'r-', label='Matter + blocked-thermal', linewidth=2)
axes[1].axhline(0, color='k', linestyle=':', alpha=0.5)
axes[1].axhline(-1, color='gray', linestyle=':', alpha=0.5, label='q=-1')
axes[1].set_xlabel('cosmic time')
axes[1].set_ylabel('deceleration parameter q')
axes[1].set_title('Deceleration parameter: transition from matter to DE')
axes[1].legend()
axes[1].grid(alpha=0.3)

plt.tight_layout()
plt.savefig('/home/claude/cosmology_comparison.png', dpi=120)
print("Saved figure: cosmology_comparison.png")

# ==============================================================
# CRITICAL QUESTIONS
# ==============================================================
print("\n" + "=" * 70)
print("CRITICAL ASSESSMENT")
print("=" * 70)
print("""
The numerical result is striking: if the blocked-radiation cooling
luminosity scales as L ~ rho^1 (which is the blackbody-like case,
since rho_rad ~ T^4 and L ~ T^4 ~ rho), then the forced-mechanical-cooling
equation H = L/(3w*rho) gives H = constant EXACTLY.

Constant H = de Sitter = cosmological constant behaviour.
Equation of state w_eff = -1 exactly.

This is a genuinely interesting result. But before getting excited,
let me check whether this is physics or a tautology.

TAUTOLOGY CHECK:
The derivation went:
  (1) Energy loss rate from blackbody-like process: dU/dt = -L*V, L ~ rho
  (2) This must equal mechanical work: dU/dt = -p*dV/dt = -w*rho*3*a^2*a_dot
  (3) Setting these equal: L*V = w*rho*3*a^2*a_dot => H = L/(3w*rho)
  (4) For L ~ rho: H = const

The logic relies on:
  - Energy conservation (OK)
  - The assumption that L ~ rho (blackbody; needs justification)
  - The assumption that NO other channel exists (bulk gravity? tension?)
  - The interpretation that mechanical work in 3D expansion drains the
    SAME energy that radiation would have drained

Issues:
  A. In standard cosmology, radiation ALSO expands adiabatically. Radiation
     doesn't "need" to be blocked from bulk radiation; it just redshifts.
     Redshift alone gives rho_rad ~ a^-4, not constant-H behaviour.
     
  B. The "forced" equation 3*w*rho*H = L is NOT a force equation. It's
     just energy conservation. It doesn't tell you what H actually IS;
     it tells you what rho does given H. To get H, you need a DYNAMICAL
     equation (Friedmann or equivalent).
     
  C. So the question is: in Jay's brane framework, what plays the role
     of the Friedmann equation? Without that, we can't actually predict H;
     we can only say IF expansion is the only cooling channel, energy
     conservation relates H to L.

This is the key issue. Let me think about this more carefully.
""")
