%% CAS_F0002_VERIFY.m — Newtonian potential and 1/r^2 limit
%  Assertion-based CAS audit block
%  Pillar: Mechanics | Chain: Poisson eq -> spherical Laplacian -> Phi = -Gm/r -> F ~ 1/r^2
%  CalRef: Mechanics_Calibration S1A (potential closure)
%
%  Structure mirrors cas_F02.txt sections A-E exactly.
%  Every derivation step produces a verifiable symbolic assertion.
%  Final output: PASS or FAIL with step-level detail.
%
%  HARDENING NOTES (v2):
%    - r declared real + assume(r>0) for redeclaration safety
%    - C1, C2 declared real to prevent complex-valued limit edge cases
%    - isAlways(..., 'Unknown', false) throughout for ambiguity safety

clear; clc;
fprintf('=== CAS AUDIT: F0002 — Newtonian potential and 1/r^2 limit ===\n\n');

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

%% ---- A. INPUTS ----
% Phi: R^3\{0} -> R  (gravitational potential)
% rho(r) = m * delta^3(r)  (point mass at origin)
% Poisson: nabla^2 Phi = 4*pi*G*rho
% g(r) := -grad(Phi)

syms r real            % radial coordinate
assume(r > 0);         % r > 0 (survives accidental redeclaration)
syms G positive       % gravitational constant
syms m positive       % source mass
syms m_test positive  % test mass
syms C1 C2 real       % integration constants (real prevents complex limit edge cases)

% Phi as symbolic function of r
syms Phi(r)

fprintf('Section A: Inputs defined.\n');
fprintf('  Poisson: nabla^2 Phi = 4*pi*G*rho\n');
fprintf('  g := -grad(Phi)\n\n');

%% ---- B. ASSUMPTIONS / DOMAINS ----
% G > 0, m > 0   (enforced by 'positive')
% Domain: r > 0   (vacuum region, source at origin only)
% Phi in C^2(R^3\{0})
% BC: lim_{r->inf} Phi(r) = 0
% Spherical symmetry: Phi = Phi(r)

fprintf('Section B: Assumptions set (G>0, m>0, r>0, Phi->0 as r->inf).\n\n');

%% ---- C. ALLOWED LEMMAS ----
% C.1: Spherical Laplacian (radial): nabla^2 f(r) = (1/r^2) d/dr(r^2 df/dr)
% C.2: Vacuum: for r>0, rho=0 => nabla^2 Phi = 0
% C.3: Gauss: integral g.dA over S_R = -4*pi*G*m
% C.4: Radial field: integral g.dA = 4*pi*R^2 * g_r(R)

fprintf('Section C: Lemmas declared.\n');
fprintf('  C.1: Spherical Laplacian (radial)\n');
fprintf('  C.2: Vacuum region (nabla^2 Phi = 0 for r>0)\n');
fprintf('  C.3: Gauss divergence theorem\n');
fprintf('  C.4: Radial field surface integral\n\n');

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

% --- Step 1: Vacuum Laplace equation in spherical coords ---
% From C.1 + C.2: (1/r^2) d/dr(r^2 dPhi/dr) = 0 for r>0
% Equivalently: d/dr(r^2 dPhi/dr) = 0
%
% Verify: for any Phi(r), compute the spherical Laplacian
% and confirm its structure.
%
% We'll work with a general Phi(r) and check that the ODE
% d/dr(r^2 * Phi'(r)) = 0 integrates correctly.

dPhi = diff(Phi, r);
laplacian_radial = (1/r^2) * diff(r^2 * dPhi, r);

% The equation is laplacian_radial = 0
% Expand: (1/r^2) * d/dr(r^2 * Phi'(r)) = 0
% => d/dr(r^2 * Phi'(r)) = 0  (since 1/r^2 != 0 for r>0)

fprintf('  Step 1  INFO  Laplace eq: (1/r^2)*d/dr(r^2*dPhi/dr) = 0\n');

% --- Step 2: Verify general solution Phi = -C1/r + C2 ---
% Claim: Phi(r) = -C1/r + C2 satisfies (1/r^2)*d/dr(r^2*dPhi/dr) = 0
Phi_general = -C1/r + C2;
dPhi_general = diff(Phi_general, r);                  % C1/r^2
inner = r^2 * dPhi_general;                           % C1
d_inner = diff(inner, r);                             % 0
laplacian_check = simplify((1/r^2) * d_inner);

total_steps = total_steps + 1;
if isAlways(laplacian_check == 0, 'Unknown', 'false')
    fprintf('  Step 2  PASS  Phi = -C1/r + C2 satisfies nabla^2 Phi = 0\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 2  FAIL  Laplacian residual: %s\n', char(laplacian_check));
    fail_count = fail_count + 1;
end

% --- Step 3: First integration check ---
% d/dr(r^2 * Phi') = 0 => r^2 * Phi' = const
% For Phi = -C1/r + C2: Phi' = C1/r^2, so r^2 * Phi' = C1
r2_dPhi = simplify(r^2 * dPhi_general);

total_steps = total_steps + 1;
if isAlways(r2_dPhi == C1, 'Unknown', 'false')
    fprintf('  Step 3  PASS  r^2 * Phi''(r) = C1 (constant)\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 3  FAIL  r^2*Phi'' = %s (expected C1)\n', char(r2_dPhi));
    fail_count = fail_count + 1;
end

% --- Step 4: Boundary condition at infinity => C2 = 0 ---
% lim_{r->inf} Phi(r) = lim_{r->inf} (-C1/r + C2) = C2
% Require C2 = 0 => Phi(r) = -C1/r
Phi_bc = limit(-C1/r + C2, r, inf);

total_steps = total_steps + 1;
if isAlways(Phi_bc == C2, 'Unknown', 'false')
    fprintf('  Step 4  PASS  lim_{r->inf} Phi = C2 (set C2=0 by BC)\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 4  FAIL  limit = %s (expected C2)\n', char(Phi_bc));
    fail_count = fail_count + 1;
end

% After BC: Phi(r) = -C1/r
Phi_after_bc = -C1/r;

% --- Step 5: Gravitational field g_r = -dPhi/dr ---
% g_r = -d/dr(-C1/r) = -C1/r^2
g_r = -diff(Phi_after_bc, r);
g_r_expected = -C1/r^2;

step5_residual = simplify(g_r - g_r_expected);

total_steps = total_steps + 1;
if isAlways(step5_residual == 0, 'Unknown', 'false')
    fprintf('  Step 5  PASS  g_r = -dPhi/dr = -C1/r^2\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 5  FAIL  g_r residual: %s\n', char(step5_residual));
    fail_count = fail_count + 1;
end

% --- Step 6: Gauss law determines C1 = Gm ---
% Surface integral: 4*pi*R^2 * g_r(R) = -4*pi*G*m  (lemmas C.3 + C.4)
% Substitute g_r(R) = -C1/R^2:
%   4*pi*R^2 * (-C1/R^2) = -4*pi*C1 = -4*pi*G*m
%   => C1 = Gm
syms R positive
flux_computed = 4*sym(pi)*R^2 * subs(g_r_expected, r, R);
flux_gauss = -4*sym(pi)*G*m;

% From flux_computed = flux_gauss, solve for C1
C1_sol = solve(flux_computed - flux_gauss, C1);

total_steps = total_steps + 1;
if isAlways(C1_sol == G*m, 'Unknown', 'false')
    fprintf('  Step 6  PASS  Gauss law => C1 = G*m\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 6  FAIL  C1 = %s (expected G*m)\n', char(C1_sol));
    fail_count = fail_count + 1;
end

% --- Step 7: Final potential Phi(r) = -Gm/r ---
Phi_final = subs(Phi_after_bc, C1, G*m);
Phi_expected = -G*m/r;

step7_residual = simplify(Phi_final - Phi_expected);

total_steps = total_steps + 1;
if isAlways(step7_residual == 0, 'Unknown', 'false')
    fprintf('  Step 7  PASS  Phi(r) = -G*m/r\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 7  FAIL  Phi residual: %s\n', char(step7_residual));
    fail_count = fail_count + 1;
end

% --- Step 8: Verify Phi = -Gm/r satisfies Laplace equation ---
% Independent cross-check: compute spherical Laplacian of -Gm/r directly
dPhi_final = diff(Phi_expected, r);
laplacian_final = simplify((1/r^2) * diff(r^2 * dPhi_final, r));

total_steps = total_steps + 1;
if isAlways(laplacian_final == 0, 'Unknown', 'false')
    fprintf('  Step 8  PASS  nabla^2(-Gm/r) = 0 for r>0 (cross-check)\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 8  FAIL  Laplacian of -Gm/r = %s\n', char(laplacian_final));
    fail_count = fail_count + 1;
end

% --- Step 9: Recover 1/r^2 force law ---
% g = -grad(Phi) = -d/dr(-Gm/r) rhat = -(Gm/r^2) rhat
% F = m_test * g => |F| = G*m*m_test/r^2
g_final = -diff(Phi_expected, r);
g_expected = -G*m/r^2;

step9a_residual = simplify(g_final - g_expected);

F_mag = m_test * abs(g_final);
F_expected = G*m*m_test/r^2;

step9b_residual = simplify(F_mag - F_expected);

total_steps = total_steps + 1;
if isAlways(step9a_residual == 0, 'Unknown', 'false')
    fprintf('  Step 9a PASS  g_r = -Gm/r^2\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 9a FAIL  g_r residual: %s\n', char(step9a_residual));
    fail_count = fail_count + 1;
end

total_steps = total_steps + 1;
if isAlways(step9b_residual == 0, 'Unknown', 'false')
    fprintf('  Step 9b PASS  |F| = G*m*m_test/r^2 (inverse-square law)\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 9b FAIL  |F| residual: %s\n', char(step9b_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) ---
fprintf('  Unit check:\n');
fprintf('    Phi = -Gm/r: [m^3/(kg*s^2)]*[kg]/[m] = [m^2/s^2] (energy/mass)\n');
fprintf('    g = -Gm/r^2: [m^2/s^2]/[m] = [m/s^2] (acceleration)\n');
fprintf('    F = m_test*g: [kg]*[m/s^2] = [N] (force)\n');
fprintf('    PASS (all units consistent)\n\n');

% --- CAS flux consistency check ---
% Verify: integral of g over sphere = -4*pi*G*m
flux_final = 4*sym(pi)*R^2 * subs(g_expected, r, R);
flux_check = simplify(flux_final - flux_gauss);

total_steps = total_steps + 1;
if isAlways(flux_check == 0, 'Unknown', 'false')
    fprintf('  Flux check: 4*pi*R^2*g_r(R) = -4*pi*G*m  PASS\n');
    pass_count = pass_count + 1;
else
    fprintf('  Flux check: FAIL (residual: %s)\n', char(flux_check));
    fail_count = fail_count + 1;
end

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

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