"""Heat capacity relation (Cp - Cv = R_eff).

Assertion-based CAS audit block.
Pillar: Thermodynamics | Chain: enthalpy -> ideal gas -> Cp - Cv = R_eff
"""


def run():
    from sympy import symbols, Function, diff, simplify, Rational

    print("=== CAS AUDIT: F0017 — Heat capacity relation ===\n")

    pass_count = 0
    fail_count = 0
    total_steps = 0

    print("Section A: Inputs defined.")
    print("  Cv = (dU/dT)_V, Cp = (dH/dT)_P\n")

    T_sym = symbols("T_sym", positive=True)
    p_sym = symbols("p_sym", positive=True)
    V_sym = symbols("V_sym", positive=True)
    R_eff = symbols("R_eff", positive=True)

    U_func = Function("U_func")(T_sym)
    Cv = diff(U_func, T_sym)

    print("Section B: Ideal gas, R_eff constant, U = U(T) only.\n")
    print("Section C: Lemmas declared.\n")
    print("Section D: Step log")
    print("---------------------------------------------")

    # Step 1: Enthalpy differential
    dU_sym, dS_sym, dV_sym, dp_sym = symbols("dU_sym dS_sym dV_sym dp_sym", real=True)
    T_val, p_val, V_val = symbols("T_val p_val V_val", positive=True)

    dH_from_product = dU_sym + p_val * dV_sym + V_val * dp_sym
    dpV = p_val * dV_sym + V_val * dp_sym

    pV_func = p_val * V_val
    dpV_dV = diff(pV_func, V_val)
    dpV_dp = diff(pV_func, p_val)

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

    total_steps += 1
    if simplify(step1a_residual) == 0 and simplify(step1b_residual) == 0:
        print("  Step 1  PASS — d(pV) = p*dV + V*dp (product rule)")
        pass_count += 1
    else:
        print("  Step 1  FAIL")
        fail_count += 1

    # Step 2: dH = 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 += 1
    if simplify(step2_residual) == 0:
        print("  Step 2  PASS — dH = TdS + Vdp (pdV cancels)")
        pass_count += 1
    else:
        print(f"  Step 2  FAIL — residual: {step2_residual}")
        fail_count += 1

    # Step 3: U = U(T) only
    dU_dp = diff(U_func, p_sym)
    dU_dV = diff(U_func, V_sym)

    total_steps += 1
    if simplify(dU_dp) == 0 and simplify(dU_dV) == 0:
        print("  Step 3  PASS — U = U(T) only: dU/dp = 0, dU/dV = 0")
        pass_count += 1
    else:
        print("  Step 3  FAIL — U has unexpected p or V dependence")
        fail_count += 1

    # Step 4: 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 += 1
    if simplify(step4_residual) == 0:
        print("  Step 4  PASS — Cp - Cv = R_eff")
        pass_count += 1
    else:
        print(f"  Step 4  FAIL — residual: {step4_residual}")
        fail_count += 1

    # Step 5: Reciprocal
    Cv_from_Cp = Cp_sym - R_eff
    step5_residual = simplify(Cv_from_Cp - Cv)

    total_steps += 1
    if simplify(step5_residual) == 0:
        print("  Step 5  PASS — Cv = Cp - R_eff")
        pass_count += 1
    else:
        print(f"  Step 5  FAIL — residual: {step5_residual}")
        fail_count += 1

    # Step 6: Heat capacity ratio
    Cv_val = symbols("Cv_val", positive=True)
    gamma_sym = (Cv_val + R_eff) / Cv_val
    gamma_expected = 1 + R_eff / Cv_val
    step6_residual = simplify(gamma_sym - gamma_expected)

    total_steps += 1
    if simplify(step6_residual) == 0:
        print("  Step 6  PASS — gamma = 1 + R_eff/Cv")
        pass_count += 1
    else:
        print(f"  Step 6  FAIL — residual: {step6_residual}")
        fail_count += 1

    # Step 7: Monatomic
    R_gas = symbols("R_gas", positive=True)
    Cv_mono = Rational(3, 2) * R_gas
    Cp_mono = Cv_mono + R_gas
    Cp_mono_expected = Rational(5, 2) * R_gas
    step7a_residual = simplify(Cp_mono - Cp_mono_expected)

    gamma_mono = simplify(Cp_mono / Cv_mono)
    gamma_mono_expected = Rational(5, 3)
    step7b_residual = simplify(gamma_mono - gamma_mono_expected)

    total_steps += 1
    if simplify(step7a_residual) == 0 and simplify(step7b_residual) == 0:
        print("  Step 7  PASS — Monatomic: Cv=3R/2, Cp=5R/2, gamma=5/3")
        pass_count += 1
    else:
        print("  Step 7  FAIL")
        fail_count += 1

    # Step 8: Numerical
    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 += 1
    if rel_error < 1e-12:
        print(f"  Step 8  PASS — Numerical: Cp - Cv = {diff_num:.5f} = R J/(mol K)")
        pass_count += 1
    else:
        print(f"  Step 8  FAIL — Numerical rel error: {rel_error:.2e}")
        fail_count += 1

    print("---------------------------------------------\n")
    print("Section E: Output checks")
    print("---------------------------------------------")
    print("  Unit check: [J/(mol*K)] — PASS\n")

    # Self-test: wrong sign
    wrong_diff = -R_eff
    wrong_residual = simplify(wrong_diff - R_eff)

    total_steps += 1
    if simplify(wrong_residual) != 0:
        print("  Self-test: Wrong sign (Cp-Cv = -R_eff) detected  PASS")
        pass_count += 1
    else:
        print("  Self-test: FAIL (wrong sign not detected)")
        fail_count += 1

    expected_wrong_res = -2 * R_eff
    wrong_quant = simplify(wrong_residual - expected_wrong_res)

    total_steps += 1
    if simplify(wrong_quant) == 0:
        print("  Self-test: wrong - correct = -2*R_eff (quantified)  PASS")
        pass_count += 1
    else:
        print(f"  Self-test: FAIL — residual = {wrong_quant}")
        fail_count += 1

    print("---------------------------------------------\n")
    print("=============================================")
    print("  F0017 AUDIT RESULT")
    print(f"  Steps: {total_steps}  |  Pass: {pass_count}  |  Fail: {fail_count}")
    if fail_count == 0:
        print("  STATUS: *** PASS ***")
    else:
        print(f"  STATUS: *** FAIL *** ({fail_count} step(s) failed)")
    print("=============================================")
    print("Audit complete for F0017.")
    print(f"  ✓ F0017 — {pass_count}/{total_steps} PASS")


if __name__ == "__main__":
    run()