%% CAS_F0013_VERIFY.m -- Wave-particle boundary continuity / IGW mode quantization
%  Assertion-based CAS audit block
%  Pillar: Particle Mechanics | Chain: Bohr-Sommerfeld -> de Broglie -> mode quantization
%  CalRef: Mathematical Bridge S5.2, Particle Mechanics preliminaries a71
%
%  Structure mirrors cas_F12.txt (= F0013) sections A-E.
%  Verifies:
%    1. Constant-p loop integral: pL = 2*pi*hbar*n
%    2. hbar substitution: hbar = h/(2*pi)
%    3. Momentum to wave number via de Broglie: k_n = 2*pi*n/L
%    4. k_n*L = 2*pi*n identity
%    5. Interface phase quantization: k_n = (2*pi*n - phi_int)/L
%    6. phi_int = 0 recovery: k_n*L = 2*pi*n
%    7. Mode spacing: Delta_k = 2*pi/L (independent of n)
%    8. Concrete numerical test
%    9. Self-test: wrong quantization (pi*n vs 2*pi*n) detected
%   10. Self-test: wrong residual quantified
%
%  HARDENING: isAlways(..., 'Unknown', 'false') throughout.

clear; clc;
fprintf('=== CAS AUDIT: F0013 -- IGW mode quantization ===\n\n');

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

%% ---- A. INPUTS ----
syms h_planck positive      % Planck constant
syms L_loop positive        % loop length
syms n_mode                 % mode number (integer, but symbolic)
assume(n_mode, 'integer')
assume(n_mode > 0)
syms phi_int real           % interface phase offset
syms p_mom positive         % momentum (on loop branch)

hbar = h_planck / (sym(2)*sym(pi));

fprintf('Section A: Inputs defined.\n');
fprintf('  Bohr-Sommerfeld: loop integral p dl = 2*pi*hbar*n\n');
fprintf('  de Broglie: lambda = h/p, k = 2*pi/lambda\n');
fprintf('  Uniform loop: p constant, length L\n');
fprintf('  Interface phase: phi_int (real parameter)\n\n');

%% ---- B. ASSUMPTIONS / DOMAINS ----
fprintf('Section B: Periodic boundary, p constant on loop, L > 0, phi_int real.\n\n');

%% ---- C. ALLOWED LEMMAS ----
fprintf('Section C: Lemmas declared.\n');
fprintf('  C.1: Constant p => loop integral p dl = p*L\n');
fprintf('  C.2: phi_geom = k*L\n');
fprintf('  C.3: Single-valuedness: k*L + phi_int = 2*pi*n\n\n');

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

% --- Step 1: Bohr-Sommerfeld with constant p ---
% loop integral p dl = 2*pi*hbar*n  and  loop integral p dl = p*L
% Therefore: p*L = 2*pi*hbar*n
% Solve for p: p = 2*pi*hbar*n / L
p_from_BS = 2*sym(pi)*hbar*n_mode / L_loop;

% Verify p*L = 2*pi*hbar*n
step1_residual = simplify(p_from_BS * L_loop - 2*sym(pi)*hbar*n_mode);

total_steps = total_steps + 1;
if isAlways(step1_residual == 0, 'Unknown', 'false')
    fprintf('  Step 1  PASS  p*L = 2*pi*hbar*n (Bohr-Sommerfeld + constant p)\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 1  FAIL  residual: %s\n', char(step1_residual));
    fail_count = fail_count + 1;
end

% --- Step 2: hbar substitution ---
% p = 2*pi*(h/(2*pi))*n/L = h*n/L
p_substituted = simplify(p_from_BS);
p_expected = h_planck * n_mode / L_loop;

step2_residual = simplify(p_substituted - p_expected);

total_steps = total_steps + 1;
if isAlways(step2_residual == 0, 'Unknown', 'false')
    fprintf('  Step 2  PASS  p = h*n/L (hbar = h/(2*pi) substituted)\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 2  FAIL  residual: %s\n', char(step2_residual));
    fail_count = fail_count + 1;
end

% --- Step 3: de Broglie to wave number ---
% k = 2*pi*p/h = 2*pi*(h*n/L)/h = 2*pi*n/L
k_n = 2*sym(pi)*p_expected / h_planck;
k_n_simplified = simplify(k_n);
k_n_expected = 2*sym(pi)*n_mode / L_loop;

step3_residual = simplify(k_n_simplified - k_n_expected);

total_steps = total_steps + 1;
if isAlways(step3_residual == 0, 'Unknown', 'false')
    fprintf('  Step 3  PASS  k_n = 2*pi*n/L (de Broglie)\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 3  FAIL  residual: %s\n', char(step3_residual));
    fail_count = fail_count + 1;
end

% --- Step 4: k_n * L = 2*pi*n identity ---
knL = simplify(k_n_expected * L_loop);
step4_residual = simplify(knL - 2*sym(pi)*n_mode);

total_steps = total_steps + 1;
if isAlways(step4_residual == 0, 'Unknown', 'false')
    fprintf('  Step 4  PASS  k_n*L = 2*pi*n\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 4  FAIL  residual: %s\n', char(step4_residual));
    fail_count = fail_count + 1;
end

% --- Step 5: Interface phase quantization ---
% k_n*L + phi_int = 2*pi*n
% => k_n = (2*pi*n - phi_int)/L
k_n_with_phase = (2*sym(pi)*n_mode - phi_int) / L_loop;

% Verify: k_n_with_phase * L + phi_int = 2*pi*n
step5_residual = simplify(k_n_with_phase * L_loop + phi_int - 2*sym(pi)*n_mode);

total_steps = total_steps + 1;
if isAlways(step5_residual == 0, 'Unknown', 'false')
    fprintf('  Step 5  PASS  k_n = (2*pi*n - phi_int)/L => k_n*L + phi_int = 2*pi*n\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 5  FAIL  residual: %s\n', char(step5_residual));
    fail_count = fail_count + 1;
end

% --- Step 6: phi_int = 0 recovery ---
% When phi_int = 0: k_n = 2*pi*n/L (matches Step 4)
k_n_no_phase = subs(k_n_with_phase, phi_int, 0);
step6_residual = simplify(k_n_no_phase - k_n_expected);

total_steps = total_steps + 1;
if isAlways(step6_residual == 0, 'Unknown', 'false')
    fprintf('  Step 6  PASS  phi_int = 0 recovers k_n = 2*pi*n/L\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 6  FAIL  residual: %s\n', char(step6_residual));
    fail_count = fail_count + 1;
end

% --- Step 7: Mode spacing ---
% Delta_k = k_{n+1} - k_n = 2*pi/L (independent of n and phi_int)
syms n_plus integer
assume(n_plus > 0)
k_n_plus1 = (2*sym(pi)*(n_mode + 1) - phi_int) / L_loop;
Delta_k = simplify(k_n_plus1 - k_n_with_phase);
Delta_k_expected = 2*sym(pi) / L_loop;

step7_residual = simplify(Delta_k - Delta_k_expected);

total_steps = total_steps + 1;
if isAlways(step7_residual == 0, 'Unknown', 'false')
    fprintf('  Step 7  PASS  Delta_k = 2*pi/L (mode spacing, n-independent)\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 7  FAIL  residual: %s\n', char(step7_residual));
    fail_count = fail_count + 1;
end

% --- Step 8: Concrete numerical test ---
% L = 1 nm = 1e-9 m, n = 3, phi_int = pi/4
% k_3 = (2*pi*3 - pi/4)/(1e-9) = (6*pi - pi/4)/(1e-9)
%      = (23*pi/4)/(1e-9) ~ 1.8064e10 /m
L_val = 1e-9;
n_val = 3;
phi_val = pi / 4;

k_computed = (2*pi*n_val - phi_val) / L_val;
k_expected_num = (23*pi/4) / L_val;  % 6*pi - pi/4 = 24*pi/4 - pi/4 = 23*pi/4

rel_error = abs(k_computed - k_expected_num) / abs(k_expected_num);

total_steps = total_steps + 1;
if rel_error < 1e-12
    fprintf('  Step 8  PASS  Numerical: k_3 = %.6e /m (L=1nm, phi=pi/4)\n', k_computed);
    pass_count = pass_count + 1;
else
    fprintf('  Step 8  FAIL  Numerical rel error: %.2e\n', rel_error);
    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('    k: [1/m], L: [m] => k*L: dimensionless\n');
fprintf('    phi_int: dimensionless (phase)\n');
fprintf('    2*pi*n: dimensionless\n');
fprintf('    PASS\n\n');

% --- Self-test: wrong quantization (pi*n instead of 2*pi*n) ---
% Standing-wave condition k_n*L = pi*n is for fixed-fixed, not periodic loop.
% The residual should be nonzero.
k_n_wrong = sym(pi)*n_mode / L_loop;
wrong_phase_residual = simplify(k_n_wrong * L_loop + phi_int - 2*sym(pi)*n_mode);
% = pi*n + phi_int - 2*pi*n = -pi*n + phi_int

total_steps = total_steps + 1;
if ~isAlways(wrong_phase_residual == 0, 'Unknown', 'false')
    fprintf('  Self-test 1: Wrong quantization (pi*n vs 2*pi*n) detected  PASS\n');
    pass_count = pass_count + 1;
else
    fprintf('  Self-test 1: FAIL (wrong quantization not detected)\n');
    fail_count = fail_count + 1;
end

% --- Self-test: quantify wrong residual ---
% wrong_phase_residual should be -pi*n + phi_int
expected_wrong_residual = -sym(pi)*n_mode + phi_int;
wrong_quant = simplify(wrong_phase_residual - expected_wrong_residual);

total_steps = total_steps + 1;
if isAlways(wrong_quant == 0, 'Unknown', 'false')
    fprintf('  Self-test 2: wrong residual = -pi*n + phi_int (quantified)  PASS\n');
    pass_count = pass_count + 1;
else
    fprintf('  Self-test 2: FAIL  residual = %s\n', char(wrong_quant));
    fail_count = fail_count + 1;
end

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

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