"""
Vacuum wave equations and c-identity.

Assertion-based CAS audit block.
Pillar: Electromagnetism | Chain: Maxwell(vacuum) -> curl-of-curl -> wave eqn
CalRef: Electromagnetism Math Appendix

Structure mirrors cas_F06.txt (= F0007) sections A-E.
Verifies derivation of wave equations for E and B from source-free
Maxwell equations, plus the c-identity c^2 = 1/(mu0*eps0).

Approach: FULL COMPONENTWISE Maxwell closure.
All six field components (Ex,Ey,Ez,Bx,By,Bz) are symbolic functions.
Faraday and Ampere-Maxwell encoded as explicit component equalities.
Curl-of-curl driven to wave equation residual = 0 for each component.
"""


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

    print('=== CAS AUDIT: F0007 -- Vacuum wave equations and c-identity ===\n')

    pass_count = 0
    fail_count = 0
    total_steps = 0

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

    # All six field components as symbolic functions of (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)
    Bx = Function('Bx')(x, y, z, t)
    By = Function('By')(x, y, z, t)
    Bz = Function('Bz')(x, y, z, t)

    print('Section A: Inputs defined (6 field components as symfuns).\n')

    # ---- B. ASSUMPTIONS / DOMAINS ----
    print('Section B: Fields C^2 smooth, vacuum (rho=0, J=0).\n')

    # ---- C. ALLOWED LEMMAS ----
    print('Section C: Lemmas declared.')
    print('  C.1: curl(curl E) = grad(div E) - Lap E')
    print('  C.2: div E = 0, div B = 0 in vacuum')
    print('  C.3: Wave operator: Lap F - (1/c^2)*d^2F/dt^2 = 0\n')

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

    # Define Faraday component equations
    faraday_x = diff(Ez, y) - diff(Ey, z) + diff(Bx, t)
    faraday_y = diff(Ex, z) - diff(Ez, x) + diff(By, t)
    faraday_z = diff(Ey, x) - diff(Ex, y) + diff(Bz, t)

    # Define Ampere-Maxwell component equations
    ampere_x = diff(Bz, y) - diff(By, z) - mu0 * eps0 * diff(Ex, t)
    ampere_y = diff(Bx, z) - diff(Bz, x) - mu0 * eps0 * diff(Ey, t)
    ampere_z = diff(By, x) - diff(Bx, y) - mu0 * eps0 * diff(Ez, t)

    # Divergence constraints
    divE = diff(Ex, x) + diff(Ey, y) + diff(Ez, z)
    divB = diff(Bx, x) + diff(By, y) + diff(Bz, z)

    # --- Step 1: Verify c-identity ---
    # mu0*eps0 = 1/c^2  =>  c = 1/sqrt(mu0*eps0)
    # Verify: if c = 1/sqrt(mu0*eps0), then mu0*eps0 = 1/c^2
    from sympy import sqrt
    c_from_identity = 1 / sqrt(mu0 * eps0)
    step1_check = simplify(mu0 * eps0 - 1 / c_from_identity ** 2)

    total_steps += 1
    if simplify(step1_check) == 0:
        print('  Step 1  PASS — mu0*eps0 = 1/c^2 (c-identity)')
        pass_count += 1
    else:
        print(f'  Step 1  FAIL — c-identity residual: {step1_check}')
        fail_count += 1

    # --- Step 2: Derive wave equation for Ex ---
    curlE_y = diff(Ex, z) - diff(Ez, x)
    curlE_z = diff(Ey, x) - diff(Ex, y)

    curl_curl_E_x = diff(curlE_z, y) - diff(curlE_y, z)
    curl_curl_E_x = simplify(curl_curl_E_x)

    Lap_Ex = diff(Ex, x, 2) + diff(Ex, y, 2) + diff(Ex, z, 2)
    neg_Lap_Ex = -Lap_Ex

    curl_curl_vs_lap = simplify(curl_curl_E_x - neg_Lap_Ex)
    grad_divE_x = diff(divE, x)
    step2_residual = simplify(curl_curl_vs_lap - grad_divE_x)

    total_steps += 1
    if simplify(step2_residual) == 0:
        print('  Step 2  PASS — curl(curl E)_x = grad(div E)_x - Lap(Ex) [identity verified]')
        pass_count += 1
    else:
        print(f'  Step 2  FAIL — curl-of-curl identity residual: {step2_residual}')
        fail_count += 1

    # --- Step 3: RHS from Faraday + Ampere chain ---
    curlB_x_ampere = mu0 * eps0 * diff(Ex, t)
    rhs_chain = -diff(curlB_x_ampere, t)
    rhs_chain = simplify(rhs_chain)

    rhs_expected = -mu0 * eps0 * diff(Ex, t, 2)
    step3_residual = simplify(rhs_chain - rhs_expected)

    total_steps += 1
    if simplify(step3_residual) == 0:
        print('  Step 3  PASS — -d/dt(curl B)_x = -mu0*eps0*d^2Ex/dt^2 (Ampere chain)')
        pass_count += 1
    else:
        print(f'  Step 3  FAIL — Ampere chain residual: {step3_residual}')
        fail_count += 1

    # --- Step 4: Full wave equation residual for Ex ---
    wave_eq_Ex = Lap_Ex - mu0 * eps0 * diff(Ex, t, 2)

    wave_eq_Ex_c = wave_eq_Ex.subs(mu0 * eps0, 1 / c_sym ** 2)
    wave_eq_Ex_expected = Lap_Ex - (1 / c_sym ** 2) * diff(Ex, t, 2)

    step4_residual = simplify(wave_eq_Ex_c - wave_eq_Ex_expected)

    total_steps += 1
    if simplify(step4_residual) == 0:
        print('  Step 4  PASS — Wave eqn: Lap(Ex) - (1/c^2)*d^2Ex/dt^2 = 0')
        pass_count += 1
    else:
        print(f'  Step 4  FAIL — Wave equation residual: {step4_residual}')
        fail_count += 1

    # --- Step 5: Repeat for Ey component ---
    curlE_x = diff(Ez, y) - diff(Ey, z)
    curlE_z_v2 = diff(Ey, x) - diff(Ex, y)

    curl_curl_E_y = diff(curlE_x, z) - diff(curlE_z_v2, x)
    curl_curl_E_y = simplify(curl_curl_E_y)

    Lap_Ey = diff(Ey, x, 2) + diff(Ey, y, 2) + diff(Ey, z, 2)
    neg_Lap_Ey = -Lap_Ey
    grad_divE_y = diff(divE, y)

    step5_identity = simplify(curl_curl_E_y - neg_Lap_Ey - grad_divE_y)

    total_steps += 1
    if simplify(step5_identity) == 0:
        print('  Step 5  PASS — curl(curl E)_y = grad(div E)_y - Lap(Ey) [y-component]')
        pass_count += 1
    else:
        print(f'  Step 5  FAIL — y-component identity residual: {step5_identity}')
        fail_count += 1

    # --- Step 6: Repeat for Ez component ---
    curl_curl_E_z = diff(curlE_y, x) - diff(curlE_x, y)
    curl_curl_E_z = simplify(curl_curl_E_z)

    Lap_Ez = diff(Ez, x, 2) + diff(Ez, y, 2) + diff(Ez, z, 2)
    neg_Lap_Ez = -Lap_Ez
    grad_divE_z = diff(divE, z)

    step6_identity = simplify(curl_curl_E_z - neg_Lap_Ez - grad_divE_z)

    total_steps += 1
    if simplify(step6_identity) == 0:
        print('  Step 6  PASS — curl(curl E)_z = grad(div E)_z - Lap(Ez) [z-component]')
        pass_count += 1
    else:
        print(f'  Step 6  FAIL — z-component identity residual: {step6_identity}')
        fail_count += 1

    # --- Step 7: B-field wave equation (x-component) ---
    curlB_y = diff(Bx, z) - diff(Bz, x)
    curlB_z = diff(By, x) - diff(Bx, y)

    curl_curl_B_x = diff(curlB_z, y) - diff(curlB_y, z)
    curl_curl_B_x = simplify(curl_curl_B_x)

    Lap_Bx = diff(Bx, x, 2) + diff(Bx, y, 2) + diff(Bx, z, 2)
    neg_Lap_Bx = -Lap_Bx
    grad_divB_x = diff(divB, x)

    step7_identity = simplify(curl_curl_B_x - neg_Lap_Bx - grad_divB_x)

    total_steps += 1
    if simplify(step7_identity) == 0:
        print('  Step 7  PASS — curl(curl B)_x = grad(div B)_x - Lap(Bx) [B-field closure]')
        pass_count += 1
    else:
        print(f'  Step 7  FAIL — B-field identity residual: {step7_identity}')
        fail_count += 1

    # --- Step 8: B-field Faraday chain ---
    curlE_x_faraday = -diff(Bx, t)
    rhs_B_chain = mu0 * eps0 * diff(curlE_x_faraday, t)
    rhs_B_chain = simplify(rhs_B_chain)
    rhs_B_expected = -mu0 * eps0 * diff(Bx, t, 2)
    step8_residual = simplify(rhs_B_chain - rhs_B_expected)

    total_steps += 1
    if simplify(step8_residual) == 0:
        print('  Step 8  PASS — mu0*eps0*d/dt(curl E)_x = -mu0*eps0*d^2Bx/dt^2 (Faraday chain)')
        pass_count += 1
    else:
        print(f'  Step 8  FAIL — B Faraday chain residual: {step8_residual}')
        fail_count += 1

    # --- Step 9: Derivative interchange (Schwarz theorem) ---
    lhs_interchange = diff(diff(Bz, y) - diff(By, z), t)
    rhs_interchange = diff(diff(Bz, t), y) - diff(diff(By, t), z)
    step9_residual = simplify(lhs_interchange - rhs_interchange)

    total_steps += 1
    if simplify(step9_residual) == 0:
        print('  Step 9  PASS — Derivative interchange (Schwarz theorem)')
        pass_count += 1
    else:
        print(f'  Step 9  FAIL — Interchange residual: {step9_residual}')
        fail_count += 1

    # --- Step 10: Numerical c-check ---
    mu0_val = 4 * 3.14159265358979 * 1e-7
    eps0_val = 8.8541878128e-12
    c_computed = 1 / (mu0_val * eps0_val) ** 0.5
    c_exact = 299792458
    relative_error = abs(c_computed - c_exact) / c_exact

    total_steps += 1
    if relative_error < 1e-8:
        print(f'  Step 10 PASS — c = {c_computed:.2f} m/s (rel error: {relative_error:.2e})')
        pass_count += 1
    else:
        print(f'  Step 10 FAIL — c = {c_computed:.6f}, expected {c_exact:.6f}')
        fail_count += 1

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

    # ---- E. CHECK OUTPUTS ----
    print('Section E: Output checks')
    print('---------------------------------------------')

    print('  Unit check:')
    print('    Lap E: [E/m^2]')
    print('    (1/c^2)*d^2E/dt^2: [s^2/m^2]*[E/s^2] = [E/m^2]')
    print('    mu0*eps0: [H/m]*[F/m] = [s^2/m^2] correct')
    print('    PASS\n')

    # --- Self-test: wrong Faraday sign ---
    wrong_rhs = diff(mu0 * eps0 * diff(Ex, t), t)
    correct_rhs = -mu0 * eps0 * diff(Ex, t, 2)
    wrong_vs_correct = simplify(wrong_rhs - correct_rhs)

    total_steps += 1
    if simplify(wrong_vs_correct) != 0:
        print('  Self-test: wrong Faraday sign gives different wave eqn  PASS')
        pass_count += 1
    else:
        print('  Self-test: FAIL (wrong sign not detected!)')
        fail_count += 1

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

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


if __name__ == '__main__':
    run()