%% CAS_F0009_VERIFY.m -- Fundamental thermodynamic relation
%  Assertion-based CAS audit block
%  Pillar: Thermodynamics | Chain: U(S,V,N) -> total differential -> dU=TdS-pdV+mudN
%  CalRef: Thermodynamics Math Appendix S1
%
%  Structure mirrors cas_F08.txt (= F0009) sections A-E.
%  Verifies:
%    1. Total differential of U(S,V,N) yields 3-term form
%    2. Thermodynamic identifications T, -p, mu from partial derivatives
%    3. Substitution produces dU = TdS - pdV + mudN
%    4. Sign structure (especially the -p term)
%    5. Dimensional consistency
%    6. Maxwell relations (cross-derivative consistency)
%
%  HARDENING: isAlways(..., 'Unknown', 'false') throughout.

clear; clc;
fprintf('=== CAS AUDIT: F0009 -- Fundamental thermodynamic relation ===\n\n');

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

%% ---- A. INPUTS ----
% U = U(S,V,N) -- internal energy as state function
% dU = (dU/dS)_{V,N} dS + (dU/dV)_{S,N} dV + (dU/dN)_{S,V} dN
% T := (dU/dS)_{V,N}
% -p := (dU/dV)_{S,N}  =>  (dU/dV) = -p
% mu := (dU/dN)_{S,V}

syms S V_vol N real        % extensive variables (V_vol to avoid clash with sym V)
assume(S > 0); assume(V_vol > 0); assume(N > 0);

syms T_therm positive      % temperature (T_therm to avoid clash)
syms p_therm positive      % pressure
syms mu_therm real          % chemical potential

% U as symbolic function of (S, V, N)
syms U(S, V_vol, N)

fprintf('Section A: Inputs defined.\n');
fprintf('  U(S,V,N), T, p, mu as thermodynamic variables\n\n');

%% ---- B. ASSUMPTIONS / DOMAINS ----
fprintf('Section B: U is C^1, homogeneous equilibrium system, single species.\n\n');

%% ---- C. ALLOWED LEMMAS ----
fprintf('Section C: Lemmas declared.\n');
fprintf('  C.1: Total differential: df = f_x dx + f_y dy + f_z dz\n');
fprintf('  C.2: Thermodynamic identifications T, -p, mu\n\n');

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

% --- Step 1: Total differential of U ---
% dU = U_S dS + U_V dV + U_N dN
% where U_S = diff(U, S), U_V = diff(U, V_vol), U_N = diff(U, N)
U_S = diff(U, S);
U_V = diff(U, V_vol);
U_N = diff(U, N);

% The total differential is the 1-form:
% dU = U_S * dS + U_V * dV + U_N * dN
% We represent coefficients symbolically.

% Verify: U has exactly 3 independent partial derivatives w.r.t. (S, V, N)
% (This is structural -- U is defined as U(S, V_vol, N))
syms dS_form dV_form dN_form real  % differential forms (symbolic placeholders)

dU_total = U_S * dS_form + U_V * dV_form + U_N * dN_form;

% Check: dU_total contains all three differential forms
vars_in_dU = symvar(dU_total);
has_dS = any(vars_in_dU == dS_form);
has_dV = any(vars_in_dU == dV_form);
has_dN = any(vars_in_dU == dN_form);

total_steps = total_steps + 1;
if has_dS && has_dV && has_dN
    fprintf('  Step 1  PASS  Total differential has 3 terms (dS, dV, dN)\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 1  FAIL  Missing differential forms\n');
    fail_count = fail_count + 1;
end

% --- Step 2: Thermodynamic identifications ---
% T = U_S, -p = U_V (so U_V = -p), mu = U_N
% After substitution:
% dU = T*dS + (-p)*dV + mu*dN = T*dS - p*dV + mu*dN

dU_thermo = T_therm * dS_form + (-p_therm) * dV_form + mu_therm * dN_form;

% Verify: the substitution U_S -> T, U_V -> -p, U_N -> mu
% gives the correct form
dU_substituted = subs(dU_total, [U_S, U_V, U_N], [T_therm, -p_therm, mu_therm]);
step2_residual = simplify(dU_substituted - dU_thermo);

total_steps = total_steps + 1;
if isAlways(step2_residual == 0, 'Unknown', 'false')
    fprintf('  Step 2  PASS  Substitution: dU = T*dS - p*dV + mu*dN\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 2  FAIL  Substitution residual: %s\n', char(step2_residual));
    fail_count = fail_count + 1;
end

% --- Step 3: Verify sign of pressure term ---
% The coefficient of dV in the fundamental relation is -p (not +p)
% Extract coefficient of dV_form from dU_thermo
coeff_dV = simplify(diff(dU_thermo, dV_form));
expected_coeff_dV = -p_therm;

step3_residual = simplify(coeff_dV - expected_coeff_dV);

total_steps = total_steps + 1;
if isAlways(step3_residual == 0, 'Unknown', 'false')
    fprintf('  Step 3  PASS  Coefficient of dV = -p (sign correct)\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 3  FAIL  dV coefficient: %s (expected -p)\n', char(coeff_dV));
    fail_count = fail_count + 1;
end

% --- Step 4: Verify coefficient of dS ---
coeff_dS = simplify(diff(dU_thermo, dS_form));
expected_coeff_dS = T_therm;

step4_residual = simplify(coeff_dS - expected_coeff_dS);

total_steps = total_steps + 1;
if isAlways(step4_residual == 0, 'Unknown', 'false')
    fprintf('  Step 4  PASS  Coefficient of dS = T\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 4  FAIL  dS coefficient: %s (expected T)\n', char(coeff_dS));
    fail_count = fail_count + 1;
end

% --- Step 5: Verify coefficient of dN ---
coeff_dN = simplify(diff(dU_thermo, dN_form));
expected_coeff_dN = mu_therm;

step5_residual = simplify(coeff_dN - expected_coeff_dN);

total_steps = total_steps + 1;
if isAlways(step5_residual == 0, 'Unknown', 'false')
    fprintf('  Step 5  PASS  Coefficient of dN = mu\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 5  FAIL  dN coefficient: %s (expected mu)\n', char(coeff_dN));
    fail_count = fail_count + 1;
end

% --- Step 6: Maxwell relation from cross-derivatives ---
% Since U is C^2, mixed partials commute:
%   d^2U/(dS dV) = d^2U/(dV dS)
%   => (dT/dV)_S = -(dp/dS)_V
%
% Verify: diff(T, V) = diff(-p, S) in structure
% Using U: diff(U_S, V) = diff(U_V, S)
% i.e., diff(diff(U, S), V) = diff(diff(U, V), S)
cross_SV = diff(U, S, V_vol);
cross_VS = diff(U, V_vol, S);

step6_residual = simplify(cross_SV - cross_VS);

total_steps = total_steps + 1;
if isAlways(step6_residual == 0, 'Unknown', 'false')
    fprintf('  Step 6  PASS  Maxwell: d^2U/dSdV = d^2U/dVdS (cross-derivatives commute)\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 6  FAIL  Cross-derivative residual: %s\n', char(step6_residual));
    fail_count = fail_count + 1;
end

% --- Step 7: Second Maxwell relation (S,N) ---
cross_SN = diff(U, S, N);
cross_NS = diff(U, N, S);

step7_residual = simplify(cross_SN - cross_NS);

total_steps = total_steps + 1;
if isAlways(step7_residual == 0, 'Unknown', 'false')
    fprintf('  Step 7  PASS  Maxwell: d^2U/dSdN = d^2U/dNdS\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 7  FAIL  (S,N) cross-derivative residual: %s\n', char(step7_residual));
    fail_count = fail_count + 1;
end

% --- Step 8: Third Maxwell relation (V,N) ---
cross_VN = diff(U, V_vol, N);
cross_NV = diff(U, N, V_vol);

step8_residual = simplify(cross_VN - cross_NV);

total_steps = total_steps + 1;
if isAlways(step8_residual == 0, 'Unknown', 'false')
    fprintf('  Step 8  PASS  Maxwell: d^2U/dVdN = d^2U/dNdV\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 8  FAIL  (V,N) cross-derivative residual: %s\n', char(step8_residual));
    fail_count = fail_count + 1;
end

% --- Step 9: Euler relation consistency (extensive scaling) ---
% For a homogeneous first-order function U(S,V,N):
%   U = T*S - p*V + mu*N  (Euler equation)
% Verify: if we define U_euler = T*S - p*V + mu*N,
% then dU_euler = T*dS + S*dT - p*dV - V*dp + mu*dN + N*dmu
% Comparing with dU = T*dS - p*dV + mu*dN:
%   S*dT - V*dp + N*dmu = 0  (Gibbs-Duhem)
%
% Verify the Euler relation structure:
syms S_val V_val N_val positive  % concrete extensive values
U_euler = T_therm * S_val - p_therm * V_val + mu_therm * N_val;

% dU_euler with all 6 differentials
syms dT_form dp_form dmu_form real
dU_euler_full = T_therm*dS_form + S_val*dT_form ...
              - p_therm*dV_form - V_val*dp_form ...
              + mu_therm*dN_form + N_val*dmu_form;

% The fundamental relation part:
dU_fundamental = T_therm*dS_form - p_therm*dV_form + mu_therm*dN_form;

% Gibbs-Duhem: the leftover terms must sum to zero
gibbs_duhem = simplify(dU_euler_full - dU_fundamental);
% = S_val*dT_form - V_val*dp_form + N_val*dmu_form
gibbs_duhem_expected = S_val*dT_form - V_val*dp_form + N_val*dmu_form;

step9_residual = simplify(gibbs_duhem - gibbs_duhem_expected);

total_steps = total_steps + 1;
if isAlways(step9_residual == 0, 'Unknown', 'false')
    fprintf('  Step 9  PASS  Gibbs-Duhem: S*dT - V*dp + N*dmu = 0 (structure verified)\n');
    pass_count = pass_count + 1;
else
    fprintf('  Step 9  FAIL  Gibbs-Duhem residual: %s\n', char(step9_residual));
    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('    T*dS: [K]*[J/K] = [J]\n');
fprintf('    p*dV: [Pa]*[m^3] = [N/m^2]*[m^3] = [J]\n');
fprintf('    mu*dN: [J/mol]*[mol] = [J]\n');
fprintf('    dU: [J] -- all terms consistent\n');
fprintf('    PASS\n\n');

% --- Self-test: wrong sign on pressure term ---
% If we incorrectly used +p instead of -p:
% dU_wrong = T*dS + p*dV + mu*dN
dU_wrong = T_therm * dS_form + p_therm * dV_form + mu_therm * dN_form;
wrong_residual = simplify(dU_wrong - dU_thermo);
% Should be nonzero (= 2*p*dV)

total_steps = total_steps + 1;
if ~isAlways(wrong_residual == 0, 'Unknown', 'false')
    fprintf('  Self-test: wrong pressure sign (+p) detected as different  PASS\n');
    pass_count = pass_count + 1;
else
    fprintf('  Self-test: FAIL (wrong sign not detected!)\n');
    fail_count = fail_count + 1;
end

% --- Self-test: verify wrong_residual = 2*p*dV ---
expected_wrong = 2 * p_therm * dV_form;
wrong_check = simplify(wrong_residual - expected_wrong);

total_steps = total_steps + 1;
if isAlways(wrong_check == 0, 'Unknown', 'false')
    fprintf('  Self-test: wrong - correct = 2*p*dV (quantified)  PASS\n');
    pass_count = pass_count + 1;
else
    fprintf('  Self-test: FAIL (wrong residual = %s, expected 2*p*dV)\n', char(wrong_residual));
    fail_count = fail_count + 1;
end

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

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