%% =========================================================
%  纯亚声速准一维喷管流
%  MacCormack 显式预测-校正格式（守恒形式）
%
%  推进变量：
%       U1 = rho*A
%       U2 = rho*V*A
%       U3 = rho*A*( T/(gamma-1) + gamma/2*V^2 )
%
%  说明：
%  本程序保持原非守恒代码的无量纲体系：V 用 a0 归一化，M=V/sqrt(T)。
%  因此动量方程中的压力项含 1/gamma，能量变量中动能项为 gamma/2*V^2。
%
%  守恒形式：
%       dU/dt + dF/dx = J
%
%  边界条件：
%       入口（亚声）：固定总温 T0=1、总压 p0=1，V 外推
%       出口（亚声）：指定出口静压 p_exit，T 和 V 外推，rho=p_exit/T
%% =========================================================

clear; clc; close all;

%% ===================== 1. 基本参数 =====================
gamma = 1.4;
CFL   = 0.5;
xL    = 0.0;
xR    = 3.0;
N     = 31;
maxIter  = 120000;
tol      = 1e-7;
plotEvery = 200;

T0 = 1.0;
p0 = 1.0;
p_exit = 0.93;

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

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

dAdx = gradient(A,dx);
[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);

U = primitiveToConservative(rho,V,T,A,gamma);
[U,rho,V,T] = applySubsonicPressureOutletBC(U,A,gamma,T0,p0,p_exit);

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

%% ===================== 4. 纯亚声速解析解 =====================
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('Conservative 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

    [rho,V,T,p,M] = conservativeToPrimitive(U,A,gamma);

    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;
    U_old   = U;

    % ---------- 7.1 预测步：前向差分 ----------
    [F,J] = computeFluxSource(U,A,dAdx,gamma);
    U_bar = U;

    for i = 2:N-1
        U_bar(i,:) = U(i,:) - dt/dx*(F(i+1,:)-F(i,:)) + dt*J(i,:);
    end

    [U_bar,~,~,~] = applySubsonicPressureOutletBC(U_bar,A,gamma,T0,p0,p_exit);

    % ---------- 7.2 校正步：后向差分 ----------
    [F_bar,J_bar] = computeFluxSource(U_bar,A,dAdx,gamma);
    U_new = U_old;

    for i = 2:N-1
        U_new(i,:) = 0.5*( ...
            U_old(i,:) + U_bar(i,:) ...
            - dt/dx*(F_bar(i,:)-F_bar(i-1,:)) ...
            + dt*J_bar(i,:) );
    end

    [U,rho,V,T] = applySubsonicPressureOutletBC(U_new,A,gamma,T0,p0,p_exit);
    p = rho.*T;
    M = V./sqrt(T);

    % ---------- 7.3 残差 ----------
    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.4 诊断量 ----------
    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));

    if any(~isfinite(U(:))) || any(rho<=0) || any(T<=0)
        warning('计算发散：iter = %d',iter);
        break;
    end

    % ---------- 7.5 动态绘图 ----------
    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( ...
            'Conservative 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

    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';

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;

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;

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;

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;

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;

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('Conservative 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;

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,h,'k-','LineWidth',1.6); plot(x,-h,'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,h,'k-','LineWidth',1.6); plot(x,-h,'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,h,'k-','LineWidth',1.6); plot(x,-h,'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,h,'k-','LineWidth',1.6); plot(x,-h,'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('Conservative Purely Subsonic Nozzle Flow: Final Contour Visualization', ...
    'FontWeight','bold','FontSize',16);

%% =========================================================
%  局部函数
%% =========================================================
function U = primitiveToConservative(rho,V,T,A,gamma)
U = zeros(length(rho),3);
U(:,1) = rho.*A;
U(:,2) = rho.*V.*A;
U(:,3) = rho.*A.*(T/(gamma-1)+gamma/2*V.^2);
end

function [rho,V,T,p,M] = conservativeToPrimitive(U,A,gamma)
rho = U(:,1)./A;
V   = U(:,2)./U(:,1);
E   = U(:,3)./U(:,1);
T   = (gamma-1)*(E-gamma/2*V.^2);
rho = max(rho,1e-10);
T   = max(T,1e-10);
p   = rho.*T;
M   = V./sqrt(T);
end

function [F,J] = computeFluxSource(U,A,dAdx,gamma)
[rho,V,T,p,~] = conservativeToPrimitive(U,A,gamma);
F = zeros(size(U));
J = zeros(size(U));
F(:,1) = U(:,2);
F(:,2) = rho.*V.^2.*A + p.*A/gamma;
F(:,3) = U(:,2).*(U(:,3)./U(:,1)+T);
J(:,2) = p.*dAdx/gamma;
end

function [U,rho,V,T] = applySubsonicPressureOutletBC(U,A,gamma,T0,p0,p_exit)
[rho,V,T,~,~] = conservativeToPrimitive(U,A,gamma);
N = length(rho);

% 入口：总温、总压固定，速度外推
V(1) = 2*V(2)-V(3);
T(1) = T0-(gamma-1)/2*V(1)^2;
T(1) = max(T(1),1e-8);
p_in = p0*(T(1)/T0)^(gamma/(gamma-1));
rho(1) = p_in/T(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);
U = primitiveToConservative(rho,V,T,A,gamma);
end

function F = areaMachFunction(M,gamma)
F = (1./M).* ...
    ((2/(gamma+1)).*(1+(gamma-1)/2.*M.^2)).^((gamma+1)/(2*(gamma-1)));
end

function M = solveAreaMach(areaRatio,gamma,branch)
fun = @(M) areaMachFunction(M,gamma)-areaRatio;
switch lower(branch)
    case 'subsonic'
        M = fzero(fun,[1e-8,1-1e-8]);
    case 'supersonic'
        M = fzero(fun,[1+1e-8,20]);
    otherwise
        error('branch 必须为 subsonic 或 supersonic');
end
end