"""
Lorentz force from covariant formulation.

Assertion-based CAS audit block.
Pillar: Electromagnetism | Chain: F^{mu nu} -> F^mu_nu -> dp^mu/dtau = qF^mu_nu u^nu
CalRef: EM field tensor -> 3-vector Lorentz force

Structure mirrors cas_F03_F0004.txt sections A-E.
Encodes the field tensor as an explicit 4x4 matrix,
performs index contraction symbolically, and verifies
that the 3-vector Lorentz force emerges.
"""


def run():
    from sympy import symbols, Matrix, simplify, diag

    print('=== CAS AUDIT: F0004 — Lorentz force (fully explicit) ===\n')

    pass_count = 0
    fail_count = 0
    total_steps = 0

    # ---- A. INPUTS ----
    q = symbols('q', real=True)
    m = symbols('m', positive=True)
    c = symbols('c', positive=True)
    gamma_L = symbols('gamma_L', positive=True)
    v1, v2, v3 = symbols('v1 v2 v3', real=True)
    Ex, Ey, Ez = symbols('Ex Ey Ez', real=True)
    Bx, By, Bz = symbols('Bx By Bz', real=True)

    # Metric tensor (Minkowski, signature -+++)
    eta = diag(-1, 1, 1, 1)

    print('Section A: Inputs defined.')
    print('  eta = diag(-1,1,1,1), F^{mu nu} from E,B, u^mu = gamma*(c,v)\n')

    # ---- B. ASSUMPTIONS / DOMAINS ----
    print('Section B: Assumptions set (m>0, q real, |v|<c).\n')

    # ---- C. ALLOWED LEMMAS ----
    print('Section C: Lemmas declared.')
    print('  C.1: Time conversion d/dtau = gamma*d/dt')
    print('  C.2: E_energy = p^0 * c')
    print('  C.3: Mixed tensor via metric\n')

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

    # --- Step 1: Build F^{mu nu} matrix ---
    F_up = Matrix([
        [0, Ex/c, Ey/c, Ez/c],
        [-Ex/c, 0, -Bz, By],
        [-Ey/c, Bz, 0, -Bx],
        [-Ez/c, -By, Bx, 0]
    ])

    # Verify antisymmetry: F^{mu nu} + F^{nu mu} = 0
    antisym_check = simplify(F_up + F_up.T)

    total_steps += 1
    if antisym_check == Matrix(4, 4, lambda i, j: 0):
        print('  Step 1  PASS — F^{mu nu} antisymmetric')
        pass_count += 1
    else:
        print('  Step 1  FAIL — F^{mu nu} not antisymmetric')
        fail_count += 1

    # --- Step 2: Compute mixed tensor F^mu_nu = eta_{nu alpha} F^{mu alpha} ---
    F_mixed = simplify(F_up * eta)

    # Check F^0_0 = 0
    total_steps += 1
    if simplify(F_mixed[0, 0]) == 0:
        print('  Step 2a PASS — F^0_0 = 0')
        pass_count += 1
    else:
        print(f'  Step 2a FAIL — F^0_0 = {F_mixed[0, 0]}')
        fail_count += 1

    # Check F^0_j = E^j/c
    E_vec = Matrix([Ex, Ey, Ez])
    F0j_expected = E_vec / c
    F0j_actual = Matrix([F_mixed[0, 1], F_mixed[0, 2], F_mixed[0, 3]])

    total_steps += 1
    if simplify(F0j_actual - F0j_expected) == Matrix(3, 1, lambda i, j: 0):
        print('  Step 2b PASS — F^0_j = E^j/c')
        pass_count += 1
    else:
        print('  Step 2b FAIL — F^0_j mismatch')
        fail_count += 1

    # Check F^i_0 = E^i/c
    Fi0_expected = E_vec / c
    Fi0_actual = Matrix([F_mixed[1, 0], F_mixed[2, 0], F_mixed[3, 0]])

    total_steps += 1
    if simplify(Fi0_actual - Fi0_expected) == Matrix(3, 1, lambda i, j: 0):
        print('  Step 2c PASS — F^i_0 = E^i/c')
        pass_count += 1
    else:
        print('  Step 2c FAIL — F^i_0 mismatch')
        fail_count += 1

    # --- Step 3: Build four-velocity u^mu ---
    u_up = Matrix([gamma_L * c, gamma_L * v1, gamma_L * v2, gamma_L * v3])

    print('  Step 3  INFO — u^mu = gamma*(c, v1, v2, v3)')

    # --- Step 4: Compute q * F^mu_nu * u^nu for mu=0 (temporal/energy) ---
    # RHS_0 is dot product of row and column vector
    RHS_0 = simplify(q * (F_mixed[0, 0] * u_up[0] + F_mixed[0, 1] * u_up[1] +
                          F_mixed[0, 2] * u_up[2] + F_mixed[0, 3] * u_up[3]))

    E_dot_v = Ex * v1 + Ey * v2 + Ez * v3
    RHS_0_expected = q * gamma_L * E_dot_v / c

    step4_residual = simplify(RHS_0 - RHS_0_expected)

    total_steps += 1
    if simplify(step4_residual) == 0:
        print('  Step 4  PASS — q*F^0_nu*u^nu = q*gamma*(E.v)/c')
        pass_count += 1
    else:
        print(f'  Step 4  FAIL — temporal RHS residual: {step4_residual}')
        fail_count += 1

    # --- Step 5: Time conversion for mu=0 ---
    dEdt_from_dp0 = c * (q * E_dot_v / c)
    dEdt_expected = q * E_dot_v

    step5_residual = simplify(dEdt_from_dp0 - dEdt_expected)

    total_steps += 1
    if simplify(step5_residual) == 0:
        print('  Step 5  PASS — dE/dt = q*(E.v) (energy equation)')
        pass_count += 1
    else:
        print(f'  Step 5  FAIL — energy eq residual: {step5_residual}')
        fail_count += 1

    # --- Step 6: Compute q * F^i_nu * u^nu for mu=i (spatial/momentum) ---
    RHS_spatial = Matrix([
        simplify(q * (F_mixed[1, 0] * u_up[0] + F_mixed[1, 1] * u_up[1] +
                      F_mixed[1, 2] * u_up[2] + F_mixed[1, 3] * u_up[3])),
        simplify(q * (F_mixed[2, 0] * u_up[0] + F_mixed[2, 1] * u_up[1] +
                      F_mixed[2, 2] * u_up[2] + F_mixed[2, 3] * u_up[3])),
        simplify(q * (F_mixed[3, 0] * u_up[0] + F_mixed[3, 1] * u_up[1] +
                      F_mixed[3, 2] * u_up[2] + F_mixed[3, 3] * u_up[3]))
    ])

    # Cross product v x B
    vxB = Matrix([
        v2 * Bz - v3 * By,
        v3 * Bx - v1 * Bz,
        v1 * By - v2 * Bx
    ])

    RHS_spatial_expected = q * gamma_L * (E_vec - vxB)

    step6_residual = simplify(RHS_spatial - RHS_spatial_expected)

    total_steps += 1
    if step6_residual == Matrix(3, 1, lambda i, j: 0):
        print('  Step 6  PASS — q*F^i_nu*u^nu = q*gamma*(E - v x B) [3 components]')
        pass_count += 1
    else:
        print(f'  Step 6  FAIL — spatial RHS residual')
        fail_count += 1

    # --- Step 7: Time conversion for spatial components ---
    force_3vec = RHS_spatial_expected / gamma_L
    force_expected = q * (E_vec - vxB)

    step7_residual = simplify(force_3vec - force_expected)

    total_steps += 1
    if step7_residual == Matrix(3, 1, lambda i, j: 0):
        print('  Step 7  PASS — dp/dt = q*(E - v x B) (Lorentz force, 3-vector)')
        pass_count += 1
    else:
        print(f'  Step 7  FAIL — Lorentz force residual')
        fail_count += 1

    # --- Step 8: Verify F_mixed = F_up * eta independently ---
    F_mixed_check = simplify(F_up * eta)
    step8_residual = simplify(F_mixed - F_mixed_check)

    total_steps += 1
    if step8_residual == Matrix(4, 4, lambda i, j: 0):
        print('  Step 8  PASS — F^mu_nu via explicit index contraction matches matrix multiply')
        pass_count += 1
    else:
        print('  Step 8  FAIL — Mixed tensor cross-check mismatch')
        fail_count += 1

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

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

    print('  Unit check:')
    print('    dp/dt: [kg*m/s]/[s] = [kg*m/s^2] = [N]')
    print('    q*E:   [C]*[V/m] = [C]*[kg*m/(A*s^3)] = [kg*m/s^2] = [N]')
    print('    q*v*B: [C]*[m/s]*[T] = [C]*[m/s]*[kg/(A*s^2)] = [kg*m/s^2] = [N]')
    print('    dE/dt: [J/s] = [W] (power)')
    print('    q*E.v: [C]*[V/m]*[m/s] = [W]')
    print('    PASS (all units consistent)\n')

    # --- Antisymmetry of F^{mu nu} cross-check ---
    trace_F = simplify(F_up.trace())
    total_steps += 1
    if simplify(trace_F) == 0:
        print('  Trace check: tr(F^{mu nu}) = 0  PASS')
        pass_count += 1
    else:
        print(f'  Trace check: FAIL (trace = {trace_F})')
        fail_count += 1

    # --- Consistency: magnetic force does no work ---
    vxB_dot_v = simplify(vxB[0] * v1 + vxB[1] * v2 + vxB[2] * v3)

    total_steps += 1
    if simplify(vxB_dot_v) == 0:
        print('  Consistency: (v x B).v = 0  PASS (magnetic force does no work)')
        pass_count += 1
    else:
        print(f'  Consistency: (v x B).v = {vxB_dot_v}  FAIL')
        fail_count += 1

    # Power consistency: F.v = q*E.v
    F_dot_v = simplify(force_expected[0] * v1 + force_expected[1] * v2 + force_expected[2] * v3)
    power_expected = q * E_dot_v
    step_power = simplify(F_dot_v - power_expected)

    total_steps += 1
    if simplify(step_power) == 0:
        print('  Power check: F.v = q*E.v  PASS')
        pass_count += 1
    else:
        print(f'  Power check: FAIL (residual: {step_power})')
        fail_count += 1

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

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


if __name__ == '__main__':
    run()