%% CAS_F0020_VERIFY.m — Hydrogen spectral fine structure (En + ΔEn)
%  Pillar : Particle Mechanics
%  CalRef : Math Appendix §5D, Calibration §5D
%  Chain  : baseline Rydberg → fine-structure correction → factorised En
%  Hardening: isAlways(..., 'Unknown', 'false') on every symbolic assertion
%            symvar-based presence checks; solver and has-function avoided

clear; clc;
fprintf('\n=== CAS F0020 : Hydrogen fine structure ===\n\n');

%% ── Symbols ──────────────────────────────────────────────────
syms R_inf h_pl c_light n real positive
syms f_val real                     % f(alpha,n,j) — dimensionless correction

%% ── A. Inputs ────────────────────────────────────────────────
%  Nonrelativistic baseline
E0_n = -R_inf * h_pl * c_light / n^2;

%  Fine-structure correction template: ΔEn = E0_n * f(α,n,j)
Delta_En = E0_n * f_val;

%  Total energy
E_n_def = E0_n + Delta_En;                  % from A.2
E_n_fac = E0_n * (1 + f_val);              % factorised form (C.1)

%% ── Step 1: Additive → factorised form ──────────────────────
%  E0 + E0*f = E0*(1+f)
res1 = simplify(E_n_def - E_n_fac);
assert(isAlways(res1 == 0, 'Unknown', 'false'), ...
    'FAIL Step 1: additive form ≠ factorised form');
fprintf('Step 1  PASS — E0 + E0·f = E0·(1+f)\n');

%% ── Step 2: Explicit factorised form with Rydberg ───────────
%  En = -(R∞ h c / n²) · [1 + f]
E_n_explicit = -R_inf * h_pl * c_light / n^2 * (1 + f_val);
res2 = simplify(E_n_fac - E_n_explicit);
assert(isAlways(res2 == 0, 'Unknown', 'false'), ...
    'FAIL Step 2: factorised form ≠ explicit Rydberg form');
fprintf('Step 2  PASS — En = -(R∞hc/n²)·[1+f]\n');

%% ── Step 3: f=0 recovers baseline ──────────────────────────
E_n_f0 = subs(E_n_fac, f_val, 0);
res3 = simplify(E_n_f0 - E0_n);
assert(isAlways(res3 == 0, 'Unknown', 'false'), ...
    'FAIL Step 3: f=0 does not recover baseline');
fprintf('Step 3  PASS — f=0 → En = E0_n (baseline recovery)\n');

%% ── Step 4: ΔEn / E0_n = f (dimensionless ratio) ───────────
ratio = simplify(Delta_En / E0_n);
res4 = simplify(ratio - f_val);
assert(isAlways(res4 == 0, 'Unknown', 'false'), ...
    'FAIL Step 4: ΔEn/E0 ≠ f');
fprintf('Step 4  PASS — ΔEn / E0_n = f (dimensionless ratio)\n');

%% ── Step 5: En / E0_n = 1 + f ──────────────────────────────
ratio_full = simplify(E_n_fac / E0_n);
res5 = simplify(ratio_full - (1 + f_val));
assert(isAlways(res5 == 0, 'Unknown', 'false'), ...
    'FAIL Step 5: En/E0 ≠ 1+f');
fprintf('Step 5  PASS — En / E0_n = 1 + f\n');

%% ── Step 6: Perturbative regime |f| << 1 ───────────────────
%  Check: if |f| << 1, then |En| < |E0_n| (correction is small)
%  Symbolically: En - E0_n = E0_n * f → |ΔEn| < |E0_n| when |f| < 1
delta_check = simplify(E_n_fac - E0_n);
res6 = simplify(delta_check - Delta_En);
assert(isAlways(res6 == 0, 'Unknown', 'false'), ...
    'FAIL Step 6: En - E0_n ≠ ΔEn');
fprintf('Step 6  PASS — En - E0_n = ΔEn (perturbative structure)\n');

%% ── Step 7: Transition energy between levels ────────────────
%  ΔE = En2 - En1 = -R∞hc·[(1+f2)/n2² - (1+f1)/n1²]
syms n1 n2 f1 f2 real positive
E_n1 = -R_inf * h_pl * c_light / n1^2 * (1 + f1);
E_n2 = -R_inf * h_pl * c_light / n2^2 * (1 + f2);
Delta_E_trans = simplify(E_n2 - E_n1);

%  When f1=f2=0, recover Rydberg transition formula
Delta_E_0 = subs(Delta_E_trans, [f1, f2], [0, 0]);
Delta_E_expected = -R_inf * h_pl * c_light * (1/n2^2 - 1/n1^2);
res7 = simplify(Delta_E_0 - Delta_E_expected);
assert(isAlways(res7 == 0, 'Unknown', 'false'), ...
    'FAIL Step 7: f=0 transition ≠ Rydberg formula');
fprintf('Step 7  PASS — f=0 transition → standard Rydberg formula\n');

%% ── Step 8: Numerical check — H Balmer-α (n=3→2) ──────────
%  R∞ = 1.0973731568539e7 m⁻¹, h = 6.62607015e-34 J·s, c = 2.99792458e8 m/s
%  f=0 (baseline): ΔE = R∞hc·(1/4 - 1/9) = R∞hc·5/36
R_val  = 1.0973731568539e7;
h_val  = 6.62607015e-34;
c_val  = 2.99792458e8;
eV_conv = 1.602176634e-19;

DE_Balmer = R_val * h_val * c_val * (1/4 - 1/9);  % positive (absorption)
DE_expected_eV = 1.889;   % ≈ 1.889 eV (656 nm Balmer-α)
DE_calc_eV = DE_Balmer / eV_conv;
assert(abs(DE_calc_eV - DE_expected_eV) < 0.01, ...
    'FAIL Step 8: Balmer-alpha numerical mismatch');
fprintf('Step 8  PASS — Balmer-α: %.4f eV (expected ≈ 1.889 eV)\n', DE_calc_eV);

%% ── Step 9: Fine-structure splitting estimate ───────────────
%  For H, α ≈ 1/137. Leading correction f ~ α²/(n·something)
%  At n=2: splitting ~ α² · R∞hc / 16 ≈ 4.5e-5 eV (0.000045 eV)
alpha_val = 7.2973525693e-3;   % fine-structure constant
%  Rough estimate: f ~ α² ≈ 5.3e-5, so ΔE_fine ~ 13.6 eV / 4 * 5.3e-5 ≈ 1.8e-4 eV
%  Actual Dirac splitting at n=2 (j=1/2 vs j=3/2) ≈ 4.5e-5 eV
f_est = alpha_val^2;   % order-of-magnitude
DE_fine_est = R_val * h_val * c_val / 4 * f_est / eV_conv;
assert(DE_fine_est > 1e-5 && DE_fine_est < 1e-3, ...
    'FAIL Step 9: fine-structure estimate out of range');
fprintf('Step 9  PASS — fine-structure scale: %.2e eV (order α²·E0)\n', DE_fine_est);

%% ── Step 10: Self-test — wrong factorisation ────────────────
%  Wrong: En = E0·(1 - f) instead of (1 + f)
E_wrong = E0_n * (1 - f_val);
res_wrong = simplify(E_wrong - E_n_fac);
%  Should be -2·E0·f ≠ 0 for general f
assert(~isAlways(res_wrong == 0, 'Unknown', 'false'), ...
    'FAIL Step 10a: wrong factorisation not detected');
fprintf('Step 10a PASS — wrong sign (1-f) detected as incorrect\n');

%  Quantify residual: E_wrong - E_n = E0·(1-f) - E0·(1+f) = -2·E0·f
res_quantify = simplify(res_wrong - (-2 * E0_n * f_val));
assert(isAlways(res_quantify == 0, 'Unknown', 'false'), ...
    'FAIL Step 10b: wrong residual ≠ -2·E0·f');
fprintf('Step 10b PASS — wrong residual = -2·E0·f (quantified)\n');

%% ── Summary ─────────────────────────────────────────────────
fprintf('\n✓ ALL 10 ASSERTIONS PASSED — F0020 audit complete.\n');
fprintf('  Chain: Rydberg baseline → fine-structure template → factorised En\n');
fprintf('  Cross-ref: standalone (f imported from EM/Dirac dictionary)\n');