"""
Relativistic dispersion + mass-energy relation.

Assertion-based CAS audit block.
Pillar: Mechanics (Particle) | Chain: Lorentz invariance -> E^2 = p^2c^2 + m^2c^4 -> E=mc^2
CalRef: Particle Mechanics (relativistic dispersion, VMS mass m = C*tau*ell)

Structure mirrors cas_F03_F0003.txt steps 1-4 + conclusion.
Every derivation step produces a verifiable symbolic assertion.
Final output: PASS or FAIL with step-level detail.
"""


def run():
    from sympy import symbols, sqrt, simplify, series

    print('=== CAS AUDIT: F0003 — Relativistic dispersion + mass-energy ===\n')

    pass_count = 0
    fail_count = 0
    total_steps = 0

    # ---- INPUTS / ASSUMPTIONS ----
    E = symbols('E', real=True)
    p = symbols('p', real=True)
    m = symbols('m', positive=True)
    c = symbols('c', positive=True)

    print('Inputs defined.')
    print('  Four-momentum norm: (E/c)^2 - p^2 = m^2*c^2')
    print('  Symbols: E (energy), p (3-momentum magnitude), m (mass), c (speed of light)\n')

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

    # --- Step 1: Invariant norm of four-momentum ---
    norm_LHS = (E / c) ** 2 - p ** 2
    norm_RHS = m ** 2 * c ** 2

    invariant_relation = norm_LHS - norm_RHS

    print('  Step 1  INFO  Invariant norm: (E/c)^2 - p^2 = m^2*c^2 (input definition)')

    # --- Step 2: Derive relativistic dispersion relation ---
    dispersion_expected = p ** 2 * c ** 2 + m ** 2 * c ** 4

    step2_residual = simplify(E ** 2 - dispersion_expected - c ** 2 * invariant_relation)

    total_steps += 1
    if simplify(step2_residual) == 0:
        print('  Step 2  PASS — E^2 = p^2*c^2 + m^2*c^4 (dispersion relation)')
        pass_count += 1
    else:
        print(f'  Step 2  FAIL — Dispersion residual: {step2_residual}')
        fail_count += 1

    # Step 2b: Cross-check by substituting dispersion into norm
    norm_substituted = norm_LHS.subs(E ** 2, dispersion_expected)
    step2b_residual = simplify(norm_substituted - norm_RHS)

    total_steps += 1
    if simplify(step2b_residual) == 0:
        print('  Step 2b PASS — Substituting dispersion back into norm gives m^2*c^2')
        pass_count += 1
    else:
        print(f'  Step 2b FAIL — Norm back-substitution residual: {step2b_residual}')
        fail_count += 1

    # --- Step 3: Rest frame limit (p = 0) ---
    E_squared_rest = dispersion_expected.subs(p, 0)
    E_squared_rest_expected = m ** 2 * c ** 4

    total_steps += 1
    if simplify(E_squared_rest - E_squared_rest_expected) == 0:
        print('  Step 3a PASS — E^2(p=0) = m^2*c^4')
        pass_count += 1
    else:
        print(f'  Step 3a FAIL — E^2(p=0) = {E_squared_rest}')
        fail_count += 1

    # E = mc^2 (positive root)
    E_rest = sqrt(E_squared_rest_expected)
    E_rest_simplified = simplify(E_rest)
    E_rest_expected = m * c ** 2

    total_steps += 1
    if simplify(E_rest_simplified - E_rest_expected) == 0:
        print('  Step 3b PASS — E(p=0) = m*c^2 (mass-energy equivalence)')
        pass_count += 1
    else:
        print(f'  Step 3b FAIL — E(p=0) = {E_rest_simplified} (expected m*c^2)')
        fail_count += 1

    # --- Step 4: Low-velocity limit (Taylor expansion) ---
    E_full = sqrt(p ** 2 * c ** 2 + m ** 2 * c ** 4)
    E_taylor = series(E_full, p, 0, n=4).removeO()  # up to O(p^4)

    E_taylor_expected = m * c ** 2 + p ** 2 / (2 * m)

    residual_taylor = simplify(E_taylor - E_taylor_expected)

    total_steps += 1
    residual_at_zero = residual_taylor.subs(p, 0)
    if simplify(residual_at_zero) == 0:
        print('  Step 4a PASS — Taylor: E ≈ mc^2 + p^2/(2m) + O(p^4)')
        pass_count += 1
    else:
        print(f'  Step 4a FAIL — Taylor residual at p=0: {residual_at_zero}')
        fail_count += 1

    # Verify kinetic energy: E_kin = E - mc^2 ≈ p^2/(2m)
    E_kin_taylor = simplify(E_taylor - m * c ** 2)
    E_kin_coeff_p2 = simplify(E_kin_taylor.diff(p, 2).subs(p, 0) / 2)

    total_steps += 1
    if simplify(E_kin_coeff_p2 - 1 / (2 * m)) == 0:
        print('  Step 4b PASS — E_kin leading term = p^2/(2m) (classical kinetic energy)')
        pass_count += 1
    else:
        print(f'  Step 4b FAIL — E_kin p^2 coefficient: {E_kin_coeff_p2} (expected 1/(2m))')
        fail_count += 1

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

    # ---- CHECK OUTPUTS ----
    print('Output checks')
    print('---------------------------------------------')

    print('  Unit check:')
    print('    E^2 = p^2*c^2 + m^2*c^4')
    print('    [J^2] = [kg*m/s]^2*[m/s]^2 + [kg]^2*[m/s]^4')
    print('           = [kg^2*m^4/s^4] + [kg^2*m^4/s^4] = [J^2]')
    print('    PASS (both sides [J^2])\n')

    # --- Consistency: dispersion relation is equivalent to norm ---
    consistency = simplify((E ** 2 - p ** 2 * c ** 2 - m ** 2 * c ** 4) - c ** 2 * invariant_relation)
    total_steps += 1
    if simplify(consistency) == 0:
        print('  Consistency: dispersion ↔ norm equivalence  PASS')
        pass_count += 1
    else:
        print(f'  Consistency: FAIL (residual: {consistency})')
        fail_count += 1

    # --- Cross-check: E(p=0) from full dispersion ---
    E_at_rest_cross = simplify(E_full.subs(p, 0))
    total_steps += 1
    if simplify(E_at_rest_cross - m * c ** 2) == 0:
        print('  Cross-check: sqrt(0 + m^2*c^4) = m*c^2  PASS')
        pass_count += 1
    else:
        print(f'  Cross-check: FAIL ({E_at_rest_cross})')
        fail_count += 1

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

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


if __name__ == '__main__':
    run()