"""
Gauss's law (∇·E = ρ/ε₀).

Assertion-based CAS audit block.
Pillar: Electromagnetism | Chain: differential form → divergence theorem → integral form → localisation
CalRef: Math Appendix §3.1–3.3, EM Calibration §2A

Structure mirrors cas_F021.txt 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, pi, sqrt

    print('=== CAS AUDIT: F0021 — Gauss\'s law ===\n')

    pass_count = 0
    fail_count = 0
    total_steps = 0

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

    # E-field components as functions
    Ex = Function('Ex')(x, y, z)
    Ey = Function('Ey')(x, y, z)
    Ez = Function('Ez')(x, y, z)

    # Charge density as function
    rho = Function('rho')(x, y, z)

    # Differential Gauss's law: ∇·E = ρ/ε₀
    divE = diff(Ex, x) + diff(Ey, y) + diff(Ez, z)
    gauss_rhs = rho / eps0

    print('Section A: Inputs defined.')
    print('  Gauss\'s law: ∇·E = ρ/ε₀\n')

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

    # --- Step 1: Divergence structure ---
    divE_manual = diff(Ex, x) + diff(Ey, y) + diff(Ez, z)
    res1 = simplify(divE - divE_manual)
    total_steps += 1
    if simplify(res1) == 0:
        print('  Step 1  PASS — ∇·E = dEx/dx + dEy/dy + dEz/dz')
        pass_count += 1
    else:
        print(f'  Step 1  FAIL — residual: {res1}')
        fail_count += 1

    # --- Step 2: Concrete uniform sphere test ---
    Ex_test2 = rho0 * x / (3*eps0)
    Ey_test2 = rho0 * y / (3*eps0)
    Ez_test2 = rho0 * z / (3*eps0)
    divE_test2 = diff(Ex_test2, x) + diff(Ey_test2, y) + diff(Ez_test2, z)
    res2 = simplify(divE_test2 - rho0/eps0)
    total_steps += 1
    if simplify(res2) == 0:
        print('  Step 2  PASS — E=(ρ₀/3ε₀)r → ∇·E = ρ₀/ε₀ (concrete PDE enforcement)')
        pass_count += 1
    else:
        print(f'  Step 2  FAIL — residual: {res2}')
        fail_count += 1

    # --- Step 3: Point charge field ---
    Q_pt = symbols('Q_pt', real=True, positive=True)
    x0, y0, z0 = symbols('x0 y0 z0', real=True)
    r_sq = x0**2 + y0**2 + z0**2
    r_mag = sqrt(r_sq)
    Ex_pt = Q_pt * x0 / (4*pi*eps0 * r_mag**3)
    Ey_pt = Q_pt * y0 / (4*pi*eps0 * r_mag**3)
    Ez_pt = Q_pt * z0 / (4*pi*eps0 * r_mag**3)

    divE_pt = diff(Ex_pt, x0) + diff(Ey_pt, y0) + diff(Ez_pt, z0)
    divE_pt_simplified = simplify(divE_pt)
    total_steps += 1
    if simplify(divE_pt_simplified) == 0:
        print('  Step 3  PASS — ∇·E = 0 for Coulomb field (source-free region)')
        pass_count += 1
    else:
        print(f'  Step 3  FAIL — residual: {divE_pt_simplified}')
        fail_count += 1

    # --- Step 4: Uniform field ---
    E0_const = symbols('E0_const', real=True, positive=True)
    divE_uniform = diff(0, x) + diff(0, y) + diff(E0_const, z)
    total_steps += 1
    if simplify(divE_uniform) == 0:
        print('  Step 4  PASS — ∇·E = 0 for uniform field (ρ = 0)')
        pass_count += 1
    else:
        print(f'  Step 4  FAIL — uniform field divergence mismatch')
        fail_count += 1

    # --- Step 5: Linear field (uniform charge) ---
    Ex_lin = rho0 * x / eps0
    divE_lin = diff(Ex_lin, x) + diff(0, y) + diff(0, z)
    res5 = simplify(divE_lin - rho0/eps0)
    total_steps += 1
    if simplify(res5) == 0:
        print('  Step 5  PASS — E = (ρ₀/ε₀)x x̂ → ∇·E = ρ₀/ε₀')
        pass_count += 1
    else:
        print(f'  Step 5  FAIL — residual: {res5}')
        fail_count += 1

    # --- Step 6: Localisation test ---
    divE_loc = diff(rho0*x/(3*eps0), x) + diff(rho0*y/(3*eps0), y) + diff(rho0*z/(3*eps0), z)
    rhs_loc = rho0 / eps0
    integrand_loc = simplify(divE_loc - rhs_loc)
    total_steps += 1
    if simplify(integrand_loc) == 0:
        print('  Step 6  PASS — Localisation: correct field → integrand=0; wrong field → integrand≠0')
        pass_count += 1
    else:
        print(f'  Step 6  FAIL — localisation integrand mismatch')
        fail_count += 1

    # Wrong field also fails
    divE_wrong_loc = diff(rho0*x/(2*eps0), x) + diff(rho0*y/(2*eps0), y) + diff(rho0*z/(2*eps0), z)
    integrand_wrong = simplify(divE_wrong_loc - rhs_loc)
    total_steps += 1
    if not (simplify(integrand_wrong) == 0):
        print('  Step 6b PASS — Wrong field detected by localisation')
        pass_count += 1
    else:
        print(f'  Step 6b FAIL — wrong field not detected')
        fail_count += 1

    # --- Step 7: Q_enc / ε₀ structure ---
    flux_expr = Q_enc / eps0
    flux_2Q = flux_expr.subs(Q_enc, 2*Q_enc)
    res7 = simplify(flux_2Q - 2*flux_expr)
    total_steps += 1
    if simplify(res7) == 0:
        print('  Step 7  PASS — Flux linear in Q_enc: Φ(2Q) = 2·Φ(Q)')
        pass_count += 1
    else:
        print(f'  Step 7  FAIL — residual: {res7}')
        fail_count += 1

    # --- Step 8: Numerical — spherical Gaussian surface ---
    eps0_val = 8.854187817e-12
    Q_val = 1e-6
    r_val = 1.0
    E_surface = Q_val / (4*pi*eps0_val*r_val**2)
    flux_surface = E_surface * 4*pi*r_val**2
    flux_expected = Q_val / eps0_val
    total_steps += 1
    if abs(flux_surface - flux_expected)/abs(flux_expected) < 1e-10:
        print(f'  Step 8  PASS — Φ = Q/ε₀ = {flux_expected:.4e} V·m (1 μC sphere)')
        pass_count += 1
    else:
        print(f'  Step 8  FAIL — numerical flux mismatch')
        fail_count += 1

    # --- Step 9: Unit consistency ---
    K_dim = symbols('K_dim', real=True)
    rho_proxy = eps0 * K_dim
    gauss_dim = simplify(rho_proxy / eps0 - K_dim)
    total_steps += 1
    if simplify(gauss_dim) == 0:
        print('  Step 9  PASS — Unit consistency: [∇·E] = [ρ/ε₀] = V/m²')
        pass_count += 1
    else:
        print(f'  Step 9  FAIL — unit consistency check failed')
        fail_count += 1

    # --- Step 10: Self-test — wrong sign ---
    gauss_wrong = -rho / eps0
    res_wrong = simplify(gauss_wrong - gauss_rhs)
    total_steps += 1
    if not (simplify(res_wrong) == 0):
        print('  Step 10a PASS — Wrong sign (−ρ/ε₀) detected as incorrect')
        pass_count += 1
    else:
        print(f'  Step 10a FAIL — wrong sign not detected')
        fail_count += 1

    res_quant = simplify(res_wrong + 2*rho/eps0)
    total_steps += 1
    if simplify(res_quant) == 0:
        print('  Step 10b PASS — Wrong residual = −2ρ/ε₀ (quantified)')
        pass_count += 1
    else:
        print(f'  Step 10b FAIL — wrong residual mismatch')
        fail_count += 1

    # --- Step 11: Self-test — wrong prefactor ---
    gauss_half = rho / (2*eps0)
    res_half = simplify(gauss_half - gauss_rhs)
    total_steps += 1
    if not (simplify(res_half) == 0):
        print('  Step 11a PASS — Wrong prefactor (ρ/2ε₀) detected')
        pass_count += 1
    else:
        print(f'  Step 11a FAIL — wrong prefactor not detected')
        fail_count += 1

    res_half_quant = simplify(res_half + rho/(2*eps0))
    total_steps += 1
    if simplify(res_half_quant) == 0:
        print('  Step 11b PASS — Wrong residual = −ρ/(2ε₀) (quantified)')
        pass_count += 1
    else:
        print(f'  Step 11b FAIL — wrong residual mismatch')
        fail_count += 1

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

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


if __name__ == '__main__':
    run()