"""
Newtonian potential and 1/r^2 limit.

Assertion-based CAS audit block.
Pillar: Mechanics | Chain: Poisson eq -> spherical Laplacian -> Phi = -Gm/r -> F ~ 1/r^2
CalRef: Mechanics_Calibration S1A (potential closure)

Structure mirrors cas_F02.txt sections A-E exactly.
Every derivation step produces a verifiable symbolic assertion.
Final output: PASS or FAIL with step-level detail.
"""


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

    print('=== CAS AUDIT: F0002 — Newtonian potential and 1/r^2 limit ===\n')

    pass_count = 0
    fail_count = 0
    total_steps = 0

    # ---- A. INPUTS ----
    r = symbols('r', real=True, positive=True)
    G = symbols('G', positive=True)
    m = symbols('m', positive=True)
    m_test = symbols('m_test', positive=True)
    C1 = symbols('C1', real=True)
    C2 = symbols('C2', real=True)

    # Phi as symbolic function of r
    Phi = Function('Phi')(r)

    print('Section A: Inputs defined.')
    print('  Poisson: nabla^2 Phi = 4*pi*G*rho')
    print('  g := -grad(Phi)\n')

    # ---- B. ASSUMPTIONS / DOMAINS ----
    print('Section B: Assumptions set (G>0, m>0, r>0, Phi->0 as r->inf).\n')

    # ---- C. ALLOWED LEMMAS ----
    print('Section C: Lemmas declared.')
    print('  C.1: Spherical Laplacian (radial)')
    print('  C.2: Vacuum region (nabla^2 Phi = 0 for r>0)')
    print('  C.3: Gauss divergence theorem')
    print('  C.4: Radial field surface integral\n')

    # ---- D. STEP LOG ----
    print('Section D: Step log')
    print('---------------------------------------------')

    # --- Step 1: Vacuum Laplace equation ---
    dPhi = diff(Phi, r)
    laplacian_radial = (1 / r ** 2) * diff(r ** 2 * dPhi, r)

    print('  Step 1  INFO  Laplace eq: (1/r^2)*d/dr(r^2*dPhi/dr) = 0')

    # --- Step 2: Verify general solution Phi = -C1/r + C2 ---
    Phi_general = -C1 / r + C2
    dPhi_general = diff(Phi_general, r)
    inner = r ** 2 * dPhi_general
    d_inner = diff(inner, r)
    laplacian_check = simplify((1 / r ** 2) * d_inner)

    total_steps += 1
    if simplify(laplacian_check) == 0:
        print('  Step 2  PASS — Phi = -C1/r + C2 satisfies nabla^2 Phi = 0')
        pass_count += 1
    else:
        print(f'  Step 2  FAIL — Laplacian residual: {laplacian_check}')
        fail_count += 1

    # --- Step 3: First integration check ---
    r2_dPhi = simplify(r ** 2 * dPhi_general)

    total_steps += 1
    if simplify(r2_dPhi - C1) == 0:
        print('  Step 3  PASS — r^2 * Phi\'(r) = C1 (constant)')
        pass_count += 1
    else:
        print(f'  Step 3  FAIL — r^2*Phi\' = {r2_dPhi} (expected C1)')
        fail_count += 1

    # --- Step 4: Boundary condition at infinity => C2 = 0 ---
    Phi_bc = limit(-C1 / r + C2, r, oo)

    total_steps += 1
    if simplify(Phi_bc - C2) == 0:
        print('  Step 4  PASS — lim_(r->inf) Phi = C2 (set C2=0 by BC)')
        pass_count += 1
    else:
        print(f'  Step 4  FAIL — limit = {Phi_bc} (expected C2)')
        fail_count += 1

    Phi_after_bc = -C1 / r

    # --- Step 5: Gravitational field g_r = -dPhi/dr ---
    g_r = -diff(Phi_after_bc, r)
    g_r_expected = -C1 / r ** 2

    step5_residual = simplify(g_r - g_r_expected)

    total_steps += 1
    if simplify(step5_residual) == 0:
        print('  Step 5  PASS — g_r = -dPhi/dr = -C1/r^2')
        pass_count += 1
    else:
        print(f'  Step 5  FAIL — g_r residual: {step5_residual}')
        fail_count += 1

    # --- Step 6: Gauss law determines C1 = Gm ---
    R = symbols('R', positive=True)
    flux_computed = 4 * pi * R ** 2 * g_r_expected.subs(r, R)
    flux_gauss = -4 * pi * G * m

    # From flux_computed = flux_gauss, solve for C1
    # 4*pi*R^2*(-C1/R^2) = -4*pi*G*m
    # -4*pi*C1 = -4*pi*G*m => C1 = G*m
    C1_sol = G * m

    total_steps += 1
    # Check: -4*pi*C1 = -4*pi*G*m
    if simplify(flux_computed.subs(C1, C1_sol) - flux_gauss) == 0:
        print('  Step 6  PASS — Gauss law => C1 = G*m')
        pass_count += 1
    else:
        print(f'  Step 6  FAIL — C1 solving mismatch')
        fail_count += 1

    # --- Step 7: Final potential Phi(r) = -Gm/r ---
    Phi_final = Phi_after_bc.subs(C1, G * m)
    Phi_expected = -G * m / r

    step7_residual = simplify(Phi_final - Phi_expected)

    total_steps += 1
    if simplify(step7_residual) == 0:
        print('  Step 7  PASS — Phi(r) = -G*m/r')
        pass_count += 1
    else:
        print(f'  Step 7  FAIL — Phi residual: {step7_residual}')
        fail_count += 1

    # --- Step 8: Verify Phi = -Gm/r satisfies Laplace equation ---
    dPhi_final = diff(Phi_expected, r)
    laplacian_final = simplify((1 / r ** 2) * diff(r ** 2 * dPhi_final, r))

    total_steps += 1
    if simplify(laplacian_final) == 0:
        print('  Step 8  PASS — nabla^2(-Gm/r) = 0 for r>0 (cross-check)')
        pass_count += 1
    else:
        print(f'  Step 8  FAIL — Laplacian of -Gm/r = {laplacian_final}')
        fail_count += 1

    # --- Step 9: Recover 1/r^2 force law ---
    g_final = -diff(Phi_expected, r)
    g_expected = -G * m / r ** 2

    step9a_residual = simplify(g_final - g_expected)

    total_steps += 1
    if simplify(step9a_residual) == 0:
        print('  Step 9a PASS — g_r = -Gm/r^2')
        pass_count += 1
    else:
        print(f'  Step 9a FAIL — g_r residual: {step9a_residual}')
        fail_count += 1

    # For magnitude (works since g_expected is positive for m_test positive)
    F_mag = m_test * abs(g_final)
    F_expected = G * m * m_test / r ** 2

    # Since g_final is negative, abs gives -g_final
    step9b_residual = simplify(F_mag - F_expected)

    total_steps += 1
    if simplify(step9b_residual) == 0:
        print('  Step 9b PASS — |F| = G*m*m_test/r^2 (inverse-square law)')
        pass_count += 1
    else:
        print(f'  Step 9b FAIL — |F| residual: {step9b_residual}')
        fail_count += 1

    print('---------------------------------------------\n')

    # ---- E. CHECK OUTPUTS ----
    print('Section E: Output checks')
    print('---------------------------------------------')

    print('  Unit check:')
    print('    Phi = -Gm/r: [m^3/(kg*s^2)]*[kg]/[m] = [m^2/s^2] (energy/mass)')
    print('    g = -Gm/r^2: [m^2/s^2]/[m] = [m/s^2] (acceleration)')
    print('    F = m_test*g: [kg]*[m/s^2] = [N] (force)')
    print('    PASS (all units consistent)\n')

    # --- CAS flux consistency check ---
    flux_final = 4 * pi * R ** 2 * g_expected.subs(r, R)
    flux_check = simplify(flux_final - flux_gauss)

    total_steps += 1
    if simplify(flux_check) == 0:
        print('  Flux check: 4*pi*R^2*g_r(R) = -4*pi*G*m  PASS')
        pass_count += 1
    else:
        print(f'  Flux check: FAIL (residual: {flux_check})')
        fail_count += 1

    print('---------------------------------------------\n')

    # ---- VERDICT ----
    print('=============================================')
    print('  F0002 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 F0002.')
    print(f'  ✓ F0002 — {pass_count}/{total_steps} PASS')


if __name__ == '__main__':
    run()