"""
Charge continuity equation (∂ρ/∂t + ∇·J = 0).

Assertion-based CAS audit block.
Pillar: Electromagnetism | Chain: div(Ampère–Maxwell) → vector identity → Gauss substitution → continuity
CalRef: Math Appendix §4.1, EM Calibration §3A
"""


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

    print('=== CAS AUDIT: F0024 — Charge continuity equation ===\n')

    pass_count = 0
    fail_count = 0
    total_steps = 0

    # ---- A. INPUTS ----
    x, y, z, t = symbols('x y z t', real=True)
    mu0, eps0 = symbols('mu0 eps0', real=True, positive=True)

    # Field components
    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)
    rho = Function('rho')(x, y, z, t)

    # Curl of B
    curlB_x = diff(Bz, y) - diff(By, z)
    curlB_y = diff(Bx, z) - diff(Bz, x)
    curlB_z = diff(By, x) - diff(Bx, y)

    # Divergence of curl(B) — should be zero
    div_curlB = diff(curlB_x, x) + diff(curlB_y, y) + diff(curlB_z, z)

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

    # --- Step 1: Vector identity ---
    res1 = simplify(div_curlB)
    total_steps += 1
    if simplify(res1) == 0:
        print('  Step 1  PASS — ∇·(∇×B) = 0 (vector identity)')
        pass_count += 1
    else:
        print('  Step 1  FAIL')
        fail_count += 1

    # --- Step 2: Divergence of Ampère–Maxwell RHS ---
    divJ = diff(Jx, x) + diff(Jy, y) + diff(Jz, z)
    divE = diff(Ex, x) + diff(Ey, y) + diff(Ez, z)
    d_divE_dt = diff(divE, t)

    div_rhs = mu0 * divJ + mu0 * eps0 * d_divE_dt

    step2_expr = simplify(div_rhs / mu0)
    step2_expected = divJ + eps0 * d_divE_dt
    res2 = simplify(step2_expr - step2_expected)
    total_steps += 1
    if simplify(res2) == 0:
        print('  Step 2  PASS — 0 = ∇·J + ε₀ ∂(∇·E)/∂t (after μ₀ division)')
        pass_count += 1
    else:
        print('  Step 2  FAIL')
        fail_count += 1

    # --- Step 3: Gauss substitution ---
    gauss_sub = rho / eps0
    eps0_times_d_gauss_dt = eps0 * diff(gauss_sub, t)
    res3 = simplify(eps0_times_d_gauss_dt - diff(rho, t))
    total_steps += 1
    if simplify(res3) == 0:
        print('  Step 3  PASS — ε₀ · ∂(ρ/ε₀)/∂t = ∂ρ/∂t (Gauss substitution)')
        pass_count += 1
    else:
        print('  Step 3  FAIL')
        fail_count += 1

    # --- Step 4: Final continuity equation ---
    continuity = diff(rho, t) + divJ

    # --- Step 5: Concrete test — uniform charge decay ---
    rho_0, tau = symbols('rho_0 tau', real=True, positive=True)
    rho_c = rho_0 * exp(-t/tau)
    Jx_c = rho_0 * x / (3*tau) * exp(-t/tau)
    Jy_c = rho_0 * y / (3*tau) * exp(-t/tau)
    Jz_c = rho_0 * z / (3*tau) * exp(-t/tau)

    drho_dt_c = diff(rho_c, t)
    divJ_c = diff(Jx_c, x) + diff(Jy_c, y) + diff(Jz_c, z)
    continuity_c = simplify(drho_dt_c + divJ_c)
    total_steps += 1
    if simplify(continuity_c) == 0:
        print('  Step 5  PASS — ρ=ρ₀e^{-t/τ}, J=(ρ₀/3τ)e^{-t/τ}r → ∂ρ/∂t+∇·J=0')
        pass_count += 1
    else:
        print('  Step 5  FAIL')
        fail_count += 1

    # --- Step 6: Static charge ---
    rho_s = symbols('rho_s', real=True, positive=True)
    rho_static = rho_s
    drho_s_dt = 0
    divJ_zero = 0
    res6 = simplify(drho_s_dt + divJ_zero)
    total_steps += 1
    if simplify(res6) == 0:
        print('  Step 6  PASS — Static ρ, J=0 → continuity trivially satisfied')
        pass_count += 1
    else:
        print('  Step 6  FAIL')
        fail_count += 1

    # --- Step 7: Steady current ---
    J0_val = symbols('J0_val', real=True, positive=True)
    divJ_uniform = diff(0, x) + diff(0, y) + diff(J0_val, z)
    res7 = simplify(divJ_uniform)
    total_steps += 1
    if simplify(res7) == 0:
        print('  Step 7  PASS — Steady uniform current: ∇·J = 0, ∂ρ/∂t = 0')
        pass_count += 1
    else:
        print('  Step 7  FAIL')
        fail_count += 1

    # --- Step 8: Cross-block — ε₀ cancellation ---
    cancel_expr = simplify(eps0 * (1/eps0) - 1)
    full_chain = mu0 * (divJ + diff(rho, t))
    reduced = simplify(full_chain / mu0 - (divJ + diff(rho, t)))
    total_steps += 1
    if simplify(cancel_expr) == 0 and simplify(reduced) == 0:
        print('  Step 8  PASS — ε₀ cancellation + μ₀ division chain verified')
        pass_count += 1
    else:
        print('  Step 8  FAIL')
        fail_count += 1

    # --- Step 9: Numerical — capacitor charging ---
    I_val = 1.0
    A_plate = 0.01
    J_plate = I_val / A_plate
    sigma_rate = I_val / A_plate
    total_steps += 1
    if abs(J_plate - sigma_rate) < 1e-15:
        print(f'  Step 9  PASS — Capacitor: J = dσ/dt = {J_plate:.0f} A/m² (I=1A, A=0.01m²)')
        pass_count += 1
    else:
        print('  Step 9  FAIL')
        fail_count += 1

    # --- Step 10: Self-test — wrong sign ---
    wrong_cont = simplify(drho_dt_c - divJ_c)
    total_steps += 1
    if not (simplify(wrong_cont) == 0):
        print('  Step 10a PASS — Wrong sign (∂ρ/∂t − ∇·J) detected as incorrect')
        pass_count += 1
    else:
        print('  Step 10a FAIL')
        fail_count += 1

    wrong_expected = -2 * rho_0 / tau * exp(-t/tau)
    res10b = simplify(wrong_cont - wrong_expected)
    total_steps += 1
    if simplify(res10b) == 0:
        print('  Step 10b PASS — Wrong residual = −2(ρ₀/τ)e^{-t/τ} (quantified)')
        pass_count += 1
    else:
        print('  Step 10b FAIL')
        fail_count += 1

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

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


if __name__ == '__main__':
    run()