%% CAS_F0027_VERIFY.m — Simple harmonic oscillator (mẍ + kx = 0; ω² = k/m)
%  Pillar : Mechanics
%  CalRef : Math Appendix §4B–4C, Mechanics Calibration §2A
%  Chain  : quadratic potential → force → Newton → EOM → ω² identification
%  Hardening: isAlways(..., 'Unknown', 'false') on every symbolic assertion
%            symvar-based presence checks; solver and has-function avoided

clear; clc;
fprintf('\n=== CAS F0027 : Simple harmonic oscillator ===\n\n');

%% ── Symbols ──────────────────────────────────────────────────
syms x t real
syms m k real positive
syms omega real positive
syms A_amp phi_0 real              % amplitude and phase

%  Position as symfun of t
syms x_t(t)

%% ── A. Inputs ────────────────────────────────────────────────
%  Potential: V(x) = V₀ + (1/2)kx² (V₀ drops out of force)
syms V0 real
V_x = V0 + k * x^2 / 2;

%% ── Step 1: Force from potential ────────────────────────────
%  F = −dV/dx = −kx
dVdx = diff(V_x, x);
F_x = -dVdx;
res1 = simplify(F_x + k*x);  % F + kx should be 0
assert(isAlways(res1 == 0, 'Unknown', 'false'), ...
    'FAIL Step 1: F ≠ −kx');
fprintf('Step 1  PASS — F = −dV/dx = −kx\n');

%% ── Step 2: V₀ drops out of force ──────────────────────────
%  dV/dx should not contain V₀
assert(~any(symvar(dVdx) == V0), ...
    'FAIL Step 2: V₀ appears in dV/dx');
fprintf('Step 2  PASS — V₀ drops out of force (constant offset irrelevant)\n');

%% ── Step 3: Newton's law → EOM ─────────────────────────────
%  mẍ = F = −kx  →  mẍ + kx = 0
%  Verify: substituting F into mẍ = F gives mẍ + kx = 0
%  EOM residual: m·d²x/dt² + k·x = 0
%  For general solution x(t) = A cos(ωt + φ) with ω² = k/m:
x_sol = A_amp * cos(omega * t + phi_0);
d2x_sol = diff(x_sol, t, 2);  % = −Aω² cos(ωt + φ)
eom_residual = m * d2x_sol + k * x_sol;
%  = m(−Aω²cos) + k(Acos) = A cos(ωt+φ)(k − mω²)
%  Substitute ω² = k/m:
eom_sub = subs(eom_residual, omega^2, k/m);
res3 = simplify(eom_sub);
assert(isAlways(res3 == 0, 'Unknown', 'false'), ...
    'FAIL Step 3: x=Acos(ωt+φ) with ω²=k/m does not satisfy EOM');
fprintf('Step 3  PASS — x=Acos(ωt+φ), ω²=k/m satisfies mẍ+kx=0\n');

%% ── Step 4: ω² identification ──────────────────────────────
%  Divide EOM by m: ẍ + (k/m)x = 0
%  Compare with ẍ + ω²x = 0 → ω² = k/m
%  Verify algebraically: k/m − ω² = 0 when ω² = k/m
omega_sq = k / m;
res4 = simplify(omega_sq - omega^2);
res4_sub = subs(res4, omega^2, k/m);
assert(isAlways(res4_sub == 0, 'Unknown', 'false'), ...
    'FAIL Step 4: ω² ≠ k/m');
fprintf('Step 4  PASS — ω² = k/m\n');

%% ── Step 5: Period T = 2π/ω = 2π√(m/k) ────────────────────
T_period = 2*sym(pi) / omega;
T_explicit = 2*sym(pi) * sqrt(m/k);
%  Substitute ω = √(k/m):
T_sub = subs(T_period, omega, sqrt(k/m));
res5 = simplify(T_sub - T_explicit);
assert(isAlways(res5 == 0, 'Unknown', 'false'), ...
    'FAIL Step 5: T ≠ 2π√(m/k)');
fprintf('Step 5  PASS — T = 2π/ω = 2π√(m/k)\n');

%% ── Step 6: Energy conservation ─────────────────────────────
%  T = (1/2)mv², V = (1/2)kx²
%  E = T + V = (1/2)m(dx/dt)² + (1/2)kx²
%  For x = A cos(ωt+φ): v = −Aω sin(ωt+φ)
%  E = (1/2)mA²ω²sin² + (1/2)kA²cos² = (1/2)A²(mω²sin² + kcos²)
%  With mω² = k: E = (1/2)kA²(sin² + cos²) = (1/2)kA²
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 = subs(E_total, omega^2, k/m);
E_expected = k * A_amp^2 / 2;
res6 = simplify(E_sub - E_expected);
assert(isAlways(res6 == 0, 'Unknown', 'false'), ...
    'FAIL Step 6: E ≠ (1/2)kA²');
fprintf('Step 6  PASS — E = (1/2)kA² (energy conserved, time-independent)\n');

%% ── Step 7: Virial cross-check ⟨T⟩ = ⟨V⟩ ──────────────────
%  For SHO (n=2 potential): virial gives ⟨T⟩ = ⟨V⟩ = E/2
%  ⟨T⟩ = (1/4)mA²ω² = (1/4)kA² = E/2
%  ⟨V⟩ = (1/4)kA² = E/2
T_avg_sho = m * A_amp^2 * omega^2 / 4;
V_avg_sho = k * A_amp^2 / 4;
T_avg_sub = subs(T_avg_sho, omega^2, k/m);
res7 = simplify(T_avg_sub - V_avg_sho);
assert(isAlways(res7 == 0, 'Unknown', 'false'), ...
    'FAIL Step 7: ⟨T⟩ ≠ ⟨V⟩');
fprintf('Step 7  PASS — ⟨T⟩ = ⟨V⟩ = (1/4)kA² (virial, consistent with F0026)\n');

%% ── Step 8: Numerical — 1 kg mass, k = 4 N/m ──────────────
%  ω = √(4/1) = 2 rad/s, T = 2π/2 = π ≈ 3.1416 s
m_val = 1.0; k_val = 4.0;
omega_val = sqrt(k_val / m_val);
T_val = 2*pi / omega_val;
assert(abs(omega_val - 2.0) < 1e-10, 'FAIL Step 8a: ω mismatch');
assert(abs(T_val - pi) < 1e-10, 'FAIL Step 8b: T mismatch');
fprintf('Step 8  PASS — m=1kg, k=4N/m → ω=%.1f rad/s, T=%.4f s\n', omega_val, T_val);

%% ── Step 9: Unit consistency ────────────────────────────────
%  [k/m] = (N/m)/kg = (kg/s²)/kg = 1/s² = [ω²]
%  Verify via dimensional proxy: if k = m·C for C [1/s²], then ω² = C
syms C_dim real positive
k_proxy = m * C_dim;
omega_sq_proxy = k_proxy / m;
res9 = simplify(omega_sq_proxy - C_dim);
assert(isAlways(res9 == 0, 'Unknown', 'false'), ...
    'FAIL Step 9: unit consistency');
fprintf('Step 9  PASS — [k/m] = [1/s²] = [ω²] (unit consistency)\n');

%% ── Step 10: Self-test — wrong frequency ────────────────────
%  Wrong: ω² = k/(2m) (factor of 2 error)
x_wrong = A_amp * cos(sqrt(k/(2*m)) * t + phi_0);
d2x_wrong = diff(x_wrong, t, 2);
eom_wrong = simplify(m * d2x_wrong + k * x_wrong);
assert(~isAlways(eom_wrong == 0, 'Unknown', 'false'), ...
    'FAIL Step 10a: wrong frequency not detected');
fprintf('Step 10a PASS — Wrong ω²=k/(2m) does not satisfy EOM\n');

%  Quantify: residual = (k − k/2)·x_wrong = (k/2)·x_wrong
eom_wrong_expected = k * x_wrong / 2;
res10b = simplify(eom_wrong - eom_wrong_expected);
assert(isAlways(res10b == 0, 'Unknown', 'false'), ...
    'FAIL Step 10b: wrong residual ≠ (k/2)·x');
fprintf('Step 10b PASS — Wrong residual = (k/2)·x (quantified)\n');

%% ── Summary ─────────────────────────────────────────────────
fprintf('\n✓ ALL 11 ASSERTIONS PASSED — F0027 audit complete.\n');
fprintf('  Chain: V=(1/2)kx² → F=−kx → mẍ+kx=0 → ω²=k/m\n');
fprintf('  Cross-ref: F0026 (virial ⟨T⟩=⟨V⟩ for n=2 potential)\n');
