"""
Hydrogen spectral fine structure (En + ΔEn).

Assertion-based CAS audit block.
Pillar: Particle Mechanics | Chain: baseline Rydberg → fine-structure correction → factorised En
CalRef: Math Appendix §5D, Calibration §5D

Structure mirrors cas_F020.txt exactly.
Every derivation step produces a verifiable symbolic assertion.
Final output: PASS or FAIL with step-level detail.
"""


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

    print('=== CAS AUDIT: F0020 — Hydrogen fine structure ===\n')

    pass_count = 0
    fail_count = 0
    total_steps = 0

    # ---- A. INPUTS ----
    R_inf, h_pl, c_light, n = symbols('R_inf h_pl c_light n', real=True, positive=True)
    f_val = symbols('f_val', real=True)

    # Nonrelativistic baseline
    E0_n = -R_inf * h_pl * c_light / n**2

    # Fine-structure correction template: ΔEn = E0_n * f(α,n,j)
    Delta_En = E0_n * f_val

    # Total energy
    E_n_def = E0_n + Delta_En
    E_n_fac = E0_n * (1 + f_val)

    print('Section A: Inputs defined.')
    print('  E0_n = -R∞hc/n²')
    print('  E_n = E0_n + E0_n·f = E0_n·(1+f)\n')

    # ---- B. ASSUMPTIONS / DOMAINS ----
    print('Section B: Assumptions set (R∞,h,c > 0, n positive integer).\n')

    # ---- C. ALLOWED LEMMAS ----
    print('Section C: Lemmas declared.')
    print('  C.1: Fine-structure template: ΔEn = E0_n·f')
    print('  C.2: Factorisation: E_n = E0_n·(1+f)\n')

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

    # --- Step 1: Additive → factorised form ---
    res1 = simplify(E_n_def - E_n_fac)
    total_steps += 1
    if simplify(res1) == 0:
        print('  Step 1  PASS — E0 + E0·f = E0·(1+f)')
        pass_count += 1
    else:
        print(f'  Step 1  FAIL — residual: {res1}')
        fail_count += 1

    # --- Step 2: Explicit factorised form with Rydberg ---
    E_n_explicit = -R_inf * h_pl * c_light / n**2 * (1 + f_val)
    res2 = simplify(E_n_fac - E_n_explicit)
    total_steps += 1
    if simplify(res2) == 0:
        print('  Step 2  PASS — En = -(R∞hc/n²)·[1+f]')
        pass_count += 1
    else:
        print(f'  Step 2  FAIL — residual: {res2}')
        fail_count += 1

    # --- Step 3: f=0 recovers baseline ---
    E_n_f0 = E_n_fac.subs(f_val, 0)
    res3 = simplify(E_n_f0 - E0_n)
    total_steps += 1
    if simplify(res3) == 0:
        print('  Step 3  PASS — f=0 → En = E0_n (baseline recovery)')
        pass_count += 1
    else:
        print(f'  Step 3  FAIL — residual: {res3}')
        fail_count += 1

    # --- Step 4: ΔEn / E0_n = f (dimensionless ratio) ---
    ratio = simplify(Delta_En / E0_n)
    res4 = simplify(ratio - f_val)
    total_steps += 1
    if simplify(res4) == 0:
        print('  Step 4  PASS — ΔEn / E0_n = f (dimensionless ratio)')
        pass_count += 1
    else:
        print(f'  Step 4  FAIL — residual: {res4}')
        fail_count += 1

    # --- Step 5: En / E0_n = 1 + f ---
    ratio_full = simplify(E_n_fac / E0_n)
    res5 = simplify(ratio_full - (1 + f_val))
    total_steps += 1
    if simplify(res5) == 0:
        print('  Step 5  PASS — En / E0_n = 1 + f')
        pass_count += 1
    else:
        print(f'  Step 5  FAIL — residual: {res5}')
        fail_count += 1

    # --- Step 6: Perturbative regime ---
    delta_check = simplify(E_n_fac - E0_n)
    res6 = simplify(delta_check - Delta_En)
    total_steps += 1
    if simplify(res6) == 0:
        print('  Step 6  PASS — En - E0_n = ΔEn (perturbative structure)')
        pass_count += 1
    else:
        print(f'  Step 6  FAIL — residual: {res6}')
        fail_count += 1

    # --- Step 7: Transition energy between levels ---
    n1, n2, f1, f2 = symbols('n1 n2 f1 f2', real=True, positive=True)
    E_n1 = -R_inf * h_pl * c_light / n1**2 * (1 + f1)
    E_n2 = -R_inf * h_pl * c_light / n2**2 * (1 + f2)
    Delta_E_trans = simplify(E_n2 - E_n1)

    Delta_E_0 = Delta_E_trans.subs([(f1, 0), (f2, 0)])
    Delta_E_expected = -R_inf * h_pl * c_light * (1/n2**2 - 1/n1**2)
    res7 = simplify(Delta_E_0 - Delta_E_expected)
    total_steps += 1
    if simplify(res7) == 0:
        print('  Step 7  PASS — f=0 transition → standard Rydberg formula')
        pass_count += 1
    else:
        print(f'  Step 7  FAIL — residual: {res7}')
        fail_count += 1

    # --- Step 8: Numerical check — H Balmer-α (n=3→2) ---
    R_val = 1.0973731568539e7
    h_val = 6.62607015e-34
    c_val = 2.99792458e8
    eV_conv = 1.602176634e-19

    DE_Balmer = R_val * h_val * c_val * (1/4 - 1/9)
    DE_expected_eV = 1.889
    DE_calc_eV = DE_Balmer / eV_conv

    total_steps += 1
    if abs(DE_calc_eV - DE_expected_eV) < 0.01:
        print(f'  Step 8  PASS — Balmer-α: {DE_calc_eV:.4f} eV (expected ≈ 1.889 eV)')
        pass_count += 1
    else:
        print(f'  Step 8  FAIL — Balmer-α numerical mismatch')
        fail_count += 1

    # --- Step 9: Fine-structure splitting estimate ---
    alpha_val = 7.2973525693e-3
    f_est = alpha_val**2
    DE_fine_est = R_val * h_val * c_val / 4 * f_est / eV_conv

    total_steps += 1
    if DE_fine_est > 1e-5 and DE_fine_est < 1e-3:
        print(f'  Step 9  PASS — fine-structure scale: {DE_fine_est:.2e} eV (order α²·E0)')
        pass_count += 1
    else:
        print(f'  Step 9  FAIL — fine-structure estimate out of range')
        fail_count += 1

    # --- Step 10: Self-test — wrong factorisation ---
    E_wrong = E0_n * (1 - f_val)
    res_wrong = simplify(E_wrong - E_n_fac)

    total_steps += 1
    if not (simplify(res_wrong) == 0):
        print('  Step 10a PASS — wrong sign (1-f) detected as incorrect')
        pass_count += 1
    else:
        print('  Step 10a FAIL — wrong sign not detected')
        fail_count += 1

    res_quantify = simplify(res_wrong - (-2 * E0_n * f_val))
    total_steps += 1
    if simplify(res_quantify) == 0:
        print('  Step 10b PASS — wrong residual = -2·E0·f (quantified)')
        pass_count += 1
    else:
        print('  Step 10b FAIL — wrong residual mismatch')
        fail_count += 1

    print('---------------------------------------------\n')

    # ---- E. CHECK OUTPUTS ----
    print('Section E: Output checks')
    print('---------------------------------------------')
    print('  Unit check:')
    print('    E0_n: [R∞]·[h]·[c]/[n²] = (1/m)·(J·s)·(m/s)/1 = J  ✓')
    print('    ΔEn: [E]·[dimensionless] = J  ✓\n')

    final_expr = simplify(E_n_def - E_n_fac)
    total_steps += 1
    if simplify(final_expr) == 0:
        print('  CAS simplification: E_n_def - E_n_fac = 0  PASS')
        pass_count += 1
    else:
        print(f'  CAS simplification: FAIL (residual: {final_expr})')
        fail_count += 1

    print('---------------------------------------------\n')

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


if __name__ == '__main__':
    run()