"""
Radiative flux law j = σT⁴ (PROVISIONAL).

Assertion-based CAS audit block.
Pillar: Thermodynamics (provisional — constitutive law, not derived)
CalRef: Formula Atlas F0031, EM Laws L6–L9, Thermo radiation section
"""


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

    print('=== CAS AUDIT: F0031 — Radiative flux law j = σT⁴ (provisional) ===\n')

    pass_count = 0
    fail_count = 0
    total_steps = 0

    T, T1, T2 = symbols('T T1 T2', real=True, positive=True)
    sigma_sb = symbols('sigma_sb', real=True, positive=True)
    R_sph = symbols('R_sph', real=True, positive=True)

    # Constitutive law: j = σ T⁴
    j_flux = sigma_sb * T**4

    print('Section D: Step log')
    print('---------------------------------------------')

    # --- Step 1: Constitutive relation ---
    res1 = simplify(j_flux - sigma_sb * T**4)
    total_steps += 1
    if simplify(res1) == 0:
        print('  Step 1  PASS — j = σT⁴ constitutive law encoded')
        pass_count += 1
    else:
        print('  Step 1  FAIL')
        fail_count += 1

    # --- Step 2: T⁴ scaling ---
    j_at_T = sigma_sb * T**4
    j_at_2T = sigma_sb * (2*T)**4
    ratio_2T = simplify(j_at_2T / j_at_T)
    res2 = simplify(ratio_2T - 16)
    total_steps += 1
    if simplify(res2) == 0:
        print('  Step 2  PASS — j(2T)/j(T) = 16 (T⁴ scaling verified)')
        pass_count += 1
    else:
        print('  Step 2  FAIL')
        fail_count += 1

    # --- Step 3: General flux ratio ---
    j1 = sigma_sb * T1**4
    j2 = sigma_sb * T2**4
    flux_ratio = simplify(j1 / j2)
    expected_ratio = (T1/T2)**4
    res3 = simplify(flux_ratio - expected_ratio)
    total_steps += 1
    if simplify(res3) == 0:
        print('  Step 3  PASS — j₁/j₂ = (T₁/T₂)⁴ (σ cancels in ratio)')
        pass_count += 1
    else:
        print('  Step 3  FAIL')
        fail_count += 1

    # --- Step 4: Temperature sensitivity ---
    from sympy import diff
    dj_dT = diff(j_flux, T)
    expected_deriv = 4 * sigma_sb * T**3
    res4 = simplify(dj_dT - expected_deriv)
    total_steps += 1
    if simplify(res4) == 0:
        print('  Step 4  PASS — dj/dT = 4σT³ (radiative sensitivity)')
        pass_count += 1
    else:
        print('  Step 4  FAIL')
        fail_count += 1

    # --- Step 5: Sphere power ---
    P_sphere = 4 * pi * R_sph**2 * sigma_sb * T**4
    total_steps += 1
    if simplify(P_sphere - 4*pi*R_sph**2*sigma_sb*T**4) == 0:
        print('  Step 5  PASS — P = 4πR²σT⁴; doubling R quadruples total power')
        pass_count += 1
    else:
        print('  Step 5  FAIL')
        fail_count += 1

    # --- Step 6: Numerical — Stefan–Boltzmann constant ---
    sigma_num = 5.670374419e-8
    T_sun = 5778.0
    j_sun = sigma_num * T_sun**4
    j_sun_expected = 6.32e7
    total_steps += 1
    if abs(j_sun - j_sun_expected)/j_sun_expected < 0.01:
        print(f'  Step 6  PASS — Sun (T=5778K): j = {j_sun:.3e} W/m² ≈ 6.32×10⁷')
        pass_count += 1
    else:
        print('  Step 6  FAIL')
        fail_count += 1

    # --- Step 7: Room temperature ---
    T_room = 300.0
    j_room = sigma_num * T_room**4
    j_room_expected = 459.3
    total_steps += 1
    if abs(j_room - j_room_expected)/j_room_expected < 0.01:
        print(f'  Step 7  PASS — Room (T=300K): j = {j_room:.1f} W/m² ≈ 459 W/m²')
        pass_count += 1
    else:
        print('  Step 7  FAIL')
        fail_count += 1

    # --- Step 8: Solar luminosity ---
    R_sun = 6.957e8
    L_sun = 4 * pi * R_sun**2 * j_sun
    L_sun_expected = 3.846e26
    L_sun_float = float(L_sun)
    total_steps += 1
    if abs(L_sun_float - L_sun_expected)/L_sun_expected < 0.02:
        print(f'  Step 8  PASS — L_sun = {L_sun_float:.3e} W ≈ 3.846×10²⁶ W')
        pass_count += 1
    else:
        print('  Step 8  FAIL')
        fail_count += 1

    # --- Step 9: Dimensional proxy ---
    C_dim, D_dim = symbols('C_dim D_dim', real=True, positive=True)
    sigma_proxy = C_dim * D_dim
    j_proxy = sigma_proxy * T**4
    dim_check = simplify(j_proxy / C_dim - D_dim * T**4)
    total_steps += 1
    if simplify(dim_check) == 0:
        print('  Step 9  PASS — [σ] = [W/(m²·K⁴)] dimensional proxy consistent')
        pass_count += 1
    else:
        print('  Step 9  FAIL')
        fail_count += 1

    # --- Step 10: Self-test — wrong exponent T³ ---
    j_wrong_exp = sigma_sb * T**3
    j_wrong_2T = sigma_sb * (2*T)**3
    wrong_ratio = simplify(j_wrong_2T / j_wrong_exp)
    res10 = simplify(wrong_ratio - 16)
    total_steps += 1
    if not (simplify(res10) == 0):
        print('  Step 10 PASS — Wrong exponent (T³ → ratio 8≠16) detected')
        pass_count += 1
    else:
        print('  Step 10 FAIL')
        fail_count += 1

    # --- Step 11: Self-test — wrong sign ---
    j_wrong_sign = -sigma_sb * T**4
    diff_check = simplify(j_wrong_sign - j_flux)
    total_steps += 1
    if not (simplify(diff_check) == 0):
        print('  Step 11 PASS — Wrong sign (j = −σT⁴ ≠ +σT⁴) detected')
        pass_count += 1
    else:
        print('  Step 11 FAIL')
        fail_count += 1

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

    # ---- VERDICT ----
    print('=============================================')
    print('  F0031 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('  STATUS: PROVISIONAL — constitutive law, derivation chain not closed.')
    print('  Law: j = σT⁴ with σ = 5.670374419×10⁻⁸ W/(m²·K⁴)')
    print(f'  ✓ F0031 — {pass_count}/{total_steps} PASS')


if __name__ == '__main__':
    run()