"""
CHANNEL 2 — SUB-HUBBLE LOW-x STRUCTURE POCKETS.
Does integrating the optical screen through CLUMPED matter (galaxies/clusters, where peculiar
acceleration -> x<1 deep-MOND, sqrt revives) push the LINE-OF-SIGHT effective gate up to
sqrt(Omega_b/Omega_m)=0.396?  Or is the LOS dominated by smooth x>1 regions so it stays ~0.24?

Gate physics recap (from derive_screen_amplitude.py / four_case_decomposition.py):
  A = 2 sqrt(alpha) * GATE.   2 sqrt(alpha)=0.17085 is the FIXED deep-MOND boost at g=cH.
  GATE interpolates with x=g/a0:
     x->0 (deep-MOND):  GATE -> sqrt(f_b)  = sqrt(Omega_b/Omega_m) = 0.3964   (the unreached ceiling)
     x->inf (Newtonian):GATE -> f_b        = Omega_b/Omega_m       = 0.1572   (the floor)
  Smooth cosmic regime x_struct(z=0)=sqrt(Om)/(2 sqrt(alpha))=3.06 -> Newtonian-leaning,
     GATE ~ f_b^p with p in [0.72,0.79] -> 0.23-0.26 -> A=0.040-0.045 (40% SHORT of 0.0686).

The honest LOS question: A = 2 sqrt(alpha) * <GATE>_LOS, where <.>_LOS is the average of the LOCAL
gate over the photon path, weighted by HOW THE SCREEN ACCUMULATES per unit path. We test three
defensible weightings and a realistic lognormal matter PDF to see if clumping lifts <GATE>.

honesty-guard: do NOT assume clumping saves it. Report exactly what the integral gives for each weighting.
"""
import numpy as np

ALPHA = 1/137.035999084
twosqrta = 2*np.sqrt(ALPHA)          # 0.17085
Ob, Och = 0.0430, 0.2306
Om = Ob + Och                         # 0.2737
fb = Ob/Om                            # 0.1572  Newtonian floor
sqrt_fb = np.sqrt(fb)                 # 0.3964  deep-MOND ceiling
A_data = 0.0686
gate_data = A_data/twosqrta           # 0.4015

print("="*88)
print("CHANNEL 2: LOS-weighted effective gate through CLUMPED structure")
print("="*88)
print(f"  2 sqrt(alpha)            = {twosqrta:.5f}")
print(f"  f_b = Ob/Om (Newt floor) = {fb:.5f}  -> A_floor   = {twosqrta*fb:.5f}")
print(f"  sqrt(f_b) (dMOND ceiling)= {sqrt_fb:.5f}  -> A_ceiling = {twosqrta*sqrt_fb:.5f}")
print(f"  DATA gate                = {gate_data:.5f}  (A_data={A_data})")
print(f"  => need <GATE>_LOS = {gate_data:.4f}, i.e. essentially the deep-MOND CEILING {sqrt_fb:.4f}.")
print(f"     The Newtonian-leaning smooth value ~0.24 must be lifted by +67% to hit it.")

# -------------------------------------------------------------------------------------
# LOCAL GATE as a function of x = g/a0.
# AQUAL/MOND mu(x)=x/(1+x). The baryon-share gate (from the task brief) interpolates the
# EXPONENT p of f_b^p between deep-MOND (p=1/2 -> sqrt) and Newtonian (p=1 -> f_b).
# A standard simple-mu mapping: nu(y)=1/2+sqrt(1/4+1/y) etc. We use the brief's own statement:
#   gate = f_b^p(x),  p(x) interpolates 1/2 (x->0) to 1 (x->inf).
# Concretely with simple-mu the boost weight gives p = 1 - 1/2 * 1/(1+x)  (p:1/2 at x=0, ->1 at inf).
def p_of_x(x):
    return 1.0 - 0.5/(1.0+x)          # x=0 ->0.5 ; x=inf ->1.0 ; matches brief's p~0.72-0.79 at x~3-6
def gate_local(x):
    return fb**p_of_x(x)

for x in [0.0,0.1,0.3,1.0,3.06,5.85,10,100]:
    print(f"     x={x:6.2f}: p={p_of_x(x):.3f}  gate=f_b^p={gate_local(x):.4f}")

# -------------------------------------------------------------------------------------
# x as a function of LOCAL overdensity delta.
# Peculiar acceleration of a structure of overdensity delta on scale R: g_pec ~ G delta rho_m R.
# In MOND units x = g/a0 where a0 = 2 sqrt(alpha) c H. For linear theory the peculiar accel field
# scales as g_pec/cH ~ f(Om) * delta * (something). The CLEAN scaling actually used in the smooth
# case is x_smooth = sqrt(Omega_m(z)) / (2 sqrt(alpha)) = sqrt(<rho>/rho_crit)/(2 sqrt a).
# Generalize to LOCAL density: x(delta) ~ sqrt( rho_local/rho_crit ) /(2 sqrt(alpha))
#   = sqrt( Om*(1+delta) ) /(2 sqrt(alpha))  = x_smooth * sqrt(1+delta).
# KEY CONSEQUENCE: overdensity (delta>0) makes x LARGER (deeper Newtonian), NOT smaller!
#   Only UNDERdensities (voids, delta -> -1) drive x -> 0 (deep-MOND).
x_smooth0 = np.sqrt(Om)/twosqrta
def x_of_delta(delta, z=0.0):
    Om_z = Om*(1+z)**3/(Om*(1+z)**3 + (1-Om))
    return np.sqrt(Om_z*(1+delta))/twosqrta
print("\n" + "-"*88)
print("  x as a function of local overdensity (x = sqrt(Om*(1+delta))/(2 sqrt alpha)):")
print(f"     smooth (delta=0):   x={x_of_delta(0):.3f}  gate={gate_local(x_of_delta(0)):.4f}")
for d in [-0.99,-0.9,-0.5,0,1,5,30,200,1000]:
    x=x_of_delta(d); print(f"     delta={d:7.1f}: x={x:6.3f}  gate={gate_local(x):.4f}")
print("  *** CRITICAL: dense clumps have HIGHER x (more Newtonian, gate->f_b=0.157, LOWER).")
print("      Deep-MOND (gate->0.396) lives in VOIDS (delta->-1), not in clusters. ***")

# -------------------------------------------------------------------------------------
# REALISTIC matter PDF along a random LOS: lognormal in (1+delta) with variance sigma_ln^2.
# Volume-weighted: most of the path is in voids/mean (1+delta<1), tiny fraction in clumps.
# Mass(density)-weighted: weights by (1+delta), pulling toward clumps.
# We test BOTH because the screen-accumulation weight is the real physics question.
def lognormal_pdf_samples(sigma_ln, N=2_000_000, seed=0):
    rng = np.random.default_rng(seed)
    # 1+delta lognormal with <1+delta>=1  => mean of normal = -sigma_ln^2/2
    g = rng.normal(-0.5*sigma_ln**2, sigma_ln, N)
    return np.exp(g)                  # = 1+delta, with ensemble mean 1

print("\n" + "-"*88)
print("  LOS average of the gate over a realistic lognormal matter field:")
print("  sigma_ln set by smoothing scale; z=0 nonlinear sigma_8~0.8 -> sigma_ln(8 Mpc)~0.8,")
print("  but the screen samples the FULL nonlinear field; we scan sigma_ln to bracket.")
print(f"  {'sigma_ln':>8s} {'<gate>_VOL':>11s} {'<gate>_MASS':>11s} {'<gate>_screen':>13s}")
for sigma_ln in [0.5, 0.8, 1.2, 1.8, 2.5, 3.5]:
    opd = lognormal_pdf_samples(sigma_ln)        # 1+delta
    delta = opd - 1.0
    x = np.sqrt(Om*opd)/twosqrta                 # local x
    g = gate_local(x)
    # VOLUME weighting (path length is uniform in physical distance = volume-weighted):
    gate_vol = g.mean()
    # MASS weighting (weight by local density 1+delta):
    gate_mass = np.average(g, weights=opd)
    # SCREEN weighting: the optical screen dpsi ~ trace coupling ~ local matter density (T^mu_mu ~ rho),
    # so the screen ACCUMULATES proportional to (1+delta) -> same as mass weighting; but if the screen
    # saturates / is sourced by the DE fraction (smooth), it is volume-weighted. Report mass as the
    # most clumping-favorable physical choice.
    gate_screen = gate_mass
    print(f"  {sigma_ln:8.2f} {gate_vol:11.4f} {gate_mass:11.4f} {gate_screen:13.4f}")

print("\n  Interpretation:")
print("   - VOLUME-weighted <gate> goes DOWN with clumping: voids (low gate? no — voids have low x")
print("     => HIGH gate ~0.396) vs the fact that mean is preserved. Let's read the actual numbers above.")

# -------------------------------------------------------------------------------------
# Decisive diagnostic: which density regions carry deep-MOND (x<1)?  Need x<1 => Om*(1+delta)<(2sqrt a)^2
x_crit_delta = (twosqrta**2)/Om - 1.0
print("\n" + "-"*88)
print(f"  x<1 (deep-MOND, gate near ceiling) requires 1+delta < (2 sqrt a)^2/Om = {(twosqrta**2)/Om:.4f}")
print(f"     => delta < {x_crit_delta:.4f}  i.e. UNDERDENSE voids only (delta < -0.97).")
print(f"  Fraction of VOLUME that deep-MOND (in lognormal sigma_ln=1.5):")
for sigma_ln in [0.8,1.5,2.5]:
    opd=lognormal_pdf_samples(sigma_ln,seed=1)
    fvol=np.mean(opd < (twosqrta**2)/Om)
    fmass=np.average(opd < (twosqrta**2)/Om, weights=opd)
    print(f"     sigma_ln={sigma_ln}: vol-frac deep-MOND={fvol:.4f}, mass-frac={fmass:.5f}")
print("  => Deep-MOND volume is the EXTREME voids; it carries ~0 mass and a modest volume slice.")
print("     Clumps (where mass is) are MORE Newtonian, gate LOWER. Clumping does NOT lift the gate.")
