%% CAS_F0015_VERIFY.m -- Momentum flux conservation (Poynting theorem)
%  Assertion-based CAS audit block
%  Pillar: Electromagnetism | Chain: Maxwell -> energy density -> Poynting theorem
%  CalRef: Electromagnetism Math Appendix S3.4, Calibration S2B
%
%  Structure mirrors cas_F14.txt (= F0015) sections A-E.
%  Strict derivation of du/dt + div(S) = 0.
%  Momentum part (dg/dt + div(T) = 0) referenced as lemma per source.
%
%  Verifies:
%    1. Time derivative of E^2: d/dt(E^2) = 2*E.dE/dt (componentwise)
%    2. Time derivative of B^2: d/dt(B^2) = 2*B.dB/dt (componentwise)
%    3. du/dt expansion via chain rule
%    4. Maxwell substitution: dE/dt and dB/dt replaced
%    5. Vector identity: div(ExB) = B.(curl E) - E.(curl B)
%    6. Poynting theorem assembly: du/dt + div(S) = 0
%    7. Poynting vector S = (1/mu0)*ExB definition consistency
%    8. Numerical: c = 1/sqrt(mu0*eps0) check
%    9. Momentum density: g = S/c^2 definition
%   10. Self-test: wrong sign in Faraday gives nonzero residual
%   11. Self-test: wrong residual quantified
%
%  HARDENING: isAlways(..., 'Unknown', 'false') throughout.

clear; clc;
fprintf('=== CAS AUDIT: F0015 -- Poynting theorem ===\n\n');

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

%% ---- A. INPUTS ----
syms x y z t real
syms eps0 mu0 positive     % vacuum permittivity, permeability

% Six EM field components as symfuns of (x,y,z,t)
syms Ex(x,y,z,t) Ey(x,y,z,t) Ez(x,y,z,t)
syms Bx(x,y,z,t) By(x,y,z,t) Bz(x,y,z,t)

fprintf('Section A: Inputs defined.\n');
fprintf('  u = (1/2)(eps0*E^2 + B^2/mu0)\n');
fprintf('  S = (1/mu0)*ExB\n');
fprintf('  Maxwell vacuum equations\n\n');

%% ---- B. ASSUMPTIONS / DOMAINS ----
fprintf('Section B: Smooth vacuum fields, no sources.\n\n');

%% ---- C. ALLOWED LEMMAS ----
fprintf('Section C: Lemmas declared.\n');
fprintf('  C.1: div(ExB) = B.(curl E) - E.(curl B)\n');
fprintf('  C.2: d/dt(E^2) = 2*E.dE/dt, d/dt(B^2) = 2*B.dB/dt\n\n');

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

% --- Step 1: Time derivative of E^2 ---
% E^2 = Ex^2 + Ey^2 + Ez^2
% d/dt(E^2) = 2*(Ex*dEx/dt + Ey*dEy/dt + Ez*dEz/dt)
E_sq = Ex^2 + Ey^2 + Ez^2;
dE_sq_dt = diff(E_sq, t);
E_dot_dEdt = 2*(Ex*diff(Ex,t) + Ey*diff(Ey,t) + Ez*diff(Ez,t));

step1_residual = simplify(dE_sq_dt - E_dot_dEdt);

total_steps = total_steps + 1;
if isAlways(step1_residual == 0, 'Unknown', 'false')
    fprintf('  Step 1  PASS  d/dt(E^2) = 2*E.dE/dt\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 1  FAIL  residual: %s\n', char(step1_residual));
    fail_count = fail_count + 1;
end

% --- Step 2: Time derivative of B^2 ---
B_sq = Bx^2 + By^2 + Bz^2;
dB_sq_dt = diff(B_sq, t);
B_dot_dBdt = 2*(Bx*diff(Bx,t) + By*diff(By,t) + Bz*diff(Bz,t));

step2_residual = simplify(dB_sq_dt - B_dot_dBdt);

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

% --- Step 3: du/dt expansion ---
% u = (1/2)*(eps0*E^2 + B^2/mu0)
% du/dt = (1/2)*(eps0*d(E^2)/dt + (1/mu0)*d(B^2)/dt)
%       = eps0*E.dE/dt + (1/mu0)*B.dB/dt
u = sym(1)/sym(2) * (eps0*E_sq + B_sq/mu0);
du_dt = diff(u, t);

% Expected form using dot products
du_dt_expected = eps0*(Ex*diff(Ex,t) + Ey*diff(Ey,t) + Ez*diff(Ez,t)) ...
               + (1/mu0)*(Bx*diff(Bx,t) + By*diff(By,t) + Bz*diff(Bz,t));

step3_residual = simplify(du_dt - du_dt_expected);

total_steps = total_steps + 1;
if isAlways(step3_residual == 0, 'Unknown', 'false')
    fprintf('  Step 3  PASS  du/dt = eps0*E.dE/dt + (1/mu0)*B.dB/dt\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 3  FAIL  residual: %s\n', char(step3_residual));
    fail_count = fail_count + 1;
end

% --- Step 4: Maxwell substitution ---
% Ampere (vacuum): curl B = mu0*eps0*dE/dt => dE/dt = (1/(mu0*eps0))*curl B
% Faraday: curl E = -dB/dt => dB/dt = -curl E
%
% Curl E components:
curl_E_x = diff(Ez, y) - diff(Ey, z);
curl_E_y = diff(Ex, z) - diff(Ez, x);
curl_E_z = diff(Ey, x) - diff(Ex, y);

% Curl B components:
curl_B_x = diff(Bz, y) - diff(By, z);
curl_B_y = diff(Bx, z) - diff(Bz, x);
curl_B_z = diff(By, x) - diff(Bx, y);

% Substitution map for dE/dt and dB/dt:
% dEx/dt = (1/(mu0*eps0))*curl_B_x, etc.
% dBx/dt = -curl_E_x, etc.
%
% du/dt after substitution:
% = eps0*(Ex*(curl_B_x/(mu0*eps0)) + Ey*(curl_B_y/(mu0*eps0)) + Ez*(curl_B_z/(mu0*eps0)))
%   + (1/mu0)*(Bx*(-curl_E_x) + By*(-curl_E_y) + Bz*(-curl_E_z))
% = (1/mu0)*(Ex*curl_B_x + Ey*curl_B_y + Ez*curl_B_z)
%   - (1/mu0)*(Bx*curl_E_x + By*curl_E_y + Bz*curl_E_z)
% = (1/mu0)*[E.(curl B) - B.(curl E)]

du_dt_maxwell = (1/mu0)*(Ex*curl_B_x + Ey*curl_B_y + Ez*curl_B_z) ...
              - (1/mu0)*(Bx*curl_E_x + By*curl_E_y + Bz*curl_E_z);

% To verify, we substitute Maxwell into du_dt_expected:
du_dt_sub = eps0*(Ex*curl_B_x/(mu0*eps0) + Ey*curl_B_y/(mu0*eps0) + Ez*curl_B_z/(mu0*eps0)) ...
          + (1/mu0)*(Bx*(-curl_E_x) + By*(-curl_E_y) + Bz*(-curl_E_z));

step4_residual = simplify(du_dt_sub - du_dt_maxwell);

total_steps = total_steps + 1;
if isAlways(step4_residual == 0, 'Unknown', 'false')
    fprintf('  Step 4  PASS  Maxwell substitution: du/dt = (1/mu0)[E.(curl B) - B.(curl E)]\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 4  FAIL  residual: %s\n', char(step4_residual));
    fail_count = fail_count + 1;
end

% --- Step 5: Vector identity div(ExB) = B.(curl E) - E.(curl B) ---
% (ExB)_x = Ey*Bz - Ez*By
% (ExB)_y = Ez*Bx - Ex*Bz
% (ExB)_z = Ex*By - Ey*Bx
cross_x = Ey*Bz - Ez*By;
cross_y = Ez*Bx - Ex*Bz;
cross_z = Ex*By - Ey*Bx;

div_ExB = diff(cross_x, x) + diff(cross_y, y) + diff(cross_z, z);

B_dot_curlE = Bx*curl_E_x + By*curl_E_y + Bz*curl_E_z;
E_dot_curlB = Ex*curl_B_x + Ey*curl_B_y + Ez*curl_B_z;

step5_residual = simplify(div_ExB - (B_dot_curlE - E_dot_curlB));

total_steps = total_steps + 1;
if isAlways(step5_residual == 0, 'Unknown', 'false')
    fprintf('  Step 5  PASS  div(ExB) = B.(curl E) - E.(curl B)\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 5  FAIL  residual: %s\n', char(step5_residual));
    fail_count = fail_count + 1;
end

% --- Step 6: Poynting theorem assembly ---
% du/dt = (1/mu0)[E.(curl B) - B.(curl E)]
%       = (1/mu0)*[-(B.(curl E) - E.(curl B))]
%       = -(1/mu0)*div(ExB)
%       = -div(S)
%
% So: du/dt + div(S) = 0
%
% Verify: du_dt_maxwell + (1/mu0)*div_ExB = 0
% du_dt_maxwell = (1/mu0)*[E_dot_curlB - B_dot_curlE]
% div_ExB = B_dot_curlE - E_dot_curlB
% Sum = (1/mu0)*[E_dot_curlB - B_dot_curlE + B_dot_curlE - E_dot_curlB] = 0

poynting_sum = du_dt_maxwell + (1/mu0)*div_ExB;
step6_residual = simplify(poynting_sum);

total_steps = total_steps + 1;
if isAlways(step6_residual == 0, 'Unknown', 'false')
    fprintf('  Step 6  PASS  du/dt + div(S) = 0 (Poynting theorem)\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 6  FAIL  residual: %s\n', char(step6_residual));
    fail_count = fail_count + 1;
end

% --- Step 7: Poynting vector definition consistency ---
% S = (1/mu0)*ExB
% S_x = (1/mu0)*cross_x, etc.
% div(S) = (1/mu0)*div(ExB)
% Already used in Step 6. Verify the factor explicitly:
Sx = cross_x / mu0;
Sy = cross_y / mu0;
Sz = cross_z / mu0;
div_S = diff(Sx, x) + diff(Sy, y) + diff(Sz, z);
div_S_expected = div_ExB / mu0;

step7_residual = simplify(div_S - div_S_expected);

total_steps = total_steps + 1;
if isAlways(step7_residual == 0, 'Unknown', 'false')
    fprintf('  Step 7  PASS  div(S) = (1/mu0)*div(ExB)\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 7  FAIL  residual: %s\n', char(step7_residual));
    fail_count = fail_count + 1;
end

% --- Step 8: Numerical c = 1/sqrt(mu0*eps0) ---
mu0_val = 4*pi*1e-7;
eps0_val = 8.854187817e-12;
c_computed = 1/sqrt(mu0_val * eps0_val);
c_expected = 299792458;
rel_error = abs(c_computed - c_expected) / c_expected;

total_steps = total_steps + 1;
if rel_error < 1e-6
    fprintf('  Step 8  PASS  c = 1/sqrt(mu0*eps0) = %.6e m/s (rel err %.2e)\n', c_computed, rel_error);
    pass_count = pass_count + 1;
else
    fprintf('  Step 8  FAIL  c rel error: %.2e\n', rel_error);
    fail_count = fail_count + 1;
end

% --- Step 9: Momentum density g = S/c^2 ---
% g_x = S_x/c^2 = (1/mu0)*cross_x / c^2
% With c^2 = 1/(mu0*eps0): g_x = eps0*cross_x
% Verify: g = eps0*(ExB)
syms c_light positive
g_x = Sx / c_light^2;
g_x_sub = subs(g_x, c_light^2, 1/(mu0*eps0));
g_x_expected = eps0 * cross_x;

step9_residual = simplify(g_x_sub - g_x_expected);

total_steps = total_steps + 1;
if isAlways(step9_residual == 0, 'Unknown', 'false')
    fprintf('  Step 9  PASS  g = S/c^2 = eps0*(ExB)\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 9  FAIL  residual: %s\n', char(step9_residual));
    fail_count = fail_count + 1;
end

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

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

% --- Unit check ---
fprintf('  Unit check:\n');
fprintf('    u: [J/m^3], S: [W/m^2]\n');
fprintf('    du/dt: [W/m^3], div(S): [W/m^3]\n');
fprintf('    du/dt + div(S): [W/m^3] = 0\n');
fprintf('    PASS\n\n');

% --- Self-test: wrong Faraday sign ---
% If dB/dt = +curl E (wrong sign), the cancellation fails.
% du/dt_wrong = (1/mu0)*[E.(curl B) + B.(curl E)]
% Poynting sum = du/dt_wrong + (1/mu0)*div(ExB)
%              = (1/mu0)*[E.curlB + B.curlE + B.curlE - E.curlB]
%              = (2/mu0)*B.curlE
du_dt_wrong = (1/mu0)*(E_dot_curlB + B_dot_curlE);
wrong_poynting = du_dt_wrong + (1/mu0)*div_ExB;
wrong_poynting_simplified = simplify(wrong_poynting);

total_steps = total_steps + 1;
if ~isAlways(wrong_poynting_simplified == 0, 'Unknown', 'false')
    fprintf('  Self-test 1: Wrong Faraday sign gives nonzero Poynting residual  PASS\n');
    pass_count = pass_count + 1;
else
    fprintf('  Self-test 1: FAIL (wrong sign not detected)\n');
    fail_count = fail_count + 1;
end

% --- Self-test: quantify wrong residual ---
% Should be (2/mu0)*B.(curl E)
expected_wrong = (2/mu0)*B_dot_curlE;
wrong_quant = simplify(wrong_poynting_simplified - expected_wrong);

total_steps = total_steps + 1;
if isAlways(wrong_quant == 0, 'Unknown', 'false')
    fprintf('  Self-test 2: wrong residual = (2/mu0)*B.(curl E) (quantified)  PASS\n');
    pass_count = pass_count + 1;
else
    fprintf('  Self-test 2: FAIL  residual = %s\n', char(wrong_quant));
    fail_count = fail_count + 1;
end

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

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