clear; clc; close all;

%% ---------------- 参数 ----------------
L = 1.0;
Nx = 400;
dx = L / Nx;

a = 1.0;
dt = 0.4 * dx / a;
T_end = 2.5;
Nt = round(T_end / dt);

x = linspace(0, L-dx, Nx);
t_all = (0:Nt)*dt;

%% =========================
% 初始条件类型
% 1: 高斯
% 2: 方波
% 3: 双峰
% 4: 正弦
% 5: 多频
% 6: 随机
% 7: 平滑随机
%% =========================
IC_type = 6;

switch IC_type
    
    case 1
        phi = exp(-200*(x-0.5).^2);
        g   = zeros(size(x));
        x0_list = 0.5;
        
    case 2
        phi = double(abs(x-0.5)<0.1);
        g   = zeros(size(x));
        x0_list = [0.4, 0.6];
        
    case 3
        phi = exp(-200*(x-0.3).^2) + exp(-200*(x-0.7).^2);
        g   = zeros(size(x));
        x0_list = [0.3, 0.7];
        
    case 4  % 单一频率
        phi = sin(2*pi*x);
        g   = zeros(size(x));
        x0_list = 0.5;
        
    case 5  % 多频叠加（色散展示神器）
        phi = sin(2*pi*x) + 0.5*sin(6*pi*x);
        g   = zeros(size(x));
        x0_list = [];
        
    case 6  % 完全随机（非光滑）
        phi = rand(1,Nx) - 0.5;
        g   = zeros(size(x));
        x0_list = [];
        
    case 7  % 平滑随机（更物理）
        phi = rand(1,Nx);
        phi = smoothdata(phi,'gaussian',20);
        phi = phi - mean(phi);
        g   = zeros(size(x));
        x0_list = [];
        
    case 8  % 非零初速度（关键教学点）
        phi = exp(-200*(x-0.5).^2);
        g   = sin(2*pi*x);   % 初始速度
        x0_list = 0.5;
        
end

%% ---------------- 初始化 ----------------
r = a*dt/dx;

u_old = phi;
u = phi + dt*g + 0.5*r^2*(circshift(phi,-1) - 2*phi + circshift(phi,1));

U = zeros(Nt+1, Nx);
U(1,:) = phi;

%% ---------------- 图像初始化 ----------------
figure('Color','w','Position',[100 100 1300 520]);

subplot(1,2,1)
h_num   = plot(x, u, 'b-', 'LineWidth', 2); hold on;
h_exact = plot(x, u, 'r--', 'LineWidth', 1.5);
h_char_dummy = plot(nan, nan, 'k--');

legend([h_num, h_exact, h_char_dummy], ...
       {'数值解','解析解','特征线'}, ...
       'Location','northeast');

xlabel('x'); ylabel('u');
grid on; box on;

ylim([min(phi)-1, max(phi)+1]);

htitle = title('Wave Equation');

subplot(1,2,2)
h_img = imagesc(t_all, x, U);
set(gca,'YDir','normal');
xlabel('x'); ylabel('t');
title('x-t 演化图');
colorbar;

%% ---------------- 特征线 ----------------
subplot(1,2,1)
nChar = numel(x0_list);
h_char = gobjects(2*nChar,1);

for k = 1:nChar
    h_char(2*k-1) = xline(x0_list(k),'k--');
    h_char(2*k)   = xline(x0_list(k),'k--');
end

%% ---------------- 时间推进 ----------------
for n = 1:Nt
    t = n*dt;

    %% 数值解
    u_new = 2*u - u_old + r^2*(circshift(u,-1) - 2*u + circshift(u,1));

    u_old = u;
    u = u_new;

    U(n+1,:) = u;

    %% 解析解（含初速度项）
    xL = mod(x - a*t, L);
    xR = mod(x + a*t, L);

    phi_L = interp1(x, phi, xL, 'linear', 0);
    phi_R = interp1(x, phi, xR, 'linear', 0);

    % 👉 完整 d'Alembert（含 g）
    g_interp = interp1(x, g, x, 'linear', 0);
    
    % 数值近似积分项
    int_term = zeros(size(x));
    if any(g ~= 0)
        for i = 1:Nx
            xi1 = mod(x(i)-a*t, L);
            xi2 = mod(x(i)+a*t, L);
            idx = (x >= min(xi1,xi2) & x <= max(xi1,xi2));
            int_term(i) = trapz(x(idx), g(idx));
        end
        int_term = int_term / (2*a);
    end

    u_exact = 0.5*(phi_L + phi_R) + int_term;

    %% 可视化
    if mod(n,5)==0

        subplot(1,2,1)
        h_num.YData   = u;
        h_exact.YData = u_exact;

        for k = 1:nChar
            h_char(2*k-1).Value = mod(x0_list(k) - a*t, L);
            h_char(2*k).Value   = mod(x0_list(k) + a*t, L);
        end

        htitle.String = sprintf('t = %.3f | IC = %d', t, IC_type);

        subplot(1,2,2)
        set(h_img, 'CData', U);

        drawnow limitrate
        pause(0.01)
    end
end