"""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

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)
"""


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

    print("=== CAS AUDIT: F0009 — Fundamental thermodynamic relation ===\n")

    pass_count = 0
    fail_count = 0
    total_steps = 0

    # ---- A. INPUTS ----
    S, V_vol, N = symbols("S V_vol N", positive=True)
    T_therm = symbols("T_therm", positive=True)
    p_therm = symbols("p_therm", positive=True)
    mu_therm = symbols("mu_therm", real=True)

    # U as symbolic function of (S, V, N)
    U = Function("U")(S, V_vol, N)

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

    # ---- B. ASSUMPTIONS / DOMAINS ----
    print("Section B: U is C^1, homogeneous equilibrium system, single species.\n")

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

    # ---- D. STEP LOG ----
    print("Section D: Step log")
    print("---------------------------------------------")

    # Total differential of U
    U_S = diff(U, S)
    U_V = diff(U, V_vol)
    U_N = diff(U, N)

    # Differential forms
    dS_form, dV_form, dN_form = symbols("dS_form dV_form dN_form", real=True)
    dU_total = U_S * dS_form + U_V * dV_form + U_N * dN_form

    # Step 1: Check structure
    total_steps += 1
    if dS_form in dU_total.free_symbols and dV_form in dU_total.free_symbols and dN_form in dU_total.free_symbols:
        print("  Step 1  PASS — Total differential has 3 terms (dS, dV, dN)")
        pass_count += 1
    else:
        print("  Step 1  FAIL — Missing differential forms")
        fail_count += 1

    # Step 2: Thermodynamic identifications
    dU_thermo = T_therm * dS_form + (-p_therm) * dV_form + mu_therm * dN_form
    dU_substituted = dU_total.subs([(U_S, T_therm), (U_V, -p_therm), (U_N, mu_therm)])
    step2_residual = simplify(dU_substituted - dU_thermo)

    total_steps += 1
    if simplify(step2_residual) == 0:
        print("  Step 2  PASS — Substitution: dU = T*dS - p*dV + mu*dN")
        pass_count += 1
    else:
        print(f"  Step 2  FAIL — Substitution residual: {step2_residual}")
        fail_count += 1

    # Step 3: Verify sign of pressure term
    coeff_dV = simplify(diff(dU_thermo, dV_form))
    expected_coeff_dV = -p_therm
    step3_residual = simplify(coeff_dV - expected_coeff_dV)

    total_steps += 1
    if simplify(step3_residual) == 0:
        print("  Step 3  PASS — Coefficient of dV = -p (sign correct)")
        pass_count += 1
    else:
        print(f"  Step 3  FAIL — dV coefficient: {coeff_dV} (expected -p)")
        fail_count += 1

    # 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 += 1
    if simplify(step4_residual) == 0:
        print("  Step 4  PASS — Coefficient of dS = T")
        pass_count += 1
    else:
        print(f"  Step 4  FAIL — dS coefficient: {coeff_dS} (expected T)")
        fail_count += 1

    # 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 += 1
    if simplify(step5_residual) == 0:
        print("  Step 5  PASS — Coefficient of dN = mu")
        pass_count += 1
    else:
        print(f"  Step 5  FAIL — dN coefficient: {coeff_dN} (expected mu)")
        fail_count += 1

    # Step 6: Maxwell relation from cross-derivatives
    cross_SV = diff(U, S, V_vol)
    cross_VS = diff(U, V_vol, S)
    step6_residual = simplify(cross_SV - cross_VS)

    total_steps += 1
    if simplify(step6_residual) == 0:
        print("  Step 6  PASS — Maxwell: d^2U/dSdV = d^2U/dVdS (cross-derivatives commute)")
        pass_count += 1
    else:
        print(f"  Step 6  FAIL — Cross-derivative residual: {step6_residual}")
        fail_count += 1

    # 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 += 1
    if simplify(step7_residual) == 0:
        print("  Step 7  PASS — Maxwell: d^2U/dSdN = d^2U/dNdS")
        pass_count += 1
    else:
        print(f"  Step 7  FAIL — (S,N) cross-derivative residual: {step7_residual}")
        fail_count += 1

    # 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 += 1
    if simplify(step8_residual) == 0:
        print("  Step 8  PASS — Maxwell: d^2U/dVdN = d^2U/dNdV")
        pass_count += 1
    else:
        print(f"  Step 8  FAIL — (V,N) cross-derivative residual: {step8_residual}")
        fail_count += 1

    # Step 9: Euler relation consistency
    S_val, V_val, N_val = symbols("S_val V_val N_val", positive=True)
    dT_form, dp_form, dmu_form = symbols("dT_form dp_form dmu_form", real=True)

    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
    )
    dU_fundamental = T_therm * dS_form - p_therm * dV_form + mu_therm * dN_form
    gibbs_duhem = simplify(dU_euler_full - dU_fundamental)
    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 += 1
    if simplify(step9_residual) == 0:
        print("  Step 9  PASS — Gibbs-Duhem: S*dT - V*dp + N*dmu = 0 (structure verified)")
        pass_count += 1
    else:
        print(f"  Step 9  FAIL — Gibbs-Duhem residual: {step9_residual}")
        fail_count += 1

    print("---------------------------------------------\n")

    # ---- E. CHECK OUTPUTS ----
    print("Section E: Output checks")
    print("---------------------------------------------")
    print("  Unit check:")
    print("    T*dS: [K]*[J/K] = [J]")
    print("    p*dV: [Pa]*[m^3] = [N/m^2]*[m^3] = [J]")
    print("    mu*dN: [J/mol]*[mol] = [J]")
    print("    dU: [J] -- all terms consistent")
    print("    PASS\n")

    # Self-test: wrong sign on pressure term
    dU_wrong = T_therm * dS_form + p_therm * dV_form + mu_therm * dN_form
    wrong_residual = simplify(dU_wrong - dU_thermo)

    total_steps += 1
    if simplify(wrong_residual) != 0:
        print("  Self-test: wrong pressure sign (+p) detected as different  PASS")
        pass_count += 1
    else:
        print("  Self-test: FAIL (wrong sign not detected!)")
        fail_count += 1

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

    total_steps += 1
    if simplify(wrong_check) == 0:
        print("  Self-test: wrong - correct = 2*p*dV (quantified)  PASS")
        pass_count += 1
    else:
        print(f"  Self-test: FAIL (wrong residual = {wrong_residual}, expected 2*p*dV)")
        fail_count += 1

    print("---------------------------------------------\n")

    # ---- VERDICT ----
    print("=============================================")
    print("  F0009 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 F0009.")
    print(f"  ✓ F0009 — {pass_count}/{total_steps} PASS")


if __name__ == "__main__":
    run()