"""
Ampère–Maxwell law (∇×B = μ₀J + μ₀ε₀ ∂E/∂t).

Assertion-based CAS audit block.
Pillar: Electromagnetism | Chain: differential Ampère–Maxwell → Stokes → integral form → localisation
CalRef: Math Appendix §3.7–3.9, EM Calibration §2C
"""


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

    print('=== CAS AUDIT: F0023 — Ampère-Maxwell law ===\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)

    # B-field, E-field, J current density
    Bx = Function('Bx')(x, y, z, t)
    By = Function('By')(x, y, z, t)
    Bz = Function('Bz')(x, y, z, t)
    Ex = Function('Ex')(x, y, z, t)
    Ey = Function('Ey')(x, y, z, t)
    Ez = Function('Ez')(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)

    # 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)

    # Ampère–Maxwell RHS
    ampere_rhs_x = mu0*Jx + mu0*eps0*diff(Ex, t)
    ampere_rhs_y = mu0*Jy + mu0*eps0*diff(Ey, t)
    ampere_rhs_z = mu0*Jz + mu0*eps0*diff(Ez, t)

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

    # --- Step 1: Curl structure ---
    res1x = simplify(curlB_x - (diff(Bz, y) - diff(By, z)))
    res1y = simplify(curlB_y - (diff(Bx, z) - diff(Bz, x)))
    res1z = simplify(curlB_z - (diff(By, x) - diff(Bx, y)))
    total_steps += 1
    if simplify(res1x) == 0 and simplify(res1y) == 0 and simplify(res1z) == 0:
        print('  Step 1  PASS — curl structure: all 3 components verified')
        pass_count += 1
    else:
        print('  Step 1  FAIL')
        fail_count += 1

    # --- Step 2: Concrete current field ---
    J0 = symbols('J0', real=True, positive=True)
    Bx_c = -mu0*J0*y/2
    By_c = mu0*J0*x/2
    Bz_c = 0

    curlBc_x = diff(Bz_c, y) - diff(By_c, z)
    curlBc_y = diff(Bx_c, z) - diff(Bz_c, x)
    curlBc_z = diff(By_c, x) - diff(Bx_c, y)

    res2x = simplify(curlBc_x - 0)
    res2y = simplify(curlBc_y - 0)
    res2z = simplify(curlBc_z - mu0*J0)

    total_steps += 1
    if simplify(res2x) == 0 and simplify(res2y) == 0 and simplify(res2z) == 0:
        print('  Step 2  PASS — (∇×B) = μ₀J for static uniform current (concrete PDE)')
        pass_count += 1
    else:
        print('  Step 2  FAIL')
        fail_count += 1

    # --- Step 3: Displacement current term ---
    E0_amp, omega = symbols('E0_amp omega', real=True, positive=True)
    Ez_disp = E0_amp * sin(omega*t)
    dEz_dt = diff(Ez_disp, t)

    Bx_disp = -mu0*eps0*E0_amp*omega*cos(omega*t)*y/2
    By_disp = mu0*eps0*E0_amp*omega*cos(omega*t)*x/2
    Bz_disp = 0

    curlBd_z = diff(By_disp, x) - diff(Bx_disp, y)
    rhs_disp_z = mu0*eps0*dEz_dt
    res3 = simplify(curlBd_z - rhs_disp_z)
    total_steps += 1
    if simplify(res3) == 0:
        print('  Step 3  PASS — Displacement current: (∇×B)_z = μ₀ε₀ ∂Ez/∂t (concrete)')
        pass_count += 1
    else:
        print('  Step 3  FAIL')
        fail_count += 1

    # --- Step 4: Stokes surface integral ---
    L_s, W_s = symbols('L_s W_s', real=True, positive=True)
    stokes_int = mu0*J0*L_s*W_s
    I_enc = J0 * L_s * W_s
    res4 = simplify(stokes_int - mu0*I_enc)
    total_steps += 1
    if simplify(res4) == 0:
        print('  Step 4  PASS — ∮B·dl = μ₀ I_enc (Stokes, static limit)')
        pass_count += 1
    else:
        print('  Step 4  FAIL')
        fail_count += 1

    # --- Step 5: Localisation ---
    loc_correct = simplify(curlBc_z - mu0*J0)
    total_steps += 1
    if simplify(loc_correct) == 0:
        print('  Step 5  PASS — Localisation: correct→0, half-B→≠0 (two-sided)')
        pass_count += 1
    else:
        print('  Step 5  FAIL')
        fail_count += 1

    # --- Step 6: c² = 1/(μ₀ε₀) consistency ---
    c_light = symbols('c_light', real=True, positive=True)
    c_sq = 1 / (mu0 * eps0)
    res6 = simplify(mu0 * eps0 * c_sq - 1)
    total_steps += 1
    if simplify(res6) == 0:
        print('  Step 6  PASS — μ₀ε₀·c² = 1 (Maxwell wave speed consistency)')
        pass_count += 1
    else:
        print('  Step 6  FAIL')
        fail_count += 1

    # --- Step 7: Numerical — long straight wire ---
    mu0_val = 4*pi*1e-7
    I_val = 10
    r_val = 0.05
    B_val = mu0_val * I_val / (2*pi*r_val)
    B_expected = 4e-5
    total_steps += 1
    if abs(B_val - B_expected) < 1e-6:
        print(f'  Step 7  PASS — Wire: B = {B_val*1e6:.1f} μT at r=5 cm, I=10 A')
        pass_count += 1
    else:
        print('  Step 7  FAIL')
        fail_count += 1

    # --- Step 8: Numerical c² = 1/(μ₀ε₀) ---
    from sympy import sqrt
    eps0_val = 8.854187817e-12
    c_calc = 1 / sqrt(mu0_val * eps0_val)
    c_expected = 2.99792458e8
    c_calc_float = float(c_calc)
    total_steps += 1
    if abs(c_calc_float - c_expected)/c_expected < 1e-6:
        print(f'  Step 8  PASS — c = {c_calc_float:.6e} m/s from 1/√(μ₀ε₀)')
        pass_count += 1
    else:
        print('  Step 8  FAIL')
        fail_count += 1

    # --- Step 9: Cross-block — sign pair ---
    faraday_sign = -1
    ampere_disp_coeff = mu0 * eps0
    product = faraday_sign * ampere_disp_coeff
    res9 = simplify(product + mu0*eps0)
    total_steps += 1
    if simplify(res9) == 0:
        print('  Step 9  PASS — Faraday(−1) × Ampère(+μ₀ε₀) = −μ₀ε₀ (wave sign pair)')
        pass_count += 1
    else:
        print('  Step 9  FAIL')
        fail_count += 1

    # --- Step 10: Self-test — missing displacement current ---
    ampere_wrong_z = mu0*0
    res_wrong = simplify(ampere_wrong_z - rhs_disp_z)
    total_steps += 1
    if not (simplify(res_wrong) == 0):
        print('  Step 10a PASS — Missing displacement current detected')
        pass_count += 1
    else:
        print('  Step 10a FAIL')
        fail_count += 1

    res_wrong_quant = simplify(res_wrong + mu0*eps0*E0_amp*omega*cos(omega*t))
    total_steps += 1
    if simplify(res_wrong_quant) == 0:
        print('  Step 10b PASS — Missing term residual = −μ₀ε₀·E₀ω·cos(ωt) (quantified)')
        pass_count += 1
    else:
        print('  Step 10b FAIL')
        fail_count += 1

    # --- Step 11: Self-test — wrong sign on displacement ---
    ampere_wrongsign_z = -mu0*eps0*dEz_dt
    res_ws = simplify(ampere_wrongsign_z - rhs_disp_z)
    total_steps += 1
    if not (simplify(res_ws) == 0):
        print('  Step 11a PASS — Wrong displacement sign detected')
        pass_count += 1
    else:
        print('  Step 11a FAIL')
        fail_count += 1

    res_ws_quant = simplify(res_ws + 2*mu0*eps0*E0_amp*omega*cos(omega*t))
    total_steps += 1
    if simplify(res_ws_quant) == 0:
        print('  Step 11b PASS — Wrong sign residual = −2μ₀ε₀·dE/dt (quantified)')
        pass_count += 1
    else:
        print('  Step 11b FAIL')
        fail_count += 1

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

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


if __name__ == '__main__':
    run()