%% CAS_F0025_VERIFY.m — Poynting theorem (∂u/∂t + ∇·S + J·E = 0)
%  Pillar : Electromagnetism
%  CalRef : Math Appendix §4.2–4.3, EM Calibration §4B
%  Chain  : u definition → ∂u/∂t → Maxwell substitution → vector identity → S → theorem
%  Hardening: isAlways(..., 'Unknown', 'false') on every symbolic assertion
%            symvar-based presence checks; solver and has-function avoided

clear; clc;
fprintf('\n=== CAS F0025 : Poynting theorem (with sources) ===\n\n');

%% ── Symbols ──────────────────────────────────────────────────
syms x y z t real
syms mu0 eps0 real positive

%  E, B, J fields 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)
syms Jx(x,y,z,t) Jy(x,y,z,t) Jz(x,y,z,t)

%% ── A. Inputs ────────────────────────────────────────────────
%  Energy density: u = (1/2)(ε₀E² + B²/μ₀)
E_sq = Ex^2 + Ey^2 + Ez^2;
B_sq = Bx^2 + By^2 + Bz^2;
u = (eps0 * E_sq + B_sq / mu0) / 2;

%  Poynting vector: S = (1/μ₀)(E × B)
Sx = (Ey*Bz - Ez*By) / mu0;
Sy = (Ez*Bx - Ex*Bz) / mu0;
Sz = (Ex*By - Ey*Bx) / mu0;

%  Ohmic dissipation: J·E
JdotE = Jx*Ex + Jy*Ey + Jz*Ez;

%% ── Step 1: ∂u/∂t chain rule ────────────────────────────────
%  ∂u/∂t = ε₀(E·∂E/∂t) + (1/μ₀)(B·∂B/∂t)
du_dt = diff(u, t);
du_dt_expected = eps0*(Ex*diff(Ex,t) + Ey*diff(Ey,t) + Ez*diff(Ez,t)) ...
              + (Bx*diff(Bx,t) + By*diff(By,t) + Bz*diff(Bz,t))/mu0;
res1 = simplify(du_dt - du_dt_expected);
assert(isAlways(res1 == 0, 'Unknown', 'false'), ...
    'FAIL Step 1: ∂u/∂t chain rule');
fprintf('Step 1  PASS — ∂u/∂t = ε₀(E·∂E/∂t) + (1/μ₀)(B·∂B/∂t)\n');

%% ── Step 2: Maxwell substitution for ∂B/∂t ─────────────────
%  Faraday: ∂B/∂t = −∇×E
%  (∂Bx/∂t) = −(∂Ez/∂y − ∂Ey/∂z)
%  (∂By/∂t) = −(∂Ex/∂z − ∂Ez/∂x)
%  (∂Bz/∂t) = −(∂Ey/∂x − ∂Ex/∂y)
dBx_dt_faraday = -(diff(Ez,y) - diff(Ey,z));
dBy_dt_faraday = -(diff(Ex,z) - diff(Ez,x));
dBz_dt_faraday = -(diff(Ey,x) - diff(Ex,y));

%  B·(∂B/∂t) with Faraday substitution
BdotdBdt_faraday = Bx*dBx_dt_faraday + By*dBy_dt_faraday + Bz*dBz_dt_faraday;

%% ── Step 3: Maxwell substitution for ∂E/∂t ─────────────────
%  Ampère–Maxwell: ∂E/∂t = (1/(μ₀ε₀))(∇×B − μ₀J)
curlB_x = diff(Bz,y) - diff(By,z);
curlB_y = diff(Bx,z) - diff(Bz,x);
curlB_z = diff(By,x) - diff(Bx,y);

dEx_dt_ampere = (curlB_x - mu0*Jx) / (mu0*eps0);
dEy_dt_ampere = (curlB_y - mu0*Jy) / (mu0*eps0);
dEz_dt_ampere = (curlB_z - mu0*Jz) / (mu0*eps0);

%  E·(∂E/∂t) with Ampère substitution
EdotdEdt_ampere = Ex*dEx_dt_ampere + Ey*dEy_dt_ampere + Ez*dEz_dt_ampere;

%% ── Step 4: du/dt after Maxwell substitution ────────────────
%  du/dt = ε₀(E·∂E/∂t) + (1/μ₀)(B·∂B/∂t)
du_dt_maxwell = eps0*EdotdEdt_ampere + BdotdBdt_faraday/mu0;

%  Expected from derivation Step 2.3:
%  du/dt = (1/μ₀)[E·(∇×B) − B·(∇×E)] − J·E
curlE_x = diff(Ez,y) - diff(Ey,z);
curlE_y = diff(Ex,z) - diff(Ez,x);
curlE_z = diff(Ey,x) - diff(Ex,y);

EdotcurlB = Ex*curlB_x + Ey*curlB_y + Ez*curlB_z;
BdotcurlE = Bx*curlE_x + By*curlE_y + Bz*curlE_z;

du_dt_step23 = (EdotcurlB - BdotcurlE)/mu0 - JdotE;
res4 = simplify(du_dt_maxwell - du_dt_step23);
assert(isAlways(res4 == 0, 'Unknown', 'false'), ...
    'FAIL Step 4: Maxwell substitution ≠ Step 2.3 form');
fprintf('Step 4  PASS — ∂u/∂t = (1/μ₀)[E·(∇×B)−B·(∇×E)] − J·E\n');

%% ── Step 5: Vector identity ∇·(E×B) ────────────────────────
%  ∇·(E×B) = B·(∇×E) − E·(∇×B)
divExB = diff(Ey*Bz - Ez*By, x) + diff(Ez*Bx - Ex*Bz, y) + diff(Ex*By - Ey*Bx, z);
identity_rhs = BdotcurlE - EdotcurlB;
res5 = simplify(divExB - identity_rhs);
assert(isAlways(res5 == 0, 'Unknown', 'false'), ...
    'FAIL Step 5: ∇·(E×B) ≠ B·(∇×E) − E·(∇×B)');
fprintf('Step 5  PASS — ∇·(E×B) = B·(∇×E) − E·(∇×B) (vector identity)\n');

%% ── Step 6: ∇·S from vector identity ───────────────────────
%  S = (1/μ₀)(E×B) → ∇·S = (1/μ₀)∇·(E×B) = (1/μ₀)[B·(∇×E) − E·(∇×B)]
divS = diff(Sx, x) + diff(Sy, y) + diff(Sz, z);
divS_expected = divExB / mu0;
res6 = simplify(divS - divS_expected);
assert(isAlways(res6 == 0, 'Unknown', 'false'), ...
    'FAIL Step 6: ∇·S ≠ (1/μ₀)∇·(E×B)');
fprintf('Step 6  PASS — ∇·S = (1/μ₀)∇·(E×B)\n');

%% ── Step 7: Poynting theorem assembly ───────────────────────
%  From Steps 4 and 5:
%  ∂u/∂t = −(1/μ₀)∇·(E×B) − J·E = −∇·S − J·E
%  → ∂u/∂t + ∇·S + J·E = 0
%  Verify: du_dt_maxwell + divS + JdotE = 0
%  Substitute divS via identity:
poynting_sum = du_dt_maxwell + divS + JdotE;
poynting_sum_simplified = simplify(poynting_sum);
assert(isAlways(poynting_sum_simplified == 0, 'Unknown', 'false'), ...
    'FAIL Step 7: ∂u/∂t + ∇·S + J·E ≠ 0');
fprintf('Step 7  PASS — ∂u/∂t + ∇·S + J·E = 0 (Poynting theorem)\n');

%% ── Step 8: Source-free limit (J=0) → F0015 consistency ────
%  When J=0: ∂u/∂t + ∇·S = 0
poynting_J0 = subs(poynting_sum, [Jx, Jy, Jz], [sym(0), sym(0), sym(0)]);
res8 = simplify(poynting_J0);
assert(isAlways(res8 == 0, 'Unknown', 'false'), ...
    'FAIL Step 8: source-free Poynting ≠ 0');
fprintf('Step 8  PASS — J=0 → ∂u/∂t + ∇·S = 0 (consistent with F0015)\n');

%% ── Step 9: Numerical — resistive dissipation ──────────────
%  σ = 5.96e7 S/m (copper), E = 0.01 V/m, J = σE = 5.96e5 A/m²
%  J·E = σE² = 5960 W/m³
sigma_Cu = 5.96e7; E_val = 0.01;
JE_val = sigma_Cu * E_val^2;
JE_expected = 5960;
assert(abs(JE_val - JE_expected) < 1, ...
    'FAIL Step 9: J·E numerical mismatch');
fprintf('Step 9  PASS — Cu dissipation: J·E = %.0f W/m³ (σ=5.96e7, E=0.01 V/m)\n', JE_val);

%% ── Step 10: Self-test — wrong Faraday sign ────────────────
%  Wrong: ∂B/∂t = +∇×E (positive instead of negative)
dBx_wrong = (diff(Ez,y) - diff(Ey,z));  % positive, missing minus
dBy_wrong = (diff(Ex,z) - diff(Ez,x));
dBz_wrong = (diff(Ey,x) - diff(Ex,y));

BdotdBdt_wrong = Bx*dBx_wrong + By*dBy_wrong + Bz*dBz_wrong;
du_dt_wrong = eps0*EdotdEdt_ampere + BdotdBdt_wrong/mu0;
poynting_wrong = du_dt_wrong + divS + JdotE;
poynting_wrong_s = simplify(poynting_wrong);
assert(~isAlways(poynting_wrong_s == 0, 'Unknown', 'false'), ...
    'FAIL Step 10: wrong Faraday sign not detected');
fprintf('Step 10 PASS — Wrong Faraday sign (+∇×E) breaks Poynting theorem\n');

%% ── Step 11: Self-test — missing J·E term ──────────────────
%  Wrong: ∂u/∂t + ∇·S = 0 (omitting dissipation when J≠0)
%  Residual should be J·E ≠ 0 for general J, E
poynting_noJ = du_dt_maxwell + divS;
poynting_noJ_s = simplify(poynting_noJ);
%  This should equal −J·E, not zero
res11_check = simplify(poynting_noJ_s + JdotE);
assert(isAlways(res11_check == 0, 'Unknown', 'false'), ...
    'FAIL Step 11: missing J·E residual not quantified');
fprintf('Step 11 PASS — Missing J·E: residual = −J·E (quantified)\n');

%% ── Summary ─────────────────────────────────────────────────
fprintf('\n✓ ALL 11 ASSERTIONS PASSED — F0025 audit complete.\n');
fprintf('  Chain: u def → du/dt → Maxwell sub → vector identity → S → Poynting\n');
fprintf('  Cross-ref: F0015 (source-free Poynting), F0022 (Faraday), F0023 (Ampere-Maxwell)\n');