%% =========================================================
%  纯亚声速准一维喷管流
%  MacCormack 显式预测-校正格式（非守恒形式）
%
%  修正版要点：
%  1. 入口：亚声速入口，固定总温 T0=1、总压 p0=1，速度外推
%  2. 出口：亚声速出口，指定出口压力 p_exit
%  3. 解析解：基于 p0=1、T0=1 的等熵关系
%  4. 与纯亚声速等熵解析解对比
%
%  无量纲状态方程：
%       p = rho*T
%
%  入口边界核心：
%       V_1 = 2*V_2 - V_3
%       T_1 = T0 - (gamma-1)/2 * V_1^2
%       rho_1 = T_1^(1/(gamma-1))
%
%  出口边界核心：
%       T_N   = 2*T_{N-1} - T_{N-2}
%       V_N   = 2*V_{N-1} - V_{N-2}
%       rho_N = p_exit / T_N
%% =========================================================

clear; clc; close all;

%% ===================== 1. 基本参数 =====================
gamma = 1.4;               % 比热比
CFL   = 0.5;               % CFL数，建议先用0.5
xL    = 0.0;
xR    = 3.0;
N     = 31;

maxIter  = 120000;         % 纯亚声速收敛较慢，迭代步数给大一些
tol      = 1e-7;
plotEvery = 200;

% 总温、总压，无量纲取 1
T0 = 1.0;
p0 = 1.0;

% ---------------------------------------------------------
% 出口压力边界
% 纯亚声速稳定算例：p_exit = 0.93
% ---------------------------------------------------------
p_exit = 0.93;

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

% ---------------------------------------------------------
% 纯亚声速算例的非对称收缩-扩张喷管
% 左侧收缩段：A/At = 1+2.2(x-1.5)^2
% 右侧扩张段：A/At = 1+0.2223(x-1.5)^2
% ---------------------------------------------------------
A = zeros(N,1);

for i = 1:N
    if x(i) <= 1.5
        A(i) = 1+2.2*(x(i)-1.5)^2;
    else
        A(i) = 1+0.2223*(x(i)-1.5)^2;
    end
end

lnA = log(A);

[At, ithroat] = min(A);

%% ===================== 3. 初始条件 =====================
% 课件中的纯亚声速初始条件
rho = 1.0-0.023*x;
T   = 1.0-0.009333*x;
V   = 0.05+0.11*x;

% 防止初始条件出现非物理值
rho = max(rho,1e-8);
T   = max(T,1e-8);

% ---------------------------------------------------------
% 对初始场也施加一次边界条件
% 入口改为总压、总温边界，而不是固定静态 rho 和 T
% ---------------------------------------------------------
V(1) = 2*V(2)-V(3);
T(1) = T0-(gamma-1)/2*V(1)^2;
T(1) = max(T(1),1e-8);
rho(1) = T(1)^(1/(gamma-1));

% 出口压力边界
T(end) = 2*T(end-1)-T(end-2);
V(end) = 2*V(end-1)-V(end-2);
T(end) = max(T(end),1e-8);
rho(end) = p_exit/T(end);

p = rho.*T;
M = V./sqrt(T);

%% ===================== 4. 纯亚声速解析解 =====================
% 思路：
% 出口压力 p_exit 已知，由等熵关系先求出口马赫数 M_e
%
%   p/p0 = [1+(gamma-1)/2*M^2]^(-gamma/(gamma-1))
%
% 然后由出口面积 Ae 和 M_e 反推出 A*
%
%   Ae/A* = F(M_e)
%
% 最后对每个 x，根据 A/A* 求亚声速分支 M(x)

Ae = A(end);

M_exit_exact = sqrt( ...
    2/(gamma-1)*((p_exit/p0)^(-(gamma-1)/gamma)-1) );

F_exit = areaMachFunction(M_exit_exact,gamma);
Astar  = Ae/F_exit;

M_exact = zeros(N,1);

for i = 1:N
    areaRatio = A(i)/Astar;
    M_exact(i) = solveAreaMach(areaRatio,gamma,'subsonic');
end

T_exact   = T0 ./ (1+(gamma-1)/2*M_exact.^2);
p_exact   = p0 .* (T_exact/T0).^(gamma/(gamma-1));
rho_exact = p_exact ./ T_exact;
V_exact   = M_exact.*sqrt(T_exact);

fprintf('\n======== 解析解信息 ========\n');
fprintf('出口压力 p_exit       = %.5f\n',p_exit);
fprintf('解析出口马赫数 M_e    = %.5f\n',M_exit_exact);
fprintf('解析喉部马赫数 M_t    = %.5f\n',M_exact(ithroat));
fprintf('几何喉部面积 A_t      = %.5f\n',At);
fprintf('等效临界面积 A_star   = %.5f\n',Astar);
fprintf('A_star/A_t            = %.5f\n',Astar/At);

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

massFlowError = zeros(maxIter,1);
entropyError  = zeros(maxIter,1);

physicalTime = 0.0;

%% ===================== 6. 动态图形窗口 =====================
fig = figure('Color','w','Position',[100 60 1400 820]);
htitle = sgtitle('Purely Subsonic Nozzle Flow', 'FontWeight','bold');

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

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

yyaxis right
hM = plot(x,M,'b-','LineWidth',2); hold on;
plot(x,M_exact,'r--','LineWidth',1.2);
yline(1.0,'k-.','LineWidth',1.0);
ylabel('Mach number M');

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

% ---------- 密度 ----------
subplot(2,3,2);
hRho = plot(x,rho,'b-','LineWidth',2); hold on;
plot(x,rho_exact,'r--','LineWidth',1.2);
xlabel('x');
ylabel('\rho');
title('Density');
legend('Numerical','Exact','Location','best');
grid on;

% ---------- 温度 ----------
subplot(2,3,3);
hT = plot(x,T,'b-','LineWidth',2); hold on;
plot(x,T_exact,'r--','LineWidth',1.2);
xlabel('x');
ylabel('T');
title('Temperature');
legend('Numerical','Exact','Location','best');
grid on;

% ---------- 速度 ----------
subplot(2,3,4);
hV = plot(x,V,'b-','LineWidth',2); hold on;
plot(x,V_exact,'r--','LineWidth',1.2);
xlabel('x');
ylabel('V');
title('Velocity');
legend('Numerical','Exact','Location','best');
grid on;

% ---------- 压力 ----------
subplot(2,3,5);
hP = plot(x,p,'b-','LineWidth',2); hold on;
plot(x,p_exact,'r--','LineWidth',1.2);
yline(p_exit,'k-.','LineWidth',1.0);
xlabel('x');
ylabel('p');
title('Pressure');
legend('Numerical','Exact','p_{exit}','Location','best');
grid on;

% ---------- 残差 ----------
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');
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+abs(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 预测值边界条件 ----------
    % 入口：亚声速入口，固定总温 T0 和总压 p0，V 外推
    V_bar(1) = 2*V_bar(2)-V_bar(3);

    T_bar(1) = T0-(gamma-1)/2*V_bar(1)^2;
    T_bar(1) = max(T_bar(1),1e-8);

    rho_bar(1) = T_bar(1)^(1/(gamma-1));

    % 出口：亚声速压力出口
    % 方法一：先外推 T，再由 p=rho*T 求 rho
    T_bar(N)   = 2*T_bar(N-1)-T_bar(N-2);
    V_bar(N)   = 2*V_bar(N-1)-V_bar(N-2);
    T_bar(N)   = max(T_bar(N),1e-8);
    rho_bar(N) = p_exit/T_bar(N);

    % 防止负值
    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 校正后边界条件 ----------
    % 入口：固定总温、总压，速度外推
    V(1) = 2*V(2)-V(3);

    T(1) = T0-(gamma-1)/2*V(1)^2;
    T(1) = max(T(1),1e-8);

    rho(1) = T(1)^(1/(gamma-1));

    % 出口：亚声速压力出口
    T(N)   = 2*T(N-1)-T(N-2);
    V(N)   = 2*V(N-1)-V(N-2);
    T(N)   = max(T(N),1e-8);
    rho(N) = p_exit/T(N);

    % 防止非物理值
    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 质量流量与等熵性诊断 ----------
    massFlow = rho.*V.*A;
    entropyLike = p./rho.^gamma;

    massFlowError(iter) = (max(massFlow)-min(massFlow))/mean(abs(massFlow));
    entropyError(iter)  = (max(entropyLike)-min(entropyLike))/mean(abs(entropyLike));

    % ---------- 7.10 非物理检测 ----------
    if any(~isfinite(rho)) || any(~isfinite(V)) || any(~isfinite(T)) ...
            || any(rho<=0) || any(T<=0)

        warning('计算发散：iter = %d',iter);
        break;
    end

    % ---------- 7.11 动态绘图 ----------
    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( ...
            'Purely Subsonic Nozzle Flow | Iter = %d | t = %.4f | max(M)=%.4f', ...
            iter,physicalTime,max(M)));

        drawnow limitrate;

        if iter==1
            waitforbuttonpress;   % 课堂演示时暂停看初始场
        end

        pause(0.01);
    end

    % ---------- 7.12 收敛判断 ----------
    if res_max < tol
        fprintf('\n计算收敛！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);
massFlowError = massFlowError(1:finalIter);
entropyError  = entropyError(1:finalIter);

fprintf('\n======== 计算结束 ========\n');
fprintf('总迭代步数              = %d\n',finalIter);
fprintf('最终物理时间            = %.6f\n',physicalTime);
fprintf('给定出口压力            = %.6f\n',p_exit);
fprintf('数值出口压力            = %.6f\n',p(end));
fprintf('数值出口马赫数          = %.6f\n',M(end));
fprintf('解析出口马赫数          = %.6f\n',M_exact(end));
fprintf('数值喉部马赫数          = %.6f\n',M(ithroat));
fprintf('解析喉部马赫数          = %.6f\n',M_exact(ithroat));
fprintf('最终最大残差            = %.3e\n',max([res_rho(end),res_V(end),res_T(end)]));
fprintf('最大马赫数 max(M)       = %.6f\n',max(M));
fprintf('质量流量相对波动        = %.3e\n',massFlowError(end));
fprintf('p/rho^gamma 相对波动    = %.3e\n',entropyError(end));

%% ===================== 9. 最终一维结果对比 =====================
figure('Color','w','Position',[120 80 1300 760]);

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

% ---------- Mach ----------
nexttile;
yyaxis left
plot(x,A,'k:','LineWidth',1.5);
ylabel('A');

yyaxis right
plot(x,M,'b-','LineWidth',2); hold on;
plot(x,M_exact,'r--','LineWidth',1.2);
yline(1.0,'k-.','LineWidth',1.0);
xlabel('x');
ylabel('M');
title('Mach Number');
legend('A','Numerical','Exact','M=1','Location','best');
grid on;

% ---------- rho ----------
nexttile;
plot(x,rho,'b-','LineWidth',2); hold on;
plot(x,rho_exact,'r--','LineWidth',1.2);
xlabel('x');
ylabel('\rho');
title('Density');
legend('Numerical','Exact','Location','best');
grid on;

% ---------- T ----------
nexttile;
plot(x,T,'b-','LineWidth',2); hold on;
plot(x,T_exact,'r--','LineWidth',1.2);
xlabel('x');
ylabel('T');
title('Temperature');
legend('Numerical','Exact','Location','best');
grid on;

% ---------- V ----------
nexttile;
plot(x,V,'b-','LineWidth',2); hold on;
plot(x,V_exact,'r--','LineWidth',1.2);
xlabel('x');
ylabel('V');
title('Velocity');
legend('Numerical','Exact','Location','best');
grid on;

% ---------- p ----------
nexttile;
plot(x,p,'b-','LineWidth',2); hold on;
plot(x,p_exact,'r--','LineWidth',1.2);
yline(p_exit,'k-.','LineWidth',1.0);
xlabel('x');
ylabel('p');
title('Pressure');
legend('Numerical','Exact','p_{exit}','Location','best');
grid on;

% ---------- Residual ----------
nexttile;
semilogy(1:finalIter,res_rho,'LineWidth',1.4); hold on;
semilogy(1:finalIter,res_V,'LineWidth',1.4);
semilogy(1:finalIter,res_T,'LineWidth',1.4);
xlabel('Iteration');
ylabel('Residual');
title('Residual History');
legend('\rho','V','T','Location','best');
grid on;

sgtitle('Purely Subsonic Quasi-1D Nozzle Flow', ...
    'FontWeight','bold','FontSize',16);

%% ===================== 10. 守恒性与等熵性诊断 =====================
figure('Color','w','Position',[160 120 1100 420]);

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

nexttile;
massFlow = rho.*V.*A;
plot(x,massFlow,'b-','LineWidth',2);
xlabel('x');
ylabel('\rho V A');
title('Mass Flow Rate');
grid on;

nexttile;
entropyLike = p./rho.^gamma;
plot(x,entropyLike,'b-','LineWidth',2);
xlabel('x');
ylabel('p/\rho^\gamma');
title('Isentropic Indicator');
grid on;

sgtitle('Steady-State Diagnostics', ...
    'FontWeight','bold','FontSize',15);

%% ===================== 11. 喷管内二维 contour 展示 =====================
% 准一维结果沿横向复制，仅用于课堂可视化展示

Ny = 101;
yScale = 0.35;
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);

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;

figure('Color','w','Position',[80 80 1400 720]);

tl3 = tiledlayout(2,2);
tl3.TileSpacing = 'compact';
tl3.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';
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';
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';
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(turbo);

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

sgtitle('Purely Subsonic Nozzle Flow: Final Contour Visualization', ...
    'FontWeight','bold','FontSize',16);

%% =========================================================
%  局部函数 1：面积-马赫数函数 A/A*
%% =========================================================
function F = areaMachFunction(M,gamma)

F = (1./M).* ...
    ( ...
    (2/(gamma+1)) .* ...
    (1+(gamma-1)/2.*M.^2) ...
    ).^((gamma+1)/(2*(gamma-1)));

end

%% =========================================================
%  局部函数 2：由 A/A* 反求 M
%% =========================================================
function M = solveAreaMach(areaRatio,gamma,branch)

fun = @(M) areaMachFunction(M,gamma)-areaRatio;

switch lower(branch)

    case 'subsonic'
        % 亚声速分支：0 < M < 1
        M_left  = 1e-8;
        M_right = 1-1e-8;

    case 'supersonic'
        % 超声速分支：M > 1
        M_left  = 1+1e-8;
        M_right = 20;

    otherwise
        error('branch 必须为 subsonic 或 supersonic');
end

M = fzero(fun,[M_left,M_right]);

end