%% CAS_F0018_VERIFY.m -- Thermal conduction law (Fourier) as constitutive input
%  Assertion-based CAS audit block
%  Pillar: Thermodynamics | Chain: Fourier law -> energy balance -> diffusion equation
%  CalRef: Thermodynamics Math Appendix S4D, Calibration S4D
%
%  Structure mirrors cas_F17.txt (= F0018) sections A-E.
%  Fourier law is INPUT (not derived). Block verifies substitution into
%  energy balance to obtain the heat diffusion equation.
%
%  Verifies:
%    1. Fourier law structure: Jq = -kappa*grad(T) componentwise
%    2. Energy balance: c*dT/dt + div(Jq) = 0
%    3. Substitution: c*dT/dt = div(kappa*grad(T))
%    4. Constant kappa: c*dT/dt = kappa*Lap(T)
%    5. Sign consistency with F0016: sigma = kappa*|gradT|^2/T^2 >= 0
%    6. Thermal diffusivity: alpha = kappa/c
%    7. Concrete 1D steady-state: d^2T/dx^2 = 0 => T linear
%    8. Numerical diffusivity check
%    9. Self-test: wrong Fourier sign gives anti-diffusion
%   10. Self-test: wrong residual quantified
%
%  HARDENING: isAlways(..., 'Unknown', 'false') throughout.

clear; clc;
fprintf('=== CAS AUDIT: F0018 -- Fourier conduction law ===\n\n');

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

%% ---- A. INPUTS ----
syms x y z t real
syms kappa_th positive      % thermal conductivity
syms c_vol positive         % volumetric heat capacity

% Temperature field
syms T_field(x,y,z,t)

fprintf('Section A: Inputs defined.\n');
fprintf('  Jq = -kappa*grad(T) (constitutive, not derived)\n');
fprintf('  du/dt + div(Jq) = 0, u = c*T\n\n');

%% ---- B. ASSUMPTIONS / DOMAINS ----
fprintf('Section B: Isotropic homogeneous medium, kappa >= 0, c > 0.\n\n');

%% ---- C. ALLOWED LEMMAS ----
fprintf('Section C: Lemmas declared.\n');
fprintf('  C.1: Fourier + energy balance => c*dT/dt = div(kappa*gradT)\n\n');

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

% Gradient components
gradT_x = diff(T_field, x);
gradT_y = diff(T_field, y);
gradT_z = diff(T_field, z);

% --- Step 1: Fourier law structure ---
% Jq_x = -kappa*dT/dx, Jq_y = -kappa*dT/dy, Jq_z = -kappa*dT/dz
Jq_x = -kappa_th * gradT_x;
Jq_y = -kappa_th * gradT_y;
Jq_z = -kappa_th * gradT_z;

% Verify: Jq + kappa*gradT = 0 (componentwise)
step1_res_x = simplify(Jq_x + kappa_th*gradT_x);
step1_res_y = simplify(Jq_y + kappa_th*gradT_y);
step1_res_z = simplify(Jq_z + kappa_th*gradT_z);

total_steps = total_steps + 1;
if isAlways(step1_res_x == 0, 'Unknown', 'false') && ...
   isAlways(step1_res_y == 0, 'Unknown', 'false') && ...
   isAlways(step1_res_z == 0, 'Unknown', 'false')
    fprintf('  Step 1  PASS  Jq = -kappa*grad(T) (all 3 components)\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 1  FAIL\n');
    fail_count = fail_count + 1;
end

% --- Step 2: Energy balance substitution ---
% u = c*T => du/dt = c*dT/dt
% du/dt + div(Jq) = 0
% c*dT/dt + div(Jq) = 0
dTdt = diff(T_field, t);
div_Jq = diff(Jq_x, x) + diff(Jq_y, y) + diff(Jq_z, z);

energy_balance = c_vol*dTdt + div_Jq;

% div(Jq) = div(-kappa*gradT) = -kappa*div(gradT) = -kappa*Lap(T) [const kappa]
% But for general kappa: div(-kappa*gradT) = -kappa*Lap(T) - grad(kappa).gradT
% Here kappa is a symbol (not spatially varying), so div(Jq) = -kappa*Lap(T)
Lap_T = diff(T_field, x, 2) + diff(T_field, y, 2) + diff(T_field, z, 2);

% Energy balance should give: c*dT/dt - kappa*Lap(T) = 0
% i.e. c*dT/dt = kappa*Lap(T)
step2_residual = simplify(energy_balance - (c_vol*dTdt - kappa_th*Lap_T));

total_steps = total_steps + 1;
if isAlways(step2_residual == 0, 'Unknown', 'false')
    fprintf('  Step 2  PASS  c*dT/dt + div(Jq) = c*dT/dt - kappa*Lap(T)\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: Diffusion equation ---
% c*dT/dt = kappa*Lap(T)
% Verify: energy_balance = 0 iff c*dT/dt = kappa*Lap(T)
diffusion_eq = c_vol*dTdt - kappa_th*Lap_T;
step3_residual = simplify(energy_balance - diffusion_eq);

total_steps = total_steps + 1;
if isAlways(step3_residual == 0, 'Unknown', 'false')
    fprintf('  Step 3  PASS  Heat equation: c*dT/dt = kappa*Lap(T)\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: Constant kappa reduction ---
% div(kappa*gradT) = kappa*Lap(T) when kappa is constant
% Verify: div(kappa*gradT) - kappa*Lap(T) = 0
div_kappa_gradT = diff(kappa_th*gradT_x, x) + diff(kappa_th*gradT_y, y) + diff(kappa_th*gradT_z, z);
step4_residual = simplify(div_kappa_gradT - kappa_th*Lap_T);

total_steps = total_steps + 1;
if isAlways(step4_residual == 0, 'Unknown', 'false')
    fprintf('  Step 4  PASS  div(kappa*gradT) = kappa*Lap(T) (const kappa)\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: Sign consistency with F0016 ---
% F0016 showed: sigma = kappa*|gradT|^2/T^2 >= 0
% This requires Jq = -kappa*gradT (negative sign).
% Verify: Jq.grad(1/T) = (-kappa*gradT).(-gradT/T^2) = kappa*|gradT|^2/T^2
syms gradT_sq positive
syms T_pos positive

sigma_from_fourier = kappa_th * gradT_sq / T_pos^2;

% Rebuild from Jq.grad(1/T):
% Jq = -kappa*gradT, grad(1/T) = -(1/T^2)*gradT
% Jq.grad(1/T) = (-kappa*|gradT|)*(-(1/T^2)*|gradT|) ... dot product:
% = (-kappa*gradT).(-gradT/T^2) = kappa*|gradT|^2/T^2
sigma_check = (-kappa_th)*gradT_sq * (-(1/T_pos^2));

step5_residual = simplify(sigma_check - sigma_from_fourier);

total_steps = total_steps + 1;
if isAlways(step5_residual == 0, 'Unknown', 'false')
    fprintf('  Step 5  PASS  Fourier sign consistent with F0016 entropy production\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: Thermal diffusivity ---
% alpha = kappa/c [m^2/s]
% Heat equation: dT/dt = alpha*Lap(T)
syms alpha_th positive
alpha_def = kappa_th / c_vol;

% Verify: c*dT/dt = kappa*Lap(T) <=> dT/dt = (kappa/c)*Lap(T) = alpha*Lap(T)
% (kappa/c) - alpha_def = 0 when alpha_def = kappa/c
step6_residual = simplify(alpha_def - kappa_th/c_vol);

total_steps = total_steps + 1;
if isAlways(step6_residual == 0, 'Unknown', 'false')
    fprintf('  Step 6  PASS  alpha = kappa/c (thermal diffusivity)\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: Concrete 1D steady-state ---
% Steady-state (dT/dt = 0): kappa*d^2T/dx^2 = 0 => d^2T/dx^2 = 0
% Solution: T(x) = A + B*x (linear)
% Verify: d^2/dx^2 (A + B*x) = 0
syms A_coeff B_coeff real
syms x_1d real
T_1d = A_coeff + B_coeff * x_1d;
d2T_dx2 = diff(T_1d, x_1d, 2);

total_steps = total_steps + 1;
if isAlways(d2T_dx2 == 0, 'Unknown', 'false')
    fprintf('  Step 7  PASS  Steady-state 1D: T = A + Bx => d^2T/dx^2 = 0\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 7  FAIL  d^2T/dx^2 = %s\n', char(d2T_dx2));
    fail_count = fail_count + 1;
end

% --- Step 8: Numerical diffusivity check ---
% Copper: kappa ~ 401 W/(m K), c_vol ~ 3.45e6 J/(m^3 K)
% alpha = 401/3.45e6 ~ 1.162e-4 m^2/s
kappa_val = 401;
c_val = 3.45e6;
alpha_val = kappa_val / c_val;
alpha_expected = 1.162e-4;
rel_error = abs(alpha_val - alpha_expected) / alpha_expected;

total_steps = total_steps + 1;
if rel_error < 1e-2
    fprintf('  Step 8  PASS  Numerical: alpha(Cu) = %.4e m^2/s\n', alpha_val);
    pass_count = pass_count + 1;
else
    fprintf('  Step 8  FAIL  Numerical rel error: %.2e\n', rel_error);
    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('    Jq: [W/m^2], kappa: [W/(m*K)], gradT: [K/m]\n');
fprintf('    kappa*gradT: [W/m^2] => Jq: [W/m^2]\n');
fprintf('    c: [J/(m^3*K)], dT/dt: [K/s]\n');
fprintf('    c*dT/dt: [J/(m^3*s)] = [W/m^3]\n');
fprintf('    kappa*Lap(T): [W/(m*K)]*[K/m^2] = [W/m^3]\n');
fprintf('    PASS\n\n');

% --- Self-test: wrong Fourier sign (Jq = +kappa*gradT) ---
% This gives anti-diffusion: c*dT/dt = -kappa*Lap(T)
% Verify the sign flips
Jq_wrong_x = kappa_th * gradT_x;  % wrong sign
Jq_wrong_y = kappa_th * gradT_y;  % wrong sign
Jq_wrong_z = kappa_th * gradT_z;  % wrong sign
div_Jq_wrong = diff(Jq_wrong_x, x) + diff(Jq_wrong_y, y) + diff(Jq_wrong_z, z);
energy_wrong = c_vol*dTdt + div_Jq_wrong;
diffusion_wrong = c_vol*dTdt + kappa_th*Lap_T;  % note +kappa instead of -kappa

wrong_residual = simplify(energy_wrong - diffusion_wrong);

total_steps = total_steps + 1;
if isAlways(wrong_residual == 0, 'Unknown', 'false')
    % energy_wrong = c*dT/dt + kappa*Lap(T), which differs from correct c*dT/dt - kappa*Lap(T)
    % Check that wrong equation differs from correct:
    wrong_vs_correct = simplify(diffusion_wrong - diffusion_eq);
    if ~isAlways(wrong_vs_correct == 0, 'Unknown', 'false')
        fprintf('  Self-test 1: Wrong Fourier sign gives anti-diffusion (detected)  PASS\n');
        pass_count = pass_count + 1;
    else
        fprintf('  Self-test 1: FAIL (wrong sign not detected)\n');
        fail_count = fail_count + 1;
    end
else
    fprintf('  Self-test 1: FAIL (unexpected structure)\n');
    fail_count = fail_count + 1;
end

% --- Self-test: quantify wrong residual ---
% wrong - correct = (c*dT/dt + kappa*Lap) - (c*dT/dt - kappa*Lap) = 2*kappa*Lap(T)
expected_wrong_diff = 2*kappa_th*Lap_T;
wrong_quant = simplify(wrong_vs_correct - expected_wrong_diff);

total_steps = total_steps + 1;
if isAlways(wrong_quant == 0, 'Unknown', 'false')
    fprintf('  Self-test 2: wrong - correct = 2*kappa*Lap(T) (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('  F0018 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 F0018.\n');