clear; clc; close all;

% 参数
L = 1.0;
Nx = 256;
dx = L / Nx;

c = 1.0;
CFL = 1.0;     % 👉 稳定性
nu = CFL;

dt = CFL * dx / c;

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

x = linspace(0, L, Nx);

% 初始条件（多频叠加）
u0 = sin(2*pi*5*x) + 0.5*sin(2*pi*20*x) + 0.2*sin(2*pi*60*x);

u = u0;

% 频率空间
k = (0:Nx-1);
theta = 2*pi*k/Nx;

% 画图窗口
figure('Color','w','Position',[10 10 1400 800]);

for n = 1:Nt
    
    % ===== Upwind推进 =====
    u_new = u;
    for i = 2:Nx
        u_new(i) = u(i) - nu*(u(i)-u(i-1));
    end
    u_new(1) = u_new(end); % 周期
    
    u = u_new;
    
    % 物理空间
    subplot(2,2,1)
    plot(x, u, 'b','LineWidth',1.5);
    title(['Physical Space (t = ', num2str(n*dt,'%.2f'), ')']);
    axis([0 1 -2 2]); grid on;
    
    % FFT 频谱
    U = fft(u);
    spectrum = abs(U)/Nx;
    
    subplot(2,2,2)
    stem(k(1:Nx/2), spectrum(1:Nx/2),'filled');
    title('Fourier Spectrum');
    xlabel('Wavenumber'); ylabel('Amplitude');
    grid on;
    
    % 放大因子 |G|
    G = 1 - nu*(1 - exp(-1i*theta));
    
    subplot(2,2,3)
    plot(theta, abs(G),'r','LineWidth',2); hold on;
    yline(1,'k--','|G|=1');
    hold off;
    
    title(['Amplification Factor |G| (CFL=',num2str(CFL),')']);
    xlabel('\theta'); ylabel('|G|');
    grid on;
    
    % 复平面稳定区
    subplot(2,2,4)
    
    % 单位圆（稳定边界）
    t_circle = linspace(0,2*pi,200);
    plot(cos(t_circle), sin(t_circle),'k--'); hold on;
    
    % G轨迹
    plot(real(G), imag(G),'r','LineWidth',1.5);
    
    hold off;
    axis equal;
    xlim([-1.5 1.5]); ylim([-1.5 1.5]);
    
    title('Stability Region in Complex Plane');
    xlabel('Re(G)'); ylabel('Im(G)');
    grid on;
    
    drawnow;
end