%% CAS_F0031_VERIFY.m — Radiative flux law j = σT⁴ (PROVISIONAL)
%  Pillar : Thermodynamics (provisional — constitutive law, not derived)
%  CalRef : Formula Atlas F0031, EM Laws L6–L9, Thermo radiation section
%  Status : Provisional — carried as calibration hook; internal derivation
%           chain not yet closed in source atlas.
%  Hardening: isAlways(..., 'Unknown', 'false') on every symbolic assertion
%            symvar-based presence checks; solver and has-function avoided

clear; clc;
fprintf('\n=== CAS F0031 : Radiative flux law j = σT⁴ (provisional) ===\n\n');

%% ── Symbols ──────────────────────────────────────────────────
syms T T1 T2 real positive          % absolute temperature
syms sigma_sb real positive          % Stefan–Boltzmann constant
syms R_sph real positive             % sphere radius

%% ── A. Inputs ────────────────────────────────────────────────
%  Constitutive law (provisional):  j = σ T⁴
%  σ is an effective radiative constant fixed by external calibration
%  This block verifies the algebraic and numerical structure of the
%  T⁴ law; it does NOT attempt spectral or geometric derivation.

%% ── Step 1: Constitutive relation — symbolic encoding ────────
%  j = σT⁴  →  j − σT⁴ = 0
j_flux = sigma_sb * T^4;
res1 = simplify(j_flux - sigma_sb * T^4);
assert(isAlways(res1 == 0, 'Unknown', 'false'), ...
    'FAIL Step 1: j ≠ σT⁴');
fprintf('Step 1  PASS — j = σT⁴ constitutive law encoded\n');

%% ── Step 2: T⁴ scaling — doubling temperature ───────────────
%  j(2T)/j(T) = (2T)⁴/T⁴ = 16
j_at_T = sigma_sb * T^4;
j_at_2T = sigma_sb * (2*T)^4;
ratio_2T = simplify(j_at_2T / j_at_T);
res2 = simplify(ratio_2T - 16);
assert(isAlways(res2 == 0, 'Unknown', 'false'), ...
    'FAIL Step 2: doubling T does not yield 16× flux');
fprintf('Step 2  PASS — j(2T)/j(T) = 16 (T⁴ scaling verified)\n');

%% ── Step 3: General flux ratio ───────────────────────────────
%  j₁/j₂ = (T₁/T₂)⁴
j1 = sigma_sb * T1^4;
j2 = sigma_sb * T2^4;
flux_ratio = simplify(j1 / j2);
expected_ratio = (T1/T2)^4;
res3 = simplify(flux_ratio - expected_ratio);
assert(isAlways(res3 == 0, 'Unknown', 'false'), ...
    'FAIL Step 3: flux ratio ≠ (T₁/T₂)⁴');
fprintf('Step 3  PASS — j₁/j₂ = (T₁/T₂)⁴ (σ cancels in ratio)\n');

%% ── Step 4: Temperature sensitivity dj/dT = 4σT³ ────────────
dj_dT = diff(j_flux, T);
expected_deriv = 4 * sigma_sb * T^3;
res4 = simplify(dj_dT - expected_deriv);
assert(isAlways(res4 == 0, 'Unknown', 'false'), ...
    'FAIL Step 4: dj/dT ≠ 4σT³');
fprintf('Step 4  PASS — dj/dT = 4σT³ (radiative sensitivity)\n');

%% ── Step 5: Total sphere power P = 4πR²σT⁴ ─────────────────
%  Isotropic emitter: P = A·j = 4πR²·σT⁴
P_sphere = 4 * sym(pi) * R_sph^2 * sigma_sb * T^4;
%  Verify structure: P contains σ, T, R
assert(isAlways(P_sphere == 4*sym(pi)*R_sph^2*sigma_sb*T^4, 'Unknown', 'false'), ...
    'FAIL Step 5: sphere power structure');
%  Verify P scales as R²T⁴:
P_double_R = subs(P_sphere, R_sph, 2*R_sph);
ratio_R = simplify(P_double_R / P_sphere);
res5 = simplify(ratio_R - 4);
assert(isAlways(res5 == 0, 'Unknown', 'false'), ...
    'FAIL Step 5b: doubling R does not quadruple power');
fprintf('Step 5  PASS — P = 4πR²σT⁴; doubling R quadruples total power\n');

%% ── Step 6: Numerical — Stefan–Boltzmann constant ────────────
%  σ = 5.670374419 × 10⁻⁸ W/(m²·K⁴)  (CODATA 2018 exact)
sigma_num = 5.670374419e-8;
%  Sun surface: T_sun = 5778 K → j ≈ 6.32 × 10⁷ W/m²
T_sun = 5778.0;
j_sun = sigma_num * T_sun^4;
j_sun_expected = 6.32e7;
assert(abs(j_sun - j_sun_expected)/j_sun_expected < 0.01, ...
    'FAIL Step 6: solar surface flux mismatch');
fprintf('Step 6  PASS — Sun (T=5778K): j = %.3e W/m² ≈ 6.32×10⁷\n', j_sun);

%% ── Step 7: Numerical — room temperature ────────────────────
%  T = 300 K → j ≈ 459.3 W/m²
T_room = 300.0;
j_room = sigma_num * T_room^4;
j_room_expected = 459.3;
assert(abs(j_room - j_room_expected)/j_room_expected < 0.01, ...
    'FAIL Step 7: room temperature flux mismatch');
fprintf('Step 7  PASS — Room (T=300K): j = %.1f W/m² ≈ 459 W/m²\n', j_room);

%% ── Step 8: Numerical — solar luminosity cross-check ─────────
%  L_sun = 4πR²_sun · σ T⁴_sun
%  R_sun = 6.957 × 10⁸ m → L_sun ≈ 3.846 × 10²⁶ W
R_sun = 6.957e8;
L_sun = 4 * pi * R_sun^2 * j_sun;
L_sun_expected = 3.846e26;
assert(abs(L_sun - L_sun_expected)/L_sun_expected < 0.02, ...
    'FAIL Step 8: solar luminosity mismatch');
fprintf('Step 8  PASS — L_sun = %.3e W ≈ 3.846×10²⁶ W\n', L_sun);

%% ── Step 9: Dimensional proxy ────────────────────────────────
%  [j] = W/m² = [σ]·[T]⁴ = [σ]·K⁴
%  → [σ] = W/(m²·K⁴)
%  Proxy: define σ = C·D where C [W/m²] and D [1/K⁴]
%  then j = C·D·T⁴ should have [W/m²] when T⁴ cancels D
syms C_dim D_dim real positive
sigma_proxy = C_dim * D_dim;  % [W/m²] · [1/K⁴]
j_proxy = sigma_proxy * T^4;
%  j_proxy / C_dim should be dimensionless (D_dim·T⁴ → K⁴/K⁴)
dim_check = simplify(j_proxy / C_dim - D_dim * T^4);
assert(isAlways(dim_check == 0, 'Unknown', 'false'), ...
    'FAIL Step 9: dimensional proxy check');
fprintf('Step 9  PASS — [σ] = [W/(m²·K⁴)] dimensional proxy consistent\n');

%% ── Step 10: Self-test — wrong exponent T³ ───────────────────
%  Wrong law: j = σT³ → j(2T)/j(T) = 8, not 16
j_wrong_exp = sigma_sb * T^3;
j_wrong_2T = sigma_sb * (2*T)^3;
wrong_ratio = simplify(j_wrong_2T / j_wrong_exp);
res10 = simplify(wrong_ratio - 16);
assert(~isAlways(res10 == 0, 'Unknown', 'false'), ...
    'FAIL Step 10: wrong exponent T³ not detected');
fprintf('Step 10 PASS — Wrong exponent (T³ → ratio 8≠16) detected\n');

%% ── Step 11: Self-test — wrong sign (negative flux) ─────────
%  Wrong law: j = −σT⁴ → negative flux (physically impossible for emission)
j_wrong_sign = -sigma_sb * T^4;
%  Sum with correct: j_correct + j_wrong should not vanish for arbitrary T
%  Actually test: wrong sign gives j_wrong < 0, which contradicts j > 0
%  Symbolically: j_wrong + j_correct = 0, so flux would vanish — wrong
sum_check = simplify(j_flux + j_wrong_sign);
assert(isAlways(sum_check == 0, 'Unknown', 'false'), ...
    'FAIL Step 11 internal: sum should be zero');
%  The real test: wrong sign means j_wrong ≠ j_correct
diff_check = simplify(j_wrong_sign - j_flux);
assert(~isAlways(diff_check == 0, 'Unknown', 'false'), ...
    'FAIL Step 11: wrong sign not detected');
fprintf('Step 11 PASS — Wrong sign (j = −σT⁴ ≠ +σT⁴) detected\n');

%% ── Summary ─────────────────────────────────────────────────
fprintf('\n✓ ALL 11 ASSERTIONS PASSED — F0031 audit complete.\n');
fprintf('  STATUS: PROVISIONAL — constitutive law, derivation chain not closed.\n');
fprintf('  Law: j = σT⁴ with σ = 5.670374419×10⁻⁸ W/(m²·K⁴)\n');
fprintf('  Cross-ref: F0030 (heat diffusion), F0029 (thermodynamic potentials)\n');
fprintf('  Note: downstream use must acknowledge provisional status.\n');
