%% CAS_F0001_VERIFY.m — Newton's law from variational closure
%  Assertion-based CAS audit block
%  Pillar: Mechanics | Chain: S1.1-S1.4 -> Euler-Lagrange -> m*a = F
%  CalRef: Mechanics_Calibration S1B (force closure)
%
%  Structure mirrors cas_F01.txt sections A-E exactly.
%  Every derivation step produces a verifiable symbolic assertion.
%  Final output: PASS or FAIL with step-level detail.
%
%  MATLAB FIX NOTES (v2):
%    - syms x(t) creates x as symfun; cannot then do syms V(x).
%      Fix: use independent variable X for V, then substitute V(X)->V(x(t)).
%    - Removed orphaned local_check function handle.

clear; clc;
fprintf('=== CAS AUDIT: F0001 — Newton''s law from variational closure ===\n\n');

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

%% ---- A. INPUTS ----
% Configuration: x : [t1,t2] -> R^3 (we work componentwise in 1D, then generalize)
% Lagrangian: L(x, xdot) = (1/2)*m*xdot^2 - V(x)
% Action: S[x] = int_{t1}^{t2} L dt
% Stationarity: delta S = 0

syms t real
syms m positive
syms x(t)

% --- MATLAB fix: use separate independent variable for V ---
% x is a symfun of t, so we define V over an independent sym X,
% then compose V(X) -> V(x(t)) by substitution.
syms X real
syms V(X)

% Compose: V evaluated along the trajectory x(t)
% Use explicit substitution (version-safe; avoids auto-composition ambiguity)
V_of_xt = subs(V(X), X, x);

L = sym(1)/2 * m * diff(x, t)^2 - V_of_xt;

fprintf('Section A: Inputs defined.\n');
fprintf('  L = (1/2)*m*xdot^2 - V(x(t))\n\n');

%% ---- B. ASSUMPTIONS / DOMAINS ----
% m > 0                          (enforced by 'positive')
% V: R -> R, differentiable      (symbolic function)
% x(t) in C^2                    (symbolic, smooth by default)
% Variations vanish at endpoints  (used implicitly in EL lemma)

fprintf('Section B: Assumptions set (m > 0, V differentiable, x in C^2).\n\n');

%% ---- C. ALLOWED LEMMAS ----
% C.1: Euler-Lagrange lemma
%      delta S = 0 for all delta x with delta x(t1) = delta x(t2) = 0
%      implies: d/dt(dL/dxdot) - dL/dx = 0
%
% C.2: Vector generalization (apply componentwise)
%      We verify in 1D; R^3 follows by identical application to each component.

fprintf('Section C: Lemmas declared.\n');
fprintf('  C.1: Euler-Lagrange (d/dt(dL/dxdot) - dL/dx = 0)\n');
fprintf('  C.2: Componentwise generalization to R^3\n\n');

%% ---- D. STEP LOG (explicit transforms with assertions) ----
fprintf('Section D: Step log\n');
fprintf('---------------------------------------------\n');

% Helper: dV/dx means d/dX V(X) evaluated at X = x(t)
dVdx = subs(diff(V, X), X, x);   % dV/dx evaluated on trajectory

% --- Step 1: Compute dL/d(xdot) ---
% Input: L = (1/2)*m*xdot^2 - V(x(t))
% Operation: partial derivative w.r.t. xdot
% Expected result: m*xdot
xdot_sym = diff(x, t);
pL_pxdot = diff(L, xdot_sym);

step1_target = m * xdot_sym;
step1_residual = simplify(pL_pxdot - step1_target);

total_steps = total_steps + 1;
if isAlways(step1_residual == 0, 'Unknown', 'false')
    fprintf('  Step 1  PASS  dL/d(xdot) = m*xdot\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 1  FAIL  dL/d(xdot) residual: %s\n', char(step1_residual));
    fail_count = fail_count + 1;
end

% --- Step 2: Time derivative of momentum ---
% Input: dL/d(xdot) = m*xdot (from step 1)
% Operation: d/dt
% Expected result: m*xddot
dt_pL_pxdot = diff(pL_pxdot, t);

step2_target = m * diff(x, t, 2);
step2_residual = simplify(dt_pL_pxdot - step2_target);

total_steps = total_steps + 1;
if isAlways(step2_residual == 0, 'Unknown', 'false')
    fprintf('  Step 2  PASS  d/dt(dL/d(xdot)) = m*xddot\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 2  FAIL  d/dt(dL/d(xdot)) residual: %s\n', char(step2_residual));
    fail_count = fail_count + 1;
end

% --- Step 3: Compute dL/dx ---
% Input: L = (1/2)*m*xdot^2 - V(x(t))
% Operation: partial derivative w.r.t. x
% Expected result: -dV/dx
%
% Note: MATLAB diff(L, x) differentiates w.r.t. x(t) directly.
% For V(x(t)), chain rule gives dV/dx * dx/dx = dV/dx (since dx/dx = 1
% for the partial derivative holding xdot fixed).
pL_px = diff(L, x);

step3_target = -dVdx;
step3_residual = simplify(pL_px - step3_target);

total_steps = total_steps + 1;
if isAlways(step3_residual == 0, 'Unknown', 'false')
    fprintf('  Step 3  PASS  dL/dx = -dV/dx\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 3  FAIL  dL/dx residual: %s\n', char(step3_residual));
    fail_count = fail_count + 1;
end

% --- Step 4: Apply Euler-Lagrange lemma (C.1) ---
% From C.1: d/dt(dL/dxdot) - dL/dx = 0
% Substituting steps 2 and 3:
%   m*xddot - (-dV/dx) = 0
%   m*xddot + dV/dx = 0
EL_expression = dt_pL_pxdot - pL_px;

step4_target = m * diff(x, t, 2) + dVdx;
step4_residual = simplify(EL_expression - step4_target);

total_steps = total_steps + 1;
if isAlways(step4_residual == 0, 'Unknown', 'false')
    fprintf('  Step 4  PASS  EL equation: m*xddot + dV/dx = 0\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 4  FAIL  EL residual: %s\n', char(step4_residual));
    fail_count = fail_count + 1;
end

% --- Step 5: Solve for equation of motion ---
% From step 4: m*xddot + dV/dx = 0  =>  m*xddot = -dV/dx
% This is Newton's second law in potential form.
Newton_LHS = m * diff(x, t, 2);
Newton_RHS = -dVdx;

% EL_expression should equal Newton_LHS - Newton_RHS = m*xddot + dV/dx
step5_residual = simplify(EL_expression - (Newton_LHS - Newton_RHS));

total_steps = total_steps + 1;
if isAlways(step5_residual == 0, 'Unknown', 'false')
    fprintf('  Step 5  PASS  m*xddot = -dV/dx (Newton''s law)\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 5  FAIL  Newton residual: %s\n', char(step5_residual));
    fail_count = fail_count + 1;
end

% --- Step 6: Define force and verify substitution ---
% F := -grad(V)  =>  m*xddot = F
% Check: substituting F = -dV/dx into m*a = F recovers step 5
F_def = -dVdx;
ma_eq = Newton_LHS;

% m*a - F = m*xddot + dV/dx = EL_expression
step6_residual = simplify(ma_eq - F_def - EL_expression);

total_steps = total_steps + 1;
if isAlways(step6_residual == 0, 'Unknown', 'false')
    fprintf('  Step 6  PASS  F := -dV/dx; m*a = F consistent with EL\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 6  FAIL  Force substitution residual: %s\n', char(step6_residual));
    fail_count = fail_count + 1;
end

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

%% ---- E. CHECK OUTPUTS ----
fprintf('Section E: Output checks\n');
fprintf('---------------------------------------------\n');

% --- Unit check (dimensional analysis) ---
% LHS: m * xddot  ->  [kg] * [m/s^2] = [kg*m/s^2] = [N]
% RHS: -dV/dx     ->  [J]/[m] = [kg*m^2/s^2]/[m] = [kg*m/s^2] = [N]
fprintf('  Unit check:\n');
fprintf('    LHS: m*xddot  -> [kg]*[m/s^2] = [N]\n');
fprintf('    RHS: -dV/dx   -> [J]/[m] = [N]\n');
fprintf('    PASS (both sides [N])\n\n');

% --- CAS simplification check ---
% The full EL expression minus (m*xddot + dV/dx) should simplify to 0
final_expr = simplify(EL_expression - step4_target);
total_steps = total_steps + 1;
if isAlways(final_expr == 0, 'Unknown', 'false')
    fprintf('  CAS simplification: EL_expr - (m*xddot + dV/dx) = 0  PASS\n');
    pass_count = pass_count + 1;
else
    fprintf('  CAS simplification: FAIL (residual: %s)\n', char(final_expr));
    fail_count = fail_count + 1;
end

% --- Consistency check: EL via functionalDerivative (independent method) ---
% MATLAB's functionalDerivative computes delta S / delta x directly
% This is a cross-check against our manual step-by-step derivation
try
    EL_auto = functionalDerivative(L, x);
    cross_check = simplify(EL_auto - EL_expression);
    total_steps = total_steps + 1;
    if isAlways(cross_check == 0, 'Unknown', 'false')
        fprintf('  Cross-check (functionalDerivative): PASS\n');
        pass_count = pass_count + 1;
    else
        fprintf('  Cross-check (functionalDerivative): FAIL (residual: %s)\n', char(cross_check));
        fail_count = fail_count + 1;
    end
catch ME
    total_steps = total_steps + 1;
    fprintf('  Cross-check (functionalDerivative): SKIPPED (%s)\n', ME.message);
end

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

%% ---- VERDICT ----
fprintf('=============================================\n');
fprintf('  F0001 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 F0001.\n');