%% CAS_F0017_VERIFY.m -- Heat capacity relation (Cp - Cv = R_eff)
%  Assertion-based CAS audit block
%  Pillar: Thermodynamics | Chain: enthalpy -> ideal gas -> Cp - Cv = R_eff
%  CalRef: Thermodynamics Math Appendix S5Y, Calibration S5Y.2
%
%  Structure mirrors cas_F16.txt (= F0017) sections A-E.
%  Verifies:
%    1. Enthalpy differential: dH = dU + pdV + Vdp
%    2. Substitution: dH = TdS + Vdp (pdV cancels)
%    3. Ideal gas U=U(T): no p,V dependence in H - R_eff*T
%    4. Cp - Cv = R_eff from dH/dT = dU/dT + R_eff
%    5. Reciprocal: Cv = Cp - R_eff
%    6. Heat capacity ratio: gamma = Cp/Cv = 1 + R_eff/Cv
%    7. Concrete: monatomic ideal gas (Cv = 3R/2, Cp = 5R/2)
%    8. Numerical R check
%    9. Self-test: wrong sign (Cp - Cv = -R_eff) detected
%   10. Self-test: wrong residual quantified
%
%  HARDENING: isAlways(..., 'Unknown', 'false') throughout.

clear; clc;
fprintf('=== CAS AUDIT: F0017 -- Heat capacity relation ===\n\n');

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

%% ---- A. INPUTS ----
syms T_sym positive         % temperature
syms p_sym positive         % pressure
syms V_sym positive         % volume
syms R_eff positive         % effective gas constant (per mole)

% U = U(T) for ideal gas: represent as generic function of T
syms U_func(T_sym)          % internal energy as function of T only

% Cv = dU/dT (definition)
Cv = diff(U_func, T_sym);

fprintf('Section A: Inputs defined.\n');
fprintf('  Cv = (dU/dT)_V, Cp = (dH/dT)_P\n');
fprintf('  dU = T*dS - p*dV\n');
fprintf('  H = U + pV, pV = R_eff*T\n\n');

%% ---- B. ASSUMPTIONS / DOMAINS ----
fprintf('Section B: Ideal gas, R_eff constant, U = U(T) only.\n\n');

%% ---- C. ALLOWED LEMMAS ----
fprintf('Section C: Lemmas declared.\n');
fprintf('  C.1: dH = dU + pdV + Vdp\n');
fprintf('  C.2: (dU/dT)_P = (dU/dT)_V = Cv (ideal gas)\n');
fprintf('  C.3: d(pV)/dT at const p = R_eff\n\n');

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

% --- Step 1: Enthalpy differential ---
% H = U + pV
% dH = dU + d(pV) = dU + pdV + Vdp
% Verify: d(pV) = pdV + Vdp (product rule)
% Using symbolic differentials as symbols:
syms dU_sym dS_sym dV_sym dp_sym real
syms T_val positive
syms p_val positive
syms V_val positive

dH_from_product = dU_sym + p_val*dV_sym + V_val*dp_sym;
% This is the product rule applied to d(pV) = p*dV + V*dp

% Step 1 is structural: dH = dU + pdV + Vdp
% Verify by checking that the differential of pV is pdV + Vdp
dpV = p_val*dV_sym + V_val*dp_sym;
% This is correct by product rule. Verify symbolically:
% d(p*V)/dV at const p = p, d(p*V)/dp at const V = V
syms pV_func
pV_func = p_val * V_val;
dpV_dV = diff(pV_func, V_val);  % = p_val
dpV_dp = diff(pV_func, p_val);  % = V_val

step1a_residual = simplify(dpV_dV - p_val);
step1b_residual = simplify(dpV_dp - V_val);

total_steps = total_steps + 1;
if isAlways(step1a_residual == 0, 'Unknown', 'false') && isAlways(step1b_residual == 0, 'Unknown', 'false')
    fprintf('  Step 1  PASS  d(pV) = p*dV + V*dp (product rule)\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 1  FAIL\n');
    fail_count = fail_count + 1;
end

% --- Step 2: dH = TdS + Vdp ---
% dU = TdS - pdV
% dH = dU + pdV + Vdp = TdS - pdV + pdV + Vdp = TdS + Vdp
dU_expanded = T_val*dS_sym - p_val*dV_sym;
dH_expanded = dU_expanded + p_val*dV_sym + V_val*dp_sym;
dH_expected = T_val*dS_sym + V_val*dp_sym;

step2_residual = simplify(dH_expanded - dH_expected);

total_steps = total_steps + 1;
if isAlways(step2_residual == 0, 'Unknown', 'false')
    fprintf('  Step 2  PASS  dH = TdS + Vdp (pdV cancels)\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: U = U(T) only (ideal gas property) ---
% For ideal gas, U depends only on T, not on p or V.
% Verify: U_func(T_sym) has no p_sym or V_sym dependence.
% Since U_func is declared as symfun of T_sym only, dU/dp = 0 and dU/dV = 0.
% This means (dU/dT)_P = (dU/dT)_V = Cv.
% Verify by checking that diff(U_func, p_sym) and diff(U_func, V_sym) are zero
% (U_func does not depend on p_sym or V_sym by construction).
dU_dp = diff(U_func, p_sym);
dU_dV = diff(U_func, V_sym);

total_steps = total_steps + 1;
if isAlways(dU_dp == 0, 'Unknown', 'false') && isAlways(dU_dV == 0, 'Unknown', 'false')
    fprintf('  Step 3  PASS  U = U(T) only: dU/dp = 0, dU/dV = 0 (ideal gas)\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 3  FAIL  U has unexpected p or V dependence\n');
    fail_count = fail_count + 1;
end

% --- Step 4: Cp - Cv = R_eff ---
% H = U(T) + pV = U(T) + R_eff*T
% Cp = dH/dT = dU/dT + R_eff = Cv + R_eff
% Therefore Cp - Cv = R_eff
H_ideal = U_func + R_eff * T_sym;
Cp_derived = diff(H_ideal, T_sym);
Cp_sym = Cv + R_eff;

diff_CpCv = simplify(Cp_derived - Cv);
step4_residual = simplify(diff_CpCv - R_eff);

total_steps = total_steps + 1;
if isAlways(step4_residual == 0, 'Unknown', 'false')
    fprintf('  Step 4  PASS  Cp - Cv = R_eff\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: Reciprocal: Cv = Cp - R_eff ---
Cv_from_Cp = Cp_sym - R_eff;
step6_residual = simplify(Cv_from_Cp - Cv);

total_steps = total_steps + 1;
if isAlways(step6_residual == 0, 'Unknown', 'false')
    fprintf('  Step 5  PASS  Cv = Cp - R_eff (reciprocal)\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 5  FAIL  residual: %s\n', char(step6_residual));
    fail_count = fail_count + 1;
end

% --- Step 6: Heat capacity ratio ---
% gamma = Cp/Cv = (Cv + R_eff)/Cv = 1 + R_eff/Cv
syms Cv_val positive  % concrete Cv value for ratio test
gamma_sym = (Cv_val + R_eff) / Cv_val;
gamma_expected = 1 + R_eff / Cv_val;

step7_residual = simplify(gamma_sym - gamma_expected);

total_steps = total_steps + 1;
if isAlways(step7_residual == 0, 'Unknown', 'false')
    fprintf('  Step 6  PASS  gamma = 1 + R_eff/Cv\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 6  FAIL  residual: %s\n', char(step7_residual));
    fail_count = fail_count + 1;
end

% --- Step 7: Concrete monatomic ideal gas ---
% Cv = (3/2)*R, Cp = (5/2)*R, gamma = 5/3
% Cp - Cv = R (R_eff = R for ideal gas)
syms R_gas positive
Cv_mono = sym(3)/sym(2) * R_gas;
Cp_mono = Cv_mono + R_gas;
Cp_mono_expected = sym(5)/sym(2) * R_gas;

step8a_residual = simplify(Cp_mono - Cp_mono_expected);

gamma_mono = simplify(Cp_mono / Cv_mono);
gamma_mono_expected = sym(5)/sym(3);

step8b_residual = simplify(gamma_mono - gamma_mono_expected);

total_steps = total_steps + 1;
if isAlways(step8a_residual == 0, 'Unknown', 'false') && isAlways(step8b_residual == 0, 'Unknown', 'false')
    fprintf('  Step 7  PASS  Monatomic: Cv=3R/2, Cp=5R/2, gamma=5/3\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 7  FAIL\n');
    fail_count = fail_count + 1;
end

% --- Step 8: Numerical R check ---
% R = 8.31446 J/(mol K)
% Monatomic: Cv = 12.4717, Cp = 20.7862
R_val = 8.31446;
Cv_num = 1.5 * R_val;
Cp_num = Cv_num + R_val;
diff_num = Cp_num - Cv_num;
rel_error = abs(diff_num - R_val) / R_val;

total_steps = total_steps + 1;
if rel_error < 1e-12
    fprintf('  Step 8  PASS  Numerical: Cp - Cv = %.5f = R J/(mol K)\n', diff_num);
    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('    Cp, Cv: [J/(mol*K)]\n');
fprintf('    R_eff: [J/(mol*K)]\n');
fprintf('    Cp - Cv = R_eff: [J/(mol*K)]\n');
fprintf('    PASS\n\n');

% --- Self-test: wrong sign (Cp - Cv = -R_eff) ---
wrong_diff = -R_eff;
wrong_residual = simplify(wrong_diff - R_eff);

total_steps = total_steps + 1;
if ~isAlways(wrong_residual == 0, 'Unknown', 'false')
    fprintf('  Self-test 1: Wrong sign (Cp-Cv = -R_eff) detected  PASS\n');
    pass_count = pass_count + 1;
else
    fprintf('  Self-test 1: FAIL (wrong sign not detected)\n');
    fail_count = fail_count + 1;
end

% --- Self-test: quantify wrong residual ---
% -R_eff - R_eff = -2*R_eff
expected_wrong_res = -2*R_eff;
wrong_quant = simplify(wrong_residual - expected_wrong_res);

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