% ============================================================
%
% 1. 二维超声 Prandtl-Meyer 膨胀波，定常 Euler 方程守恒形式
% 2. MacCormack predictor-corrector 空间推进
% 3. 贴体坐标 eta=(y-ys)/h
% 4. 人工粘性
% 5. 简化 Abbett 壁面边界条件
% 6. 数值结果与 Prandtl-Meyer 解析解对比
%
% ============================================================

clear; clc; close all;

% ===================== 1. 基本参数 =====================

gamma = 1.4;
R     = 287.0;

% 来流条件
M1   = 2.0;
p1   = 1.01e5;        % N/m^2
rho1 = 1.23;          % kg/m^3
T1   = 286.1;         % K

a1 = sqrt(gamma*R*T1);
u1 = M1*a1;
v1 = 0.0;

% 几何参数
E     = 10.0;                 % 膨胀角位置
H     = 40.0;                 % 上游通道高度
theta = deg2rad(5.352);       % 壁面向下转角
xEnd  = 66.23;                % 计算终点

% 网格参数
Nj   = 41;                    % eta方向网格点数
eta  = linspace(0,1,Nj);
deta = eta(2)-eta(1);

% 数值参数
CFL = 0.5;
Cy  = 0.6;

% 需要保存和绘图的站位
targetX = [0, 16.17, 32.31, 48.99, 66.23];

% 速度剖面绘图缩放
profileScale = 0.035;

% ===================== 2. 解析解：波后状态 =====================

nu1 = PMFunction(M1,gamma);
nu2 = nu1 + theta;
M2  = inversePM(nu2,gamma);

T0 = T1*(1+(gamma-1)/2*M1^2);
p0 = p1*(1+(gamma-1)/2*M1^2)^(gamma/(gamma-1));

T2   = T0/(1+(gamma-1)/2*M2^2);
p2   = p0/(1+(gamma-1)/2*M2^2)^(gamma/(gamma-1));
rho2 = p2/(R*T2);
a2   = sqrt(gamma*R*T2);
V2   = M2*a2;

u2 = V2*cos(theta);
v2 = -V2*sin(theta);

mu1 = asin(1/M1);
mu2 = asin(1/M2);

fprintf('\n===== Prandtl-Meyer 解析波后状态 =====\n');
fprintf('theta = %.4f deg\n',rad2deg(theta));
fprintf('M1    = %.4f\n',M1);
fprintf('M2    = %.4f\n',M2);
fprintf('p2    = %.4e N/m^2\n',p2);
fprintf('rho2  = %.4f kg/m^3\n',rho2);
fprintf('T2    = %.2f K\n',T2);
fprintf('u2    = %.2f m/s\n',u2);
fprintf('v2    = %.2f m/s\n',v2);

% ===================== 3. 初始条件 =====================

x = 0.0;

rho = rho1*ones(1,Nj);
u   = u1*ones(1,Nj);
v   = v1*ones(1,Nj);
p   = p1*ones(1,Nj);
T   = T1*ones(1,Nj);

F = primitiveToF(rho,u,v,p,gamma);

% 保存剖面
profileData = struct([]);
profileData(1).x   = x;
profileData(1).eta = eta;
profileData(1).rho = rho;
profileData(1).u   = u;
profileData(1).v   = v;
profileData(1).p   = p;
profileData(1).T   = T;
profileData(1).M   = M1*ones(1,Nj);

nextTargetID = 2;

xHist = x;
dxiHist = [];

% ===================== 4. 下游推进 =====================

iter = 0;
maxIter = 10000;

while x < xEnd-1e-10 && iter < maxIter

    iter = iter + 1;

    % 当前站位几何
    [ys,h] = geometryPM(x,E,H,theta);
    etax = metricEtaX(x,eta,E,H,theta);

    % 当前原始变量
    [rho,u,v,p,T,M] = FToPrimitive(F,gamma,R);

    % 当前 G 通量
    G = primitiveToG(rho,u,v,p,gamma);

    % ========== 4.1 空间推进步长 Delta xi ==========

    flowAngle = atan2(v,u);
    M_safe = max(M,1.000001);
    mu = asin(1./M_safe);

    tanPlus  = abs(tan(flowAngle+mu));
    tanMinus = abs(tan(flowAngle-mu));
    tanMax = max([tanPlus,tanMinus],[],'all');

    % 当前站位的物理横向网格间距
    dY  = h*deta;
    dxi = CFL*dY/tanMax;

    % 为了正好保存目标站位，必要时缩短当前步长
    if nextTargetID <= numel(targetX)
        if x < targetX(nextTargetID) && x+dxi > targetX(nextTargetID)
            dxi = targetX(nextTargetID)-x;
        end
    end

    if x+dxi > xEnd
        dxi = xEnd-x;
    end

    % ========== 4.2 Predictor step ==========

    dFdxi = zeros(size(F));

    for j = 2:Nj-1
        dFdxi(:,j) = etax(j)*(F(:,j)-F(:,j+1))/deta ...
                   + (1/h)*(G(:,j)-G(:,j+1))/deta;
    end

    % 人工粘性，使用当前站位变量
    SF = artificialViscosity(F,p,Cy);

    Fbar = F;

    for j = 2:Nj-1
        Fbar(:,j) = F(:,j) + dFdxi(:,j)*dxi + SF(:,j);
    end

    xNew = x + dxi;

    % 预测值的边界处理
    % 上边界：本算例中膨胀波尚未穿出上边界，取来流状态
    Fbar(:,Nj) = primitiveToF(rho1,u1,v1,p1,gamma);

    % 壁面：由近壁内点构造壁面预测值
    Fbar(:,1) = wallBoundaryCondition(Fbar(:,2),xNew,E,theta,gamma,R);

    % 由预测通量恢复预测原始变量
    [rhob,ub,vb,pb,Tb,Mb] = FToPrimitive(Fbar,gamma,R);
    Gbar = primitiveToG(rhob,ub,vb,pb,gamma);

    % ========== 4.3 Corrector step ==========

    dFbardxi = zeros(size(F));

    % 教材中校正步使用预测值和后向差分
    % 这里的 h 和 etax 仍使用当前站位 i 的值
    for j = 2:Nj-1
        dFbardxi(:,j) = etax(j)*(Fbar(:,j-1)-Fbar(:,j))/deta ...
                      + (1/h)*(Gbar(:,j-1)-Gbar(:,j))/deta;
    end

    dFavg = 0.5*(dFdxi + dFbardxi);

    % 校正步人工粘性，使用预测压力
    SFbar = artificialViscosity(Fbar,pb,Cy);

    Fnew = F;

    for j = 2:Nj-1
        Fnew(:,j) = F(:,j) + dFavg(:,j)*dxi + SFbar(:,j);
    end

    % ========== 4.4 边界条件 ==========

    % 上边界：均匀来流
    Fnew(:,Nj) = primitiveToF(rho1,u1,v1,p1,gamma);

    % 壁面边界：速度与壁面相切
    Fnew(:,1) = wallBoundaryCondition(Fnew(:,2),xNew,E,theta,gamma,R);

    % ========== 4.5 数值保护 ==========

    [rhoN,uN,vN,pN,TN,MN] = FToPrimitive(Fnew,gamma,R);

    % 若出现异常，停止并提示
    if any(~isfinite(rhoN)) || any(~isfinite(pN)) || any(rhoN <= 0) || any(pN <= 0)
        warning('数值解出现非物理状态，计算停止。iter=%d, x=%.4f',iter,xNew);
        break;
    end

    % ========== 4.6 更新站位 ==========

    x = xNew;
    F = Fnew;

    xHist(end+1) = x; %#ok<SAGROW>
    dxiHist(end+1) = dxi; %#ok<SAGROW>

    % ========== 4.7 保存目标站位 ==========

    if nextTargetID <= numel(targetX)
        if abs(x-targetX(nextTargetID)) < 1e-8
            [rho,u,v,p,T,M] = FToPrimitive(F,gamma,R);

            profileData(nextTargetID).x   = x;
            profileData(nextTargetID).eta = eta;
            profileData(nextTargetID).rho = rho;
            profileData(nextTargetID).u   = u;
            profileData(nextTargetID).v   = v;
            profileData(nextTargetID).p   = p;
            profileData(nextTargetID).T   = T;
            profileData(nextTargetID).M   = M;

            nextTargetID = nextTargetID + 1;
        end
    end
end

fprintf('\n===== 数值计算结束 =====\n');
fprintf('迭代步数 iter = %d\n',iter);
fprintf('最终 x = %.4f m\n',x);
fprintf('平均 dxi = %.4f m\n',mean(dxiHist));
fprintf('最小 dxi = %.4f m\n',min(dxiHist));
fprintf('最大 dxi = %.4f m\n',max(dxiHist));

% ===================== 5. 最终站位误差估计 =====================

last = profileData(end);
xLast = last.x;
[ysLast,hLast] = geometryPM(xLast,E,H,theta);
yLast = ysLast + last.eta*hLast;

% 选取波后区域：低于后缘线且避开壁面点
dxLast = xLast - E;
yTrail = dxLast*tan(mu2-theta);
idxDownstream = find(yLast < yTrail & last.eta > 0.05 & last.eta < 0.55);

if ~isempty(idxDownstream)
    uMean   = mean(last.u(idxDownstream));
    vMean   = mean(last.v(idxDownstream));
    pMean   = mean(last.p(idxDownstream));
    rhoMean = mean(last.rho(idxDownstream));
    TMean   = mean(last.T(idxDownstream));
    MMean   = mean(last.M(idxDownstream));

    fprintf('\n===== 最下游波后区平均误差估计 =====\n');
    fprintf('M  error = %.3f %%\n',abs(MMean-M2)/M2*100);
    fprintf('p  error = %.3f %%\n',abs(pMean-p2)/p2*100);
    fprintf('rho error = %.3f %%\n',abs(rhoMean-rho2)/rho2*100);
    fprintf('T  error = %.3f %%\n',abs(TMean-T2)/T2*100);
    fprintf('u  error = %.3f %%\n',abs(uMean-u2)/u2*100);
    fprintf('v  error = %.3f %%\n',abs(vMean-v2)/abs(v2)*100);
end

% ===================== 6. 绘图：速度剖面对比 =====================

figure('Color','w','Position',[80 80 1450 640]);
hold on; box on;

for k = 1:numel(profileData)

    xk = profileData(k).x;
    etak = profileData(k).eta;

    [ysk,hk] = geometryPM(xk,E,H,theta);
    yk = ysk + etak*hk;

    ukNum = profileData(k).u;
    ukAna = analyticUProfile(xk,yk,M1,p1,T1,rho1,gamma,R,E,theta);

    % 数值解：实线
    plot(xk + profileScale*(ukNum-u1), yk, ...
        'k-', 'LineWidth', 1.6);

    % 解析解：黑点
    plot(xk + profileScale*(ukAna-u1), yk, ...
        'ko', 'MarkerFaceColor','k', 'MarkerSize',4);

    % 站位线
    plot([xk xk],[min(yk) max(yk)],'k:','LineWidth',0.6);

    text(xk,min(yk)-2.0,sprintf('x=%.2f',xk), ...
        'HorizontalAlignment','center','FontSize',10);
end

% 壁面
xWall = linspace(0,xEnd,500);
yWall = zeros(size(xWall));

for n = 1:numel(xWall)
    [ysn,~] = geometryPM(xWall(n),E,H,theta);
    yWall(n) = ysn;
end

plot(xWall,yWall,'k-','LineWidth',2.0);

% 膨胀波前缘与后缘
xFan = linspace(E,xEnd,300);
yLead  = (xFan-E)*tan(mu1);
yTrail = (xFan-E)*tan(mu2-theta);

plot(xFan,yLead,'k--','LineWidth',1.3);
plot(xFan,yTrail,'k--','LineWidth',1.3);

text(E+16,(16)*tan(mu1)+1.2, ...
    'Leading edge of expansion wave', ...
    'Rotation',rad2deg(atan(tan(mu1))), ...
    'FontSize',10);

text(E+16,(16)*tan(mu2-theta)-1.2, ...
    'Trailing edge of expansion wave', ...
    'Rotation',rad2deg(atan(tan(mu2-theta))), ...
    'FontSize',10);

xlabel('x + scale \times (u-u_\infty)');
ylabel('y');
title('Prandtl-Meyer 膨胀波：数值解与解析解速度剖面对比');

legend({'数值解','解析解'},'Location','northwest');

grid on;
xlim([-2 xEnd+4]);
ylim([-8 44]);

% ===================== 7. 绘图：最终站位变量分布 =====================

last = profileData(end);
xk = last.x;
[ysk,hk] = geometryPM(xk,E,H,theta);
yk = ysk + last.eta*hk;

figure('Color','w','Position',[120 90 1150 780]);
tiledlayout(2,2,'Padding','compact','TileSpacing','compact');

nexttile;
plot(last.u,yk,'k-o','LineWidth',1.3,'MarkerSize',4);
xlabel('u (m/s)');
ylabel('y (m)');
title(sprintf('x = %.2f m: u',xk));
grid on;

nexttile;
plot(last.v,yk,'k-o','LineWidth',1.3,'MarkerSize',4);
xlabel('v (m/s)');
ylabel('y (m)');
title(sprintf('x = %.2f m: v',xk));
grid on;

nexttile;
plot(last.p,yk,'k-o','LineWidth',1.3,'MarkerSize',4);
xlabel('p (N/m^2)');
ylabel('y (m)');
title(sprintf('x = %.2f m: p',xk));
grid on;

nexttile;
plot(last.M,yk,'k-o','LineWidth',1.3,'MarkerSize',4);
xlabel('M');
ylabel('y (m)');
title(sprintf('x = %.2f m: Mach number',xk));
grid on;


fprintf('\n===== 观察点 =====\n');
fprintf('1. x=0 处应为均匀来流。\n');
fprintf('2. 膨胀角附近误差最大，这是几何奇异性和网格分辨率共同造成的。\n');
fprintf('3. 下游越远，膨胀波越宽，数值解与解析解通常更接近。\n');
fprintf('4. 壁面点可能存在数值速度层，不应只用壁面点判断整体精度。\n');
fprintf('5. 增大 Cy 会增强数值耗散；Cy=0 时可能出现振荡。\n');
fprintf('6. 增大 Nj 可检查网格收敛趋势。\n');











%% ============================================================
%                           局部函数
% ============================================================

function F = primitiveToF(rho,u,v,p,gamma)
    % x方向守恒通量 F
    F = zeros(4,numel(rho));

    F(1,:) = rho.*u;
    F(2,:) = rho.*u.^2 + p;
    F(3,:) = rho.*u.*v;
    F(4,:) = gamma/(gamma-1)*p.*u ...
           + rho.*u.*(u.^2+v.^2)/2;
end

function G = primitiveToG(rho,u,v,p,gamma)
    % y方向守恒通量 G
    G = zeros(4,numel(rho));

    G(1,:) = rho.*v;
    G(2,:) = rho.*u.*v;
    G(3,:) = rho.*v.^2 + p;
    G(4,:) = gamma/(gamma-1)*p.*v ...
           + rho.*v.*(u.^2+v.^2)/2;
end

function [rho,u,v,p,T,M] = FToPrimitive(F,gamma,R)
    % 由 F1-F4 反解原始变量
    %
    % 教材 Eq. (8.8)-(8.12)
    %
    % 注意关键公式：
    % A = F3^2/(2F1) - F4
    % 不是 F3^2/(2F1^2) - F4

    F1 = F(1,:);
    F2 = F(2,:);
    F3 = F(3,:);
    F4 = F(4,:);

    A = F3.^2./(2*F1) - F4;
    B = gamma/(gamma-1)*F1.*F2;
    C = -(gamma+1)/(2*(gamma-1))*F1.^3;

    disc = B.^2 - 4*A.*C;

    % 数值保护
    disc = max(disc,0);

    rho = (-B + sqrt(disc))./(2*A);

    rho = max(rho,1e-10);

    u = F1./rho;
    v = F3./F1;
    p = F2 - F1.*u;

    p = max(p,1e-8);

    T = p./(rho*R);
    T = max(T,1e-8);

    a = sqrt(gamma*R*T);
    M = sqrt(u.^2+v.^2)./a;
end

function [ys,h,dhdx,dysdx] = geometryPM(x,E,H,theta)
    % 物理平面几何
    %
    % x <= E:
    %   ys = 0
    %   h  = H
    %
    % x > E:
    %   ys = -(x-E)tan(theta)
    %   h  = H + (x-E)tan(theta)

    if x <= E
        ys = 0.0;
        h  = H;
        dysdx = 0.0;
        dhdx  = 0.0;
    else
        ys = -(x-E)*tan(theta);
        h  = H + (x-E)*tan(theta);
        dysdx = -tan(theta);
        dhdx  =  tan(theta);
    end
end

function etax = metricEtaX(x,eta,E,H,theta)
    % eta_x = partial eta / partial x
    %
    % x <= E:
    %   eta_x = 0
    %
    % x > E:
    %   eta_x = (1-eta)*tan(theta)/h

    [~,h] = geometryPM(x,E,H,theta);

    if x <= E
        etax = zeros(size(eta));
    else
        etax = (1-eta)*tan(theta)/h;
    end
end

function SF = artificialViscosity(F,p,Cy)
    % Jameson 型压力传感器人工粘性，沿 eta 方向
    %
    % SF_j = sensor_j * (F_{j+1} - 2F_j + F_{j-1})
    %
    % sensor_j = Cy*|p_{j+1}-2p_j+p_{j-1}|/(p_{j+1}+2p_j+p_{j-1})

    [~,Nj] = size(F);
    SF = zeros(size(F));

    for j = 2:Nj-1
        denom = p(j+1)+2*p(j)+p(j-1);
        denom = max(denom,1e-12);

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

        for k = 1:4
            SF(k,j) = sensor*(F(k,j+1)-2*F(k,j)+F(k,j-1));
        end
    end
end

function Fwall = wallBoundaryCondition(Fnear,x,E,theta,gamma,R)
    % Abbott 壁面边界条件
    %
    % 目的：
    %   使壁面处速度与壁面相切，即 V dot n = 0
    %
    % 处理：
    %   1. 用近壁内点估计壁面试算状态
    %   2. 计算试算速度方向与壁面方向的偏差
    %   3. 用局部 PM 修正得到 Mact
    %   4. 用 Mact,Tact 得到速度大小
    %   5. 强制速度沿壁面切向

    [~,uCal,vCal,pCal,TCal,MCal] = FToPrimitive(Fnear,gamma,R);

    if x <= E
        wallAngle = 0.0;
    else
        wallAngle = -theta;
    end

    alphaCal = atan2(vCal,uCal);

    % 需要旋转的角度
    phi = alphaCal - wallAngle;

    nuCal = PMFunction(MCal,gamma);
    nuAct = nuCal + phi;

    % 防止落到 M<1
    nuAct = max(nuAct,1e-10);

    MAct = inversePM(nuAct,gamma);

    ratioT = ...
        (1+(gamma-1)/2*MCal^2) ...
        /(1+(gamma-1)/2*MAct^2);

    TAct   = TCal*ratioT;
    pAct   = pCal*ratioT^(gamma/(gamma-1));
    rhoAct = pAct/(R*TAct);

    aAct = sqrt(gamma*R*TAct);
    VAct = MAct*aAct;

    uAct = VAct*cos(wallAngle);
    vAct = VAct*sin(wallAngle);

    Fwall = primitiveToF(rhoAct,uAct,vAct,pAct,gamma);
end

function nu = PMFunction(M,gamma)
    % Prandtl-Meyer 函数，单位 rad

    M = max(M,1.0000001);

    nu = sqrt((gamma+1)/(gamma-1))*atan( ...
        sqrt((gamma-1)/(gamma+1)*(M.^2-1)) ) ...
        - atan(sqrt(M.^2-1));
end

function M = inversePM(nuTarget,gamma)
    % 由 Prandtl-Meyer 函数反求 Mach 数
    % 简单二分法，教学演示

    nuTarget = max(nuTarget,0);

    Mlow  = 1.000001;
    Mhigh = 30.0;

    for k = 1:100
        Mmid = 0.5*(Mlow+Mhigh);
        nuMid = PMFunction(Mmid,gamma);

        if nuMid < nuTarget
            Mlow = Mmid;
        else
            Mhigh = Mmid;
        end
    end

    M = 0.5*(Mlow+Mhigh);
end

function uExact = analyticUProfile(x,y,M1,p1,T1,rho1,gamma,R,E,theta)
    % 计算给定 x 站位上不同 y 点的 PM 解析 u 分布
    %
    % 区域判断：
    %   beta > mu1             : 波前，仍为来流
    %   beta < mu2 - theta     : 波后，完全转过 theta
    %   其余                  : 膨胀波内部

    %#ok<INUSD> rho1 not used explicitly, kept for interface completeness

    a1 = sqrt(gamma*R*T1);
    V1 = M1*a1;

    nu1 = PMFunction(M1,gamma);
    nu2 = nu1 + theta;
    M2  = inversePM(nu2,gamma);

    mu1 = asin(1/M1);
    mu2 = asin(1/M2);

    T0 = T1*(1+(gamma-1)/2*M1^2);

    uExact = zeros(size(y));

    if x <= E
        uExact(:) = V1;
        return;
    end

    dx = x - E;

    betaLead  = mu1;
    betaTrail = mu2 - theta;

    for n = 1:numel(y)

        beta = atan2(y(n),dx);

        if beta >= betaLead
            % 膨胀波前
            M = M1;
            delta = 0.0;

        elseif beta <= betaTrail
            % 膨胀波后
            M = M2;
            delta = theta;

        else
            % 膨胀波内部
            % 几何关系：
            % beta = mu(M) - delta
            % delta = nu(M) - nu1

            Mlo = M1;
            Mhi = M2;

            for k = 1:80
                Mm = 0.5*(Mlo+Mhi);

                nuM = PMFunction(Mm,gamma);
                muM = asin(1/Mm);
                deltaM = nuM - nu1;

                betaM = muM - deltaM;

                if betaM > beta
                    Mlo = Mm;
                else
                    Mhi = Mm;
                end
            end

            M = 0.5*(Mlo+Mhi);
            delta = PMFunction(M,gamma)-nu1;
        end

        T = T0/(1+(gamma-1)/2*M^2);
        a = sqrt(gamma*R*T);
        V = M*a;

        uExact(n) = V*cos(delta);
    end
end