#!/usr/bin/env python3
"""
煤炭燃烧热辐射微积分计算演示
Coal Combustion Heat Radiation Calculation using Calculus
"""

import math

# 物理常数
STEFAN_BOLTZMANN = 5.67e-8  # W/(m²·K⁴)
COAL_HEAT_VALUE = 25e6      # J/kg (煤炭热值)
EMISSIVITY = 0.8            # 煤炭发射率

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

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

def simpson_integration(func, a, b, n=1000):
    """辛普森积分法数值积分"""
    if n % 2 == 1:
        n += 1
    h = (b - a) / n
    x = a
    sum_val = func(x)
    
    for i in range(1, n):
        x += h
        if i % 2 == 0:
            sum_val += 2 * func(x)
        else:
            sum_val += 4 * func(x)
    
    sum_val += func(b)
    return sum_val * h / 3

def total_radiated_energy(t_max, T0=1500, A=1.0):
    """积分计算总辐射能量: E = ∫₀ᵗ P(t) dt"""
    def integrand(t):
        T = temperature_profile(t, T0)
        return radiation_power(T, A)
    
    result = simpson_integration(integrand, 0, t_max)
    return result

def main():
    print("🔥 煤炭燃烧热辐射微积分计算")
    print("=" * 40)
    
    # 参数设置
    T0 = 1500  # 初始温度 K
    A = 0.1    # 表面积 m²
    t_max = 100  # 计算时间 s
    
    # 微积分计算总能量
    total_energy = total_radiated_energy(t_max, T0, A)
    
    print(f"初始温度: {T0} K")
    print(f"表面积: {A} m²")
    print(f"计算时间: {t_max} s")
    print(f"总辐射能量: {total_energy:.2e} J")
    
    # 瞬时功率在不同时间点
    time_points = [0, 10, 30, 60, 100]
    print(f"\n📊 不同时间点的温度和辐射功率:")
    print("时间(s) | 温度(K) | 功率(W)")
    print("-" * 30)
    
    for t in time_points:
        T_instant = temperature_profile(t, T0)
        P_instant = radiation_power(T_instant, A)
        print(f"{t:7.0f} | {T_instant:7.1f} | {P_instant:8.2f}")
    
    # 微积分原理说明
    print(f"\n🧮 微积分计算原理:")
    print("1. 温度函数: T(t) = T₀ × e^(-kt)")
    print("2. 辐射功率: P(t) = εσAT⁴(t)")
    print("3. 总能量积分: E = ∫₀ᵗ P(t) dt")
    print("4. 使用辛普森积分法数值求解")
    
    # 解析解对比 (指数衰减的四次方积分)
    k = 0.1
    analytical = (EMISSIVITY * STEFAN_BOLTZMANN * A * T0**4 / k) * (1 - math.exp(-k * t_max))
    print(f"\n🔬 解析解对比:")
    print(f"数值积分: {total_energy:.2e} J")
    print(f"解析解: {analytical:.2e} J")
    print(f"相对误差: {abs(total_energy - analytical) / analytical * 100:.2f}%")

if __name__ == "__main__":
    main()