import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate
from matplotlib.animation import FuncAnimation

plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei', 'WenQuanYi Micro Hei']
plt.rcParams['axes.unicode_minus'] = False

def f(x):
    return np.pi * np.cos(np.exp(x)) / (x ** np.e)

a, b = 1.0, 4.0
true_value, _ = integrate.quad(f, a, b)

def animate_trapezoidal(func, a, b, n=20, filename='trapezoidal.gif'):
    x = np.linspace(a, b, n + 1)
    y = func(x)
    h = (b - a) / n

    fig, ax = plt.subplots(figsize=(10, 5))
    xs = np.linspace(a, b, 500)
    ax.plot(xs, func(xs), 'k-', lw=2, label=r'$f(x)$')

    def update(i):
        ax.cla()
        ax.plot(xs, func(xs), 'k-', lw=2, label=r'$f(x)$')
        for j in range(0, i+1):
            xi = [x[j], x[j+1], x[j+1], x[j]]
            yi = [0, 0, y[j+1], y[j]]
            ax.fill(xi, yi, color='skyblue', alpha=0.6, edgecolor='blue', linewidth=0.8)
        approx = (h / 2) * (y[0] + 2 * np.sum(y[1:i+1]) + y[i+1])
        error = abs(approx - true_value)
        ax.set_title(f'梯形法 动画演示\n第{i+1}次迭代, 积分 ≈ {approx:.6f}, 误差 = {error:.2e}', fontsize=14)
        ax.set_xlabel('x')
        ax.set_ylabel('f(x)')
        ax.legend()
        ax.grid(True, linestyle='--', alpha=0.5)

    ani = FuncAnimation(fig, update, frames=n, interval=1000, repeat=False)
    ani.save(filename, writer='pillow')
    plt.show()

def animate_simpson(func, a, b, n=20, filename='simpson.gif'):
    if n % 2 == 1:
        n += 1
    x = np.linspace(a, b, n + 1)
    y = func(x)
    h = (b - a) / n

    fig, ax = plt.subplots(figsize=(10, 5))
    xs = np.linspace(a, b, 500)
    ax.plot(xs, func(xs), 'k-', lw=2, label=r'$f(x)$')

    def update(i):
        ax.cla()
        ax.plot(xs, func(xs), 'k-', lw=2, label=r'$f(x)$')
        for j in range(0, min((i+1)*2, n), 2):
            x0, x1, x2 = x[j], x[j+1], x[j+2]
            y0, y1, y2 = y[j], y[j+1], y[j+2]

            xx = np.linspace(x0, x2, 50)
            yy = ((xx-x1)*(xx-x2))/((x0-x1)*(x0-x2))*y0 + ((xx-x0)*(xx-x2))/((x1-x0)*(x1-x2))*y1 + ((xx-x0)*(xx-x1))/((x2-x0)*(x2-x1))*y2
            ax.fill_between(xx, 0, yy, color='lightcoral', alpha=0.6, edgecolor='red', linewidth=0.8)
        ax.set_title(f'辛普森法 动画演示\n第{i+1}次迭代', fontsize=14)
        ax.set_xlabel('x')
        ax.set_ylabel('f(x)')
        ax.legend()
        ax.grid(True, linestyle='--', alpha=0.5)

    ani = FuncAnimation(fig, update, frames=n//2, interval=1000, repeat=False)
    ani.save(filename, writer='pillow')
    plt.show()

def simpson_basic(func, a, b):
    c = (a + b) / 2.0
    return (b - a) / 6.0 * (func(a) + 4.0 * func(c) + func(b))

def adaptive_simpson_with_log(func, a, b, tol, max_depth=15, log=None, depth=0):
    if log isNone:
        log = []
    if depth > max_depth:
        log.append((a, b, True, depth))
        return simpson_basic(func, a, b)
    c = (a + b) / 2.0
    whole = simpson_basic(func, a, b)
    left = simpson_basic(func, a, c)
    right = simpson_basic(func, c, b)
    error = abs(left + right - whole)
    if error < 15 * tol:
        log.append((a, b, True, depth))
        return left + right
    else:
        log.append((a, b, False, depth))
        left_val = adaptive_simpson_with_log(func, a, c, tol/2, max_depth, log, depth+1)
        right_val = adaptive_simpson_with_log(func, c, b, tol/2, max_depth, log, depth+1)
        return left_val + right_val

# ✅ 修改后的自适应辛普森动画函数：使用抛物线填充，不是直线！
def plot_adaptive_simpson(func, a, b, tol=1e-5, filename='adaptive_simpson.gif'):
    log = []
    approx = adaptive_simpson_with_log(func, a, b, tol, max_depth=12, log=log)
    error = abs(approx - true_value)

    fig, ax = plt.subplots(figsize=(10, 5))
    xs = np.linspace(a, b, 500)
    ax.plot(xs, func(xs), 'k-', lw=2, label=r'$f(x)$')

    cmap = plt.get_cmap('viridis')
    max_depth = max(d for (_, _, _, d) in log)

    def update(i):
        ax.cla()
        ax.plot(xs, func(xs), 'k-', lw=2, label=r'$f(x)$')
        for j in range(i + 1):
            a_, b_, accepted, depth = log[j]
            if accepted:
                # 构造三点：左端点、中点、右端点
                x0, x2 = a_, b_
                x1 = (x0 + x2) / 2.0
                y0, y1, y2 = func(x0), func(x1), func(x2)

                # 在 [x0, x2] 上生成密集点用于绘制抛物线
                xx = np.linspace(x0, x2, 50)
                # 使用拉格朗日二次插值（与辛普森一致）
                yy = ((xx - x1) * (xx - x2)) / ((x0 - x1) * (x0 - x2)) * y0 + \
                     ((xx - x0) * (xx - x2)) / ((x1 - x0) * (x1 - x2)) * y1 + \
                     ((xx - x0) * (xx - x1)) / ((x2 - x0) * (x2 - x1)) * y2

                color = cmap(depth / max(max_depth, 1))
                ax.fill_between(xx, 0, yy, color=color, alpha=0.6, edgecolor='purple', linewidth=0.7)
        ax.set_title(f'自适应辛普森法 动画演示\n第{i+1}次迭代', fontsize=14)
        ax.set_xlabel('x')
        ax.set_ylabel('f(x)')
        ax.legend()
        ax.grid(True, linestyle='--', alpha=0.5)

    ani = FuncAnimation(fig, update, frames=len(log), interval=1000, repeat=False)
    ani.save(filename, writer='pillow')
    plt.show()

if __name__ == "__main__":
    print("参考真值（SciPy quad）:", f"{true_value:.10f}")

    # Generate animations and save them as GIF files.
    animate_trapezoidal(f, a, b, n=12, filename='trapezoidal.gif')
    animate_simpson(f, a, b, n=12, filename='simpson.gif')
    plot_adaptive_simpson(f, a, b, tol=1e-4, filename='adaptive_simpson.gif')