"""Ideal-gas calibration (PV=nRT, Z definition).

Assertion-based CAS audit block.
Pillar: Thermodynamics | Chain: ideal gas -> R_eff -> Z -> virial
CalRef: Thermodynamics Calibration S5Y.1-5Y.3

Verifies:
  1. R_eff = PV/(nT) definition
  2. Ideal-gas substitution: PV=nRT => R_eff=R
  3. Z = PV/(nRT) = R_eff/R
  4. Ideal limit: Z=1
  5. Virial expansion structure: Z = 1 + B*n/V + C*(n/V)^2 + ...
  6. Dimensional consistency
  7. Non-ideal detection (Z != 1 when B != 0)
"""


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

    print("=== CAS AUDIT: F0010 — Ideal-gas calibration (PV=nRT, Z) ===\n")

    pass_count = 0
    fail_count = 0
    total_steps = 0

    # ---- A. INPUTS ----
    P, V_gas, n, T_gas = symbols("P V_gas n T_gas", positive=True)
    R_gas = symbols("R_gas", positive=True)
    Z_comp = symbols("Z_comp", positive=True)
    B_vir, C_vir = symbols("B_vir C_vir", real=True)

    print("Section A: Inputs defined.")
    print("  P, V, n, T, R, R_eff, Z, virial coefficients B, C\n")

    # ---- B. ASSUMPTIONS / DOMAINS ----
    print("Section B: Single-component, classical, low density.\n")

    # ---- C. ALLOWED LEMMAS ----
    print("Section C: Lemmas declared.")
    print("  C.1: PV=nRT => R_eff=R, Z=1")
    print("  C.2: Z = R_eff/R")
    print("  C.3: Virial: Z = 1 + B*n/V + ...\n")

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

    # Step 1: R_eff definition
    R_eff_def = P * V_gas / (n * T_gas)

    total_steps += 1
    if R_gas not in R_eff_def.free_symbols:
        print("  Step 1  PASS — R_eff = PV/(nT) is independent of R (pure measurement)")
        pass_count += 1
    else:
        print("  Step 1  FAIL — R_eff definition contains R")
        fail_count += 1

    # Step 2: Ideal-gas substitution => R_eff = R
    P_ideal = n * R_gas * T_gas / V_gas
    R_eff_ideal = simplify(R_eff_def.subs(P, P_ideal))
    step2_residual = simplify(R_eff_ideal - R_gas)

    total_steps += 1
    if simplify(step2_residual) == 0:
        print("  Step 2  PASS — PV=nRT => R_eff = R")
        pass_count += 1
    else:
        print(f"  Step 2  FAIL — R_eff_ideal = {R_eff_ideal} (expected R)")
        fail_count += 1

    # Step 3: Z definition
    Z_def = P * V_gas / (n * R_gas * T_gas)
    Z_from_Reff = R_eff_def / R_gas
    step3_residual = simplify(Z_def - Z_from_Reff)

    total_steps += 1
    if simplify(step3_residual) == 0:
        print("  Step 3  PASS — Z = PV/(nRT) = R_eff/R")
        pass_count += 1
    else:
        print(f"  Step 3  FAIL — Z identity residual: {step3_residual}")
        fail_count += 1

    # Step 4: Ideal limit Z = 1
    Z_ideal = simplify(Z_def.subs(P, P_ideal))
    step4_residual = simplify(Z_ideal - 1)

    total_steps += 1
    if simplify(step4_residual) == 0:
        print("  Step 4  PASS — Ideal gas: Z = 1")
        pass_count += 1
    else:
        print(f"  Step 4  FAIL — Z_ideal = {Z_ideal} (expected 1)")
        fail_count += 1

    # Step 5: R_eff(ideal) = R (algebraic confirmation)
    step5_residual = simplify(simplify(R_eff_ideal) - R_gas)

    total_steps += 1
    if simplify(step5_residual) == 0:
        print("  Step 5  PASS — R_eff(ideal) = R (algebraic confirmation, all state vars cancel)")
        pass_count += 1
    else:
        print(f"  Step 5  FAIL — R_eff(ideal) - R = {step5_residual}")
        fail_count += 1

    # Step 6: Virial expansion structure
    rho = symbols("rho", positive=True)
    Z_virial = 1 + B_vir * rho + C_vir * rho**2
    Z_virial_ideal = Z_virial.subs(rho, 0)
    step6_residual = simplify(Z_virial_ideal - 1)

    total_steps += 1
    if simplify(step6_residual) == 0:
        print("  Step 6  PASS — Virial at rho=0: Z = 1 (ideal limit recovered)")
        pass_count += 1
    else:
        print(f"  Step 6  FAIL — Z_virial(rho=0) = {Z_virial_ideal}")
        fail_count += 1

    # Step 7: Virial leading correction
    from sympy import diff
    dZ_drho = diff(Z_virial, rho)
    dZ_at_zero = dZ_drho.subs(rho, 0)
    step7_residual = simplify(dZ_at_zero - B_vir)

    total_steps += 1
    if simplify(step7_residual) == 0:
        print("  Step 7  PASS — dZ/drho|_{rho=0} = B (leading virial coefficient)")
        pass_count += 1
    else:
        print(f"  Step 7  FAIL — dZ/drho(0) = {dZ_at_zero} (expected B)")
        fail_count += 1

    # Step 8: Second derivative gives C
    d2Z_drho2 = diff(Z_virial, rho, 2)
    d2Z_at_zero = d2Z_drho2.subs(rho, 0)
    step8_residual = simplify(d2Z_at_zero / 2 - C_vir)

    total_steps += 1
    if simplify(step8_residual) == 0:
        print("  Step 8  PASS — (1/2)*d^2Z/drho^2|_{rho=0} = C (third virial)")
        pass_count += 1
    else:
        print(f"  Step 8  FAIL — (1/2)*d2Z(0) = {d2Z_at_zero/2} (expected C)")
        fail_count += 1

    # Step 9: Virial substitution rho = n/V
    Z_virial_nV = Z_virial.subs(rho, n / V_gas)
    Z_virial_nV = simplify(Z_virial_nV)
    Z_virial_expected = 1 + B_vir * n / V_gas + C_vir * (n / V_gas) ** 2
    step9_residual = simplify(Z_virial_nV - Z_virial_expected)

    total_steps += 1
    if simplify(step9_residual) == 0:
        print("  Step 9  PASS — Z(n/V) = 1 + B*n/V + C*(n/V)^2")
        pass_count += 1
    else:
        print(f"  Step 9  FAIL — Substitution residual: {step9_residual}")
        fail_count += 1

    # Step 10: Numerical check with R = 8.314462618 J/(mol*K)
    k_B_val = 1.380649e-23
    N_A_val = 6.02214076e23
    R_computed = k_B_val * N_A_val
    R_CODATA = 8.314462618

    R_rel_error = abs(R_computed - R_CODATA) / R_CODATA

    total_steps += 1
    if R_rel_error < 1e-9:
        print(f"  Step 10 PASS — R = k_B*N_A = {R_computed:.9f} J/(mol*K) (rel err: {R_rel_error:.2e})")
        pass_count += 1
    else:
        print(f"  Step 10 FAIL — R = {R_computed:.9f} (rel error: {R_rel_error:.2e})")
        fail_count += 1

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

    # ---- E. CHECK OUTPUTS ----
    print("Section E: Output checks")
    print("---------------------------------------------")
    print("  Unit check:")
    print("    R_eff = PV/(nT): [Pa*m^3/(mol*K)] = [J/(mol*K)]")
    print("    Z = PV/(nRT): dimensionless")
    print("    Virial: rho = n/V [mol/m^3], B [m^3/mol], C [m^6/mol^2]")
    print("    PASS\n")

    # Self-test: non-ideal gas (B != 0) gives Z != 1
    Z_nonideal = Z_virial.subs([(B_vir, Rational(-1, 100)), (C_vir, 0), (rho, Rational(1, 10))])
    Z_nonideal = simplify(Z_nonideal)

    total_steps += 1
    if simplify(Z_nonideal) != 1:
        print(f"  Self-test: B=-0.01, rho=0.1 => Z = {Z_nonideal} != 1 (non-ideal detected)  PASS")
        pass_count += 1
    else:
        print("  Self-test: FAIL (non-ideal not detected)")
        fail_count += 1

    # Self-test: reciprocal consistency Z*R = R_eff
    Reff_from_Z = Z_def * R_gas
    Reff_from_Z = simplify(Reff_from_Z)
    step_recip = simplify(Reff_from_Z - R_eff_def)

    total_steps += 1
    if simplify(step_recip) == 0:
        print("  Self-test: Z*R = R_eff (reciprocal consistency)  PASS")
        pass_count += 1
    else:
        print("  Self-test: FAIL (Z*R != R_eff)")
        fail_count += 1

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

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


if __name__ == "__main__":
    run()