%% CAS_F0008_VERIFY.m -- Hydrogen Balmer H-alpha dictionary hook
%  Assertion-based CAS audit block
%  Pillar: Spectroscopy Anchors | Chain: Rydberg -> Balmer -> H-alpha
%  CalRef: Spectroscopy Anchors section
%
%  Structure mirrors cas_F07.txt (= F0008) sections A-E.
%  Verifies:
%    1. Rydberg wavelength formula derivation (symbolic)
%    2. Balmer H-alpha specialization (n=3->2)
%    3. H-beta (4->2) and H-gamma (5->2)
%    4. Wavelength ratio checks (R_inf cancels)
%    5. Numerical comparison with NIST ASD values
%    6. Rydberg constant from CODATA fundamentals
%
%  HARDENING: isAlways(..., 'Unknown', 'false') throughout.

clear; clc;
fprintf('=== CAS AUDIT: F0008 -- Hydrogen Balmer H-alpha dictionary hook ===\n\n');

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

%% ---- A. INPUTS ----
% R_inf = m_e * e^4 / (8 * eps0^2 * h^3 * c)
% E_n = -h*c*R_inf / n^2
% E_gamma = h*c / lambda
% 1/lambda = R_inf * (1/nf^2 - 1/ni^2)  for ni > nf
% Balmer: nf=2; H-alpha: ni=3

syms R_inf positive       % Rydberg constant [1/m]
syms h_planck positive    % Planck constant
syms c_light positive     % speed of light
syms n_i n_f positive integer  % quantum numbers
syms lambda_sym positive  % wavelength

% Enforce emission ordering: n_i > n_f (required for positive wavelength)
assume(n_i > n_f);

fprintf('Section A: Inputs defined.\n');
fprintf('  R_inf, h, c, quantum numbers n_i, n_f\n\n');

%% ---- B. ASSUMPTIONS / DOMAINS ----
fprintf('Section B: Nonrelativistic H, vacuum, infinite proton mass.\n');
fprintf('  CODATA 2022 values. Jc tolerance <= 1e-6 (ppm).\n\n');

%% ---- C. ALLOWED LEMMAS ----
fprintf('Section C: Lemmas declared.\n');
fprintf('  C.1: E_gamma = hc/lambda = E_nf - E_ni\n');
fprintf('  C.2: Balmer specialization nf=2\n');
fprintf('  C.3: Wavelength ratios: R_inf cancels\n\n');

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

% --- Step 1: Derive Rydberg wavelength formula symbolically ---
% E_n = -h*c*R_inf / n^2
% Delta E = E_nf - E_ni = -h*c*R_inf*(1/nf^2 - 1/ni^2)
% For emission (ni > nf), Delta E < 0 and photon energy = |Delta E|
% E_gamma = h*c/lambda = h*c*R_inf*(1/nf^2 - 1/ni^2)
%   [note: 1/nf^2 - 1/ni^2 > 0 when nf < ni]
% Cancel h*c: 1/lambda = R_inf*(1/nf^2 - 1/ni^2)

E_nf = -h_planck * c_light * R_inf / n_f^2;
E_ni = -h_planck * c_light * R_inf / n_i^2;

delta_E = E_nf - E_ni;
delta_E = simplify(delta_E);

% delta_E = -h*c*R_inf*(1/nf^2 - 1/ni^2)
% Photon: E_gamma = h*c/lambda => lambda = h*c / E_gamma
% Since emission: |delta_E| = h*c*R_inf*(1/nf^2 - 1/ni^2)
% => 1/lambda = R_inf*(1/nf^2 - 1/ni^2)

% Verify: delta_E = -h*c*R_inf*(1/nf^2 - 1/ni^2)
expected_delta_E = -h_planck * c_light * R_inf * (1/n_f^2 - 1/n_i^2);
step1_residual = simplify(delta_E - expected_delta_E);

total_steps = total_steps + 1;
if isAlways(step1_residual == 0, 'Unknown', 'false')
    fprintf('  Step 1  PASS  Delta E = -hcR_inf*(1/nf^2 - 1/ni^2)\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 1  FAIL  Delta E residual: %s\n', char(step1_residual));
    fail_count = fail_count + 1;
end

% --- Step 2: Cancel h*c to get 1/lambda ---
% 1/lambda = R_inf*(1/nf^2 - 1/ni^2)
% Equivalently: lambda = 1 / [R_inf*(1/nf^2 - 1/ni^2)]
inv_lambda = R_inf * (1/n_f^2 - 1/n_i^2);

% Verify: h*c cancels from delta_E / (h*c) = -R_inf*(1/nf^2 - 1/ni^2)
ratio_check = simplify(-delta_E / (h_planck * c_light));
step2_residual = simplify(ratio_check - inv_lambda);

total_steps = total_steps + 1;
if isAlways(step2_residual == 0, 'Unknown', 'false')
    fprintf('  Step 2  PASS  h*c cancels: 1/lambda = R_inf*(1/nf^2 - 1/ni^2)\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 2  FAIL  Cancellation residual: %s\n', char(step2_residual));
    fail_count = fail_count + 1;
end

% --- Step 3: h and c independence ---
% The 1/lambda expression should not contain h or c
inv_lambda_vars = symvar(inv_lambda);
has_h = any(inv_lambda_vars == h_planck);
has_c = any(inv_lambda_vars == c_light);

total_steps = total_steps + 1;
if ~has_h && ~has_c
    fprintf('  Step 3  PASS  1/lambda independent of h and c\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 3  FAIL  1/lambda still contains h=%d, c=%d\n', has_h, has_c);
    fail_count = fail_count + 1;
end

% --- Step 4: Balmer H-alpha specialization (3->2) ---
% 1/lambda_Ha = R_inf*(1/4 - 1/9) = R_inf*(5/36)
% lambda_Ha = 36/(5*R_inf)
inv_lambda_Ha = subs(inv_lambda, [n_f, n_i], [sym(2), sym(3)]);
inv_lambda_Ha = simplify(inv_lambda_Ha);
inv_lambda_Ha_expected = R_inf * sym(5) / sym(36);

step4_residual = simplify(inv_lambda_Ha - inv_lambda_Ha_expected);

total_steps = total_steps + 1;
if isAlways(step4_residual == 0, 'Unknown', 'false')
    fprintf('  Step 4  PASS  H-alpha: 1/lambda = (5/36)*R_inf\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 4  FAIL  H-alpha residual: %s\n', char(step4_residual));
    fail_count = fail_count + 1;
end

% --- Step 5: H-beta specialization (4->2) ---
% 1/lambda_Hb = R_inf*(1/4 - 1/16) = R_inf*(3/16)
inv_lambda_Hb = subs(inv_lambda, [n_f, n_i], [sym(2), sym(4)]);
inv_lambda_Hb = simplify(inv_lambda_Hb);
inv_lambda_Hb_expected = R_inf * sym(3) / sym(16);

step5_residual = simplify(inv_lambda_Hb - inv_lambda_Hb_expected);

total_steps = total_steps + 1;
if isAlways(step5_residual == 0, 'Unknown', 'false')
    fprintf('  Step 5  PASS  H-beta: 1/lambda = (3/16)*R_inf\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 5  FAIL  H-beta residual: %s\n', char(step5_residual));
    fail_count = fail_count + 1;
end

% --- Step 6: H-gamma specialization (5->2) ---
% 1/lambda_Hg = R_inf*(1/4 - 1/25) = R_inf*(21/100)
inv_lambda_Hg = subs(inv_lambda, [n_f, n_i], [sym(2), sym(5)]);
inv_lambda_Hg = simplify(inv_lambda_Hg);
inv_lambda_Hg_expected = R_inf * sym(21) / sym(100);

step6_residual = simplify(inv_lambda_Hg - inv_lambda_Hg_expected);

total_steps = total_steps + 1;
if isAlways(step6_residual == 0, 'Unknown', 'false')
    fprintf('  Step 6  PASS  H-gamma: 1/lambda = (21/100)*R_inf\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 6  FAIL  H-gamma residual: %s\n', char(step6_residual));
    fail_count = fail_count + 1;
end

% --- Step 7: Wavelength ratio H-alpha/H-beta (R_inf cancels) ---
% lambda_Ha / lambda_Hb = (1/inv_lambda_Ha) / (1/inv_lambda_Hb)
%                       = inv_lambda_Hb / inv_lambda_Ha
%                       = (3/16) / (5/36) = (3/16)*(36/5) = 108/80 = 27/20
ratio_Ha_Hb = simplify(inv_lambda_Hb_expected / inv_lambda_Ha_expected);
expected_ratio = sym(27) / sym(20);  % = 1.35

step7_residual = simplify(ratio_Ha_Hb - expected_ratio);

total_steps = total_steps + 1;
if isAlways(step7_residual == 0, 'Unknown', 'false')
    fprintf('  Step 7  PASS  lambda_Ha/lambda_Hb = 27/20 = 1.35 (R_inf cancels)\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 7  FAIL  Ratio residual: %s\n', char(step7_residual));
    fail_count = fail_count + 1;
end

% --- Step 8: R_inf independence of wavelength ratios ---
ratio_symbolic = simplify(inv_lambda_Hb / inv_lambda_Ha);
has_R = any(symvar(ratio_symbolic) == R_inf);

total_steps = total_steps + 1;
if ~has_R
    fprintf('  Step 8  PASS  Wavelength ratio independent of R_inf\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 8  FAIL  Ratio still contains R_inf: %s\n', char(ratio_symbolic));
    fail_count = fail_count + 1;
end

% --- Step 9: Emission positivity check ---
% For emission: 1/nf^2 - 1/ni^2 > 0 requires nf < ni
% Verify with concrete Balmer values: 1/4 - 1/9 = 5/36 > 0
emission_factor_Ha = sym(1)/sym(4) - sym(1)/sym(9);
emission_factor_Hb = sym(1)/sym(4) - sym(1)/sym(16);
emission_factor_Hg = sym(1)/sym(4) - sym(1)/sym(25);

total_steps = total_steps + 1;
if isAlways(emission_factor_Ha > 0, 'Unknown', 'false') && ...
   isAlways(emission_factor_Hb > 0, 'Unknown', 'false') && ...
   isAlways(emission_factor_Hg > 0, 'Unknown', 'false')
    fprintf('  Step 9  PASS  Emission positivity: 1/nf^2 - 1/ni^2 > 0 for Ha,Hb,Hg\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 9  FAIL  Emission factor not positive\n');
    fail_count = fail_count + 1;
end

% --- Step 10: Numerical R_inf from CODATA fundamentals ---
% R_inf = m_e * e^4 / (8 * eps0^2 * h^3 * c)
% CODATA 2022:
m_e_val    = 9.1093837139e-31;   % kg
e_val      = 1.602176634e-19;    % C (exact)
eps0_val   = 8.8541878128e-12;   % F/m
h_val      = 6.62607015e-34;     % J*s (exact)
c_val      = 299792458;          % m/s (exact)

R_inf_computed = m_e_val * e_val^4 / (8 * eps0_val^2 * h_val^3 * c_val);
R_inf_CODATA = 10973731.568157;  % 1/m (CODATA 2022)

R_inf_rel_error = abs(R_inf_computed - R_inf_CODATA) / R_inf_CODATA;

total_steps = total_steps + 1;
if R_inf_rel_error < 1e-6
    fprintf('  Step 10 PASS  R_inf = %.6f 1/m (rel error: %.2e)\n', ...
            R_inf_computed, R_inf_rel_error);
    pass_count = pass_count + 1;
else
    fprintf('  Step 10 FAIL  R_inf = %.6f, expected %.6f (rel error: %.2e)\n', ...
            R_inf_computed, R_inf_CODATA, R_inf_rel_error);
    fail_count = fail_count + 1;
end

% --- Step 11: Numerical H-alpha wavelength ---
lambda_Ha_computed = 36 / (5 * R_inf_CODATA);  % in meters
lambda_Ha_nm = lambda_Ha_computed * 1e9;
lambda_Ha_NIST = 656.281;  % nm (NIST ASD v5.11 2023, vacuum)

lambda_Ha_rel_error = abs(lambda_Ha_nm - lambda_Ha_NIST) / lambda_Ha_NIST;

total_steps = total_steps + 1;
if lambda_Ha_rel_error < 1e-3  % 0.1% tolerance (infinite proton mass, no reduced-mass/QED correction)
    fprintf('  Step 11 PASS  lambda_Ha = %.3f nm (NIST: %.3f nm, rel err: %.2e)\n', ...
            lambda_Ha_nm, lambda_Ha_NIST, lambda_Ha_rel_error);
    pass_count = pass_count + 1;
else
    fprintf('  Step 11 FAIL  lambda_Ha = %.3f nm (rel error: %.2e)\n', ...
            lambda_Ha_nm, lambda_Ha_rel_error);
    fail_count = fail_count + 1;
end

% --- Step 12: Numerical ratio check vs NIST ---
% lambda_Ha/lambda_Hb from NIST values
lambda_Hb_NIST = 486.133;  % nm
lambda_Hg_NIST = 434.047;  % nm

ratio_exp = lambda_Ha_NIST / lambda_Hb_NIST;
ratio_theory = 27/20;  % = 1.35

ratio_ppm = abs(ratio_exp - ratio_theory) / ratio_theory * 1e6;

total_steps = total_steps + 1;
if ratio_ppm < 10  % within 10 ppm
    fprintf('  Step 12 PASS  Ha/Hb ratio: theory=%.6f, NIST=%.6f (%.1f ppm)\n', ...
            ratio_theory, ratio_exp, ratio_ppm);
    pass_count = pass_count + 1;
else
    fprintf('  Step 12 FAIL  Ha/Hb ratio discrepancy: %.1f ppm\n', ratio_ppm);
    fail_count = fail_count + 1;
end

% --- Step 13: H-beta and H-gamma numerical wavelengths ---
lambda_Hb_computed = 16 / (3 * R_inf_CODATA) * 1e9;  % nm
lambda_Hg_computed = 100 / (21 * R_inf_CODATA) * 1e9; % nm

lambda_Hb_rel = abs(lambda_Hb_computed - lambda_Hb_NIST) / lambda_Hb_NIST;
lambda_Hg_rel = abs(lambda_Hg_computed - lambda_Hg_NIST) / lambda_Hg_NIST;

total_steps = total_steps + 1;
if lambda_Hb_rel < 1e-3 && lambda_Hg_rel < 1e-3  % 0.1% (infinite proton mass)
    fprintf('  Step 13 PASS  Hb=%.3f nm (NIST:%.3f, err:%.2e), Hg=%.3f nm (NIST:%.3f, err:%.2e)\n', ...
            lambda_Hb_computed, lambda_Hb_NIST, lambda_Hb_rel, ...
            lambda_Hg_computed, lambda_Hg_NIST, lambda_Hg_rel);
    pass_count = pass_count + 1;
else
    fprintf('  Step 13 FAIL  Hb rel err: %.2e, Hg rel err: %.2e\n', lambda_Hb_rel, lambda_Hg_rel);
    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('    R_inf: [1/m]\n');
fprintf('    1/lambda = R_inf*(..): [1/m] => lambda: [m]\n');
fprintf('    Ratios: [1/m]/[1/m] = [1] (dimensionless)\n');
fprintf('    PASS\n\n');

% --- Self-test: wrong quantum numbers should give wrong wavelength ---
% If we accidentally use n_i=2, n_f=3 (swapped), 1/lambda becomes negative
inv_lambda_wrong = subs(inv_lambda, [n_f, n_i], [sym(3), sym(2)]);
inv_lambda_wrong = simplify(inv_lambda_wrong);
% This should be negative (unphysical for emission)
inv_lambda_wrong_val = double(subs(inv_lambda_wrong, R_inf, R_inf_CODATA));

total_steps = total_steps + 1;
if inv_lambda_wrong_val < 0
    fprintf('  Self-test: swapped n_i/n_f gives 1/lambda < 0 (unphysical)  PASS\n');
    pass_count = pass_count + 1;
else
    fprintf('  Self-test: FAIL (swapped quantum numbers not detected)\n');
    fail_count = fail_count + 1;
end

fprintf('---------------------------------------------\n\n');

%% ---- VERDICT ----
fprintf('=============================================\n');
fprintf('  F0008 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 F0008.\n');