"""
Heat diffusion equation (∂T/∂t = κ∇²T).

Assertion-based CAS audit block.
Pillar: Thermodynamics | Chain: energy conservation + Fourier's law + u=ρcT → diffusion equation
CalRef: Math Appendix §4D, Thermodynamics Calibration §4D
"""


def run():
    from sympy import symbols, diff, simplify, sin, cos, exp, pi

    print('=== CAS AUDIT: F0030 — Heat diffusion equation ===\n')

    pass_count = 0
    fail_count = 0
    total_steps = 0

    x, t = symbols('x t', real=True)
    rho, c_sp, k_th, kappa = symbols('rho c_sp k_th kappa', real=True, positive=True)
    alpha_w = symbols('alpha_w', real=True, positive=True)

    # Test solution
    T_sol = exp(-kappa * alpha_w**2 * t) * sin(alpha_w * x)
    dT_dt = diff(T_sol, t)
    d2T_dx2 = diff(T_sol, x, 2)

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

    # --- Step 1: Sinusoidal mode solution ---
    res1 = simplify(dT_dt - kappa * d2T_dx2)
    total_steps += 1
    if simplify(res1) == 0:
        print('  Step 1  PASS — T=e^{−κα²t}sin(αx) satisfies ∂T/∂t = κ∇²T')
        pass_count += 1
    else:
        print('  Step 1  FAIL')
        fail_count += 1

    # --- Step 2: Fourier substitution chain ---
    conservation_residual = rho * c_sp * dT_dt - k_th * d2T_dx2
    conservation_sub = conservation_residual.subs(kappa, k_th/(rho*c_sp))
    res2 = simplify(conservation_sub)
    total_steps += 1
    if simplify(res2) == 0:
        print('  Step 2  PASS — ρc·∂T/∂t − k·∂²T/∂x² = 0 when κ=k/(ρc)')
        pass_count += 1
    else:
        print('  Step 2  FAIL')
        fail_count += 1

    # --- Step 3: Thermal diffusivity definition ---
    res3_sub = simplify(k_th - rho * c_sp * kappa).subs(kappa, k_th/(rho*c_sp))
    total_steps += 1
    if simplify(res3_sub) == 0:
        print('  Step 3  PASS — κ = k/(ρc) (thermal diffusivity definition)')
        pass_count += 1
    else:
        print('  Step 3  FAIL')
        fail_count += 1

    # --- Step 4: Cosine mode ---
    T_cos = exp(-kappa * alpha_w**2 * t) * cos(alpha_w * x)
    dT_cos_dt = diff(T_cos, t)
    d2T_cos_dx2 = diff(T_cos, x, 2)
    res4 = simplify(dT_cos_dt - kappa * d2T_cos_dx2)
    total_steps += 1
    if simplify(res4) == 0:
        print('  Step 4  PASS — T=e^{−κα²t}cos(αx) also satisfies ∂T/∂t = κ∇²T')
        pass_count += 1
    else:
        print('  Step 4  FAIL')
        fail_count += 1

    # --- Step 5: Superposition (linearity) ---
    a_coeff, b_coeff = symbols('a_coeff b_coeff', real=True)
    T_super = a_coeff * T_sol + b_coeff * T_cos
    dT_super_dt = diff(T_super, t)
    d2T_super_dx2 = diff(T_super, x, 2)
    res5 = simplify(dT_super_dt - kappa * d2T_super_dx2)
    total_steps += 1
    if simplify(res5) == 0:
        print('  Step 5  PASS — Superposition preserved (linearity of diffusion eq)')
        pass_count += 1
    else:
        print('  Step 5  FAIL')
        fail_count += 1

    # --- Step 6: Two-sided localisation ---
    wrong_residual = simplify(dT_dt - kappa * diff(T_sol, x))
    total_steps += 1
    if not (simplify(wrong_residual) == 0):
        print('  Step 6  PASS — Wrong PDE (∂T/∂t = κ∂T/∂x) correctly rejected')
        pass_count += 1
    else:
        print('  Step 6  FAIL')
        fail_count += 1

    # --- Step 7: Steady-state limit ---
    A_lin, B_lin = symbols('A_lin B_lin', real=True)
    T_steady = A_lin * x + B_lin
    d2T_steady = diff(T_steady, x, 2)
    res7 = simplify(d2T_steady)
    total_steps += 1
    if simplify(res7) == 0:
        print('  Step 7  PASS — Steady state: T=Ax+B satisfies ∇²T=0 (Laplace limit)')
        pass_count += 1
    else:
        print('  Step 7  FAIL')
        fail_count += 1

    # --- Step 8: Decay rate depends on wave number ---
    alpha1, alpha2 = symbols('alpha1 alpha2', real=True, positive=True)
    rate1 = kappa * alpha1**2
    rate2 = kappa * alpha2**2
    ratio = simplify(rate2.subs(alpha2, 2*alpha1) / rate1)
    res8 = simplify(ratio - 4)
    total_steps += 1
    if simplify(res8) == 0:
        print('  Step 8  PASS — Decay rate ∝ α²: doubling α quadruples decay')
        pass_count += 1
    else:
        print('  Step 8  FAIL')
        fail_count += 1

    # --- Step 9: Numerical — iron bar diffusion ---
    k_num = 80.0
    rho_num = 7874.0
    c_num = 449.0
    kappa_num = k_num / (rho_num * c_num)
    L_num = 0.1
    alpha_num = pi / L_num
    tau_num = 1.0 / (kappa_num * alpha_num**2)
    kappa_expected = 2.263e-5
    tau_expected = 44.8
    total_steps += 1
    if abs(kappa_num - kappa_expected)/kappa_expected < 0.01 and abs(tau_num - tau_expected)/tau_expected < 0.01:
        print(f'  Step 9  PASS — Iron: κ={kappa_num:.3e} m²/s, τ={tau_num:.1f} s (L=0.1m fundamental mode)')
        pass_count += 1
    else:
        print('  Step 9  FAIL')
        fail_count += 1

    # --- Step 10: Unit consistency ---
    D_dim = symbols('D_dim', real=True, positive=True)
    k_proxy = rho * c_sp * D_dim
    kappa_proxy = k_proxy / (rho * c_sp)
    res10 = simplify(kappa_proxy - D_dim)
    total_steps += 1
    if simplify(res10) == 0:
        print('  Step 10 PASS — [κ] = [m²/s], [∂T/∂t] = [κ∇²T] = [K/s]')
        pass_count += 1
    else:
        print('  Step 10 FAIL')
        fail_count += 1

    # --- Step 11: Self-test — wrong sign in Fourier's law ---
    wrong_antidiff = simplify(dT_dt + kappa * d2T_dx2)
    total_steps += 1
    if not (simplify(wrong_antidiff) == 0):
        print('  Step 11 PASS — Wrong Fourier sign → anti-diffusion detected')
        pass_count += 1
    else:
        print('  Step 11 FAIL')
        fail_count += 1

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

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


if __name__ == '__main__':
    run()