"""Reproduce the dimensional checks in the black hole universe study.

Python 3 standard library only. No observational data are fitted.
The bounce is an explicitly assumed homogeneous effective toy model.
"""
from pathlib import Path
import json
import math

G = 6.67430e-11                 # m^3 kg^-1 s^-2
c = 299792458.0                # m s^-1
M_sun = 1.98847e30             # kg, adopted conventional value
pc = 3.085677581491367e16      # m
year = 365.25 * 86400          # s, Julian year
H0_km_s_Mpc = 67.4             # illustrative Planck base-LCDM benchmark
H0 = H0_km_s_Mpc * 1000 / (1e6 * pc)

def schwarzschild(mass_solar):
    M = mass_solar * M_sun
    rs = 2 * G * M / c**2
    density = 3 * M / (4 * math.pi * rs**3)
    t_e1 = 4 * G * M / (3 * c**3)
    t_bound = math.pi * G * M / c**3
    assert math.isclose(density, 3*c**6/(32*math.pi*G**3*M**2), rel_tol=1e-12)
    return dict(mass_solar=mass_solar, rs_m=rs, reference_density_kg_m3=density,
                e1_proper_time_s=t_e1, schwarzschild_upper_bound_s=t_bound)

def bounce(x, w=1.0):
    """x=t/t_b; default massless-scalar/stiff-fluid w=1."""
    if w <= -1:
        raise ValueError('This analytic branch requires w > -1')
    n = 3 * (1+w)
    rho_ratio = 1 / (1+x*x)
    a_ratio = (1+x*x)**(1/n)
    H_tb = 2*x/(n*(1+x*x))
    # (8 pi G rho_b/3)*t_b^2 = 4/n^2 by definition of t_b.
    lhs = H_tb**2
    rhs = 4/n**2 * rho_ratio*(1-rho_ratio)
    assert math.isclose(lhs, rhs, rel_tol=1e-12, abs_tol=1e-14)
    return dict(t_over_tb=x, a_over_ab=a_ratio, rho_over_rhob=rho_ratio, H_times_tb=H_tb)

RH = c/H0
rhoc = 3*H0**2/(8*math.pi*G)
MH = 4*math.pi/3 * rhoc * RH**3
assert math.isclose(2*G*MH/c**2, RH, rel_tol=1e-12)
results = {
    'constants':dict(G_SI=G, c_SI=c, solar_mass_kg=M_sun, parsec_m=pc, julian_year_s=year),
    'schwarzschild_examples':[schwarzschild(m) for m in [10,1e6,1e9]],
    'flat_friedmann_hubble_identity':dict(H0_km_s_Mpc=H0_km_s_Mpc,H0_s_inv=H0,
        hubble_radius_m=RH,hubble_radius_Gly=RH/(1e9*c*year),critical_density_kg_m3=rhoc,
        hubble_energy_equivalent_mass_kg=MH,hubble_mass_solar=MH/M_sun,
        schwarzschild_ratio=2*G*MH/(c*c*RH),matter_only_ratio=0.315),
    'stiff_fluid_bounce_samples':[bounce(x) for x in [-3,-1,0,1,3]],
    'notes':['Schwarzschild proper times are idealized test-clock quantities, not daughter-universe ages.',
        'Reference density uses Euclidean volume at the areal radius; it is not a local GR density.',
        'Hubble mass identity assumes spatially flat FLRW and total mass-equivalent energy density.',
        'Bounce equation is assumed, not derived from collapse or observations.']
}

if __name__ == '__main__':
    out = Path(__file__).parent/'output'
    out.mkdir(exist_ok=True)
    (out/'calculation_results.json').write_text(json.dumps(results, indent=2)+'\n')
    print(json.dumps(results, indent=2))