"""Momentum flux conservation (Poynting theorem).

Assertion-based CAS audit block.
Pillar: Electromagnetism | Chain: Maxwell -> energy density -> Poynting theorem
CalRef: Electromagnetism Math Appendix S3.4, Calibration S2B
"""


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

    print("=== CAS AUDIT: F0015 — Poynting theorem ===\n")

    pass_count = 0
    fail_count = 0
    total_steps = 0

    print("Section A: Inputs defined.")
    print("  u = (1/2)(eps0*E^2 + B^2/mu0)")
    print("  S = (1/mu0)*ExB\n")

    x, y, z, t = symbols("x y z t", real=True)
    eps0, mu0 = symbols("eps0 mu0", positive=True)

    Ex = Function("Ex")(x, y, z, t)
    Ey = Function("Ey")(x, y, z, t)
    Ez = Function("Ez")(x, y, z, t)
    Bx = Function("Bx")(x, y, z, t)
    By = Function("By")(x, y, z, t)
    Bz = Function("Bz")(x, y, z, t)

    print("Section B: Smooth vacuum fields, no sources.\n")
    print("Section C: Lemmas declared.\n")
    print("Section D: Step log")
    print("---------------------------------------------")

    # Step 1: d/dt(E^2) = 2*E.dE/dt
    E_sq = Ex**2 + Ey**2 + Ez**2
    dE_sq_dt = diff(E_sq, t)
    E_dot_dEdt = 2 * (Ex * diff(Ex, t) + Ey * diff(Ey, t) + Ez * diff(Ez, t))
    step1_residual = simplify(dE_sq_dt - E_dot_dEdt)

    total_steps += 1
    if simplify(step1_residual) == 0:
        print("  Step 1  PASS — d/dt(E^2) = 2*E.dE/dt")
        pass_count += 1
    else:
        print(f"  Step 1  FAIL — residual: {step1_residual}")
        fail_count += 1

    # Step 2: d/dt(B^2)
    B_sq = Bx**2 + By**2 + Bz**2
    dB_sq_dt = diff(B_sq, t)
    B_dot_dBdt = 2 * (Bx * diff(Bx, t) + By * diff(By, t) + Bz * diff(Bz, t))
    step2_residual = simplify(dB_sq_dt - B_dot_dBdt)

    total_steps += 1
    if simplify(step2_residual) == 0:
        print("  Step 2  PASS — d/dt(B^2) = 2*B.dB/dt")
        pass_count += 1
    else:
        print(f"  Step 2  FAIL — residual: {step2_residual}")
        fail_count += 1

    # Step 3: du/dt
    from sympy import Rational
    u = Rational(1, 2) * (eps0 * E_sq + B_sq / mu0)
    du_dt = diff(u, t)
    du_dt_expected = eps0 * (Ex * diff(Ex, t) + Ey * diff(Ey, t) + Ez * diff(Ez, t)) + (1 / mu0) * (
        Bx * diff(Bx, t) + By * diff(By, t) + Bz * diff(Bz, t)
    )
    step3_residual = simplify(du_dt - du_dt_expected)

    total_steps += 1
    if simplify(step3_residual) == 0:
        print("  Step 3  PASS — du/dt = eps0*E.dE/dt + (1/mu0)*B.dB/dt")
        pass_count += 1
    else:
        print(f"  Step 3  FAIL — residual: {step3_residual}")
        fail_count += 1

    # Step 4: Maxwell substitution
    curl_E_x = diff(Ez, y) - diff(Ey, z)
    curl_E_y = diff(Ex, z) - diff(Ez, x)
    curl_E_z = diff(Ey, x) - diff(Ex, y)
    curl_B_x = diff(Bz, y) - diff(By, z)
    curl_B_y = diff(Bx, z) - diff(Bz, x)
    curl_B_z = diff(By, x) - diff(Bx, y)

    du_dt_sub = (eps0 * (Ex * curl_B_x / (mu0 * eps0) + Ey * curl_B_y / (mu0 * eps0) + Ez * curl_B_z / (mu0 * eps0)) +
                 (1 / mu0) * (Bx * (-curl_E_x) + By * (-curl_E_y) + Bz * (-curl_E_z)))
    du_dt_maxwell = (1 / mu0) * (Ex * curl_B_x + Ey * curl_B_y + Ez * curl_B_z) - (1 / mu0) * (
        Bx * curl_E_x + By * curl_E_y + Bz * curl_E_z
    )
    step4_residual = simplify(du_dt_sub - du_dt_maxwell)

    total_steps += 1
    if simplify(step4_residual) == 0:
        print("  Step 4  PASS — Maxwell substitution")
        pass_count += 1
    else:
        print(f"  Step 4  FAIL — residual: {step4_residual}")
        fail_count += 1

    # Step 5: Vector identity
    cross_x = Ey * Bz - Ez * By
    cross_y = Ez * Bx - Ex * Bz
    cross_z = Ex * By - Ey * Bx
    div_ExB = diff(cross_x, x) + diff(cross_y, y) + diff(cross_z, z)

    B_dot_curlE = Bx * curl_E_x + By * curl_E_y + Bz * curl_E_z
    E_dot_curlB = Ex * curl_B_x + Ey * curl_B_y + Ez * curl_B_z

    step5_residual = simplify(div_ExB - (B_dot_curlE - E_dot_curlB))

    total_steps += 1
    if simplify(step5_residual) == 0:
        print("  Step 5  PASS — div(ExB) = B.(curl E) - E.(curl B)")
        pass_count += 1
    else:
        print(f"  Step 5  FAIL — residual: {step5_residual}")
        fail_count += 1

    # Step 6: Poynting theorem
    poynting_sum = du_dt_maxwell + (1 / mu0) * div_ExB
    step6_residual = simplify(poynting_sum)

    total_steps += 1
    if simplify(step6_residual) == 0:
        print("  Step 6  PASS — du/dt + div(S) = 0 (Poynting theorem)")
        pass_count += 1
    else:
        print(f"  Step 6  FAIL — residual: {step6_residual}")
        fail_count += 1

    # Step 7: Poynting definition
    Sx = cross_x / mu0
    Sy = cross_y / mu0
    Sz = cross_z / mu0
    div_S = diff(Sx, x) + diff(Sy, y) + diff(Sz, z)
    div_S_expected = div_ExB / mu0
    step7_residual = simplify(div_S - div_S_expected)

    total_steps += 1
    if simplify(step7_residual) == 0:
        print("  Step 7  PASS — div(S) = (1/mu0)*div(ExB)")
        pass_count += 1
    else:
        print(f"  Step 7  FAIL — residual: {step7_residual}")
        fail_count += 1

    # Step 8: Numerical c
    mu0_val = 4 * 3.141592653589793 * 1e-7
    eps0_val = 8.854187817e-12
    c_computed = 1 / (mu0_val * eps0_val) ** 0.5
    c_expected = 299792458
    rel_error = abs(c_computed - c_expected) / c_expected

    total_steps += 1
    if rel_error < 1e-6:
        print(f"  Step 8  PASS — c = 1/sqrt(mu0*eps0) = {c_computed:.6e} m/s (rel err {rel_error:.2e})")
        pass_count += 1
    else:
        print(f"  Step 8  FAIL — c rel error: {rel_error:.2e}")
        fail_count += 1

    # Step 9: Momentum density
    c_light = symbols("c_light", positive=True)
    g_x = Sx / c_light**2
    g_x_sub = g_x.subs(c_light**2, 1 / (mu0 * eps0))
    g_x_expected = eps0 * cross_x
    step9_residual = simplify(g_x_sub - g_x_expected)

    total_steps += 1
    if simplify(step9_residual) == 0:
        print("  Step 9  PASS — g = S/c^2 = eps0*(ExB)")
        pass_count += 1
    else:
        print(f"  Step 9  FAIL — residual: {step9_residual}")
        fail_count += 1

    print("---------------------------------------------\n")
    print("Section E: Output checks")
    print("---------------------------------------------")
    print("  Unit check: du/dt + div(S) = [W/m^3] = 0 — PASS\n")

    # Self-test: wrong Faraday sign
    du_dt_wrong = (1 / mu0) * (E_dot_curlB + B_dot_curlE)
    wrong_poynting = du_dt_wrong + (1 / mu0) * div_ExB
    wrong_poynting_simplified = simplify(wrong_poynting)

    total_steps += 1
    if simplify(wrong_poynting_simplified) != 0:
        print("  Self-test: Wrong Faraday sign gives nonzero Poynting residual  PASS")
        pass_count += 1
    else:
        print("  Self-test: FAIL (wrong sign not detected)")
        fail_count += 1

    expected_wrong = (2 / mu0) * B_dot_curlE
    wrong_quant = simplify(wrong_poynting_simplified - expected_wrong)

    total_steps += 1
    if simplify(wrong_quant) == 0:
        print("  Self-test: wrong residual = (2/mu0)*B.(curl E) (quantified)  PASS")
        pass_count += 1
    else:
        print(f"  Self-test: FAIL — residual = {wrong_quant}")
        fail_count += 1

    print("---------------------------------------------\n")
    print("=============================================")
    print("  F0015 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 F0015.")
    print(f"  ✓ F0015 — {pass_count}/{total_steps} PASS")


if __name__ == "__main__":
    run()