%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%數值模擬二維收縮噴管流場
%顯式時間推進、有限體積法、結構化網格，2階MUSCL格式，超限插值網格
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function nozzle()                     
    Data=setdata;                  
    U=init(Data);                  
    [U,para]=getpara(U,Data);     
    
    fp=fopen('progress.txt','w');
    fclose(fp);
    
    N=para.N-1;
    M=para.M-1;
    F(1,:,:)=para.p;壓力
    F(2,:,:)=para.rho; 密度
    F(3,:,:)=para.u;
    F(4,:,:)=para.v;
    F(5,:,:)=para.a;
    F(6,:,:)=para.E0;
    F(7,:,:)=para.M0;
    fp=fopen('result.plt', 'w');
    for j=2:N
        for i=2:M
            fprintf(fp,'%12.6f  %12.6f  %12.6f  %12.6f  %12.6f  %12.6f  %12.6f  %12.6f  %12.6f\r\n', para.X(j, i), para.Y(j, i), F(1, j, i), F(2, j, i), F(3, j, i), F(4, j, i), F(5, j, i), F(6, j, i), F(7, j, i));
        end
        fprintf(fp,'\r\n');  格心坐標 對應參數值
    end
    fclose(fp);
    
    run=1;                         
    t=0;
    sumt=1;
    step=0;
    sumstep=1000;
    error_constant=1e-4;
    while run
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       fp=fopen('key.txt','r');
        exit=fscanf(fp,'%d');.
        fclose(fp);
        if exit==0
            break;    
        end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
        oldU=U(1,:,:);
        [U,dt]=Update(U,para);
        [U,para]=getpara(U,Data);   
        maxerror=max(max(abs(oldU(1,:,:)-U(1,:,:))));
       
        t=t+dt;
        step=step+1;
        
        run=stopcalculate(t,sumt,step,sumstep,error_constant,maxerror);    
   
        fp1=fopen('progress.txt','a+');
        str1=num2str(t);
        str2=num2str(step);
        str3=num2str(maxerror);
        fprintf(fp1,'\r\nsumtime：%d,sumstep：%d,error_constant=%f',sumt,sumstep,error_constant);
        fprintf(fp1,'\r\ncurrent_t=%s\r\ncurrent_step=%s',str1,str2);
        fprintf(fp1,'\r\ncurrent_maxerror=%s\r\n',str3);
        fclose(fp1);
   
        N=para.N-1;
        M=para.M-1;
        F(1,:,:)=para.p;
        F(2,:,:)=para.rho;
        F(3,:,:)=para.u;
        F(4,:,:)=para.v;
        F(5,:,:)=para.a;
        F(6,:,:)=para.E0;
        F(7,:,:)=para.M0;
        
        fp=fopen('result.plt', 'w');
        for j=2:N
            for i=2:M
                fprintf(fp,'%12.6f  %12.6f  %12.6f  %12.6f  %12.6f  %12.6f  %12.6f  %12.6f  %12.6f\r\n', para.X(j, i), para.Y(j, i), F(1, j, i), F(2, j, i), F(3, j, i), F(4, j, i), F(5, j, i), F(6, j, i), F(7, j, i));
            end
            fprintf(fp,'\r\n');
        end
        fclose(fp);
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
function Data=setdata()
    cfl=0.8;                      
    gam=1.4; 
    [PointX,PointY,Nx,Ny]=setgrid2D; %     劃分網格 網格數據節點
    [CyPointX,CyPointY,CxPointX,CxPointY]=getcentrepoint(PointX,PointY,Nx,Ny);% 包圍節點的  矩形四個點坐標
    [UPointX,UPointY,PointX,PointY]=getunitpoint(PointX,PointY,Nx,Ny);  % 格心坐標
    method='char';                  
    Data=struct('gam',gam,'cfl',cfl,'Nx',Nx,'Ny',Ny,'PX',PointX,'PY',PointY,'CyPX',CyPointX,'CyPY',CyPointY,'CxPX',CxPointX,'CxPY',CxPointY,'UPX',UPointX,'UPY',UPointY,'method',method);       
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
function [PointX,PointY,Nx,Ny]=setgrid2D() %應用超限插值法構造網格單元
clear all;
R1=1;   %噴管入口半徑                       
R2=20;  %候道的彎曲半徑
Cta=15; %候道的圓周半徑
S=0.1;  %單位長度上網格的大小                  
M2=ceil(R2*Cta*pi/180)/S;   %x方向的網格點數
N2=R1*2/S;    %y方向的網格點數
Nx=M2;
Ny=N2;

%Px(i,j)--網格(i,j)點的x坐標
%Py(i,j)--網格(i,j)點的y坐標
    Px(1,1)=0;
    Px(M2,1)=R2*sin(Cta/360*pi)*2;
    Px(1,N2)=0;
    Px(M2,N2)=R2*sin(Cta/360*pi)*2;

    for j=1:M2   %構造Lagrange 插值型函數
        qx1(j)=1-(j-1)/(M2-1);
        qx2(j)=(j-1)/(M2-1);   
    end
    cta=Cta*pi/180; %角度轉換成弧度
    dcta=(Cta*pi/180)/(M2-1);   %M2-1等分
    l=2*R2*sin(dcta/2);    %每等份的玄長
    
    a1=R2*sin(cta/2);
    b1=-R2*cos(cta/2);
    a2=a1;
    b2=2*R1+R2*cos(cta/2);    %圓心偏移點
    
    x(1,1)=Px(1,1);
    x(1,N2)=Px(1,1);
    y(1,1)=0;
    y(1,N2)=R1*2;
    for i=2:M2
        x(i,1)=x(i-1,1)+l*cos(cta/2-qx2(i)*dcta/2);
        x(i,N2)=x(i,1);        
        y(i,1)=sqrt(R2^2-(x(i,1)-a1)^2)+b1;
        y(i,N2)=-sqrt(R2^2-(x(i,N2)-a2)^2)+b2;
    end

    Py(1,1)=0;
    Py(M2,1)=0; 
    Py(1,N2)=R1*2;
    Py(M2,N2)=R1*2;
    for j=1:N2
        qy1(j)=1-(j-1)/(N2-1);
        qy2(j)=(j-1)/(N2-1);
    end 
    for j=1:N2
        x(1,j)=0;
        x(M2,j)=R2*sin(Cta/360*pi)*2;
        y(1,j)=qy1(j)*Py(1,1)+qy2(j)*Py(1,N2);
        y(M2,j)=qy1(j)*Py(M2,1)+qy2(j)*Py(M2,N2);
    end

    for i=1:M2
       for j=1:N2
           x(i,j)=qx1(i)*x(1,j)+qx2(i)*x(M2,j);
           y(i,j)=qy1(j)*y(i,1)+qy2(j)*y(i,N2);
       end
    end
    hold on
    plot(x,y,'b')
    plot(x',y','b')
    hold off
   
PointX=x.';
PointY=y.';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
function [CyPointX,CyPointY,CxPointX,CxPointY]=getcentrepoint(PointX,PointY,Nx,Ny)
    for iindex=1:Ny-1
        for jindex=1:Nx
            CyPointX(iindex,jindex)=( PointX(iindex+1, jindex) + PointX(iindex, jindex) )/2;
            CyPointY(iindex,jindex)=( PointY(iindex+1, jindex) + PointY(iindex, jindex) )/2;
        end
    end
    for iindex=1:Ny
        for jindex=1:Nx-1
            CxPointX(iindex,jindex)=( PointX(iindex, jindex+1) + PointX(iindex, jindex) )/2;
            CxPointY(iindex,jindex)=( PointY(iindex, jindex+1) + PointY(iindex, jindex) )/2;
        end
    end
    
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
function [UPointX,UPointY,PointX,PointY]=getunitpoint(PointX,PointY,Nx,Ny)%            
    UPointX=zeros(Ny+1,Nx+1);
    UPointY=zeros(Ny+1,Nx+1);
     for iindex=1:Ny-1
        for jindex=1:Nx-1
            UPointX(iindex+1,jindex+1)=( PointX(iindex,jindex) + PointX(iindex+1,jindex) + PointX(iindex,jindex+1) + PointX(iindex+1,jindex+1))/4;
            UPointY(iindex+1,jindex+1)=( PointY(iindex,jindex) + PointY(iindex+1,jindex) + PointY(iindex,jindex+1) + PointY(iindex+1,jindex+1))/4;
        end
    end
    for iindex=1:Ny-1
        UPointX(iindex+1,1) = PointX(iindex,1) + PointX(iindex+1,1) - UPointX(iindex+1,2);
        UPointY(iindex+1,1) = PointY(iindex,1) + PointY(iindex+1,1) - UPointY(iindex+1,2);
        UPointX(iindex+1,Nx+1) = PointX(iindex,Nx) + PointX(iindex+1,Nx) - UPointX(iindex+1,Nx);
        UPointY(iindex+1,Nx+1) = PointY(iindex,Nx) + PointY(iindex+1,Nx) - UPointY(iindex+1,Nx);
    end
    for jindex=1:Nx-1
        UPointX(1,jindex+1) = PointX(1,jindex) + PointX(1,jindex+1) - UPointX(2,jindex+1);
        UPointY(1,jindex+1) = PointY(1,jindex) + PointY(1,jindex+1) - UPointY(2,jindex+1);
        UPointX(Ny+1,jindex+1) = PointX(Ny,jindex+1) + PointX(Ny,jindex) - UPointX(Ny,jindex+1);
        UPointY(Ny+1,jindex+1) = PointY(Ny,jindex+1) + PointY(Ny,jindex) - UPointY(Ny,jindex+1);
    end
    UPointX(1,1) = 2*PointX(1,1) - UPointX(2,2);
    UPointY(1,1) = 2*PointY(1,1) - UPointY(2,2);
    UPointX(1,Nx+1) = 2*PointX(1,Nx) - UPointX(2,Nx);
    UPointY(1,Nx+1) = 2*PointY(1,Nx) - UPointY(2,Nx);
    UPointX(Ny+1,1) = 2*PointX(Ny,1) - UPointX(Ny,2);
    UPointY(Ny+1,1) = 2*PointY(Ny,1) - UPointY(Ny,2);
    UPointX(Ny+1,Nx+1) = 2*PointX(Ny,Nx) - UPointX(Ny,Nx);
    UPointY(Ny+1,Nx+1) = 2*PointY(Ny,Nx) - UPointY(Ny,Nx);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [X,Y]=changerow(X,Y,n)
    for iindex=1:n
        tempX(iindex,:)=X(n-iindex+1,:);
        tempY(iindex,:)=Y(n-iindex+1,:);
    end
    X=tempX;
    Y=tempY;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
function U=init(Data)      
    gam=Data.gam;
    M=Data.Nx+1;       N=Data.Ny+1;
    PX=Data.PX;        PY=Data.PY;    % 生成的網格節點坐標
    CyX=Data.CyPX;     CyY=Data.CyPY;% 包圍節點的四個點坐標
    CxX=Data.CxPX;     CxY=Data.CxPY;
    UX=Data.UPX;       UY=Data.UPY;% 格心  ，包括 虛格心 m X n
    U=zeros(4,N,M);
    
   
    P_in=130000;%入口靜壓Pa
    T_in=300.0;%入口靜溫K
    Vt=200.0;%入口平均速度m
    afa=0;%入口角
    R=287.;%氣體常數
    rho_in=P_in/R/T_in;
   
    %入口初始化
    uVelocity=Vt*cos(afa);
    vVelocity=Vt*sin(afa);
    a_in=sqrt(gam*R*T_in);
    Ma=Vt/a_in;
    T0_in=T_in*(1+Ma^2*(gam-1)/2);
    P0_in=P_in*(T0_in/T_in)^(gam/(gam-1));
    rho0_in=rho_in*(T0_in/T_in)^(1/(gam-1));
    
    for iindex=1:2
        for jindex=1:N
            U(1,jindex,iindex)=rho_in;
            U(2,jindex,iindex)=rho_in*uVelocity;  % x方向的速度×密度
            U(3,jindex,iindex)=rho_in*vVelocity; % y方向的速度×密度
            U(4,jindex,iindex)=P_in/(gam-1)+rho_in*Vt^2/2;
        end
    end
    %出口初始化，初步假設出口氣流角一致afa
    P_out=101325;
    P0_out=P0_in;
    T0_out=T0_in;
    rho0_out=rho0_in;
    Ma=sqrt(2*((P0_out/P_out)^((gam-1)/gam)-1)/(gam-1));
    T_out=T0_out/(1+(gam-1)*Ma^2/2);
    rho_out=rho0_out/((1+(gam-1)*Ma^2/2)^(1/(gam-1)));
    Vt=sqrt(2*gam*R*(T0_out-T_out)/(gam-1));
    
    tg_afa=(UY(N-2,M-1)-UY(N-2,M-2))/(UX(N-2,M-1)-UX(N-2,M-2));
    cos_afa=1/sqrt(tg_afa^2+1);
    sin_afa=sqrt(1-cos_afa^2);
    uVelocity=Vt*cos_afa;
    vVelocity=Vt*sin_afa;

    for iindex=M-1:M
        for jindex=1:N
            U(1,jindex,iindex)=rho_out;
            U(2,jindex,iindex)=rho_out*uVelocity;
            U(3,jindex,iindex)=rho_out*vVelocity;