"""Photoelectric relation (K_max = h*nu - phi).

Assertion-based CAS audit block.
Pillar: Electromagnetism | Chain: photon energy -> energy balance -> threshold
CalRef: Electromagnetism Math Appendix S5C, EM Calibration S5C
"""


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

    print("=== CAS AUDIT: F0019 — Photoelectric relation ===\n")

    pass_count = 0
    fail_count = 0
    total_steps = 0

    print("Section A: Inputs defined.")
    print("  E_gamma = h*nu, phi > 0\n")

    h_planck = symbols("h_planck", positive=True)
    nu = symbols("nu", positive=True)
    phi_work = symbols("phi_work", positive=True)
    K_kin = symbols("K_kin", real=True)
    e_charge = symbols("e_charge", positive=True)

    print("Section B: Single-photon single-electron, no additional loss.\n")
    print("Section C: Lemmas declared.\n")
    print("Section D: Step log")
    print("---------------------------------------------")

    # Step 1: Energy balance
    K_from_balance = h_planck * nu - phi_work
    step1_residual = simplify((K_from_balance + phi_work) - h_planck * nu)

    total_steps += 1
    if simplify(step1_residual) == 0:
        print("  Step 1  PASS — h*nu = phi + K => K = h*nu - phi")
        pass_count += 1
    else:
        print(f"  Step 1  FAIL — residual: {step1_residual}")
        fail_count += 1

    # Step 2: Threshold frequency
    nu_0 = phi_work / h_planck
    step2_residual = simplify(h_planck * nu_0 - phi_work)

    total_steps += 1
    if simplify(step2_residual) == 0:
        print("  Step 2  PASS — nu_0 = phi/h (threshold frequency)")
        pass_count += 1
    else:
        print(f"  Step 2  FAIL — residual: {step2_residual}")
        fail_count += 1

    # Step 3: K_max
    K_max = h_planck * nu - phi_work
    step3_residual = simplify(K_max - K_from_balance)

    total_steps += 1
    if simplify(step3_residual) == 0:
        print("  Step 3  PASS — K_max = h*nu - phi (no-loss model)")
        pass_count += 1
    else:
        print(f"  Step 3  FAIL — residual: {step3_residual}")
        fail_count += 1

    # Step 4: Stopping potential
    V0_expr = K_max / e_charge
    step4_residual = simplify(e_charge * V0_expr - K_max)

    total_steps += 1
    if simplify(step4_residual) == 0:
        print("  Step 4  PASS — e*V0 = K_max => V0 = (h*nu - phi)/e")
        pass_count += 1
    else:
        print(f"  Step 4  FAIL — residual: {step4_residual}")
        fail_count += 1

    # Step 5: Linear form
    from sympy import diff
    slope = diff(K_max, nu)
    intercept = K_max.subs(nu, 0)

    step5a_residual = simplify(slope - h_planck)
    step5b_residual = simplify(intercept - (-phi_work))

    total_steps += 1
    if simplify(step5a_residual) == 0 and simplify(step5b_residual) == 0:
        print("  Step 5  PASS — K_max vs nu: slope = h, intercept = -phi")
        pass_count += 1
    else:
        print("  Step 5  FAIL")
        fail_count += 1

    # Step 6: At threshold
    K_at_threshold = K_max.subs(nu, nu_0)
    step6_residual = simplify(K_at_threshold)

    total_steps += 1
    if simplify(step6_residual) == 0:
        print("  Step 6  PASS — K_max(nu_0) = 0 (threshold)")
        pass_count += 1
    else:
        print(f"  Step 6  FAIL — residual: {step6_residual}")
        fail_count += 1

    # Step 7: Below threshold
    nu_half = nu_0 / 2
    K_below = K_max.subs(nu, nu_half)
    K_below_simplified = simplify(K_below)
    expected_below = -phi_work / 2
    step7_residual = simplify(K_below_simplified - expected_below)

    total_steps += 1
    if simplify(step7_residual) == 0 and simplify(K_below_simplified) < 0:
        print("  Step 7  PASS — Below threshold: K_max(nu_0/2) = -phi/2 < 0")
        pass_count += 1
    else:
        print("  Step 7  FAIL")
        fail_count += 1

    # Step 8: Concrete numerical
    h_val = 6.626e-34
    eV_to_J = 1.602e-19
    phi_val = 2.1 * eV_to_J
    nu_val = 1.0e15

    K_max_num = h_val * nu_val - phi_val
    K_max_eV = K_max_num / eV_to_J
    K_expected_eV = h_val * nu_val / eV_to_J - 2.1
    rel_error = abs(K_max_eV - K_expected_eV) / abs(K_expected_eV)

    total_steps += 1
    if rel_error < 1e-10 and K_max_num > 0:
        print(f"  Step 8  PASS — Numerical: K_max(Cs, 1e15 Hz) = {K_max_eV:.3f} eV")
        pass_count += 1
    else:
        print(f"  Step 8  FAIL — rel error: {rel_error:.2e}, K = {K_max_num:.4e} J")
        fail_count += 1

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

    # Self-test: wrong balance
    K_wrong = phi_work - h_planck * nu
    wrong_residual = simplify(K_wrong - K_max)

    total_steps += 1
    if simplify(wrong_residual) != 0:
        print("  Self-test 1: Wrong balance (h*nu = phi - K) detected  PASS")
        pass_count += 1
    else:
        print("  Self-test 1: FAIL (wrong balance not detected)")
        fail_count += 1

    # Self-test: quantify
    expected_wrong = -2 * K_max
    wrong_quant = simplify(wrong_residual - expected_wrong)

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

    print("---------------------------------------------\n")
    print("=============================================")
    print("  F0019 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 F0019.")
    print(f"  ✓ F0019 — {pass_count}/{total_steps} PASS")


if __name__ == "__main__":
    run()