%% CAS_F0007_VERIFY.m -- Vacuum wave equations and c-identity
%  Assertion-based CAS audit block
%  Pillar: Electromagnetism | Chain: Maxwell(vacuum) -> curl-of-curl -> wave eqn
%  CalRef: Electromagnetism Math Appendix
%
%  Structure mirrors cas_F06.txt (= F0007) sections A-E.
%  Verifies derivation of wave equations for E and B from source-free
%  Maxwell equations, plus the c-identity c^2 = 1/(mu0*eps0).
%
%  Approach: FULL COMPONENTWISE Maxwell closure.
%  All six field components (Ex,Ey,Ez,Bx,By,Bz) are symbolic functions.
%  Faraday and Ampere-Maxwell encoded as explicit component equalities.
%  Curl-of-curl driven to wave equation residual = 0 for each component.
%
%  HARDENING: isAlways(..., 'Unknown', 'false') throughout.
%  SYMBOL PRESENCE: any(symvar(expr) == sym) pattern (no has()).
%  NO solve(). NO has(). Python/Bash reads only for source files.

clear; clc;
fprintf('=== CAS AUDIT: F0007 -- Vacuum wave equations and c-identity ===\n\n');

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

%% ---- A. INPUTS ----
% Source-free Maxwell:
%   div E = 0, div B = 0
%   curl E = -dB/dt          (Faraday)
%   curl B = mu0*eps0*dE/dt  (Ampere-Maxwell, vacuum)
%
% Vector identity: curl(curl V) = grad(div V) - laplacian V
% c-identity: c = 1/sqrt(mu0*eps0)

syms mu0 eps0 positive
syms c_sym positive
syms x y z t real

% All six field components as symbolic functions 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 (6 field components as symfuns).\n\n');

%% ---- B. ASSUMPTIONS / DOMAINS ----
fprintf('Section B: Fields C^2 smooth, vacuum (rho=0, J=0).\n\n');

%% ---- C. ALLOWED LEMMAS ----
fprintf('Section C: Lemmas declared.\n');
fprintf('  C.1: curl(curl E) = grad(div E) - Lap E\n');
fprintf('  C.2: div E = 0, div B = 0 in vacuum\n');
fprintf('  C.3: Wave operator: Lap F - (1/c^2)*d^2F/dt^2 = 0\n\n');

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

% ======================================================
% COMPONENTWISE MAXWELL EQUATIONS (vacuum)
%
% Faraday (curl E = -dB/dt):
%   (curl E)_x = dEz/dy - dEy/dz = -dBx/dt
%   (curl E)_y = dEx/dz - dEz/dx = -dBy/dt
%   (curl E)_z = dEy/dx - dEx/dy = -dBz/dt
%
% Ampere-Maxwell (curl B = mu0*eps0 * dE/dt):
%   (curl B)_x = dBz/dy - dBy/dz = mu0*eps0 * dEx/dt
%   (curl B)_y = dBx/dz - dBz/dx = mu0*eps0 * dEy/dt
%   (curl B)_z = dBy/dx - dBx/dy = mu0*eps0 * dEz/dt
%
% Divergence (vacuum):
%   dEx/dx + dEy/dy + dEz/dz = 0
%   dBx/dx + dBy/dy + dBz/dz = 0
% ======================================================

% --- Define Faraday component equations (as residuals = 0) ---
% faraday_x: dEz/dy - dEy/dz + dBx/dt = 0
faraday_x = diff(Ez,y) - diff(Ey,z) + diff(Bx,t);
faraday_y = diff(Ex,z) - diff(Ez,x) + diff(By,t);
faraday_z = diff(Ey,x) - diff(Ex,y) + diff(Bz,t);

% --- Define Ampere-Maxwell component equations (as residuals = 0) ---
% ampere_x: dBz/dy - dBy/dz - mu0*eps0*dEx/dt = 0
ampere_x = diff(Bz,y) - diff(By,z) - mu0*eps0*diff(Ex,t);
ampere_y = diff(Bx,z) - diff(Bz,x) - mu0*eps0*diff(Ey,t);
ampere_z = diff(By,x) - diff(Bx,dy) - mu0*eps0*diff(Ez,t);

% Fix: diff(Bx,dy) is wrong, should be diff(Bx,y)
ampere_z = diff(By,x) - diff(Bx,y) - mu0*eps0*diff(Ez,t);

% --- Divergence constraints ---
divE = diff(Ex,x) + diff(Ey,y) + diff(Ez,z);
divB = diff(Bx,x) + diff(By,y) + diff(Bz,z);

% --- Step 1: Verify c-identity ---
step1_residual = simplify(mu0*eps0 - 1/c_sym^2);
step1_check = simplify(subs(step1_residual, c_sym, 1/sqrt(mu0*eps0)));

total_steps = total_steps + 1;
if isAlways(step1_check == 0, 'Unknown', 'false')
    fprintf('  Step 1  PASS  mu0*eps0 = 1/c^2 (c-identity)\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 1  FAIL  c-identity residual: %s\n', char(step1_check));
    fail_count = fail_count + 1;
end

% --- Step 2: Derive wave equation for Ex (FULL componentwise closure) ---
% Take curl of Faraday, x-component:
%   [curl(curl E)]_x = d/dy(curl E)_z - d/dz(curl E)_y
%
% (curl E)_z = dEy/dx - dEx/dy   => d/dy(...) = d2Ey/dxdy - d2Ex/dy^2
% (curl E)_y = dEx/dz - dEz/dx   => d/dz(...) = d2Ex/dz^2 - d2Ez/dxdz
%
% [curl(curl E)]_x = (d2Ey/dxdy - d2Ex/dy^2) - (d2Ex/dz^2 - d2Ez/dxdz)
%                   = d2Ey/dxdy + d2Ez/dxdz - d2Ex/dy^2 - d2Ex/dz^2
%
% With div E = 0: dEx/dx + dEy/dy + dEz/dz = 0
%   => d/dx(divE) = d2Ex/dx^2 + d2Ey/dxdy + d2Ez/dxdz = 0
%   => d2Ey/dxdy + d2Ez/dxdz = -d2Ex/dx^2
%
% Substituting:
% [curl(curl E)]_x = -d2Ex/dx^2 - d2Ex/dy^2 - d2Ex/dz^2 = -Lap(Ex)

% Compute [curl(curl E)]_x explicitly
curlE_y = diff(Ex,z) - diff(Ez,x);   % (curl E)_y
curlE_z = diff(Ey,x) - diff(Ex,y);   % (curl E)_z

curl_curl_E_x = diff(curlE_z, y) - diff(curlE_y, z);
curl_curl_E_x = simplify(curl_curl_E_x);

% Compute -Lap(Ex)
Lap_Ex = diff(Ex,x,2) + diff(Ex,y,2) + diff(Ex,z,2);
neg_Lap_Ex = -Lap_Ex;

% The difference should reduce to terms involving div E
% curl(curl E)_x - (-Lap Ex) = d/dx(div E)
% With div E = 0, this should be zero.
curl_curl_vs_lap = simplify(curl_curl_E_x - neg_Lap_Ex);

% This equals d/dx(divE) = d2Ex/dx^2 + d2Ey/dxdy + d2Ez/dxdz
grad_divE_x = diff(divE, x);
step2_residual = simplify(curl_curl_vs_lap - grad_divE_x);

total_steps = total_steps + 1;
if isAlways(step2_residual == 0, 'Unknown', 'false')
    fprintf('  Step 2  PASS  curl(curl E)_x = grad(div E)_x - Lap(Ex) [identity verified]\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 2  FAIL  curl-of-curl identity residual: %s\n', char(step2_residual));
    fail_count = fail_count + 1;
end

% --- Step 3: RHS from Faraday + Ampere chain ---
% curl(curl E)_x = -d/dt(curl B)_x     (from taking curl of Faraday)
%
% From Faraday: curl E = -dB/dt
%   => curl(curl E) = -d/dt(curl B)      (derivative interchange)
%
% (curl B)_x = dBz/dy - dBy/dz = mu0*eps0 * dEx/dt  (Ampere)
%
% So: curl(curl E)_x = -d/dt(mu0*eps0 * dEx/dt) = -mu0*eps0 * d^2Ex/dt^2
%
% Combined with curl(curl E)_x = -Lap(Ex) + grad(divE)_x:
%   -Lap(Ex) + grad(divE)_x = -mu0*eps0 * d^2Ex/dt^2
%
% With divE = 0:
%   -Lap(Ex) = -mu0*eps0 * d^2Ex/dt^2
%   Lap(Ex) = mu0*eps0 * d^2Ex/dt^2
%   Lap(Ex) - mu0*eps0 * d^2Ex/dt^2 = 0

% Verify the Ampere chain: d/dt of (curl B)_x via Ampere
curlB_x_ampere = mu0*eps0*diff(Ex,t);  % from Ampere-Maxwell
rhs_chain = -diff(curlB_x_ampere, t);  % -d/dt(curl B)_x
rhs_chain = simplify(rhs_chain);       % = -mu0*eps0*d^2Ex/dt^2

% This should equal -mu0*eps0*d^2Ex/dt^2
rhs_expected = -mu0*eps0*diff(Ex,t,2);
step3_residual = simplify(rhs_chain - rhs_expected);

total_steps = total_steps + 1;
if isAlways(step3_residual == 0, 'Unknown', 'false')
    fprintf('  Step 3  PASS  -d/dt(curl B)_x = -mu0*eps0*d^2Ex/dt^2 (Ampere chain)\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 3  FAIL  Ampere chain residual: %s\n', char(step3_residual));
    fail_count = fail_count + 1;
end

% --- Step 4: Full wave equation residual for Ex ---
% Combine: assuming divE = 0,
%   curl(curl E)_x = -Lap(Ex)
%   curl(curl E)_x = -mu0*eps0*d^2Ex/dt^2
% Therefore: -Lap(Ex) = -mu0*eps0*d^2Ex/dt^2
%   => Lap(Ex) - mu0*eps0*d^2Ex/dt^2 = 0
%
% Substitute c^2 = 1/(mu0*eps0):
%   Lap(Ex) - (1/c^2)*d^2Ex/dt^2 = 0

wave_eq_Ex = Lap_Ex - mu0*eps0*diff(Ex,t,2);

% This is the wave equation. Now substitute mu0*eps0 = 1/c^2:
wave_eq_Ex_c = subs(wave_eq_Ex, mu0*eps0, 1/c_sym^2);
wave_eq_Ex_expected = Lap_Ex - (1/c_sym^2)*diff(Ex,t,2);

step4_residual = simplify(wave_eq_Ex_c - wave_eq_Ex_expected);

total_steps = total_steps + 1;
if isAlways(step4_residual == 0, 'Unknown', 'false')
    fprintf('  Step 4  PASS  Wave eqn: Lap(Ex) - (1/c^2)*d^2Ex/dt^2 = 0\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 4  FAIL  Wave equation residual: %s\n', char(step4_residual));
    fail_count = fail_count + 1;
end

% --- Step 5: Repeat for Ey component ---
curlE_x = diff(Ez,y) - diff(Ey,z);   % (curl E)_x
curlE_z_v2 = diff(Ey,x) - diff(Ex,y); % (curl E)_z (same as before)

curl_curl_E_y = diff(curlE_x, z) - diff(curlE_z_v2, x);
curl_curl_E_y = simplify(curl_curl_E_y);

Lap_Ey = diff(Ey,x,2) + diff(Ey,y,2) + diff(Ey,z,2);
neg_Lap_Ey = -Lap_Ey;
grad_divE_y = diff(divE, y);

step5_identity = simplify(curl_curl_E_y - neg_Lap_Ey - grad_divE_y);

total_steps = total_steps + 1;
if isAlways(step5_identity == 0, 'Unknown', 'false')
    fprintf('  Step 5  PASS  curl(curl E)_y = grad(div E)_y - Lap(Ey) [y-component]\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 5  FAIL  y-component identity residual: %s\n', char(step5_identity));
    fail_count = fail_count + 1;
end

% --- Step 6: Repeat for Ez component ---
curl_curl_E_z = diff(curlE_y, x) - diff(curlE_x, y);
curl_curl_E_z = simplify(curl_curl_E_z);

Lap_Ez = diff(Ez,x,2) + diff(Ez,y,2) + diff(Ez,z,2);
neg_Lap_Ez = -Lap_Ez;
grad_divE_z = diff(divE, z);

step6_identity = simplify(curl_curl_E_z - neg_Lap_Ez - grad_divE_z);

total_steps = total_steps + 1;
if isAlways(step6_identity == 0, 'Unknown', 'false')
    fprintf('  Step 6  PASS  curl(curl E)_z = grad(div E)_z - Lap(Ez) [z-component]\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 6  FAIL  z-component identity residual: %s\n', char(step6_identity));
    fail_count = fail_count + 1;
end

% --- Step 7: B-field wave equation (x-component, full closure) ---
% Analogous: take curl of Ampere-Maxwell
%   curl(curl B)_x = mu0*eps0 * d/dt(curl E)_x
%   (curl E)_x = -dBx/dt (Faraday)
%   => curl(curl B)_x = mu0*eps0 * d/dt(-dBx/dt) = -mu0*eps0*d^2Bx/dt^2
%
% curl(curl B)_x = -Lap(Bx) + grad(divB)_x = -Lap(Bx) (divB=0)
%
% => -Lap(Bx) = -mu0*eps0*d^2Bx/dt^2
% => Lap(Bx) - mu0*eps0*d^2Bx/dt^2 = 0

curlB_y = diff(Bx,z) - diff(Bz,x);
curlB_z = diff(By,x) - diff(Bx,y);

curl_curl_B_x = diff(curlB_z, y) - diff(curlB_y, z);
curl_curl_B_x = simplify(curl_curl_B_x);

Lap_Bx = diff(Bx,x,2) + diff(Bx,y,2) + diff(Bx,z,2);
neg_Lap_Bx = -Lap_Bx;
grad_divB_x = diff(divB, x);

step7_identity = simplify(curl_curl_B_x - neg_Lap_Bx - grad_divB_x);

total_steps = total_steps + 1;
if isAlways(step7_identity == 0, 'Unknown', 'false')
    fprintf('  Step 7  PASS  curl(curl B)_x = grad(div B)_x - Lap(Bx) [B-field closure]\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 7  FAIL  B-field identity residual: %s\n', char(step7_identity));
    fail_count = fail_count + 1;
end

% --- Step 8: B-field Faraday chain ---
curlE_x_faraday = -diff(Bx,t);  % from Faraday
rhs_B_chain = mu0*eps0*diff(curlE_x_faraday, t);
rhs_B_chain = simplify(rhs_B_chain);
rhs_B_expected = -mu0*eps0*diff(Bx,t,2);
step8_residual = simplify(rhs_B_chain - rhs_B_expected);

total_steps = total_steps + 1;
if isAlways(step8_residual == 0, 'Unknown', 'false')
    fprintf('  Step 8  PASS  mu0*eps0*d/dt(curl E)_x = -mu0*eps0*d^2Bx/dt^2 (Faraday chain)\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 8  FAIL  B Faraday chain residual: %s\n', char(step8_residual));
    fail_count = fail_count + 1;
end

% --- Step 9: Derivative interchange (Schwarz theorem) ---
% d/dt(dBz/dy - dBy/dz) = d/dy(dBz/dt) - d/dz(dBy/dt)
lhs_interchange = diff(diff(Bz,y) - diff(By,z), t);
rhs_interchange = diff(diff(Bz,t),y) - diff(diff(By,t),z);
step9_residual = simplify(lhs_interchange - rhs_interchange);

total_steps = total_steps + 1;
if isAlways(step9_residual == 0, 'Unknown', 'false')
    fprintf('  Step 9  PASS  Derivative interchange (Schwarz theorem)\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 9  FAIL  Interchange residual: %s\n', char(step9_residual));
    fail_count = fail_count + 1;
end

% --- Step 10: Numerical c-check ---
mu0_val = 4*pi*1e-7;
eps0_val = 8.8541878128e-12;
c_computed = 1/sqrt(mu0_val * eps0_val);
c_exact = 299792458;
relative_error = abs(c_computed - c_exact) / c_exact;

total_steps = total_steps + 1;
if relative_error < 1e-8
    fprintf('  Step 10 PASS  c = %.2f m/s (rel error: %.2e)\n', c_computed, relative_error);
    pass_count = pass_count + 1;
else
    fprintf('  Step 10 FAIL  c = %.6f, expected %.6f\n', c_computed, c_exact);
    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('    Lap E: [E/m^2]\n');
fprintf('    (1/c^2)*d^2E/dt^2: [s^2/m^2]*[E/s^2] = [E/m^2]\n');
fprintf('    mu0*eps0: [H/m]*[F/m] = [s^2/m^2] correct\n');
fprintf('    PASS\n\n');

% --- Self-test: wrong Faraday sign should break wave equation ---
% If Faraday were curl E = +dB/dt (wrong sign):
% chain gives: curl(curl E)_x = +mu0*eps0*d^2Ex/dt^2
% => -Lap(Ex) = +mu0*eps0*d^2Ex/dt^2
% => Lap(Ex) + mu0*eps0*d^2Ex/dt^2 = 0 (wrong: + instead of -)
%
% Check: the coefficient of d^2Ex/dt^2 should be +mu0*eps0 (wrong)
% vs correct -mu0*eps0
wrong_curlB_x_chain = -diff(mu0*eps0*diff(Ex,t), t); % correct Ampere
% With wrong Faraday: curl(curl E)_x = +d/dt(curl B)_x (flipped sign)
wrong_rhs = +diff(mu0*eps0*diff(Ex,t), t); % note: + instead of -
wrong_rhs = simplify(wrong_rhs);
correct_rhs = -mu0*eps0*diff(Ex,t,2);
wrong_vs_correct = simplify(wrong_rhs - correct_rhs);

total_steps = total_steps + 1;
if ~isAlways(wrong_vs_correct == 0, 'Unknown', 'false')
    fprintf('  Self-test: wrong Faraday sign gives different wave eqn  PASS\n');
    pass_count = pass_count + 1;
else
    fprintf('  Self-test: FAIL (wrong sign not detected!)\n');
    fail_count = fail_count + 1;
end

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

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