%% CAS_F0030_VERIFY.m — Heat diffusion equation (∂T/∂t = κ∇²T)
%  Pillar : Thermodynamics
%  CalRef : Math Appendix §4D, Thermodynamics Calibration §4D
%  Chain  : energy conservation + Fourier's law + u=ρcT → diffusion equation
%  Hardening: isAlways(..., 'Unknown', 'false') on every symbolic assertion
%            symvar-based presence checks; solver and has-function avoided

clear; clc;
fprintf('\n=== CAS F0030 : Heat diffusion equation ===\n\n');

%% ── Symbols ──────────────────────────────────────────────────
syms x t real
syms rho c_sp k_th kappa real positive   % density, specific heat, conductivity, diffusivity
syms alpha_w real positive               % wave number for test solutions

%% ── A. Inputs ────────────────────────────────────────────────
%  1. Local energy conservation: ∂u/∂t + ∇·J_q = 0
%  2. Fourier's law: J_q = −k∇T
%  3. Linear relation: u = ρcT
%  Combine to get: ∂T/∂t = κ∇²T  with κ = k/(ρc)

%% ── Step 1: Fourier substitution into conservation (1D) ──────
%  ∂u/∂t + ∂J_q/∂x = 0, J_q = −k ∂T/∂x
%  → ∂u/∂t − k ∂²T/∂x² = 0
%  With u = ρcT → ρc ∂T/∂t − k ∂²T/∂x² = 0
%  Test with concrete T(x,t): sinusoidal diffusion mode
%  T(x,t) = exp(−κα²t)·sin(αx)
T_sol = exp(-kappa * alpha_w^2 * t) * sin(alpha_w * x);
dT_dt = diff(T_sol, t);
d2T_dx2 = diff(T_sol, x, 2);

%  Verify ∂T/∂t = κ ∂²T/∂x²
res1 = simplify(dT_dt - kappa * d2T_dx2);
assert(isAlways(res1 == 0, 'Unknown', 'false'), ...
    'FAIL Step 1: sinusoidal mode does not satisfy diffusion equation');
fprintf('Step 1  PASS — T=e^{−κα²t}sin(αx) satisfies ∂T/∂t = κ∇²T\n');

%% ── Step 2: Fourier substitution — algebraic chain ──────────
%  Starting from conservation ∂u/∂t + ∇·J_q = 0 with J_q = −k∇T:
%  ∂u/∂t − k∇²T = 0
%  With u = ρcT: ρc ∂T/∂t − k∇²T = 0
%  Verify: ρc·(∂T/∂t) − k·(∂²T/∂x²) = 0 for our test solution
%  when κ = k/(ρc)
conservation_residual = rho * c_sp * dT_dt - k_th * d2T_dx2;
%  Substitute κ = k/(ρc):
conservation_sub = subs(conservation_residual, kappa, k_th/(rho*c_sp));
res2 = simplify(conservation_sub);
assert(isAlways(res2 == 0, 'Unknown', 'false'), ...
    'FAIL Step 2: conservation chain ρc∂T/∂t − k∂²T/∂x² ≠ 0');
fprintf('Step 2  PASS — ρc·∂T/∂t − k·∂²T/∂x² = 0 when κ=k/(ρc)\n');

%% ── Step 3: Thermal diffusivity definition ───────────────────
%  κ = k/(ρc) → k − ρcκ = 0
res3 = simplify(k_th - rho * c_sp * kappa);
%  This is the defining relation; substitute κ = k/(ρc):
res3_sub = subs(res3, kappa, k_th/(rho*c_sp));
assert(isAlways(res3_sub == 0, 'Unknown', 'false'), ...
    'FAIL Step 3: κ ≠ k/(ρc)');
fprintf('Step 3  PASS — κ = k/(ρc) (thermal diffusivity definition)\n');

%% ── Step 4: Cosine mode also satisfies diffusion ─────────────
%  T(x,t) = exp(−κα²t)·cos(αx) should also be a solution
T_cos = exp(-kappa * alpha_w^2 * t) * cos(alpha_w * x);
dT_cos_dt = diff(T_cos, t);
d2T_cos_dx2 = diff(T_cos, x, 2);
res4 = simplify(dT_cos_dt - kappa * d2T_cos_dx2);
assert(isAlways(res4 == 0, 'Unknown', 'false'), ...
    'FAIL Step 4: cosine mode does not satisfy diffusion equation');
fprintf('Step 4  PASS — T=e^{−κα²t}cos(αx) also satisfies ∂T/∂t = κ∇²T\n');

%% ── Step 5: Superposition (linearity check) ─────────────────
%  If T₁ and T₂ are solutions, aT₁ + bT₂ is also a solution
syms a_coeff b_coeff real
T_super = a_coeff * T_sol + b_coeff * T_cos;
dT_super_dt = diff(T_super, t);
d2T_super_dx2 = diff(T_super, x, 2);
res5 = simplify(dT_super_dt - kappa * d2T_super_dx2);
assert(isAlways(res5 == 0, 'Unknown', 'false'), ...
    'FAIL Step 5: superposition violates diffusion equation');
fprintf('Step 5  PASS — Superposition preserved (linearity of diffusion eq)\n');

%% ── Step 6: Two-sided localisation ───────────────────────────
%  Correct solution satisfies PDE (shown above)
%  Wrong equation: ∂T/∂t = κ ∂T/∂x (first derivative, not second)
%  Check that our sinusoidal mode does NOT satisfy the wrong equation
wrong_residual = simplify(dT_dt - kappa * diff(T_sol, x));
assert(~isAlways(wrong_residual == 0, 'Unknown', 'false'), ...
    'FAIL Step 6: wrong equation ∂T/∂t = κ∂T/∂x not rejected');
fprintf('Step 6  PASS — Wrong PDE (∂T/∂t = κ∂T/∂x) correctly rejected\n');

%% ── Step 7: Steady-state limit ───────────────────────────────
%  When ∂T/∂t = 0: κ∇²T = 0 → ∇²T = 0 (Laplace equation)
%  In 1D: ∂²T/∂x² = 0 → T(x) = Ax + B (linear profile)
%  Verify: linear profile has zero second derivative
syms A_lin B_lin real
T_steady = A_lin * x + B_lin;
d2T_steady = diff(T_steady, x, 2);
res7 = simplify(d2T_steady);
assert(isAlways(res7 == 0, 'Unknown', 'false'), ...
    'FAIL Step 7: linear profile does not satisfy ∇²T=0');
fprintf('Step 7  PASS — Steady state: T=Ax+B satisfies ∇²T=0 (Laplace limit)\n');

%% ── Step 8: Decay rate depends on wave number ────────────────
%  For T = e^{−κα²t}sin(αx), the decay rate is κα²
%  Doubling α quadruples decay rate (shorter wavelengths decay faster)
syms alpha1 alpha2 real positive
rate1 = kappa * alpha1^2;
rate2 = kappa * alpha2^2;
%  When α₂ = 2α₁: rate₂/rate₁ = 4
ratio = simplify(subs(rate2, alpha2, 2*alpha1) / rate1);
res8 = simplify(ratio - 4);
assert(isAlways(res8 == 0, 'Unknown', 'false'), ...
    'FAIL Step 8: doubling wave number does not quadruple decay rate');
fprintf('Step 8  PASS — Decay rate ∝ α²: doubling α quadruples decay\n');

%% ── Step 9: Numerical — iron bar diffusion ───────────────────
%  Iron: k = 80 W/(m·K), ρ = 7874 kg/m³, c = 449 J/(kg·K)
%  κ = 80/(7874×449) = 2.263×10⁻⁵ m²/s
%  Mode α = π/L with L = 0.1 m: α = 31.416 1/m
%  Decay time τ = 1/(κα²) ≈ 44.8 s
k_num = 80.0;
rho_num = 7874.0;
c_num = 449.0;
kappa_num = k_num / (rho_num * c_num);
L_num = 0.1;
alpha_num = pi / L_num;
tau_num = 1.0 / (kappa_num * alpha_num^2);
kappa_expected = 2.263e-5;
tau_expected = 44.8;
assert(abs(kappa_num - kappa_expected)/kappa_expected < 0.01, ...
    'FAIL Step 9a: κ for iron mismatch');
assert(abs(tau_num - tau_expected)/tau_expected < 0.01, ...
    'FAIL Step 9b: decay time mismatch');
fprintf('Step 9  PASS — Iron: κ=%.3e m²/s, τ=%.1f s (L=0.1m fundamental mode)\n', kappa_num, tau_num);

%% ── Step 10: Unit consistency ────────────────────────────────
%  [κ] = [k/(ρc)] = W/(m·K) / (kg/m³ · J/(kg·K))
%       = W/(m·K) · m³·K / J = W·m² / J = (J/s)·m² / J = m²/s
%  LHS: [∂T/∂t] = K/s
%  RHS: [κ∇²T] = (m²/s)·(K/m²) = K/s  ✓
%  Verify via dimensional proxy: κ = k/(ρc) with k = ρ·c·D for D [m²/s]
syms D_dim real positive
k_proxy = rho * c_sp * D_dim;
kappa_proxy = k_proxy / (rho * c_sp);
res10 = simplify(kappa_proxy - D_dim);
assert(isAlways(res10 == 0, 'Unknown', 'false'), ...
    'FAIL Step 10: unit consistency proxy');
fprintf('Step 10 PASS — [κ] = [m²/s], [∂T/∂t] = [κ∇²T] = [K/s]\n');

%% ── Step 11: Self-test — wrong sign in Fourier's law ─────────
%  Wrong: J_q = +k∇T (positive instead of negative)
%  Then: ∂u/∂t + ∇·(k∇T) = 0 → ρc∂T/∂t + k∂²T/∂x² = 0
%  → ∂T/∂t = −κ∇²T (anti-diffusion, physically unstable)
%  Our sinusoidal mode should NOT satisfy anti-diffusion
wrong_antidiff = simplify(dT_dt + kappa * d2T_dx2);
%  dT_dt = −κα²·T, d2T_dx2 = −α²·T
%  wrong = −κα²T + κ(−α²T) = −2κα²T ≠ 0
assert(~isAlways(wrong_antidiff == 0, 'Unknown', 'false'), ...
    'FAIL Step 11: anti-diffusion not detected');
fprintf('Step 11 PASS — Wrong Fourier sign → anti-diffusion detected\n');

%% ── Summary ─────────────────────────────────────────────────
fprintf('\n✓ ALL 11 ASSERTIONS PASSED — F0030 audit complete.\n');
fprintf('  Chain: ∂u/∂t+∇·J_q=0, J_q=−k∇T, u=ρcT → ∂T/∂t=κ∇²T\n');
fprintf('  Cross-ref: F0029 (thermodynamic potentials), F0015 (energy conservation)\n');