%% CAS_F0003_VERIFY.m — 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.
%
%  NOTE: F0003 source does not follow the full A-E format (no explicit
%  lemma section). The audit verifies the algebraic derivation chain as given.

clear; clc;
fprintf('=== CAS AUDIT: F0003 — Relativistic dispersion + mass-energy ===\n\n');

pass_count = 0;
fail_count = 0;
total_steps = 0;

%% ---- INPUTS / ASSUMPTIONS ----
% Four-momentum: p^mu = (E/c, p_vec)
% Minkowski norm: p^mu p_mu = (E/c)^2 - |p_vec|^2
% Invariant rest mass: p^mu p_mu = m^2 c^2
% Lorentz invariance requires this be frame-independent.

syms E real            % total energy
syms p real            % magnitude of 3-momentum |p_vec| (scalar, >= 0)
syms m positive        % rest mass, m > 0
syms c positive        % speed of light, c > 0

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

%% ---- D. STEP LOG ----
fprintf('Step log\n');
fprintf('---------------------------------------------\n');

% --- Step 1: Invariant norm of four-momentum ---
% p^mu p_mu = (E/c)^2 - p^2 = m^2 c^2
% This is the DEFINITION (input). We encode it as an identity.
norm_LHS = (E/c)^2 - p^2;
norm_RHS = m^2 * c^2;

% The invariant mass relation:
% (E/c)^2 - p^2 - m^2*c^2 = 0
invariant_relation = norm_LHS - norm_RHS;

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

% --- Step 2: Derive relativistic dispersion relation ---
% From (E/c)^2 - p^2 = m^2*c^2
% Multiply both sides by c^2:
%   E^2 - p^2*c^2 = m^2*c^4
% Rearrange:
%   E^2 = p^2*c^2 + m^2*c^4
%
% Verify algebraically (no solve() — avoids version-sensitive return types):
% Multiply invariant relation by c^2:
%   c^2 * ((E/c)^2 - p^2 - m^2*c^2) = E^2 - p^2*c^2 - m^2*c^4
% So E^2 = p^2*c^2 + m^2*c^4  iff  c^2 * invariant_relation = E^2 - dispersion

dispersion_expected = p^2*c^2 + m^2*c^4;

% Step 2: Direct residual check
% (E^2 - dispersion_expected) should equal c^2 * invariant_relation identically
step2_residual = simplify(E^2 - dispersion_expected - c^2 * invariant_relation);

total_steps = total_steps + 1;
if isAlways(step2_residual == 0, 'Unknown', 'false')
    fprintf('  Step 2  PASS  E^2 = p^2*c^2 + m^2*c^4 (dispersion relation)\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 2  FAIL  Dispersion residual: %s\n', char(step2_residual));
    fail_count = fail_count + 1;
end

% Step 2b: Cross-check by substituting dispersion into norm
% If E^2 = p^2*c^2 + m^2*c^4, then (E/c)^2 - p^2 = m^2*c^2
% Substitute E^2 = dispersion_expected into norm_LHS:
norm_substituted = subs(norm_LHS, E^2, dispersion_expected);
step2b_residual = simplify(norm_substituted - norm_RHS);

total_steps = total_steps + 1;
if isAlways(step2b_residual == 0, 'Unknown', 'false')
    fprintf('  Step 2b PASS  Substituting dispersion back into norm gives m^2*c^2\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 2b FAIL  Norm back-substitution residual: %s\n', char(step2b_residual));
    fail_count = fail_count + 1;
end

% --- Step 3: Rest frame limit (p = 0) ---
% Set p = 0 in E^2 = p^2*c^2 + m^2*c^4
% => E^2 = m^2*c^4
% => E = m*c^2 (taking positive root for energy)
E_squared_rest = subs(dispersion_expected, p, 0);
E_squared_rest_expected = m^2 * c^4;

total_steps = total_steps + 1;
if isAlways(simplify(E_squared_rest - E_squared_rest_expected) == 0, 'Unknown', 'false')
    fprintf('  Step 3a PASS  E^2(p=0) = m^2*c^4\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 3a FAIL  E^2(p=0) = %s\n', char(E_squared_rest));
    fail_count = fail_count + 1;
end

% 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 = total_steps + 1;
if isAlways(simplify(E_rest_simplified - E_rest_expected) == 0, 'Unknown', 'false')
    fprintf('  Step 3b PASS  E(p=0) = m*c^2 (mass-energy equivalence)\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 3b FAIL  E(p=0) = %s (expected m*c^2)\n', char(E_rest_simplified));
    fail_count = fail_count + 1;
end

% --- Step 4: Low-velocity limit (Taylor expansion) ---
% E = sqrt(p^2*c^2 + m^2*c^4)
% For p << mc, expand to second order in p:
%   E ≈ mc^2 + p^2/(2m) + O(p^4)
% Kinetic energy: E_kin = E - mc^2 ≈ p^2/(2m)
%
% Taylor expand sqrt(p^2*c^2 + m^2*c^4) around p = 0 to order p^2
E_full = sqrt(p^2*c^2 + m^2*c^4);
E_taylor = taylor(E_full, p, 0, 'Order', 4);  % up to O(p^4)

% Expected: m*c^2 + p^2/(2*m) [plus O(p^4) terms]
% Extract the p^0 and p^2 coefficients
E_taylor_expected = m*c^2 + p^2/(2*m);

% The Taylor expansion should match to order p^2
% Subtract and check that remaining terms are O(p^4) or higher
residual_taylor = simplify(E_taylor - E_taylor_expected);

total_steps = total_steps + 1;
% The residual should contain only p^4 and higher terms, or be exactly 0
% For a clean check: evaluate residual at p=0, should be 0
residual_at_zero = subs(residual_taylor, p, 0);
if isAlways(residual_at_zero == 0, 'Unknown', 'false')
    fprintf('  Step 4a PASS  Taylor: E ≈ mc^2 + p^2/(2m) + O(p^4)\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 4a FAIL  Taylor residual at p=0: %s\n', char(residual_at_zero));
    fail_count = fail_count + 1;
end

% Verify kinetic energy: E_kin = E - mc^2 ≈ p^2/(2m)
E_kin_taylor = simplify(E_taylor - m*c^2);
% Check coefficient of p^2
E_kin_coeff_p2 = subs(diff(E_kin_taylor, p, 2), p, 0) / 2;  % Taylor coeff

total_steps = total_steps + 1;
if isAlways(simplify(E_kin_coeff_p2 - 1/(2*m)) == 0, 'Unknown', 'false')
    fprintf('  Step 4b PASS  E_kin leading term = p^2/(2m) (classical kinetic energy)\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 4b FAIL  E_kin p^2 coefficient: %s (expected 1/(2m))\n', char(E_kin_coeff_p2));
    fail_count = fail_count + 1;
end

fprintf('---------------------------------------------\n\n');

%% ---- CHECK OUTPUTS ----
fprintf('Output checks\n');
fprintf('---------------------------------------------\n');

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

% --- Consistency: dispersion relation is equivalent to norm ---
% E^2 - p^2*c^2 - m^2*c^4 should equal c^2 * ((E/c)^2 - p^2 - m^2*c^2)
consistency = simplify((E^2 - p^2*c^2 - m^2*c^4) - c^2*invariant_relation);
total_steps = total_steps + 1;
if isAlways(consistency == 0, 'Unknown', 'false')
    fprintf('  Consistency: dispersion ↔ norm equivalence  PASS\n');
    pass_count = pass_count + 1;
else
    fprintf('  Consistency: FAIL (residual: %s)\n', char(consistency));
    fail_count = fail_count + 1;
end

% --- Cross-check: E(p=0) from full dispersion ---
E_at_rest_cross = simplify(subs(E_full, p, 0));
total_steps = total_steps + 1;
if isAlways(E_at_rest_cross == m*c^2, 'Unknown', 'false')
    fprintf('  Cross-check: sqrt(0 + m^2*c^4) = m*c^2  PASS\n');
    pass_count = pass_count + 1;
else
    fprintf('  Cross-check: FAIL (%s)\n', char(E_at_rest_cross));
    fail_count = fail_count + 1;
end

fprintf('---------------------------------------------\n\n');

%% ---- VERDICT ----
fprintf('=============================================\n');
fprintf('  F0003 AUDIT RESULT\n');
fprintf('  Steps: %d  |  Pass: %d  |  Fail: %d\n', total_steps, pass_count, fail_count);
if fail_count == 0
    fprintf('  STATUS: *** PASS ***\n');
else
    fprintf('  STATUS: *** FAIL *** (%d step(s) failed)\n', fail_count);
end
fprintf('=============================================\n');
fprintf('Audit complete for F0003.\n');