%% CAS_F0019_VERIFY.m -- Photoelectric relation (K_max = h*nu - phi)
%  Assertion-based CAS audit block
%  Pillar: Electromagnetism | Chain: photon energy -> energy balance -> threshold
%  CalRef: Electromagnetism Math Appendix S5C, EM Calibration S5C
%
%  Structure mirrors cas_F18.txt (= F0019) sections A-E.
%  Verifies:
%    1. Energy balance: h*nu = phi + K => K = h*nu - phi
%    2. Threshold frequency: nu_0 = phi/h
%    3. K_max = h*nu - phi (no-loss model)
%    4. Stopping potential: e*V0 = K_max
%    5. Linear form: K_max vs nu has slope h, intercept -phi
%    6. At threshold: K_max(nu_0) = 0
%    7. Below threshold: K_max < 0 (no emission)
%    8. Concrete numerical: cesium (phi ~ 2.1 eV, UV light)
%    9. Self-test: wrong balance (h*nu = phi - K) detected
%   10. Self-test: wrong residual quantified
%
%  HARDENING: isAlways(..., 'Unknown', 'false') throughout.

clear; clc;
fprintf('=== CAS AUDIT: F0019 -- Photoelectric relation ===\n\n');

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

%% ---- A. INPUTS ----
syms h_planck positive      % Planck constant
syms nu positive            % photon frequency
syms phi_work positive      % work function
syms K_kin real             % kinetic energy (can be negative for sub-threshold test)
syms e_charge positive      % elementary charge

fprintf('Section A: Inputs defined.\n');
fprintf('  E_gamma = h*nu, phi > 0\n');
fprintf('  h*nu = phi + K (eventwise balance, no-loss model)\n\n');

%% ---- B. ASSUMPTIONS / DOMAINS ----
fprintf('Section B: Single-photon single-electron, no additional loss channels.\n\n');

%% ---- C. ALLOWED LEMMAS ----
fprintf('Section C: Lemmas declared.\n');
fprintf('  C.1: K >= 0\n');
fprintf('  C.2: K = h*nu - phi\n\n');

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

% --- Step 1: Energy balance ---
% h*nu = phi + K => K = h*nu - phi
K_from_balance = h_planck*nu - phi_work;

% Verify: K + phi = h*nu
step1_residual = simplify((K_from_balance + phi_work) - h_planck*nu);

total_steps = total_steps + 1;
if isAlways(step1_residual == 0, 'Unknown', 'false')
    fprintf('  Step 1  PASS  h*nu = phi + K => K = h*nu - phi\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: Threshold frequency ---
% K >= 0 => h*nu - phi >= 0 => nu >= phi/h =: nu_0
nu_0 = phi_work / h_planck;

% Verify: h*nu_0 - phi = 0
step2_residual = simplify(h_planck*nu_0 - phi_work);

total_steps = total_steps + 1;
if isAlways(step2_residual == 0, 'Unknown', 'false')
    fprintf('  Step 2  PASS  nu_0 = phi/h (threshold frequency)\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: K_max in no-loss model ---
% K_max = h*nu - phi (equality, not inequality, in ideal model)
K_max = h_planck*nu - phi_work;

% Verify K_max = K_from_balance
step3_residual = simplify(K_max - K_from_balance);

total_steps = total_steps + 1;
if isAlways(step3_residual == 0, 'Unknown', 'false')
    fprintf('  Step 3  PASS  K_max = h*nu - phi (no-loss model)\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: Stopping potential ---
% e*V0 = K_max => V0 = K_max/e = (h*nu - phi)/e
syms V0_stop
V0_expr = K_max / e_charge;

% Verify: e*V0_expr = K_max
step4_residual = simplify(e_charge * V0_expr - K_max);

total_steps = total_steps + 1;
if isAlways(step4_residual == 0, 'Unknown', 'false')
    fprintf('  Step 4  PASS  e*V0 = K_max => V0 = (h*nu - phi)/e\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: Linear form ---
% K_max = h*nu - phi is linear in nu
% Slope = dK_max/d(nu) = h
% Intercept (at nu=0, extrapolated) = -phi
slope = diff(K_max, nu);
intercept = subs(K_max, nu, 0);
% Note: subs(nu, 0) gives -phi_work, but nu is positive.
% This is the extrapolated y-intercept, not a physical point.

step5a_residual = simplify(slope - h_planck);
step5b_residual = simplify(intercept - (-phi_work));

total_steps = total_steps + 1;
if isAlways(step5a_residual == 0, 'Unknown', 'false') && isAlways(step5b_residual == 0, 'Unknown', 'false')
    fprintf('  Step 5  PASS  K_max vs nu: slope = h, intercept = -phi\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 5  FAIL\n');
    fail_count = fail_count + 1;
end

% --- Step 6: At threshold ---
% K_max(nu_0) = h*(phi/h) - phi = phi - phi = 0
K_at_threshold = subs(K_max, nu, nu_0);
step6_residual = simplify(K_at_threshold);

total_steps = total_steps + 1;
if isAlways(step6_residual == 0, 'Unknown', 'false')
    fprintf('  Step 6  PASS  K_max(nu_0) = 0 (threshold)\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: Below threshold ---
% nu < nu_0 => K_max < 0 (no emission physically)
% Test with nu = nu_0/2:
syms nu_half positive
nu_half = nu_0 / 2;
K_below = subs(K_max, nu, nu_half);
K_below_simplified = simplify(K_below);
% K_below = h*(phi/(2h)) - phi = phi/2 - phi = -phi/2 < 0

expected_below = -phi_work / 2;
step7_residual = simplify(K_below_simplified - expected_below);

total_steps = total_steps + 1;
if isAlways(step7_residual == 0, 'Unknown', 'false') && isAlways(K_below_simplified < 0, 'Unknown', 'false')
    fprintf('  Step 7  PASS  Below threshold: K_max(nu_0/2) = -phi/2 < 0\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 7  FAIL\n');
    fail_count = fail_count + 1;
end

% --- Step 8: Concrete numerical ---
% Cesium: phi ~ 2.1 eV = 2.1 * 1.602e-19 J
% UV light: nu = 1.0e15 Hz
% h = 6.626e-34 J*s
% K_max = h*nu - phi
h_val = 6.626e-34;
eV_to_J = 1.602e-19;
phi_val = 2.1 * eV_to_J;  % J
nu_val = 1.0e15;           % Hz

K_max_num = h_val * nu_val - phi_val;
K_max_eV = K_max_num / eV_to_J;

% K_max should be ~ (6.626e-34*1e15)/(1.602e-19) - 2.1 eV
%                 ~ 4.137 - 2.1 = 2.037 eV
K_expected_eV = h_val * nu_val / eV_to_J - 2.1;
rel_error = abs(K_max_eV - K_expected_eV) / abs(K_expected_eV);

total_steps = total_steps + 1;
if rel_error < 1e-10 && K_max_num > 0
    fprintf('  Step 8  PASS  Numerical: K_max(Cs, 1e15 Hz) = %.3f eV\n', K_max_eV);
    pass_count = pass_count + 1;
else
    fprintf('  Step 8  FAIL  rel error: %.2e, K = %.4e J\n', rel_error, K_max_num);
    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('    h*nu: [J*s]*[1/s] = [J]\n');
fprintf('    phi: [J]\n');
fprintf('    K_max: [J]\n');
fprintf('    e*V0: [C]*[V] = [J]\n');
fprintf('    PASS\n\n');

% --- Self-test: wrong balance (h*nu = phi - K) ---
% This gives K_wrong = phi - h*nu (sign flipped)
K_wrong = phi_work - h_planck*nu;
wrong_residual = simplify(K_wrong - K_max);

total_steps = total_steps + 1;
if ~isAlways(wrong_residual == 0, 'Unknown', 'false')
    fprintf('  Self-test 1: Wrong balance (h*nu = phi - K) detected  PASS\n');
    pass_count = pass_count + 1;
else
    fprintf('  Self-test 1: FAIL (wrong balance not detected)\n');
    fail_count = fail_count + 1;
end

% --- Self-test: quantify wrong residual ---
% K_wrong - K_max = (phi - h*nu) - (h*nu - phi) = 2*(phi - h*nu) = -2*K_max
expected_wrong = -2*K_max;
wrong_quant = simplify(wrong_residual - expected_wrong);

total_steps = total_steps + 1;
if isAlways(wrong_quant == 0, 'Unknown', 'false')
    fprintf('  Self-test 2: wrong - correct = -2*K_max = 2*(phi - h*nu) (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('  F0019 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 F0019.\n');