clear; clc; close all;

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

c = 1.0;
CFL = 0.5;
dt = CFL * dx / c;

T_end = 5.0;
Nt = round(T_end / dt);

x = linspace(0, L, Nx);
t = (0:Nt)*dt;

%% =========================
% 初始条件选择
% 1: Gaussian
% 2: 方波
% 3: 正弦
% 4: 多峰
% 5: 完全随机
% 6: 平滑随机
%% =========================
IC_type = 5;

switch IC_type
    
    case 1  % Gaussian
        u = exp(-100 * (x - 0.3).^2);
        
    case 2  % 方波
        u = zeros(size(x));
        u(x > 0.2 & x < 0.5) = 1;
        
    case 3  % 正弦
        u = sin(2*pi*x);
        
    case 4  % 多峰
        u = exp(-80*(x-0.25).^2) + 0.8*exp(-120*(x-0.7).^2);
        
    case 5  % 完全随机（不连续）
        u = rand(1, Nx);
        
    case 6  % 平滑随机（更物理）
        u = rand(1, Nx);
        u = smoothdata(u, 'gaussian', 15);  % 平滑
        u = u / max(u);                     % 归一化
        
end

u0 = u;

%% ---------------- 存储 ----------------
U = zeros(Nt+1, Nx);
U(1,:) = u;

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

subplot(1,2,1)
h_num   = plot(x, u, 'b-', 'LineWidth',2); hold on;
h_exact = plot(x, u, 'r--', 'LineWidth',1.5);
legend('数值解','精确解');
xlabel('x'); ylabel('u');
grid on;

ylim_auto = [min(u)-0.2, max(u)+0.2];
axis([0 L ylim_auto]);

htitle = title('初始化');

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

%% ---------------- 时间推进 ----------------
for n = 1:Nt
    
    % ========= Upwind =========
    u_new = u;
    u_new(2:Nx) = u(2:Nx) - CFL * (u(2:Nx) - u(1:Nx-1));
    u_new(1) = u_new(Nx);   % 周期边界
    
    u = u_new;
    
    % ========= 存储 =========
    U(n+1,:) = u;
    
    % ========= 精确解 =========
    shift = x - c*n*dt;
    shift = shift - floor(shift/L)*L;
    
    % 👉 非随机才有严格精确解
    if IC_type <= 4
        switch IC_type
            case 1
                u_exact = exp(-100 * (shift - 0.3).^2);
            case 2
                u_exact = zeros(size(x));
                u_exact(shift > 0.2 & shift < 0.5) = 1;
            case 3
                u_exact = sin(2*pi*shift);
            case 4
                u_exact = exp(-80*(shift-0.25).^2) + ...
                          0.8*exp(-120*(shift-0.7).^2);
        end
    else
        u_exact = nan(size(x)); % 随机无解析解
    end
    
    % ========= 可视化 =========
    if mod(n,5)==0
        
        subplot(1,2,1)
        h_num.YData = u;
        h_exact.YData = u_exact;
        
        htitle.String = ['t = ', num2str(n*dt,'%.2f'), ...
                         ' | Step = ', num2str(n)];
        
        subplot(1,2,2)
        set(h_img, 'CData', U);
        
        drawnow limitrate
        pause(0.01)
    end
end