%% =========================================================
% 1D Euler + 特征线
%% =========================================================
clear; clc; close all;

%% =========================
% 基本参数
%% =========================
gamma = 1.4;
nx = 800;
x = linspace(0,1,nx);
dx = x(2) - x(1);

CFL = 0.5;
tEnd = 0.25;
t = 0;

%% =========================
% 初始条件（Sod 激波管）
%% =========================
rho = ones(1,nx);
u   = zeros(1,nx);
p   = ones(1,nx);

rho(x < 0.5) = 1.0;
p(x < 0.5)   = 1.0;

rho(x >= 0.5) = 0.125;
p(x >= 0.5)   = 0.1;

E = p/(gamma-1) + 0.5*rho.*u.^2;

U = [rho;
     rho.*u;
     E];

%% =========================
% 图像初始化（两个subplot）
%% =========================
figure('Color','w','Position',[100 100 1200 600]);

% --- 左：流场 ---
subplot(1,2,1);
hold on;
h1 = plot(x, rho, 'LineWidth',2);
xlabel('x'); ylabel('\rho');
title('Flow field');
grid on;

% --- 右：特征线 ---
subplot(1,2,2);
hold on;
xlabel('x'); ylabel('t');
title('Characteristic lines');
grid on;

%% =========================
% 初始化多条特征线
%% =========================
nChar = 15;
x0 = linspace(0.3,0.7,nChar);

x_char1 = x0;
x_char2 = x0;
x_char3 = x0;

t_hist = [];
x1_hist = [];
x2_hist = [];
x3_hist = [];

%% =========================
% 时间推进
%% =========================
iter = 0;

while t < tEnd
    
    rho = U(1,:);
    u   = U(2,:) ./ rho;
    E   = U(3,:);
    p   = (gamma-1)*(E - 0.5*rho.*u.^2);
    
    a = sqrt(gamma*p./rho);
    
    maxSpeed = max(abs(u) + a);
    dt = CFL * dx / maxSpeed;
    
    if t + dt > tEnd
        dt = tEnd - t;
    end
    
    %% -------------------------
    % 通量
    %% -------------------------
    F = [rho.*u;
         rho.*u.^2 + p;
         u.*(E + p)];
    
    %% -------------------------
    % Lax-Friedrichs
    %% -------------------------
    U_new = U;
    
    for i = 2:nx-1
        U_new(:,i) = 0.5*(U(:,i+1) + U(:,i-1)) ...
                   - dt/(2*dx)*(F(:,i+1) - F(:,i-1));
    end
    
    U_new(:,1)   = U_new(:,2);
    U_new(:,end) = U_new(:,end-1);
    
    U = U_new;
    
    %% =====================================================
    % ⭐ 特征线推进
    %% =====================================================
    rho = U(1,:);
    u   = U(2,:) ./ rho;
    E   = U(3,:);
    p   = (gamma-1)*(E - 0.5*rho.*u.^2);
    a   = sqrt(gamma*p./rho);
    
    for k = 1:nChar
        
        uk = interp1(x,u,x_char1(k),'linear','extrap');
        ak = interp1(x,a,x_char1(k),'linear','extrap');
        x_char1(k) = x_char1(k) + (uk - ak)*dt;
        
        uk = interp1(x,u,x_char2(k),'linear','extrap');
        x_char2(k) = x_char2(k) + uk*dt;
        
        uk = interp1(x,u,x_char3(k),'linear','extrap');
        ak = interp1(x,a,x_char3(k),'linear','extrap');
        x_char3(k) = x_char3(k) + (uk + ak)*dt;
        
    end
    
    % 存储轨迹
    t_hist(end+1) = t;
    x1_hist(:,end+1) = x_char1';
    x2_hist(:,end+1) = x_char2';
    x3_hist(:,end+1) = x_char3';
    
    %% =====================================================
    % 激波检测
    %% =====================================================
    drho_dx = abs(diff(rho)/dx);
    [~,shock_idx] = max(drho_dx);
    x_shock = x(shock_idx);
    
    %% 更新时间
    t = t + dt;
    iter = iter + 1;
    
    %% =====================================================
    % 更新显示
    %% =====================================================
    if mod(iter,2)==0
        
        % ===== 左图：流场 =====
        subplot(1,2,1);
        set(h1,'YData',rho);
        title(['Flow field, t = ', num2str(t,'%.4f')]);
        
        % ===== 右图：特征线 =====
        subplot(1,2,2);
        cla; hold on;
        
        for k = 1:nChar
            plot(x1_hist(k,:), t_hist, 'r-');
            plot(x2_hist(k,:), t_hist, 'k--');
            plot(x3_hist(k,:), t_hist, 'b-');
        end
        
        % 激波位置
        plot(x_shock, t, 'ko', 'MarkerFaceColor','y','MarkerSize',6);
        
        xlabel('x'); ylabel('t');
        title('Characteristic convergence → shock');
        grid on;
        xlim([0 1])
        ylim([0 tEnd])
        drawnow limitrate;
        pause(0.01)
    end
    
end