%% CAS_F0026_VERIFY.m — Virial theorem (2⟨T⟩ = ⟨r·∇V⟩)
%  Pillar : Mechanics
%  CalRef : Math Appendix §4A, Mechanics Calibration §3B
%  Chain  : G = Σ p·r → dG/dt = 2T + Σ F·r → F = −∇V → time average → virial
%  Hardening: isAlways(..., 'Unknown', 'false') on every symbolic assertion
%            symvar-based presence checks; solver and has-function avoided

clear; clc;
fprintf('\n=== CAS F0026 : Virial theorem ===\n\n');

%% ── Symbols ──────────────────────────────────────────────────
syms t real
syms m real positive              % single-particle mass (1-body demonstration)

%  Position and velocity as symfuns of t
syms rx(t) ry(t) rz(t)
vx = diff(rx, t); vy = diff(ry, t); vz = diff(rz, t);
ax = diff(vx, t); ay = diff(vy, t); az = diff(vz, t);

%% ── A. Inputs ────────────────────────────────────────────────
%  Kinetic energy: T = (1/2) m v²
T_kin = m * (vx^2 + vy^2 + vz^2) / 2;

%  Virial: G = m v·r = m(vx·rx + vy·ry + vz·rz)
G_virial = m * (vx*rx + vy*ry + vz*rz);

%% ── Step 1: dG/dt product rule ──────────────────────────────
%  dG/dt = m(a·r + v·v) = m(a·r) + m v² = m(a·r) + 2T
dG_dt = diff(G_virial, t);

%  Expected: m(ax·rx + ay·ry + az·rz) + m(vx² + vy² + vz²)
dG_expected = m*(ax*rx + ay*ry + az*rz) + m*(vx^2 + vy^2 + vz^2);
res1 = simplify(dG_dt - dG_expected);
assert(isAlways(res1 == 0, 'Unknown', 'false'), ...
    'FAIL Step 1: dG/dt product rule');
fprintf('Step 1  PASS — dG/dt = m(a·r) + m v²\n');

%% ── Step 2: Identify 2T in dG/dt ───────────────────────────
%  m v² = 2T
%  So dG/dt = m(a·r) + 2T
twoT = 2 * T_kin;
res2 = simplify(m*(vx^2 + vy^2 + vz^2) - twoT);
assert(isAlways(res2 == 0, 'Unknown', 'false'), ...
    'FAIL Step 2: m v² ≠ 2T');
fprintf('Step 2  PASS — m v² = 2T\n');

%% ── Step 3: Newton substitution m·a = F = −∇V ──────────────
%  For single particle: m a_i = F_i = −∂V/∂r_i
%  So m(a·r) = F·r = −(∇V)·r = −r·∇V
%  dG/dt = 2T − r·∇V
%  Use symbolic force components
syms Fx(t) Fy(t) Fz(t)

%  m a·r = F·r (Newton substitution)
ma_dot_r = m*(ax*rx + ay*ry + az*rz);
F_dot_r = Fx*rx + Fy*ry + Fz*rz;

%  After substitution m a_i → F_i, these are equal
%  Verify the structural replacement: substitute ax → Fx/m etc.
ma_dot_r_sub = subs(ma_dot_r, [diff(vx,t), diff(vy,t), diff(vz,t)], ...
                               [Fx/m, Fy/m, Fz/m]);
res3 = simplify(ma_dot_r_sub - F_dot_r);
assert(isAlways(res3 == 0, 'Unknown', 'false'), ...
    'FAIL Step 3: Newton substitution m·a·r ≠ F·r');
fprintf('Step 3  PASS — m(a·r) → F·r (Newton substitution)\n');

%% ── Step 4: F = −∇V gives dG/dt = 2T − r·∇V ──────────────
%  With F_i = −∂V/∂r_i:
%  F·r = −(∇V)·r
%  dG/dt = 2T + F·r = 2T − r·∇V
syms gradV_x(t) gradV_y(t) gradV_z(t)
F_from_V_x = -gradV_x;
F_from_V_y = -gradV_y;
F_from_V_z = -gradV_z;

FdotR_fromV = F_from_V_x*rx + F_from_V_y*ry + F_from_V_z*rz;
rdotgradV = gradV_x*rx + gradV_y*ry + gradV_z*rz;

res4 = simplify(FdotR_fromV + rdotgradV);  % should be 0: F·r = −r·∇V
assert(isAlways(res4 == 0, 'Unknown', 'false'), ...
    'FAIL Step 4: F·r ≠ −r·∇V');
fprintf('Step 4  PASS — F·r = −r·∇V (potential substitution)\n');

%% ── Step 5: Virial theorem structure ────────────────────────
%  ⟨dG/dt⟩ = 0 for bound systems
%  0 = 2⟨T⟩ − ⟨r·∇V⟩
%  → 2⟨T⟩ = ⟨r·∇V⟩
%  Verify algebraic rearrangement: if dG/dt = 2T − r·∇V and ⟨dG/dt⟩ = 0,
%  then 2T = r·∇V (at the time-average level)
syms T_avg rdgV_avg real   % time-averaged quantities
%  ⟨dG/dt⟩ = 2⟨T⟩ − ⟨r·∇V⟩ = 0
%  → 2⟨T⟩ = ⟨r·∇V⟩
virial_lhs = 2*T_avg;
virial_rhs_val = rdgV_avg;
%  From ⟨dG/dt⟩ = 0: 2T_avg − rdgV_avg = 0 → T_avg = rdgV_avg/2
virial_residual = simplify(virial_lhs - virial_rhs_val);
%  This is 2*T_avg - rdgV_avg. Under the constraint this = 0, substitute:
virial_sub = subs(virial_residual, rdgV_avg, 2*T_avg);
assert(isAlways(virial_sub == 0, 'Unknown', 'false'), ...
    'FAIL Step 5: virial theorem algebraic structure');
fprintf('Step 5  PASS — 2⟨T⟩ = ⟨r·∇V⟩ (virial theorem)\n');

%% ── Step 6: Power-law test — V ∝ rⁿ ────────────────────────
%  For V = C rⁿ (homogeneous potential of degree n):
%  r·∇V = n·V (Euler's theorem)
%  Virial: 2⟨T⟩ = n⟨V⟩
%  Gravity: n = −1 → 2⟨T⟩ = −⟨V⟩ → ⟨E⟩ = ⟨T⟩ + ⟨V⟩ = −⟨T⟩
syms r_var C_coeff real positive
syms n_pow real
V_power = C_coeff * r_var^n_pow;
rdgV_power = r_var * diff(V_power, r_var);
euler_check = simplify(rdgV_power - n_pow * V_power);
assert(isAlways(euler_check == 0, 'Unknown', 'false'), ...
    'FAIL Step 6: Euler theorem r·dV/dr ≠ n·V');
fprintf('Step 6  PASS — Euler''s theorem: r·dV/dr = n·V for V=Crⁿ\n');

%% ── Step 7: Gravity test n = −1 ────────────────────────────
%  V = −GMm/r → r·∇V = −V → 2⟨T⟩ = −⟨V⟩
%  Total energy ⟨E⟩ = ⟨T⟩ + ⟨V⟩ = −⟨T⟩
syms T_grav V_grav real
%  2T = −V → V = −2T → E = T + V = T − 2T = −T
V_from_virial = -2*T_grav;
E_total = T_grav + V_from_virial;
res7 = simplify(E_total + T_grav);
assert(isAlways(res7 == 0, 'Unknown', 'false'), ...
    'FAIL Step 7: gravity E ≠ −T');
fprintf('Step 7  PASS — Gravity (n=−1): ⟨E⟩ = −⟨T⟩\n');

%% ── Step 8: Harmonic oscillator test n = 2 ─────────────────
%  V = (1/2)kx² → r·∇V = 2V → 2⟨T⟩ = 2⟨V⟩ → ⟨T⟩ = ⟨V⟩
%  Numerical: k=1, A=1, ω=1 → ⟨T⟩ = ⟨V⟩ = E/2
syms k_spr A_osc omega_osc real positive
%  x(t) = A cos(ωt), v = −Aω sin(ωt)
%  ⟨T⟩ = (1/2)m⟨v²⟩ = (1/4)mA²ω²
%  ⟨V⟩ = (1/2)k⟨x²⟩ = (1/4)kA²
%  With ω² = k/m: ⟨T⟩ = (1/4)kA² = ⟨V⟩
T_avg_sho = m * A_osc^2 * omega_osc^2 / 4;
V_avg_sho = k_spr * A_osc^2 / 4;
%  Substitute ω² = k/m
T_sub = subs(T_avg_sho, omega_osc^2, k_spr/m);
res8 = simplify(T_sub - V_avg_sho);
assert(isAlways(res8 == 0, 'Unknown', 'false'), ...
    'FAIL Step 8: SHO ⟨T⟩ ≠ ⟨V⟩');
fprintf('Step 8  PASS — SHO (n=2): ⟨T⟩ = ⟨V⟩ (energy equipartition)\n');

%% ── Step 9: Numerical — Earth orbit ────────────────────────
%  Circular orbit: v = √(GM/r)
%  T = (1/2)mv² = GMm/(2r), V = −GMm/r
%  2T = GMm/r = −V ✓ (virial for n=−1)
%  E = T + V = −GMm/(2r) = −T
G_grav = 6.674e-11; M_sun = 1.989e30; m_earth = 5.972e24;
r_orbit = 1.496e11;  % 1 AU
T_earth = 0.5 * m_earth * G_grav * M_sun / r_orbit;
V_earth = -G_grav * M_sun * m_earth / r_orbit;
virial_check = 2*T_earth + V_earth;  % should ≈ 0
assert(abs(virial_check) / abs(V_earth) < 1e-10, ...
    'FAIL Step 9: Earth orbit virial mismatch');
fprintf('Step 9  PASS — Earth orbit: 2T/|V| = %.6f (expected 1.0)\n', 2*T_earth/abs(V_earth));

%% ── Step 10: Self-test — wrong virial coefficient ───────────
%  Wrong: T = r·∇V (missing factor of 2)
%  For gravity: T = −V instead of 2T = −V → E = T+V = 0 (wrong)
syms T_w V_w real
V_wrong = -T_w;           % wrong: T = −V
E_wrong = T_w + V_wrong;  % = 0
%  Correct: V = −2T → E = −T ≠ 0
assert(~isAlways(E_wrong + T_w == 0, 'Unknown', 'false'), ...
    'FAIL Step 10a: wrong virial coefficient not detected');
fprintf('Step 10a PASS — Wrong coefficient (T=−V) gives E=0 ≠ −T\n');

%  Quantify: correct E = −T, wrong E = 0, difference = T
E_correct = -T_w;
res10b = simplify(E_wrong - E_correct - T_w);
assert(isAlways(res10b == 0, 'Unknown', 'false'), ...
    'FAIL Step 10b: wrong residual ≠ T');
fprintf('Step 10b PASS — Wrong residual = T (quantified)\n');

%% ── Summary ─────────────────────────────────────────────────
fprintf('\n✓ ALL 11 ASSERTIONS PASSED — F0026 audit complete.\n');
fprintf('  Chain: G=p·r → dG/dt=2T+F·r → F=−∇V → ⟨dG/dt⟩=0 → 2⟨T⟩=⟨r·∇V⟩\n');
fprintf('  Cross-ref: standalone (Mechanics pillar)\n');