%% ============================================================
%  Shock-Capturing in a Quasi-One-Dimensional Nozzle
%
%  功能：
%  1. 使用 Anderson 教材中的准一维守恒变量
%  2. MacCormack 显式预测-校正格式
%  3. 加入压力传感器人工粘性
%  4. 比较无人工粘性与有人工粘性
%  5. 加入含正激波的准一维解析解对比
%  6. 最后给出二维 contour 可视化
%
%  注意：
%  二维 contour 只是把准一维结果沿横向复制，用于课堂展示；
%  不是二维 Euler 方程求解结果。
%% ============================================================

clear; clc; close all;

%% ---------------- 基本参数 ----------------
gamma = 1.4;

Nx = 61;
x  = linspace(0,3,Nx);
dx = x(2)-x(1);

A = 1 + 2.2*(x-1.5).^2;
dAdx = gradient(A,dx);

[At,ithroat] = min(A);
x_throat = x(ithroat);

CFL = 0.50;

% 指定出口静压
p_exit = 0.6784;

% 迭代步数
Nt_noAV = 800;      % 无人工粘性：只做短时间演示
Nt_AV   = 3000;     % 有人工粘性：用于收敛计算

Cx_noAV = 0.0;
Cx_AV   = 0.20;

doAnimate = true;

%% ---------------- 解析解 ----------------
% 含内激波的准一维解析解：
% 喉部前等熵亚声速；
% 喉部后激波前等熵超声速；
% 激波后总压损失，再按新的总压等熵亚声速流动。
exact = exact_nozzle_shock_solution(x,A,gamma,p_exit);

%% ---------------- 初始条件 ----------------
% 采用教材中类似的分段初始条件，使初始场接近含激波解
rho0 = zeros(1,Nx);
T0   = zeros(1,Nx);

for i = 1:Nx
    xi = x(i);

    if xi <= 1.5
        rho0(i) = 1.0 - 0.366*xi;
        T0(i)   = 1.0 - 0.167*xi;

    elseif xi <= 2.1
        rho0(i) = 0.634  - 0.7020*(xi-1.5);
        T0(i)   = 0.833  - 0.4908*(xi-1.5);

    else
        rho0(i) = 0.5892 + 0.10228*(xi-2.1);
        T0(i)   = 0.93968 + 0.0622*(xi-2.1);
    end
end

% 防止初始条件中出现极小值
rho0 = max(rho0,1e-6);
T0   = max(T0,1e-6);

% 初始质量流量
mdot0 = 0.59;
V0 = mdot0./(rho0.*A);

% Anderson 无量纲守恒变量
U0 = zeros(3,Nx);
U0(1,:) = rho0.*A;
U0(2,:) = rho0.*V0.*A;
U0(3,:) = rho0.*A.*(T0/(gamma-1)+gamma/2*V0.^2);

%% ---------------- 计算 ----------------
result_noAV = run_nozzle_solver( ...
    U0,A,dAdx,x,dx,gamma,CFL,Nt_noAV,p_exit,Cx_noAV, ...
    doAnimate,'No artificial viscosity');

result_AV = run_nozzle_solver( ...
    U0,A,dAdx,x,dx,gamma,CFL,Nt_AV,p_exit,Cx_AV, ...
    doAnimate,'With artificial viscosity');

%% ---------------- 最终一维对比图：数值解 vs 解析解 ----------------
figure('Color','w','Position',[80 80 1450 820]);
tiledlayout(2,3,'TileSpacing','compact','Padding','compact');

% ---------- Pressure ----------
nexttile;
plot(x,result_noAV.p,'--','LineWidth',1.2); hold on;
plot(x,result_AV.p,'-','LineWidth',2.0);
plot(x,exact.p,'k-.','LineWidth',1.8);
xline(result_AV.x_shock,'b:','LineWidth',1.5);
xline(exact.x_shock,'k:','LineWidth',1.8);
grid on;
xlabel('x/L');
ylabel('p/p_0');
title('Pressure');
legend('No AV','With AV','Exact','Numerical shock','Exact shock','Location','best');

% ---------- Mach Number ----------
nexttile;
plot(x,result_noAV.M,'--','LineWidth',1.2); hold on;
plot(x,result_AV.M,'-','LineWidth',2.0);
plot(x,exact.M,'k-.','LineWidth',1.8);
yline(1,'k--','LineWidth',1.0);
xline(result_AV.x_shock,'b:','LineWidth',1.5);
xline(exact.x_shock,'k:','LineWidth',1.8);
grid on;
xlabel('x/L');
ylabel('M');
title('Mach Number');
legend('No AV','With AV','Exact','M=1','Numerical shock','Exact shock','Location','best');

% ---------- Density ----------
nexttile;
plot(x,result_noAV.rho,'--','LineWidth',1.2); hold on;
plot(x,result_AV.rho,'-','LineWidth',2.0);
plot(x,exact.rho,'k-.','LineWidth',1.8);
xline(result_AV.x_shock,'b:','LineWidth',1.5);
xline(exact.x_shock,'k:','LineWidth',1.8);
grid on;
xlabel('x/L');
ylabel('\rho/\rho_0');
title('Density');
legend('No AV','With AV','Exact','Numerical shock','Exact shock','Location','best');

% ---------- Temperature ----------
nexttile;
plot(x,result_noAV.T,'--','LineWidth',1.2); hold on;
plot(x,result_AV.T,'-','LineWidth',2.0);
plot(x,exact.T,'k-.','LineWidth',1.8);
xline(result_AV.x_shock,'b:','LineWidth',1.5);
xline(exact.x_shock,'k:','LineWidth',1.8);
grid on;
xlabel('x/L');
ylabel('T/T_0');
title('Temperature');
legend('No AV','With AV','Exact','Numerical shock','Exact shock','Location','best');

% ---------- Velocity ----------
nexttile;
plot(x,result_noAV.V,'--','LineWidth',1.2); hold on;
plot(x,result_AV.V,'-','LineWidth',2.0);
plot(x,exact.V,'k-.','LineWidth',1.8);
xline(result_AV.x_shock,'b:','LineWidth',1.5);
xline(exact.x_shock,'k:','LineWidth',1.8);
grid on;
xlabel('x/L');
ylabel('V/a_0');
title('Velocity');
legend('No AV','With AV','Exact','Numerical shock','Exact shock','Location','best');

% ---------- Residual ----------
nexttile;
semilogy(result_noAV.res+eps,'--','LineWidth',1.2); hold on;
semilogy(result_AV.res+eps,'-','LineWidth',1.8);
grid on;
xlabel('Iteration');
ylabel('Residual');
title('Residual History');
legend('No AV','With AV','Location','best');

sgtitle('Shock-Capturing Nozzle Flow: Numerical Solution vs Exact Solution', ...
    'FontWeight','bold','FontSize',16);

%% ---------------- 误差统计 ----------------
err_M   = sqrt(mean((result_AV.M   - exact.M).^2));
err_p   = sqrt(mean((result_AV.p   - exact.p).^2));
err_rho = sqrt(mean((result_AV.rho - exact.rho).^2));
err_T   = sqrt(mean((result_AV.T   - exact.T).^2));
err_V   = sqrt(mean((result_AV.V   - exact.V).^2));

fprintf('\n========== Final summary ==========\n');
fprintf('No AV:  iter = %d, shock x/L ≈ %.4f, exit p = %.4f\n', ...
    result_noAV.iter_end,result_noAV.x_shock,result_noAV.p(end));
fprintf('AV:     iter = %d, shock x/L ≈ %.4f, exit p = %.4f\n', ...
    result_AV.iter_end,result_AV.x_shock,result_AV.p(end));

fprintf('\n========== Numerical vs Exact ==========\n');
fprintf('解析激波位置 x_s        = %.6f\n',exact.x_shock);
fprintf('数值激波位置 x_s        = %.6f\n',result_AV.x_shock);
fprintf('激波位置误差            = %.6e\n',abs(result_AV.x_shock-exact.x_shock));
fprintf('Mach L2 error           = %.6e\n',err_M);
fprintf('Pressure L2 error       = %.6e\n',err_p);
fprintf('Density L2 error        = %.6e\n',err_rho);
fprintf('Temperature L2 error    = %.6e\n',err_T);
fprintf('Velocity L2 error       = %.6e\n',err_V);

fprintf('\n========== Exact shock information ==========\n');
fprintf('Exact shock x/L         = %.6f\n',exact.x_shock);
fprintf('Exact shock area        = %.6f\n',exact.A_shock);
fprintf('M1 before shock         = %.6f\n',exact.M1);
fprintf('M2 after shock          = %.6f\n',exact.M2);
fprintf('p2/p1                  = %.6f\n',exact.p2_p1);
fprintf('p02/p01                = %.6f\n',exact.p02_p01);

%% ===================== 二维 contour 展示 =====================
% 准一维结果沿横向复制，仅用于课堂可视化展示
% 注意：这不是二维流场求解结果，而是准一维结果的二维几何展示

% 选择用于展示的结果：有人工粘性的稳定结果
rho_c = result_AV.rho;
T_c   = result_AV.T;
p_c   = result_AV.p;
M_c   = result_AV.M;
V_c   = result_AV.V;
mdot_c = result_AV.mdot;

x_shock_num   = result_AV.x_shock;
x_shock_exact = exact.x_shock;

% 保证都是行向量
x      = x(:).';
A      = A(:).';
rho_c  = rho_c(:).';
T_c    = T_c(:).';
p_c    = p_c(:).';
M_c    = M_c(:).';
V_c    = V_c(:).';
mdot_c = mdot_c(:).';

% 横向网格
Ny = 101;
eta = linspace(-1,1,Ny);

% 喷管半高，用面积 A 构造二维外形
% yScale 只影响显示厚度，不改变物理计算
yScale = 0.35;
h = yScale*A/2;

[X2D,ETA2D] = meshgrid(x,eta);

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

% 将准一维变量复制到二维截面
rho2D  = repmat(rho_c,Ny,1);
T2D    = repmat(T_c,Ny,1);
p2D    = repmat(p_c,Ny,1);
M2D    = repmat(M_c,Ny,1);
V2D    = repmat(V_c,Ny,1);
mdot2D = repmat(mdot_c,Ny,1);

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

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

x_wall = x;
y_top  = h;
y_bot  = -h;

figure('Color','w','Position',[80 80 1450 780]);

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

nLevel = 60;

% ---------- Pressure ----------
ax1 = 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);
xline(x_shock_num,'k--','LineWidth',1.4);
xline(x_shock_exact,'w--','LineWidth',1.4);
xlabel('x/L');
ylabel('y');
title('Pressure');
cb1 = colorbar;
cb1.Label.String = 'p/p_0';
grid on;
box on;
axis tight;
pbaspect([3.2 1 1]);

% ---------- Mach Number ----------
ax2 = 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_shock_num,'k--','LineWidth',1.4);
xline(x_shock_exact,'w--','LineWidth',1.4);
xline(x_throat,'w:','LineWidth',1.4);
xlabel('x/L');
ylabel('y');
title('Mach Number');
cb2 = colorbar;
cb2.Label.String = 'M';
grid on;
box on;
axis tight;
pbaspect([3.2 1 1]);

% ---------- Density ----------
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);
xline(x_shock_num,'k--','LineWidth',1.4);
xline(x_shock_exact,'w--','LineWidth',1.4);
xlabel('x/L');
ylabel('y');
title('Density');
cb3 = colorbar;
cb3.Label.String = '\rho/\rho_0';
grid on;
box on;
axis tight;
pbaspect([3.2 1 1]);

% ---------- Temperature ----------
ax4 = 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);
xline(x_shock_num,'k--','LineWidth',1.4);
xline(x_shock_exact,'w--','LineWidth',1.4);
xlabel('x/L');
ylabel('y');
title('Temperature');
cb4 = colorbar;
cb4.Label.String = 'T/T_0';
grid on;
box on;
axis tight;
pbaspect([3.2 1 1]);

% ---------- Velocity ----------
ax5 = nexttile;
contourf(X2D,Y2D,V2D,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_shock_num,'k--','LineWidth',1.4);
xline(x_shock_exact,'w--','LineWidth',1.4);
xlabel('x/L');
ylabel('y');
title('Velocity');
cb5 = colorbar;
cb5.Label.String = 'V/a_0';
grid on;
box on;
axis tight;
pbaspect([3.2 1 1]);

% ---------- Mass Flow ----------
ax6 = nexttile;
contourf(X2D,Y2D,mdot2D,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_shock_num,'k--','LineWidth',1.4);
xline(x_shock_exact,'w--','LineWidth',1.4);
xlabel('x/L');
ylabel('y');
title('Mass Flow');
cb6 = colorbar;
cb6.Label.String = '\rho V A';
grid on;
box on;
axis tight;
pbaspect([3.2 1 1]);

% ---------- 统一外观 ----------
colormap(turbo);

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

sgtitle('Shock-Capturing Nozzle Flow: Final Contour Visualization', ...
    'FontWeight','bold','FontSize',16);

fprintf('\n二维 contour 中：黑色虚线为数值激波位置，白色虚线为解析激波位置。\n');

%% ============================================================
%  主求解函数
%% ============================================================
function result = run_nozzle_solver(U,A,dAdx,x,dx,gamma,CFL,Nt,p_exit,Cx,doAnimate,caseName)

    Nx = length(x);
    res = nan(Nt,1);

    if doAnimate
        figure('Color','w','Position',[100 80 1320 720]);
        tiledlayout(2,2,'TileSpacing','compact','Padding','compact');

        ax1 = nexttile;
        hP = plot(x,nan(size(x)),'LineWidth',2); hold on;
        hPs = xline(0,'k:','LineWidth',1.4);
        grid on;
        xlabel('x/L');
        ylabel('p/p_0');
        title(['Pressure: ',caseName]);

        ax2 = nexttile;
        hM = plot(x,nan(size(x)),'LineWidth',2); hold on;
        yline(1,'k:','LineWidth',1.2);
        hMs = xline(0,'k:','LineWidth',1.4);
        grid on;
        xlabel('x/L');
        ylabel('M');
        title('Mach number');

        ax3 = nexttile;
        hMass = plot(x,nan(size(x)),'LineWidth',2);
        grid on;
        xlabel('x/L');
        ylabel('\rho V A');
        title('Mass flow');

        ax4 = nexttile;
        hRes = semilogy(nan,nan,'LineWidth',1.8);
        grid on;
        xlabel('Iteration');
        ylabel('Residual');
        title('Residual history');
    end

    iter_end = Nt;

    for n = 1:Nt

        %% ---------------- 解码原始变量 ----------------
        [rho,V,T,p,M,mdot] = primitive_from_U(U,A,gamma);

        if any(~isfinite(rho)) || any(~isfinite(T)) || any(rho <= 0) || any(T <= 0)
            warning('%s stopped at iter %d because rho/T became nonphysical.',caseName,n);
            iter_end = n-1;
            break;
        end

        a = sqrt(T);
        dt = CFL*dx/max(abs(V)+a);

        %% ====================================================
        %  Predictor step
        %% ====================================================
        F = flux_from_U(U,gamma);
        J = source_from_pressure(p,dAdx,gamma);

        U_bar = U;

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

        S = artificial_viscosity(U,p,Cx);

        for i = 2:Nx-1
            U_bar(:,i) = U_bar(:,i) + S(:,i);
        end

        U_bar = apply_boundary_conditions(U_bar,A,gamma,p_exit);

        [rho_b,V_b,T_b,p_b,M_b,mdot_b] = primitive_from_U(U_bar,A,gamma);

        if any(~isfinite(rho_b)) || any(~isfinite(T_b)) || any(rho_b <= 0) || any(T_b <= 0)
            warning('%s stopped at predictor iter %d because rho/T became nonphysical.',caseName,n);
            iter_end = n-1;
            break;
        end

        %% ====================================================
        %  Corrector step
        %% ====================================================
        F_bar = flux_from_U(U_bar,gamma);
        J_bar = source_from_pressure(p_b,dAdx,gamma);

        U_new = U;

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

        S_bar = artificial_viscosity(U_bar,p_b,Cx);

        for i = 2:Nx-1
            U_new(:,i) = U_new(:,i) + S_bar(:,i);
        end

        U_new = apply_boundary_conditions(U_new,A,gamma,p_exit);

        [rho_new,V_new,T_new,p_new,M_new,mdot_new] = primitive_from_U(U_new,A,gamma);

        if any(~isfinite(rho_new)) || any(~isfinite(T_new)) || any(rho_new <= 0) || any(T_new <= 0)
            warning('%s stopped at corrector iter %d because rho/T became nonphysical.',caseName,n);
            iter_end = n-1;
            break;
        end

        %% ---------------- 残差 ----------------
        res(n) = max(abs(U_new(:)-U(:)));
        U = U_new;

        %% ---------------- 动画 ----------------
        if doAnimate && (mod(n,5)==0 || n==1 || n==Nt)

            [rho,V,T,p,M,mdot] = primitive_from_U(U,A,gamma);
            x_shock = estimate_shock_location(x,M,p);

            set(hP,'YData',p);
            set(hPs,'Value',x_shock);

            set(hM,'YData',M);
            set(hMs,'Value',x_shock);

            set(hMass,'YData',mdot);

            set(hRes,'XData',1:n,'YData',res(1:n)+eps);

            title(ax1,sprintf('Pressure: %s, iter=%d, C_x=%.2f',caseName,n,Cx));
            title(ax2,sprintf('Mach number, shock x/L \\approx %.3f',x_shock));

            drawnow limitrate;
            pause(0.01)
        end
    end

    res = res(1:max(iter_end,1));

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

    result.U = U;
    result.rho = rho;
    result.V = V;
    result.T = T;
    result.p = p;
    result.M = M;
    result.mdot = mdot;
    result.res = res;
    result.iter_end = iter_end;
    result.x_shock = estimate_shock_location(x,M,p);
end

%% ============================================================
%  Anderson 无量纲守恒变量解码
%% ============================================================
function [rho,V,T,p,M,mdot] = primitive_from_U(U,A,gamma)

    rho = U(1,:)./A;
    V   = U(2,:)./U(1,:);

    % U3/U1 = T/(gamma-1)+gamma/2*V^2
    T = (gamma-1)*(U(3,:)./U(1,:) - gamma/2*V.^2);

    p = rho.*T;

    M = V./sqrt(T);

    mdot = rho.*V.*A;
end

%% ============================================================
%  Anderson 教材形式的通量
%% ============================================================
function F = flux_from_U(U,gamma)

    F = zeros(size(U));

    U1 = U(1,:);
    U2 = U(2,:);
    U3 = U(3,:);

    F(1,:) = U2;

    F(2,:) = U2.^2./U1 ...
        + (gamma-1)/gamma*(U3 - gamma/2*U2.^2./U1);

    F(3,:) = gamma*U2.*U3./U1 ...
        - gamma*(gamma-1)/2*U2.^3./U1.^2;
end

%% ============================================================
%  几何源项
%% ============================================================
function J = source_from_pressure(p,dAdx,gamma)

    J = zeros(3,length(p));

    % Anderson 无量纲形式中，动量方程源项为：
    % J2 = (1/gamma)*p*dA/dx
    J(2,:) = p.*dAdx/gamma;
end

%% ============================================================
%  人工粘性
%% ============================================================
function S = artificial_viscosity(U,p,Cx)

    [nvar,Nx] = size(U);
    S = zeros(nvar,Nx);

    if Cx == 0
        return;
    end

    for i = 2:Nx-1

        denom = p(i+1)+2*p(i)+p(i-1);
        denom = max(denom,1e-12);

        sensor = abs(p(i+1)-2*p(i)+p(i-1))/denom;

        S(:,i) = Cx*sensor*(U(:,i+1)-2*U(:,i)+U(:,i-1));
    end
end

%% ============================================================
%  边界条件
%% ============================================================
function U = apply_boundary_conditions(U,A,gamma,p_exit)

    Nx = length(A);

    %% ---------------- 入口边界 ----------------
    % 亚声速入口：
    % 固定入口总量的简化教学处理，rho=T=1；
    % 速度从内部外推。

    V2 = U(2,2)/U(1,2);
    V3 = U(2,3)/U(1,3);

    V_in = 2*V2 - V3;
    V_in = max(V_in,0.02);
    V_in = min(V_in,0.99);

    rho_in = 1.0;
    T_in   = 1.0;

    U(1,1) = rho_in*A(1);
    U(2,1) = rho_in*V_in*A(1);
    U(3,1) = rho_in*A(1)*(T_in/(gamma-1)+gamma/2*V_in^2);

    %% ---------------- 出口边界 ----------------
    % 亚声速出口：
    % U1、U2外推，p_exit指定，由p_exit修正U3。

    U(1,Nx) = 2*U(1,Nx-1) - U(1,Nx-2);
    U(2,Nx) = 2*U(2,Nx-1) - U(2,Nx-2);

    rho_e = U(1,Nx)/A(Nx);
    rho_e = max(rho_e,1e-8);

    V_e = U(2,Nx)/U(1,Nx);

    T_e = p_exit/rho_e;
    T_e = max(T_e,1e-8);

    U(3,Nx) = rho_e*A(Nx)*(T_e/(gamma-1)+gamma/2*V_e^2);
end

%% ============================================================
%  激波位置估计
%% ============================================================
function x_shock = estimate_shock_location(x,M,p)

    idx = find(M(1:end-1)>1 & M(2:end)<1,1,'first');

    if ~isempty(idx)
        x1 = x(idx);
        x2 = x(idx+1);
        M1 = M(idx);
        M2 = M(idx+1);

        x_shock = x1 + (1-M1)*(x2-x1)/(M2-M1);
    else
        [~,idx] = max(abs(gradient(p,x)));
        x_shock = x(idx);
    end
end

%% ============================================================
%  含正激波的准一维喷管理想解析解
%% ============================================================
function exact = exact_nozzle_shock_solution(x,A,gamma,p_exit)

    x = x(:).';
    A = A(:).';

    [At,ithroat] = min(A);
    x_throat = x(ithroat);

    % 喉部阻塞，入口总压总温无量纲取 1
    Astar1 = At;
    p01_ratio = 1.0;
    T01_ratio = 1.0;

    % 搜索激波位置，使出口压力等于指定 p_exit
    xs_left  = x_throat + 1e-6;
    xs_right = x(end) - 1e-6;

    shock_fun = @(xs) predicted_exit_pressure(xs,x,A,gamma,Astar1,p01_ratio) - p_exit;

    xs_scan = linspace(xs_left,xs_right,400);
    f_scan = zeros(size(xs_scan));

    for k = 1:length(xs_scan)
        f_scan(k) = shock_fun(xs_scan(k));
    end

    idx = find(f_scan(1:end-1).*f_scan(2:end) <= 0,1,'first');

    if isempty(idx)
        warning('未找到与出口压力匹配的内激波位置，改用压力差最小的位置。');
        [~,idx_min] = min(abs(f_scan));
        x_shock = xs_scan(idx_min);
    else
        x_shock = fzero(shock_fun,[xs_scan(idx),xs_scan(idx+1)]);
    end

    A_shock = interp1(x,A,x_shock,'linear');

    % 激波前马赫数：超声速分支
    M1 = solveAreaMach(A_shock/Astar1,gamma,'supersonic');

    % 正激波关系
    M2 = normalShockM2(M1,gamma);
    p2_p1 = normalShockPressureRatio(M1,gamma);
    p02_p01 = normalShockTotalPressureRatio(M1,gamma);

    % 激波后新的临界面积
    Astar2 = A_shock/areaMachFunction(M2,gamma);

    p02_ratio = p01_ratio*p02_p01;
    T02_ratio = T01_ratio;

    % 初始化解析解
    M   = zeros(size(x));
    T   = zeros(size(x));
    p   = zeros(size(x));
    rho = zeros(size(x));
    V   = zeros(size(x));

    for i = 1:length(x)

        if x(i) <= x_throat
            % 喉部前：亚声速等熵
            M(i) = solveAreaMach(A(i)/Astar1,gamma,'subsonic');
            T(i) = T01_ratio/(1+(gamma-1)/2*M(i)^2);
            p(i) = p01_ratio*(T(i)/T01_ratio)^(gamma/(gamma-1));

        elseif x(i) < x_shock
            % 激波前：超声速等熵
            M(i) = solveAreaMach(A(i)/Astar1,gamma,'supersonic');
            T(i) = T01_ratio/(1+(gamma-1)/2*M(i)^2);
            p(i) = p01_ratio*(T(i)/T01_ratio)^(gamma/(gamma-1));

        else
            % 激波后：以损失后的总压为基准，亚声速等熵
            M(i) = solveAreaMach(A(i)/Astar2,gamma,'subsonic');
            T(i) = T02_ratio/(1+(gamma-1)/2*M(i)^2);
            p(i) = p02_ratio*(T(i)/T02_ratio)^(gamma/(gamma-1));
        end

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

    mdot = rho.*V.*A;

    exact.x_shock = x_shock;
    exact.A_shock = A_shock;
    exact.M1 = M1;
    exact.M2 = M2;
    exact.p2_p1 = p2_p1;
    exact.p02_p01 = p02_p01;
    exact.Astar1 = Astar1;
    exact.Astar2 = Astar2;

    exact.M = M;
    exact.T = T;
    exact.p = p;
    exact.rho = rho;
    exact.V = V;
    exact.mdot = mdot;
end

%% ============================================================
%  给定激波位置时，预测出口压力
%% ============================================================
function p_exit_pred = predicted_exit_pressure(xs,x,A,gamma,Astar1,p01_ratio)

    A_shock = interp1(x,A,xs,'linear');

    % 激波前超声速马赫数
    M1 = solveAreaMach(A_shock/Astar1,gamma,'supersonic');

    % 激波后马赫数与总压损失
    M2 = normalShockM2(M1,gamma);
    p02_p01 = normalShockTotalPressureRatio(M1,gamma);

    % 激波后临界面积
    Astar2 = A_shock/areaMachFunction(M2,gamma);

    % 出口马赫数：激波后亚声速分支
    Ae = A(end);
    Me = solveAreaMach(Ae/Astar2,gamma,'subsonic');

    % 出口压力，相对于入口总压 p01
    p02_ratio = p01_ratio*p02_p01;

    Te_T02 = 1/(1+(gamma-1)/2*Me^2);
    pe_p02 = Te_T02^(gamma/(gamma-1));

    p_exit_pred = p02_ratio*pe_p02;
end

%% ============================================================
%  面积-马赫数函数 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

%% ============================================================
%  由 A/A* 反求 M
%% ============================================================
function M = solveAreaMach(areaRatio,gamma,branch)

    areaRatio = max(areaRatio,1.0+1e-12);

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

    switch lower(branch)

        case 'subsonic'
            M_left  = 1e-8;
            M_right = 1-1e-8;

        case 'supersonic'
            M_left  = 1+1e-8;
            M_right = 20;

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

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

%% ============================================================
%  正激波后马赫数
%% ============================================================
function M2 = normalShockM2(M1,gamma)

    M2 = sqrt( ...
        (1+(gamma-1)/2*M1^2) / ...
        (gamma*M1^2-(gamma-1)/2) );
end

%% ============================================================
%  正激波静压比 p2/p1
%% ============================================================
function p2_p1 = normalShockPressureRatio(M1,gamma)

    p2_p1 = 1 + 2*gamma/(gamma+1)*(M1^2-1);
end

%% ============================================================
%  正激波总压比 p02/p01
%% ============================================================
function p02_p01 = normalShockTotalPressureRatio(M1,gamma)

    term1 = ((gamma+1)*M1^2) / ((gamma-1)*M1^2+2);
    term2 = (gamma+1) / (2*gamma*M1^2-(gamma-1));

    p02_p01 = term1^(gamma/(gamma-1)) * term2^(1/(gamma-1));
end