#!/usr/bin/env python3
"""
煤炭燃烧热辐射微积分计算演示 - 精确版本
Coal Combustion Heat Radiation - Accurate Calculus Implementation
"""

import math

# 物理常数
STEFAN_BOLTZMANN = 5.67e-8  # W/(m²·K⁴)
EMISSIVITY = 0.8            # 煤炭发射率

def temperature_profile(t, T0=1500, k=0.1):
    """温度随时间变化: T(t) = T₀ × e^(-kt)"""
    return T0 * math.exp(-k * t)

def radiation_power(T, A=1.0):
    """斯特藩-玻尔兹曼定律: P = εσAT⁴"""
    return EMISSIVITY * STEFAN_BOLTZMANN * A * T**4

def analytical_solution(t_max, T0=1500, A=1.0, k=0.1):
    """
    解析解: ∫₀ᵗ εσA(T₀e^(-kt))⁴ dt
    = εσAT₀⁴ ∫₀ᵗ e^(-4kt) dt
    = εσAT₀⁴ × [-1/(4k) × e^(-4kt)]₀ᵗ
    = εσAT₀⁴/(4k) × (1 - e^(-4kt))
    """
    coefficient = EMISSIVITY * STEFAN_BOLTZMANN * A * T0**4
    return coefficient / (4 * k) * (1 - math.exp(-4 * k * t_max))

def trapezoidal_integration(func, a, b, n=10000):
    """梯形积分法"""
    h = (b - a) / n
    result = 0.5 * (func(a) + func(b))
    for i in range(1, n):
        result += func(a + i * h)
    return result * h

def main():
    print("🔥 煤炭燃烧热辐射微积分计算 - 精确版本")
    print("=" * 50)
    
    # 参数设置
    T0 = 1500   # 初始温度 K
    A = 0.1     # 表面积 m²
    k = 0.1     # 衰减常数 s⁻¹
    t_max = 50  # 计算时间 s (减少以提高精度)
    
    print(f"🔧 物理参数:")
    print(f"  初始温度 T₀ = {T0} K")
    print(f"  表面积 A = {A} m²")
    print(f"  衰减常数 k = {k} s⁻¹")
    print(f"  发射率 ε = {EMISSIVITY}")
    print(f"  斯特藩-玻尔兹曼常数 σ = {STEFAN_BOLTZMANN} W/(m²·K⁴)")
    
    # 定义被积函数
    def integrand(t):
        T = temperature_profile(t, T0, k)
        return radiation_power(T, A)
    
    # 数值积分
    numerical_result = trapezoidal_integration(integrand, 0, t_max, 10000)
    
    # 解析解
    analytical_result = analytical_solution(t_max, T0, A, k)
    
    print(f"\n📊 积分计算结果 (t = 0 到 {t_max} s):")
    print(f"  数值积分 (梯形法): {numerical_result:.6e} J")
    print(f"  解析解:           {analytical_result:.6e} J")
    print(f"  相对误差:         {abs(numerical_result - analytical_result) / analytical_result * 100:.4f}%")
    
    # 微分计算 - 功率变化率
    print(f"\n🧮 微分分析:")
    t_points = [0, 5, 10, 20]
    print("时间(s) | 温度(K) | 功率(W) | dP/dt(W/s)")
    print("-" * 45)
    
    for t in t_points:
        T = temperature_profile(t, T0, k)
        P = radiation_power(T, A)
        
        # 功率对时间的导数: dP/dt = d/dt[εσAT₀⁴e^(-4kt)] = -4kεσAT₀⁴e^(-4kt)
        dP_dt = -4 * k * EMISSIVITY * STEFAN_BOLTZMANN * A * T0**4 * math.exp(-4 * k * t)
        
        print(f"{t:7.0f} | {T:7.1f} | {P:7.2f} | {dP_dt:9.2f}")
    
    # 积分的物理意义
    print(f"\n🔬 微积分物理意义:")
    print("1. 温度衰减: T(t) = T₀e^(-kt)")
    print("2. 辐射功率: P(t) = εσAT⁴(t) = εσAT₀⁴e^(-4kt)")
    print("3. 功率导数: dP/dt = -4kεσAT₀⁴e^(-4kt)")
    print("4. 总辐射能: E = ∫₀ᵗ P(τ)dτ = εσAT₀⁴/(4k)[1-e^(-4kt)]")
    
    # 能量效率分析
    initial_power = radiation_power(T0, A)
    efficiency_50 = (1 - math.exp(-4 * k * t_max)) * 100
    
    print(f"\n⚡ 能量分析:")
    print(f"  初始辐射功率: {initial_power:.2f} W")
    print(f"  {t_max}s内辐射效率: {efficiency_50:.1f}%")
    print(f"  平均功率: {analytical_result/t_max:.2f} W")

if __name__ == "__main__":
    main()