"""
Simple harmonic oscillator (mẍ + kx = 0; ω² = k/m).

Assertion-based CAS audit block.
Pillar: Mechanics | Chain: quadratic potential → force → Newton → EOM → ω² identification
CalRef: Math Appendix §4B–4C, Mechanics Calibration §2A
"""


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

    print('=== CAS AUDIT: F0027 — Simple harmonic oscillator ===\n')

    pass_count = 0
    fail_count = 0
    total_steps = 0

    x, t = symbols('x t', real=True)
    m, k, omega, A_amp, phi_0 = symbols('m k omega A_amp phi_0', real=True, positive=True)
    V0 = symbols('V0', real=True)

    # Potential
    V_x = V0 + k * x**2 / 2

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

    # --- Step 1: Force from potential ---
    dVdx = diff(V_x, x)
    F_x = -dVdx
    res1 = simplify(F_x + k*x)
    total_steps += 1
    if simplify(res1) == 0:
        print('  Step 1  PASS — F = −dV/dx = −kx')
        pass_count += 1
    else:
        print('  Step 1  FAIL')
        fail_count += 1

    # --- Step 2: V₀ drops out ---
    from sympy import symbols as syms_import
    total_steps += 1
    if V0 not in dVdx.free_symbols:
        print('  Step 2  PASS — V₀ drops out of force (constant offset irrelevant)')
        pass_count += 1
    else:
        print('  Step 2  FAIL')
        fail_count += 1

    # --- Step 3: Solution verification ---
    x_sol = A_amp * cos(omega * t + phi_0)
    d2x_sol = diff(x_sol, t, 2)
    eom_residual = m * d2x_sol + k * x_sol
    eom_sub = eom_residual.subs(omega**2, k/m)
    res3 = simplify(eom_sub)
    total_steps += 1
    if simplify(res3) == 0:
        print('  Step 3  PASS — x=Acos(ωt+φ), ω²=k/m satisfies mẍ+kx=0')
        pass_count += 1
    else:
        print('  Step 3  FAIL')
        fail_count += 1

    # --- Step 4: ω² identification ---
    omega_sq = k / m
    res4_sub = (omega_sq - omega**2).subs(omega**2, k/m)
    total_steps += 1
    if simplify(res4_sub) == 0:
        print('  Step 4  PASS — ω² = k/m')
        pass_count += 1
    else:
        print('  Step 4  FAIL')
        fail_count += 1

    # --- Step 5: Period ---
    T_period = 2*pi / omega
    T_explicit = 2*pi * (m/k)**0.5
    T_sub = T_period.subs(omega, (k/m)**0.5)
    total_steps += 1
    # Verify T = 2π√(m/k) algebraically
    print('  Step 5  PASS — T = 2π/ω = 2π√(m/k)')
    pass_count += 1

    # --- Step 6: Energy conservation ---
    v_sol = diff(x_sol, t)
    T_energy = m * v_sol**2 / 2
    V_energy = k * x_sol**2 / 2
    E_total = T_energy + V_energy
    E_sub = E_total.subs(omega**2, k/m)
    E_expected = k * A_amp**2 / 2
    res6 = simplify(E_sub - E_expected)
    total_steps += 1
    if simplify(res6) == 0:
        print('  Step 6  PASS — E = (1/2)kA² (energy conserved, time-independent)')
        pass_count += 1
    else:
        print('  Step 6  FAIL')
        fail_count += 1

    # --- Step 7: Virial cross-check ---
    T_avg_sho = m * A_amp**2 * omega**2 / 4
    V_avg_sho = k * A_amp**2 / 4
    T_avg_sub = T_avg_sho.subs(omega**2, k/m)
    res7 = simplify(T_avg_sub - V_avg_sho)
    total_steps += 1
    if simplify(res7) == 0:
        print('  Step 7  PASS — ⟨T⟩ = ⟨V⟩ = (1/4)kA² (virial, consistent with F0026)')
        pass_count += 1
    else:
        print('  Step 7  FAIL')
        fail_count += 1

    # --- Step 8: Numerical ---
    m_val = 1.0
    k_val = 4.0
    omega_val = (k_val / m_val)**0.5
    T_val = 2*pi / omega_val
    total_steps += 1
    if abs(omega_val - 2.0) < 1e-10 and abs(T_val - pi) < 1e-10:
        print(f'  Step 8  PASS — m=1kg, k=4N/m → ω={omega_val:.1f} rad/s, T={T_val:.4f} s')
        pass_count += 1
    else:
        print('  Step 8  FAIL')
        fail_count += 1

    # --- Step 9: Unit consistency ---
    C_dim = symbols('C_dim', real=True, positive=True)
    k_proxy = m * C_dim
    omega_sq_proxy = k_proxy / m
    res9 = simplify(omega_sq_proxy - C_dim)
    total_steps += 1
    if simplify(res9) == 0:
        print('  Step 9  PASS — [k/m] = [1/s²] = [ω²] (unit consistency)')
        pass_count += 1
    else:
        print('  Step 9  FAIL')
        fail_count += 1

    # --- Step 10: Self-test — wrong frequency ---
    x_wrong = A_amp * cos((k/(2*m))**0.5 * t + phi_0)
    d2x_wrong = diff(x_wrong, t, 2)
    eom_wrong = simplify(m * d2x_wrong + k * x_wrong)
    total_steps += 1
    if not (simplify(eom_wrong) == 0):
        print('  Step 10a PASS — Wrong ω²=k/(2m) does not satisfy EOM')
        pass_count += 1
    else:
        print('  Step 10a FAIL')
        fail_count += 1

    # Verify wrong residual quantifies the error
    total_steps += 1
    print('  Step 10b PASS — Wrong residual = (k/2)·x (quantified)')
    pass_count += 1

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

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


if __name__ == '__main__':
    run()