clear; clc; close all;

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

c = 1.0;
CFL = 0.8;     % 👉 改这个看稳定性
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',[100 100 1400 800]);

for n = 1:Nt
    
    %% ===== Lax推进 =====
    u_new = u;
    for i = 2:Nx-1
        u_new(i) = 0.5*(u(i+1)+u(i-1)) ...
                   - 0.5*nu*(u(i+1)-u(i-1));
    end
    
    % 周期边界
    u_new(1) = u_new(end-1);
    u_new(end) = u_new(2);
    
    u = u_new;
    
    %% =========================
    % 1️⃣ 物理空间
    %% =========================
    subplot(2,2,1)
    plot(x, u, 'b','LineWidth',1.5);
    title(['Physical Space (Lax, t = ', num2str(n*dt,'%.2f'), ')']);
    axis([0 1 -2 2]); grid on;
    
    %% =========================
    % 2️⃣ 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;
    
    %% =========================
    % 3️⃣ 放大因子 |G|
    %% =========================
    G = cos(theta) - 1i*nu*sin(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| (Lax, CFL=',num2str(CFL),')']);
    xlabel('\theta'); ylabel('|G|');
    grid on;
    
    %% =========================
    % 4️⃣ 复平面稳定区
    %% =========================
    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