%% 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');