"""
Poynting theorem (∂u/∂t + ∇·S + J·E = 0).

Assertion-based CAS audit block.
Pillar: Electromagnetism | Chain: u definition → ∂u/∂t → Maxwell substitution → vector identity → S → theorem
CalRef: Math Appendix §4.2–4.3, EM Calibration §4B
"""


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

    print('=== CAS AUDIT: F0025 — Poynting theorem ===\n')

    pass_count = 0
    fail_count = 0
    total_steps = 0

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

    # E, B, J fields
    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)
    Jx = Function('Jx')(x, y, z, t)
    Jy = Function('Jy')(x, y, z, t)
    Jz = Function('Jz')(x, y, z, t)

    # Energy density
    E_sq = Ex**2 + Ey**2 + Ez**2
    B_sq = Bx**2 + By**2 + Bz**2
    u = (eps0 * E_sq + B_sq / mu0) / 2

    # Poynting vector
    Sx = (Ey*Bz - Ez*By) / mu0
    Sy = (Ez*Bx - Ex*Bz) / mu0
    Sz = (Ex*By - Ey*Bx) / mu0

    # Ohmic dissipation
    JdotE = Jx*Ex + Jy*Ey + Jz*Ez

    print('Section D: Step log')
    print('---------------------------------------------')

    # --- Step 1: ∂u/∂t chain rule ---
    du_dt = diff(u, t)
    du_dt_expected = eps0*(Ex*diff(Ex,t) + Ey*diff(Ey,t) + Ez*diff(Ez,t)) + (Bx*diff(Bx,t) + By*diff(By,t) + Bz*diff(Bz,t))/mu0
    res1 = simplify(du_dt - du_dt_expected)
    total_steps += 1
    if simplify(res1) == 0:
        print('  Step 1  PASS — ∂u/∂t = ε₀(E·∂E/∂t) + (1/μ₀)(B·∂B/∂t)')
        pass_count += 1
    else:
        print('  Step 1  FAIL')
        fail_count += 1

    # --- Step 2-3: Maxwell substitutions are verified structurally ---
    # Faraday: ∂B/∂t = −∇×E
    dBx_dt_faraday = -(diff(Ez, y) - diff(Ey, z))
    dBy_dt_faraday = -(diff(Ex, z) - diff(Ez, x))
    dBz_dt_faraday = -(diff(Ey, x) - diff(Ex, y))

    BdotdBdt_faraday = Bx*dBx_dt_faraday + By*dBy_dt_faraday + Bz*dBz_dt_faraday

    # Ampère–Maxwell: ∂E/∂t = (1/(μ₀ε₀))(∇×B − μ₀J)
    curlB_x = diff(Bz,y) - diff(By,z)
    curlB_y = diff(Bx,z) - diff(Bz,x)
    curlB_z = diff(By,x) - diff(Bx,y)

    dEx_dt_ampere = (curlB_x - mu0*Jx) / (mu0*eps0)
    dEy_dt_ampere = (curlB_y - mu0*Jy) / (mu0*eps0)
    dEz_dt_ampere = (curlB_z - mu0*Jz) / (mu0*eps0)

    EdotdEdt_ampere = Ex*dEx_dt_ampere + Ey*dEy_dt_ampere + Ez*dEz_dt_ampere

    # --- Step 4: du/dt after Maxwell substitution ---
    du_dt_maxwell = eps0*EdotdEdt_ampere + BdotdBdt_faraday/mu0

    curlE_x = diff(Ez,y) - diff(Ey,z)
    curlE_y = diff(Ex,z) - diff(Ez,x)
    curlE_z = diff(Ey,x) - diff(Ex,y)

    EdotcurlB = Ex*curlB_x + Ey*curlB_y + Ez*curlB_z
    BdotcurlE = Bx*curlE_x + By*curlE_y + Bz*curlE_z

    du_dt_step23 = (EdotcurlB - BdotcurlE)/mu0 - JdotE
    res4 = simplify(du_dt_maxwell - du_dt_step23)
    total_steps += 1
    if simplify(res4) == 0:
        print('  Step 4  PASS — ∂u/∂t = (1/μ₀)[E·(∇×B)−B·(∇×E)] − J·E')
        pass_count += 1
    else:
        print('  Step 4  FAIL')
        fail_count += 1

    # --- Step 5: Vector identity ---
    divExB = diff(Ey*Bz - Ez*By, x) + diff(Ez*Bx - Ex*Bz, y) + diff(Ex*By - Ey*Bx, z)
    identity_rhs = BdotcurlE - EdotcurlB
    res5 = simplify(divExB - identity_rhs)
    total_steps += 1
    if simplify(res5) == 0:
        print('  Step 5  PASS — ∇·(E×B) = B·(∇×E) − E·(∇×B) (vector identity)')
        pass_count += 1
    else:
        print('  Step 5  FAIL')
        fail_count += 1

    # --- Step 6: ∇·S from vector identity ---
    divS = diff(Sx, x) + diff(Sy, y) + diff(Sz, z)
    divS_expected = divExB / mu0
    res6 = simplify(divS - divS_expected)
    total_steps += 1
    if simplify(res6) == 0:
        print('  Step 6  PASS — ∇·S = (1/μ₀)∇·(E×B)')
        pass_count += 1
    else:
        print('  Step 6  FAIL')
        fail_count += 1

    # --- Step 7: Poynting theorem assembly ---
    poynting_sum = du_dt_maxwell + divS + JdotE
    poynting_sum_simplified = simplify(poynting_sum)
    total_steps += 1
    if simplify(poynting_sum_simplified) == 0:
        print('  Step 7  PASS — ∂u/∂t + ∇·S + J·E = 0 (Poynting theorem)')
        pass_count += 1
    else:
        print('  Step 7  FAIL')
        fail_count += 1

    # --- Step 8: Source-free limit ---
    poynting_J0 = poynting_sum.subs([(Jx, 0), (Jy, 0), (Jz, 0)])
    res8 = simplify(poynting_J0)
    total_steps += 1
    if simplify(res8) == 0:
        print('  Step 8  PASS — J=0 → ∂u/∂t + ∇·S = 0 (consistent with F0015)')
        pass_count += 1
    else:
        print('  Step 8  FAIL')
        fail_count += 1

    # --- Step 9: Numerical — resistive dissipation ---
    sigma_Cu = 5.96e7
    E_val = 0.01
    JE_val = sigma_Cu * E_val**2
    JE_expected = 5960
    total_steps += 1
    if abs(JE_val - JE_expected) < 1:
        print(f'  Step 9  PASS — Cu dissipation: J·E = {JE_val:.0f} W/m³ (σ=5.96e7, E=0.01 V/m)')
        pass_count += 1
    else:
        print('  Step 9  FAIL')
        fail_count += 1

    # --- Step 10: Self-test — wrong Faraday sign ---
    dBx_wrong = (diff(Ez, y) - diff(Ey, z))
    dBy_wrong = (diff(Ex, z) - diff(Ez, x))
    dBz_wrong = (diff(Ey, x) - diff(Ex, y))

    BdotdBdt_wrong = Bx*dBx_wrong + By*dBy_wrong + Bz*dBz_wrong
    du_dt_wrong = eps0*EdotdEdt_ampere + BdotdBdt_wrong/mu0
    poynting_wrong = du_dt_wrong + divS + JdotE
    poynting_wrong_s = simplify(poynting_wrong)
    total_steps += 1
    if not (simplify(poynting_wrong_s) == 0):
        print('  Step 10 PASS — Wrong Faraday sign (+∇×E) breaks Poynting theorem')
        pass_count += 1
    else:
        print('  Step 10 FAIL')
        fail_count += 1

    # --- Step 11: Self-test — missing J·E term ---
    poynting_noJ = du_dt_maxwell + divS
    res11_check = simplify(poynting_noJ + JdotE)
    total_steps += 1
    if simplify(res11_check) == 0:
        print('  Step 11 PASS — Missing J·E: residual = −J·E (quantified)')
        pass_count += 1
    else:
        print('  Step 11 FAIL')
        fail_count += 1

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

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


if __name__ == '__main__':
    run()