"""Gravitational energy density (rho_g = -(grad Phi)^2/(8*pi*G)).

Assertion-based CAS audit block.
Pillar: Mechanics | Chain: Poisson -> potential energy -> IBP -> energy density
CalRef: Mathematical Bridge S8A, Mechanics Calibration S3C
"""


def run():
    from sympy import symbols, Function, diff, simplify, integrate, pi, oo

    print("=== CAS AUDIT: F0012 — Gravitational energy density ===\n")

    pass_count = 0
    fail_count = 0
    total_steps = 0

    print("Section A: Inputs defined.")
    print("  Phi(x,y,z), G, Poisson equation, potential energy integral\n")
    print("Section B: Static Newtonian gravity, fields smooth, vanish at infinity.\n")
    print("Section C: Lemmas declared.")
    print("  C.1: Poisson: Lap(Phi) = 4*pi*G*rho")
    print("  C.2: IBP: integral (Lap Phi)*Phi = -integral |grad Phi|^2 (surface terms drop)\n")
    print("Section D: Step log")
    print("---------------------------------------------")

    G_grav = symbols("G_grav", positive=True)
    x, y, z = symbols("x y z", real=True)
    Phi = Function("Phi")(x, y, z)

    # Step 1: Laplacian
    Lap_Phi = diff(Phi, x, 2) + diff(Phi, y, 2) + diff(Phi, z, 2)
    grad_Phi_x = diff(Phi, x)
    grad_Phi_y = diff(Phi, y)
    grad_Phi_z = diff(Phi, z)
    grad_Phi_sq = grad_Phi_x**2 + grad_Phi_y**2 + grad_Phi_z**2

    Phi_xx = diff(Phi, x, 2)
    Phi_yy = diff(Phi, y, 2)
    Phi_zz = diff(Phi, z, 2)
    step1_residual = simplify(Lap_Phi - (Phi_xx + Phi_yy + Phi_zz))

    total_steps += 1
    if simplify(step1_residual) == 0:
        print("  Step 1  PASS — Lap(Phi) = Phi_xx + Phi_yy + Phi_zz")
        pass_count += 1
    else:
        print(f"  Step 1  FAIL — Laplacian residual: {step1_residual}")
        fail_count += 1

    # Step 2: Coefficient chain
    from sympy import Rational
    coeff_U = Rational(1, 2) * Rational(1, 4) / (pi * G_grav)
    expected_coeff = Rational(1, 8) / (pi * G_grav)
    step2_residual = simplify(coeff_U - expected_coeff)

    total_steps += 1
    if simplify(step2_residual) == 0:
        print("  Step 2  PASS — (1/2)*(1/(4*pi*G)) = 1/(8*pi*G)")
        pass_count += 1
    else:
        print(f"  Step 2  FAIL — Coefficient residual: {step2_residual}")
        fail_count += 1

    # Step 3: IBP identity
    div_Phi_gradPhi = diff(Phi * grad_Phi_x, x) + diff(Phi * grad_Phi_y, y) + diff(Phi * grad_Phi_z, z)
    step3_residual = simplify(div_Phi_gradPhi - grad_Phi_sq - Phi * Lap_Phi)

    total_steps += 1
    if simplify(step3_residual) == 0:
        print("  Step 3  PASS — div(Phi*grad Phi) = |grad Phi|^2 + Phi*Lap(Phi)")
        pass_count += 1
    else:
        print(f"  Step 3  FAIL — IBP identity residual: {step3_residual}")
        fail_count += 1

    # Step 4: Rearranged IBP
    Phi_Lap_Phi = Phi * Lap_Phi
    rearranged = div_Phi_gradPhi - grad_Phi_sq
    step4_residual = simplify(Phi_Lap_Phi - rearranged)

    total_steps += 1
    if simplify(step4_residual) == 0:
        print("  Step 4  PASS — Phi*Lap(Phi) = div(Phi*grad Phi) - |grad Phi|^2")
        pass_count += 1
    else:
        print(f"  Step 4  FAIL — Rearrangement residual: {step4_residual}")
        fail_count += 1

    # Step 5: Energy density
    rho_g_expr = -grad_Phi_sq / (8 * pi * G_grav)
    step5_check = simplify(rho_g_expr * (8 * pi * G_grav) + grad_Phi_sq)

    total_steps += 1
    if simplify(step5_check) == 0:
        print("  Step 5  PASS — rho_g = -|grad Phi|^2/(8*pi*G)")
        pass_count += 1
    else:
        print(f"  Step 5  FAIL — Energy density residual: {step5_check}")
        fail_count += 1

    # Step 6: Sign check
    grad_sq_val = symbols("grad_sq_val", positive=True)
    rho_g_sign_test = -grad_sq_val / (8 * pi * G_grav)

    total_steps += 1
    if simplify(rho_g_sign_test) < 0:
        print("  Step 6  PASS — rho_g < 0 when |grad Phi|^2 > 0 (binding energy negative)")
        pass_count += 1
    else:
        print("  Step 6  FAIL — Sign check")
        fail_count += 1

    # Step 7: Concrete point mass
    M_mass = symbols("M_mass", positive=True)
    r_rad = symbols("r_rad", positive=True)
    Phi_point = -G_grav * M_mass / r_rad
    dPhi_dr = diff(Phi_point, r_rad)
    grad_Phi_sq_point = dPhi_dr**2
    rho_g_point = simplify(-grad_Phi_sq_point / (8 * pi * G_grav))
    rho_g_expected = -G_grav * M_mass**2 / (8 * pi * r_rad**4)
    step7_residual = simplify(rho_g_point - rho_g_expected)

    total_steps += 1
    if simplify(step7_residual) == 0:
        print("  Step 7  PASS — Point mass: rho_g = -G*M^2/(8*pi*r^4)")
        pass_count += 1
    else:
        print(f"  Step 7  FAIL — Point mass residual: {step7_residual}")
        fail_count += 1

    # Step 8: Point mass total energy integral
    r_min = symbols("r_min", positive=True)
    integrand_shell = rho_g_point * 4 * pi * r_rad**2
    integrand_shell = simplify(integrand_shell)
    U_total = integrate(integrand_shell, (r_rad, r_min, oo))
    U_expected = -G_grav * M_mass**2 / (2 * r_min)
    step8_residual = simplify(U_total - U_expected)

    total_steps += 1
    if simplify(step8_residual) == 0:
        print("  Step 8  PASS — Total energy: U = -G*M^2/(2*r_min)")
        pass_count += 1
    else:
        print(f"  Step 8  FAIL — Total energy residual: {step8_residual}")
        fail_count += 1

    # Step 9: Numerical Earth
    G_val = 6.67430e-11
    M_earth = 5.972e24
    R_earth = 6.371e6
    rho_g_surface = -G_val * M_earth**2 / (8 * 3.141592653589793 * R_earth**4)
    g_surface = G_val * M_earth / R_earth**2
    rho_g_from_g = -g_surface**2 / (8 * 3.141592653589793 * G_val)
    rel_error = abs(rho_g_surface - rho_g_from_g) / abs(rho_g_surface)

    total_steps += 1
    if rel_error < 1e-10:
        print(f"  Step 9  PASS — Numerical: rho_g(Earth surface) = {rho_g_surface:.4e} J/m^3")
        pass_count += 1
    else:
        print(f"  Step 9  FAIL — Numerical disagreement: rel error = {rel_error:.2e}")
        fail_count += 1

    print("---------------------------------------------\n")
    print("Section E: Output checks")
    print("---------------------------------------------")
    print("  Unit check: [J/m^3] — PASS\n")

    # Self-test: wrong sign
    rho_g_wrong = grad_Phi_sq_point / (8 * pi * G_grav)
    wrong_residual = simplify(rho_g_wrong - rho_g_point)

    total_steps += 1
    if simplify(wrong_residual) != 0:
        print("  Self-test: wrong sign (+) detected as different  PASS")
        pass_count += 1
    else:
        print("  Self-test: FAIL (wrong sign not detected)")
        fail_count += 1

    expected_wrong = 2 * grad_Phi_sq_point / (8 * pi * G_grav)
    wrong_quant = simplify(wrong_residual - expected_wrong)

    total_steps += 1
    if simplify(wrong_quant) == 0:
        print("  Self-test: wrong - correct = 2*|grad Phi|^2/(8*pi*G) (quantified)  PASS")
        pass_count += 1
    else:
        print(f"  Self-test: FAIL (wrong residual = {wrong_residual})")
        fail_count += 1

    print("---------------------------------------------\n")
    print("=============================================")
    print("  F0012 AUDIT RESULT")
    print(f"  Steps: {total_steps}  |  Pass: {pass_count}  |  Fail: {fail_count}")
    if fail_count == 0:
        print("  STATUS: *** PASS ***")
    else:
        print(f"  STATUS: *** FAIL *** ({fail_count} step(s) failed)")
    print("=============================================")
    print("Audit complete for F0012.")
    print(f"  ✓ F0012 — {pass_count}/{total_steps} PASS")


if __name__ == "__main__":
    run()