close all; clear; clc;

% 准一维喷管等熵流动
% MacCormack 显式预测-校正格式（非守恒形式）
%
% 控制方程（无量纲、非守恒形式）：
% drho/dt + V*drho/dx + rho*dV/dx + rho*V*d(lnA)/dx = 0
% dV/dt   + V*dV/dx = -(1/gamma)*(dT/dx + (T/rho)*drho/dx)
% dT/dt   + V*dT/dx = -(gamma-1)*T*(dV/dx + V*d(lnA)/dx)
%
% 主要功能：
% 1) 展示从初始场到稳态场的推进过程
% 2) 展示残差收敛
% 3) 与解析等熵解进行对比
%
% 边界条件：
% 入口（亚声）：rho=1, T=1, V外推
% 出口（超声）：rho, T, V全部外推


% 一：只看现象。
% 观察马赫数、压力、密度、温度如何从初始场逐渐靠近解析解。
% 
% 二：主循环。
% 7.1 到 7.6：CFL 时间步、预测步、边界条件、校正步、平均更新。
% 
% 三：公式对应。
% 把连续方程、动量方程、能量方程分别和 drho_dt、dV_dt、dT_dt 对应起来。
% 
% 四：讨论网格收敛性、CFL数。
% ------------------------------------------------------------

clear; clc; close all;

% ===================== 1. 基本参数 =====================
gamma = 1.4;          % 比热比
CFL   = 0.5;          % CFL 数 0.5/1.0/1.1/1.2
xL    = 0.0;          % 计算域左端
xR    = 3.0;          % 计算域右端
N     = 101;           % 网格点数（可改为 41, 61, 81 观察差异）
maxIter = 5000;       % 最大迭代步数
tol     = 1e-6;       % 收敛阈值（残差）
plotEvery = 10;       % 每隔多少步刷新一次图像

% ===================== 2. 网格与喷管几何 ================
x  = linspace(xL, xR, N)';
dx = x(2) - x(1);

% 喷管面积分布：A = 1 + 2.2 (x-1.5)^2
A = 1 + 2.2*(x - 1.5).^2;
lnA = log(A);

% 喉部位置与面积
[Ast, ithroat] = min(A); 

% ===================== 3. 初始条件 =====================
% 这里与课件中的"聪明初始化"一致，同学们可自行更改
rho = 1 - 0.3146*x;
T   = 1 - 0.2314*x;
V   = (0.1 + 1.09*x).*sqrt(T);

% 可由状态方程得到压力
p = rho .* T;

% ===================== 4. 解析解 =======================
% 对于该喷管，A* 是喉部面积，设计工况为喉部 M=1 的跨声速等熵解
M_exact = zeros(N,1);
for i = 1:N
    areaRatio = A(i)/min(A);
    if i < ithroat
        % 喉部左侧：亚声支
        M_exact(i) = solveAreaMach(areaRatio, gamma, 'subsonic');
    elseif i > ithroat
        % 喉部右侧：超声支
        M_exact(i) = solveAreaMach(areaRatio, gamma, 'supersonic');
    elseif i == ithroat
        % 喉部：M=1
        M_exact(i) = 1;
    end
end

% 等熵关系获得精确解：注意无量纲体系
T_exact   = 1 ./ (1 + (gamma-1)/2 .* M_exact.^2);
rho_exact = T_exact.^(1/(gamma-1));
p_exact   = T_exact.^(gamma/(gamma-1));
V_exact   = M_exact .* sqrt(T_exact);   % 无量纲声速 a = sqrt(T)


% ===================== 5. 预分配残差与时间记录 ==============
res_rho = zeros(maxIter,1);
res_V   = zeros(maxIter,1);
res_T   = zeros(maxIter,1);
time_history = zeros(maxIter,1);

physicalTime = 0.0;

% ===================== 6. 建立图形窗口 =====================
fig = figure('Color','w','Position',[10 10 1400 820]);

htitle = sgtitle('初始化', 'FontWeight','bold');

% (1) 马赫数 + 面积 
subplot(2,3,1); 

yyaxis left
hA = plot(x, A, 'k:', 'LineWidth', 1.5); hold on;
ylabel('Area A(x)')

yyaxis right
hM = plot(x, zeros(size(x)), 'b-', 'LineWidth', 2); hold on;
hM_exact = plot(x, M_exact, 'r--', 'LineWidth', 1.2);
ylabel('Mach number M')

xlabel('x')
title('Area & Mach Number')
legend('A(x)', 'M (numerical)', 'M (exact)', 'Location','best')
grid on;


% (2) 密度
subplot(2,3,2);

hRho = plot(x, rho, 'b-', 'LineWidth', 2); hold on;
hRho_exact = plot(x, rho_exact, 'r--');

xlabel('x')
ylabel('\rho / \rho_0')
title('Density Distribution')

legend('Numerical', 'Exact', 'Location','best')
grid on;


% (3) 温度 
subplot(2,3,3);

hT = plot(x, T, 'b-', 'LineWidth', 2); hold on;
hT_exact = plot(x, T_exact, 'r--');

xlabel('x')
ylabel('T / T_0')
title('Temperature Distribution')

legend('Numerical', 'Exact', 'Location','best')
grid on;


% (4) 速度
subplot(2,3,4);

hV = plot(x, V, 'b-', 'LineWidth', 2); hold on;
hV_exact = plot(x, V_exact, 'r--');

xlabel('x')
ylabel('Velocity V / a_0')
title('Velocity Distribution')

legend('Numerical', 'Exact', 'Location','best')
grid on;


% (5) 压力
subplot(2,3,5);

hP = plot(x, p, 'b-', 'LineWidth', 2); hold on;
hP_exact = plot(x, p_exact, 'r--');

xlabel('x')
ylabel('p / p_0')
title('Pressure Distribution')

legend('Numerical', 'Exact', 'Location','best')
grid on;


% (6) 残差
subplot(2,3,6);

hRes1 = semilogy(1,1,'LineWidth',1.5); hold on;
hRes2 = semilogy(1,1,'LineWidth',1.5);
hRes3 = semilogy(1,1,'LineWidth',1.5);

xlabel('Iteration')
ylabel('Residual (log scale)')
title('Convergence History')

legend('\rho','V','T','Location','best')
grid on;

% ===================== 7. 主循环：MacCormack推进 ===========
for iter = 1:maxIter

    % ---------- 7.1 根据 CFL 计算全局时间步 ----------
    a  = sqrt(T);  % 无量纲局部声速
    dt_local = CFL * dx ./ (a + V);
    dt = min(dt_local);
    physicalTime = physicalTime + dt;
    time_history(iter) = physicalTime;

    % 保存旧值，用于计算残差
    rho_old = rho;
    V_old   = V;
    T_old   = T;

    % ---------- 7.2 预测步（forward difference） ----------
    drho_dt = zeros(N,1);
    dV_dt   = zeros(N,1);
    dT_dt   = zeros(N,1);

    for i = 2:N-1
        dlnAdx_f = (lnA(i+1) - lnA(i)) / dx;
        dVdx_f   = (V(i+1)   - V(i))   / dx;
        drhodx_f = (rho(i+1) - rho(i)) / dx;
        dTdx_f   = (T(i+1)   - T(i))   / dx;

        drho_dt(i) = - rho(i)*dVdx_f ...
                     - rho(i)*V(i)*dlnAdx_f ...
                     - V(i)*drhodx_f;

        dV_dt(i)   = - V(i)*dVdx_f ...
                     - (1/gamma)*(dTdx_f + (T(i)/rho(i))*drhodx_f);

        dT_dt(i)   = - V(i)*dTdx_f ...
                     - (gamma-1)*T(i)*(dVdx_f + V(i)*dlnAdx_f);
    end

    % 预测值
    rho_bar = rho + drho_dt*dt;
    V_bar   = V   + dV_dt*dt;
    T_bar   = T   + dT_dt*dt;

    % ---------- 7.3 对预测值施加边界条件 ----------
    % 入口：rho=1, T=1, V外推
    rho_bar(1) = 1.0;
    T_bar(1)   = 1.0;
    V_bar(1)   = 2*V_bar(2) - V_bar(3);

    % 出口：全部线性外推（超声出口）
    rho_bar(N) = 2*rho_bar(N-1) - rho_bar(N-2);
    T_bar(N)   = 2*T_bar(N-1)   - T_bar(N-2);
    V_bar(N)   = 2*V_bar(N-1)   - V_bar(N-2);

    % 防止数值误差导致负值
    rho_bar = max(rho_bar, 1e-8);
    T_bar   = max(T_bar,   1e-8);

    % ---------- 7.4 校正步（backward difference） ----------
    drho_bar_dt = zeros(N,1);
    dV_bar_dt   = zeros(N,1);
    dT_bar_dt   = zeros(N,1);

    for i = 2:N-1
        dlnAdx_b = (lnA(i)     - lnA(i-1))     / dx;
        dVdx_b   = (V_bar(i)   - V_bar(i-1))   / dx;
        drhodx_b = (rho_bar(i) - rho_bar(i-1)) / dx;
        dTdx_b   = (T_bar(i)   - T_bar(i-1))   / dx;

        drho_bar_dt(i) = - rho_bar(i)*dVdx_b ...
                         - rho_bar(i)*V_bar(i)*dlnAdx_b ...
                         - V_bar(i)*drhodx_b;

        dV_bar_dt(i)   = - V_bar(i)*dVdx_b ...
                         - (1/gamma)*(dTdx_b + (T_bar(i)/rho_bar(i))*drhodx_b);

        dT_bar_dt(i)   = - V_bar(i)*dTdx_b ...
                         - (gamma-1)*T_bar(i)*(dVdx_b + V_bar(i)*dlnAdx_b);
    end

    % ---------- 7.5 平均导数并更新 ----------
    drho_dt_av = 0.5*(drho_dt + drho_bar_dt);
    dV_dt_av   = 0.5*(dV_dt   + dV_bar_dt);
    dT_dt_av   = 0.5*(dT_dt   + dT_bar_dt);

    rho(2:N-1) = rho_old(2:N-1) + drho_dt_av(2:N-1)*dt;
    V(2:N-1)   = V_old(2:N-1)   + dV_dt_av(2:N-1)*dt;
    T(2:N-1)   = T_old(2:N-1)   + dT_dt_av(2:N-1)*dt;

    % ---------- 7.6 对校正后的解施加边界条件 ----------
    % 入口：指定 rho, T，V 由内部外推
    rho(1) = 1.0;
    T(1)   = 1.0;
    V(1)   = 2*V(2) - V(3);

    % 出口：超声外推
    rho(N) = 2*rho(N-1) - rho(N-2);
    T(N)   = 2*T(N-1)   - T(N-2);
    V(N)   = 2*V(N-1)   - V(N-2);

    % 防止负值
    rho = max(rho, 1e-8);
    T   = max(T,   1e-8);

    % ---------- 7.7 计算辅助物理量 ----------
    p = rho .* T;
    M = V ./ sqrt(T);

    % ---------- 7.8 定义残差 ----------
    % 这里用每一步更新量的最大绝对值作为简洁残差
    res_rho(iter) = max(abs(rho - rho_old));
    res_V(iter)   = max(abs(V   - V_old));
    res_T(iter)   = max(abs(T   - T_old));

    res_max = max([res_rho(iter), res_V(iter), res_T(iter)]);

    % ---------- 7.9 动态画图 ----------
    if mod(iter, plotEvery) == 0 || iter == 1

        % 更新数据（核心！）
        set(hM, 'YData', M);
        set(hRho, 'YData', rho);
        set(hT, 'YData', T);
        set(hV, 'YData', V);
        set(hP, 'YData', p);
    
        % 更新残差
        set(hRes1, 'XData', 1:iter, 'YData', res_rho(1:iter));
        set(hRes2, 'XData', 1:iter, 'YData', res_V(1:iter));
        set(hRes3, 'XData', 1:iter, 'YData', res_T(1:iter));
    
        set(htitle, 'String', sprintf('Quasi-1D Nozzle Flow | Iter = %d | t = %.4f', iter, physicalTime));

        drawnow limitrate  
        if iter ==1
            waitforbuttonpress; % 这里用于课堂演示，初始条件
        end
        pause(0.01) % 这里用于课堂演示，实际计算中没有必要
    end

    % ---------- 7.10 收敛判据 ----------
    if res_max < tol
        fprintf('计算收敛！iter = %d, time = %.6f, max residual = %.3e\n', ...
            iter, physicalTime, res_max);
        break;
    end

end

% ===================== 8. 输出最终结果 =====================
finalIter = iter;
res_rho = res_rho(1:finalIter);
res_V   = res_V(1:finalIter);
res_T   = res_T(1:finalIter);
time_history = time_history(1:finalIter);

fprintf('\n======== 计算结束 ========\n');
fprintf('总迭代步数      = %d\n', finalIter);
fprintf('最终物理时间    = %.6f\n', physicalTime);
fprintf('出口马赫数      = %.4f\n', M(end));
fprintf('喉部马赫数      = %.4f\n', M(ithroat));
fprintf('最终最大残差    = %.3e\n', max([res_rho(end), res_V(end), res_T(end)]));

%% ===================== 9. 喷管内二维 contour 展示 =====================
% 说明：
% 准一维喷管计算得到的是沿 x 方向的截面平均量。
% 为了教学展示，这里将一维结果沿喷管横向复制，
% 构造一个二维可视化场，用于显示喷管内温度、压力、密度和马赫数分布。

Ny = 101;                         % 横向网格点数

% ------------------------------------------------------------
% 为了让图像长宽比更自然，这里只对可视化几何做纵向压缩
% 不改变一维计算结果，只改变显示效果
% ------------------------------------------------------------
yScale = 0.35;                     % 纵向显示缩放因子，可调为 0.25~0.5
h = yScale * A / 2;                % 喷管半高，用于绘图显示

eta = linspace(-1, 1, Ny);         % 归一化横向坐标
[X2D, ETA2D] = meshgrid(x, eta);   % 二维显示网格

H2D = repmat(h', Ny, 1);
Y2D = ETA2D .* H2D;

% 将一维准一维结果扩展为二维显示场
rho2D = repmat(rho', Ny, 1);
T2D   = repmat(T',   Ny, 1);
p2D   = repmat(p',   Ny, 1);
M2D   = repmat(M',   Ny, 1);

% 喷管外部区域置为 NaN
mask = abs(Y2D) > H2D;

rho2D(mask) = NaN;
T2D(mask)   = NaN;
p2D(mask)   = NaN;
M2D(mask)   = NaN;

% 喷管壁面
x_wall = x;
y_top  = h;
y_bot  = -h;

% ===================== 紧密 panel 布局 =====================
figure('Color','w','Position',[10 10 1400 720]);

tl = tiledlayout(2,2);
tl.TileSpacing = 'compact';
tl.Padding     = 'compact';

% 统一等值线数量
nLevel = 50;

% ---------- 温度 ----------
ax1 = nexttile;
contourf(X2D, Y2D, T2D, nLevel, 'LineColor', 'none');
hold on;
plot(x_wall, y_top, 'k-', 'LineWidth', 1.6);
plot(x_wall, y_bot, 'k-', 'LineWidth', 1.6);
xlabel('x');
ylabel('y');
title('Temperature');
cb1 = colorbar;
cb1.Label.String = 'T / T_0';
grid on;
box on;
axis tight;
pbaspect([3.2 1 1]);

% ---------- 压力 ----------
ax2 = nexttile;
contourf(X2D, Y2D, p2D, nLevel, 'LineColor', 'none');
hold on;
plot(x_wall, y_top, 'k-', 'LineWidth', 1.6);
plot(x_wall, y_bot, 'k-', 'LineWidth', 1.6);
xlabel('x');
ylabel('y');
title('Pressure');
cb2 = colorbar;
cb2.Label.String = 'p / p_0';
grid on;
box on;
axis tight;
pbaspect([3.2 1 1]);

% ---------- 密度 ----------
ax3 = nexttile;
contourf(X2D, Y2D, rho2D, nLevel, 'LineColor', 'none');
hold on;
plot(x_wall, y_top, 'k-', 'LineWidth', 1.6);
plot(x_wall, y_bot, 'k-', 'LineWidth', 1.6);
xlabel('x');
ylabel('y');
title('Density');
cb3 = colorbar;
cb3.Label.String = '\rho / \rho_0';
grid on;
box on;
axis tight;
pbaspect([3.2 1 1]);

% ---------- 马赫数 ----------
ax4 = nexttile;
contourf(X2D, Y2D, M2D, nLevel, 'LineColor', 'none');
hold on;
plot(x_wall, y_top, 'k-', 'LineWidth', 1.6);
plot(x_wall, y_bot, 'k-', 'LineWidth', 1.6);
xline(x(ithroat), 'k--', 'LineWidth', 1.2);
xlabel('x');
ylabel('y');
title('Mach Number');
cb4 = colorbar;
cb4.Label.String = 'M';
grid on;
box on;
axis tight;
pbaspect([3.2 1 1]);

% 统一 colormap 和字体
colormap(turbo);

set([ax1 ax2 ax3 ax4], ...
    'FontName','Times New Roman', ...
    'FontSize',12, ...
    'LineWidth',1.0, ...
    'Layer','top');

sgtitle('Quasi-1D Nozzle Flow: Final Contour Visualization', ...
    'FontWeight','bold', ...
    'FontSize',16);