"""
Hydrogen Balmer H-alpha dictionary hook.

Assertion-based CAS audit block.
Pillar: Spectroscopy Anchors | Chain: Rydberg -> Balmer -> H-alpha
CalRef: Spectroscopy Anchors section

Structure mirrors cas_F07.txt (= F0008) sections A-E.
Verifies:
  1. Rydberg wavelength formula derivation (symbolic)
  2. Balmer H-alpha specialization (n=3->2)
  3. H-beta (4->2) and H-gamma (5->2)
  4. Wavelength ratio checks (R_inf cancels)
  5. Numerical comparison with NIST ASD values
  6. Rydberg constant from CODATA fundamentals
"""


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

    print('=== CAS AUDIT: F0008 -- Hydrogen Balmer H-alpha dictionary hook ===\n')

    pass_count = 0
    fail_count = 0
    total_steps = 0

    # ---- A. INPUTS ----
    R_inf = symbols('R_inf', positive=True)
    h_planck = symbols('h_planck', positive=True)
    c_light = symbols('c_light', positive=True)
    n_i = symbols('n_i', positive=True, integer=True)
    n_f = symbols('n_f', positive=True, integer=True)
    lambda_sym = symbols('lambda_sym', positive=True)

    print('Section A: Inputs defined.')
    print('  R_inf, h, c, quantum numbers n_i, n_f\n')

    # ---- B. ASSUMPTIONS / DOMAINS ----
    print('Section B: Nonrelativistic H, vacuum, infinite proton mass.')
    print('  CODATA 2022 values. Jc tolerance <= 1e-6 (ppm).\n')

    # ---- C. ALLOWED LEMMAS ----
    print('Section C: Lemmas declared.')
    print('  C.1: E_gamma = hc/lambda = E_nf - E_ni')
    print('  C.2: Balmer specialization nf=2')
    print('  C.3: Wavelength ratios: R_inf cancels\n')

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

    # --- Step 1: Derive Rydberg wavelength formula symbolically ---
    E_nf = -h_planck * c_light * R_inf / n_f ** 2
    E_ni = -h_planck * c_light * R_inf / n_i ** 2

    delta_E = E_nf - E_ni
    delta_E = simplify(delta_E)

    expected_delta_E = -h_planck * c_light * R_inf * (1 / n_f ** 2 - 1 / n_i ** 2)
    step1_residual = simplify(delta_E - expected_delta_E)

    total_steps += 1
    if simplify(step1_residual) == 0:
        print('  Step 1  PASS — Delta E = -hcR_inf*(1/nf^2 - 1/ni^2)')
        pass_count += 1
    else:
        print(f'  Step 1  FAIL — Delta E residual: {step1_residual}')
        fail_count += 1

    # --- Step 2: Cancel h*c to get 1/lambda ---
    inv_lambda = R_inf * (1 / n_f ** 2 - 1 / n_i ** 2)

    ratio_check = simplify(-delta_E / (h_planck * c_light))
    step2_residual = simplify(ratio_check - inv_lambda)

    total_steps += 1
    if simplify(step2_residual) == 0:
        print('  Step 2  PASS — h*c cancels: 1/lambda = R_inf*(1/nf^2 - 1/ni^2)')
        pass_count += 1
    else:
        print(f'  Step 2  FAIL — Cancellation residual: {step2_residual}')
        fail_count += 1

    # --- Step 3: h and c independence ---
    inv_lambda_vars = inv_lambda.free_symbols
    has_h = h_planck in inv_lambda_vars
    has_c = c_light in inv_lambda_vars

    total_steps += 1
    if not has_h and not has_c:
        print('  Step 3  PASS — 1/lambda independent of h and c')
        pass_count += 1
    else:
        print(f'  Step 3  FAIL — 1/lambda still contains h={has_h}, c={has_c}')
        fail_count += 1

    # --- Step 4: Balmer H-alpha specialization (3->2) ---
    inv_lambda_Ha = inv_lambda.subs([(n_f, 2), (n_i, 3)])
    inv_lambda_Ha = simplify(inv_lambda_Ha)
    inv_lambda_Ha_expected = R_inf * Rational(5, 36)

    step4_residual = simplify(inv_lambda_Ha - inv_lambda_Ha_expected)

    total_steps += 1
    if simplify(step4_residual) == 0:
        print('  Step 4  PASS — H-alpha: 1/lambda = (5/36)*R_inf')
        pass_count += 1
    else:
        print(f'  Step 4  FAIL — H-alpha residual: {step4_residual}')
        fail_count += 1

    # --- Step 5: H-beta specialization (4->2) ---
    inv_lambda_Hb = inv_lambda.subs([(n_f, 2), (n_i, 4)])
    inv_lambda_Hb = simplify(inv_lambda_Hb)
    inv_lambda_Hb_expected = R_inf * Rational(3, 16)

    step5_residual = simplify(inv_lambda_Hb - inv_lambda_Hb_expected)

    total_steps += 1
    if simplify(step5_residual) == 0:
        print('  Step 5  PASS — H-beta: 1/lambda = (3/16)*R_inf')
        pass_count += 1
    else:
        print(f'  Step 5  FAIL — H-beta residual: {step5_residual}')
        fail_count += 1

    # --- Step 6: H-gamma specialization (5->2) ---
    inv_lambda_Hg = inv_lambda.subs([(n_f, 2), (n_i, 5)])
    inv_lambda_Hg = simplify(inv_lambda_Hg)
    inv_lambda_Hg_expected = R_inf * Rational(21, 100)

    step6_residual = simplify(inv_lambda_Hg - inv_lambda_Hg_expected)

    total_steps += 1
    if simplify(step6_residual) == 0:
        print('  Step 6  PASS — H-gamma: 1/lambda = (21/100)*R_inf')
        pass_count += 1
    else:
        print(f'  Step 6  FAIL — H-gamma residual: {step6_residual}')
        fail_count += 1

    # --- Step 7: Wavelength ratio H-alpha/H-beta (R_inf cancels) ---
    ratio_Ha_Hb = simplify(inv_lambda_Hb_expected / inv_lambda_Ha_expected)
    expected_ratio = Rational(27, 20)

    step7_residual = simplify(ratio_Ha_Hb - expected_ratio)

    total_steps += 1
    if simplify(step7_residual) == 0:
        print('  Step 7  PASS — lambda_Ha/lambda_Hb = 27/20 = 1.35 (R_inf cancels)')
        pass_count += 1
    else:
        print(f'  Step 7  FAIL — Ratio residual: {step7_residual}')
        fail_count += 1

    # --- Step 8: R_inf independence of wavelength ratios ---
    ratio_symbolic = simplify(inv_lambda_Hb / inv_lambda_Ha)
    has_R = R_inf in ratio_symbolic.free_symbols

    total_steps += 1
    if not has_R:
        print('  Step 8  PASS — Wavelength ratio independent of R_inf')
        pass_count += 1
    else:
        print(f'  Step 8  FAIL — Ratio still contains R_inf: {ratio_symbolic}')
        fail_count += 1

    # --- Step 9: Emission positivity check ---
    emission_factor_Ha = Rational(1, 4) - Rational(1, 9)
    emission_factor_Hb = Rational(1, 4) - Rational(1, 16)
    emission_factor_Hg = Rational(1, 4) - Rational(1, 25)

    total_steps += 1
    if emission_factor_Ha > 0 and emission_factor_Hb > 0 and emission_factor_Hg > 0:
        print('  Step 9  PASS — Emission positivity: 1/nf^2 - 1/ni^2 > 0 for Ha,Hb,Hg')
        pass_count += 1
    else:
        print('  Step 9  FAIL — Emission factor not positive')
        fail_count += 1

    # --- Step 10: Numerical R_inf from CODATA fundamentals ---
    m_e_val = 9.1093837139e-31
    e_val = 1.602176634e-19
    eps0_val = 8.8541878128e-12
    h_val = 6.62607015e-34
    c_val = 299792458

    R_inf_computed = m_e_val * e_val ** 4 / (8 * eps0_val ** 2 * h_val ** 3 * c_val)
    R_inf_CODATA = 10973731.568157

    R_inf_rel_error = abs(R_inf_computed - R_inf_CODATA) / R_inf_CODATA

    total_steps += 1
    if R_inf_rel_error < 1e-6:
        print(f'  Step 10 PASS — R_inf = {R_inf_computed:.6f} 1/m (rel error: {R_inf_rel_error:.2e})')
        pass_count += 1
    else:
        print(f'  Step 10 FAIL — R_inf = {R_inf_computed:.6f}, expected {R_inf_CODATA:.6f} (rel error: {R_inf_rel_error:.2e})')
        fail_count += 1

    # --- Step 11: Numerical H-alpha wavelength ---
    lambda_Ha_computed = 36 / (5 * R_inf_CODATA)
    lambda_Ha_nm = lambda_Ha_computed * 1e9
    lambda_Ha_NIST = 656.281

    lambda_Ha_rel_error = abs(lambda_Ha_nm - lambda_Ha_NIST) / lambda_Ha_NIST

    total_steps += 1
    if lambda_Ha_rel_error < 1e-3:
        print(f'  Step 11 PASS — lambda_Ha = {lambda_Ha_nm:.3f} nm (NIST: {lambda_Ha_NIST:.3f} nm, rel err: {lambda_Ha_rel_error:.2e})')
        pass_count += 1
    else:
        print(f'  Step 11 FAIL — lambda_Ha = {lambda_Ha_nm:.3f} nm (rel error: {lambda_Ha_rel_error:.2e})')
        fail_count += 1

    # --- Step 12: Numerical ratio check vs NIST ---
    lambda_Hb_NIST = 486.133
    lambda_Hg_NIST = 434.047

    ratio_exp = lambda_Ha_NIST / lambda_Hb_NIST
    ratio_theory = 27 / 20

    ratio_ppm = abs(ratio_exp - ratio_theory) / ratio_theory * 1e6

    total_steps += 1
    if ratio_ppm < 10:
        print(f'  Step 12 PASS — Ha/Hb ratio: theory={ratio_theory:.6f}, NIST={ratio_exp:.6f} ({ratio_ppm:.1f} ppm)')
        pass_count += 1
    else:
        print(f'  Step 12 FAIL — Ha/Hb ratio discrepancy: {ratio_ppm:.1f} ppm')
        fail_count += 1

    # --- Step 13: H-beta and H-gamma numerical wavelengths ---
    lambda_Hb_computed = 16 / (3 * R_inf_CODATA) * 1e9
    lambda_Hg_computed = 100 / (21 * R_inf_CODATA) * 1e9

    lambda_Hb_rel = abs(lambda_Hb_computed - lambda_Hb_NIST) / lambda_Hb_NIST
    lambda_Hg_rel = abs(lambda_Hg_computed - lambda_Hg_NIST) / lambda_Hg_NIST

    total_steps += 1
    if lambda_Hb_rel < 1e-3 and lambda_Hg_rel < 1e-3:
        print(f'  Step 13 PASS — Hb={lambda_Hb_computed:.3f} nm (NIST:{lambda_Hb_NIST:.3f}, err:{lambda_Hb_rel:.2e}), Hg={lambda_Hg_computed:.3f} nm (NIST:{lambda_Hg_NIST:.3f}, err:{lambda_Hg_rel:.2e})')
        pass_count += 1
    else:
        print(f'  Step 13 FAIL — Hb rel err: {lambda_Hb_rel:.2e}, Hg rel err: {lambda_Hg_rel:.2e}')
        fail_count += 1

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

    # ---- E. CHECK OUTPUTS ----
    print('Section E: Output checks')
    print('---------------------------------------------')

    print('  Unit check:')
    print('    R_inf: [1/m]')
    print('    1/lambda = R_inf*(..): [1/m] => lambda: [m]')
    print('    Ratios: [1/m]/[1/m] = [1] (dimensionless)')
    print('    PASS\n')

    # --- Self-test: wrong quantum numbers ---
    inv_lambda_wrong = inv_lambda.subs([(n_f, 3), (n_i, 2)])
    inv_lambda_wrong = simplify(inv_lambda_wrong)
    inv_lambda_wrong_val = float(inv_lambda_wrong.subs(R_inf, R_inf_CODATA))

    total_steps += 1
    if inv_lambda_wrong_val < 0:
        print('  Self-test: swapped n_i/n_f gives 1/lambda < 0 (unphysical)  PASS')
        pass_count += 1
    else:
        print('  Self-test: FAIL (swapped quantum numbers not detected)')
        fail_count += 1

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

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


if __name__ == '__main__':
    run()