U(4,jindex,iindex)=P_in/(gam-1)+rho_in*Vt^2/2;
        end
    end
    
    U_in = U(:,1,1);
    U_out = U(:,1,M);
    
   
    for jindex=1:N
        for iindex=3:M-2
            U(:,jindex,iindex) = U_in*((M - 2 - iindex)/(M-4)) + U_out*((iindex - 2)/(M-4));
        end
    end    
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
function [U,para]=getpara(U,Data)   
    gam=Data.gam;
    cfl=Data.cfl;
    X=Data.CyPX;
    Y=Data.CyPY;
    M=Data.Nx+1;
    N=Data.Ny+1;

    for iindex=1:N        % N 為行數，M為列數 第2列為進口 M－1 列為出口
        for jindex=1:M
            rho0(iindex,jindex)=U(1,iindex,jindex);  %初始時刻的流場參數值
            u0(iindex,jindex)=U(2,iindex,jindex)/rho0(iindex,jindex);
            v0(iindex,jindex)=U(3,iindex,jindex)/rho0(iindex,jindex);
            E0(iindex,jindex)=U(4,iindex,jindex)/rho0(iindex,jindex);
            p0(iindex,jindex)=(gam-1)*rho0(iindex,jindex)*(E0(iindex,jindex)-0.5*((u0(iindex,jindex))^2+(v0(iindex,jindex))^2));
            a0(iindex,jindex)=sqrt(gam*p0(iindex,jindex)/rho0(iindex,jindex));
            M0(iindex,jindex)=sqrt((u0(iindex,jindex))^2+(v0(iindex,jindex))^2)/a0(iindex,jindex);
        end
    end
    para=struct('gam',gam,'cfl',cfl,'rho',rho0,'p',p0,'u',u0,'v',v0,'a',a0,'N',N,'M',M,'LRPX',Data.CyPX,'LRPY',Data.CyPY,'BTPX',Data.CxPX,'BTPY',Data.CxPY,'X',Data.UPX,'Y',Data.UPY,'PX',Data.PX,'PY',Data.PY,'E0',E0,'M0',M0);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
function [U,dt]=Update(U,para)    
    k=5;
    RK_ceof(1)= 0.1067;  龍格－褲它 格式 各項時間多部系數 5 步
    RK_ceof(2)= 0.1979;
    RK_ceof(3)= 0.3232;
    RK_ceof(4)= 0.5201;
    RK_ceof(5)= 1.0000;
    while k
        [Uh_L,Uh_R,Uh_B,Uh_T]=GetU_L_R_B_T_Vanleer(U,para);   
        dt=tstep(Uh_L,Uh_R,Uh_B,Uh_T,para); 
        [FhX,GhX,FhY,GhY]=GetFh_MUSCL(U,Uh_L,Uh_R,Uh_B,Uh_T,para);
        U=Runge_Kutta(RK_ceof(6-k),U,FhX,GhX,FhY,GhY,dt,para);
        k=k-1;
    end
    
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5   
function U=Runge_Kutta(afa,U,FhX,GhX,FhY,GhY,dt,para)
    gam=para.gam;
    M=para.M;
    N=para.N;
    X=para.X;
    Y=para.Y;
    PX=para.PX;
    PY=para.PY;
    
    for i=2:N-1
        for j=2:M-1
            dx12=PX(i-1,j-1)-PX(i,j-1);
            dy12=PY(i-1,j-1)-PY(i,j-1);
        
            dx34=PX(i,j)-PX(i-1,j);
            dy34=PY(i,j)-PY(i-1,j);
           
            dx23=PX(i-1,j)-PX(i-1,j-1);
            dy23=PY(i-1,j)-PY(i-1,j-1);
            
            dx41=PX(i,j-1)-PX(i,j);
            dy41=PY(i,j-1)-PY(i,j);
            
            S=0.5*abs((PX(i,j)-PX(i-1,j-1))*(PY(i-1,j)-PY(i,j-1))-(PY(i,j)-PY(i-1,j-1))*(PX(i-1,j)-PX(i,j-1)));
            
            U(:,i,j)=U(:,i,j)-(dt/S)*afa*( (FhX(:,i,j-1)*dy12 - GhX(:,i,j-1)*dx12) + (FhX(:,i,j)*dy34 - GhX(:,i,j)*dx34) + (FhY(:,i-1,j)*dy23 - GhY(:,i-1,j)*dx23) + (FhY(:,i,j)*dy41 - GhY(:,i,j)*dx41) );
        end
    end
    
    %在這里進行邊界條件的設置
    %進口,出口邊界
    P_in=130000;%入口靜壓Pa
    P_out=101325;%出口反壓
    Vt=200.0;%入口平均速度m
    afa=0;%入口角
    for iindex=2:N-1
        for jindex=2:M-1
            rho0(iindex,jindex)=U(1,iindex,jindex);
            u0(iindex,jindex)=U(2,iindex,jindex)/rho0(iindex,jindex);
            v0(iindex,jindex)=U(3,iindex,jindex)/rho0(iindex,jindex);
            E0(iindex,jindex)=U(4,iindex,jindex)/rho0(iindex,jindex);
            p0(iindex,jindex)=(gam-1)*rho0(iindex,jindex)*(E0(iindex,jindex)-0.5*((u0(iindex,jindex))^2+(v0(iindex,jindex))^2));
            a0(iindex,jindex)=sqrt(gam*p0(iindex,jindex)/rho0(iindex,jindex));
            M0(iindex,jindex)=sqrt((u0(iindex,jindex))^2+(v0(iindex,jindex))^2)/a0(iindex,jindex);
        end

        %進口處理
        rho0(iindex,1)=1.579;
        u0(iindex,1)=Vt*cos(afa);
        v0(iindex,1)=Vt*sin(afa);
        E0(iindex,1)=2*E0(iindex,2)-E0(iindex,3);
        p0(iindex,1)=(gam-1)*rho0(iindex,1)*(E0(iindex,1)-0.5*((u0(iindex,1))^2+(v0(iindex,1))^2));
        a0(iindex,1)=sqrt(gam*p0(iindex,1)/rho0(iindex,1));
        M0(iindex,1)=sqrt((u0(iindex,1))^2+(v0(iindex,1))^2)/a0(iindex,1);
        
        U(1,iindex,1)=rho0(iindex,1);
        U(2,iindex,1)=rho0(iindex,1)*u0(iindex,1);
        U(3,iindex,1)=rho0(iindex,1)*v0(iindex,1);
        U(4,iindex,1)=p0(iindex,1)/(gam-1)+rho0(iindex,1)*((u0(iindex,1))^2+(v0(iindex,1))^2)/2;
        
	    rho0(iindex,2)=(rho0(iindex,1)+rho0(iindex,3))/2;
        u0(iindex,2)=(u0(iindex,1)+u0(iindex,3))/2;
        v0(iindex,2)=(v0(iindex,1)+v0(iindex,3))/2;

        U(1,iindex,2)=rho0(iindex,2);
        U(2,iindex,2)=rho0(iindex,2)*u0(iindex,2);
        U(3,iindex,2)=rho0(iindex,2)*v0(iindex,2);
        U(4,iindex,2)=p0(iindex,2)/(gam-1)+rho0(iindex,2)*((u0(iindex,2))^2+(v0(iindex,2))^2)/2;

        %出口處理
        u0(iindex,M)=2*u0(iindex,M-1)-u0(iindex,M-2);
        v0(iindex,M)=2*v0(iindex,M-1)-v0(iindex,M-2);
        rho0(iindex,M)=2*rho0(iindex,M-1)-rho0(iindex,M-2);
        p0(iindex,M)=P_out;
        E0(iindex,M)=p0(iindex,M)/(gam-1)/rho0(iindex,M)+(u0(iindex,M))^2+(v0(iindex,M))^2;
        a0(iindex,M)=sqrt(gam*p0(iindex,M)/rho0(iindex,M));
        M0(iindex,M)=sqrt((u0(iindex,M))^2+(v0(iindex,M))^2)/a0(iindex,M);
    
        U(1,iindex,M)=rho0(iindex,M);
        U(2,iindex,M)=rho0(iindex,M)*u0(iindex,M);
        U(3,iindex,M)=rho0(iindex,M)*v0(iindex,M);
        U(4,iindex,M)=p0(iindex,M)/(gam-1)+rho0(iindex,M)*((u0(iindex,M))^2+(v0(iindex,M))^2)/2;

        p0(iindex,M-1)=(p0(iindex,M)+p0(iindex,M-2))/2;    
        U(1,iindex,M-1)=rho0(iindex,M-1);
        U(2,iindex,M-1)=rho0(iindex,M-1)*u0(iindex,M-1);
        U(3,iindex,M-1)=rho0(iindex,M-1)*v0(iindex,M-1);
        U(4,iindex,M-1)=p0(iindex,M-1)/(gam-1)+rho0(iindex,M-1)*((u0(iindex,M-1))^2+(v0(iindex,M-1))^2)/2;
    end

   %壁面邊界    
   for iindex=2:M-1
       temp_rho0=U(1,2,iindex);
       temp_u0=U(2,2,iindex)/temp_rho0;
       temp_v0=U(3,2,iindex)/temp_rho0;
       temp_E0=U(4,2,iindex)/temp_rho0;
       temp_p0=(gam-1)*temp_rho0*(temp_E0-0.5*(temp_u0^2+temp_v0^2));
       
       %法向速度置為0
        V = sqrt(temp_u0^2+temp_v0^2);
        tg_a1=temp_v0/temp_u0;
        a1=atan(tg_a1);
        tg_a2=(Y(N-1,iindex)-Y(N-1,iindex-1))/(X(N-1,iindex)-X(N-1,iindex-1));
        a2=atan(tg_a2);
        Vt=V*cos(a1-a2);
        temp_u0=Vt*cos(a2);
        temp_v0=-Vt*sin(a2);
         
        temp_rho02=U(1,3,iindex);
        temp_u02=U(2,3,iindex)/temp_rho02;
        temp_v02=U(3,3,iindex)/temp_rho02;
        temp_E02=U(4,3,iindex)/temp_rho02;
        temp_p02=(gam-1)*temp_rho02*(temp_E02-0.5*(temp_u02^2+temp_v02^2));
	    temp_p0=temp_p02;
        temp_E0=temp_E02;
        
        %再根據新值來確定U
        U(1,2,iindex) = temp_rho0;
        U(2,2,iindex) = temp_rho0*temp_u0;
        U(3,2,iindex) = temp_rho0*temp_v0;
        U(4,2,iindex) = temp_rho0*temp_E0;
        
        temp_rho0=U(1,N-1,iindex);
        temp_u0=U(2,N-1,iindex)/temp_rho0;
        temp_v0=U(3,N-1,iindex)/temp_rho0;
        temp_E0=U(4,N-1,iindex)/temp_rho0;
        temp_p0=(gam-1)*temp_rho0*(temp_E0-0.5*(temp_u0^2+temp_v0^2));
       
       %法向速度置為0
        V = sqrt(temp_u0^2+temp_v0^2);
        tg_a1=temp_v0/temp_u0;
        a1=atan(tg_a1);
        tg_a2=(Y(N-1,iindex)-Y(N-1,iindex-1))/(X(N-1,iindex)-X(N-1,iindex-1));
        a2=atan(tg_a2);
        Vt=V*cos(a1-a2);
        temp_u0=Vt*cos(a2);
        temp_v0=Vt*sin(a2);
        
        temp_rho02=U(1,N-2,iindex);
        temp_u02=U(2,N-2,iindex)/temp_rho02;
        temp_v02=U(3,N-2,iindex)/temp_rho02;
        temp_E02=U(4,N-2,iindex)/temp_rho02;
        temp_p02=(gam-1)*temp_rho02*(temp_E02-0.5*(temp_u02^2+temp_v02^2));
	    temp_p0=temp_p02;
        temp_E0=temp_E02;

        %再根據新值來確定U
        U(1,N-1,iindex) = temp_rho0;
        U(2,N-1,iindex) = temp_rho0*temp_u0;
        U(3,N-1,iindex) = temp_rho0*temp_v0;
        U(4,N-1,iindex) = temp_rho0*temp_E0;  
        
        U(:,1,iindex) = 2*U(:,2,iindex) - U(:,3,iindex);
        U(:,N,iindex) = 2*U(:,N-1,iindex) - U(:,N-2,iindex);
    end
            
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
function [Uh_L,Uh_R,Uh_B,Uh_T]=GetU_L_R_B_T_Vanleer(U,para);   %  構造 MUSCL 格式

    gam=para.gam;
    M=para.M;
    N=para.N;
    X=para.X;
    Y=para.Y;
    temp_e=0.00000001;
    temp_g=1/3;
    
    %單方向分步求解
    for i=1:N
        for j=2:M-1
            for k=1:4
                deltaUp_X(k)=U(k,i,j+1)-U(k,i,j);
                deltaDown_X(k)=U(k,i,j)-U(k,i,j-1);
                q_X(k)=2* deltaUp_X(k) * deltaDown_X(k)/(temp_e + deltaUp_X(k)^2 + deltaDown_X(k)^2);
                S_L(k)=0.5*q_X(k)*((1-q_X(k)*temp_g)*deltaUp_X(k)+(1+q_X(k)*temp_g)*deltaDown_X(k));
                S_R(k)=0.5*q_X(k)*((1+q_X(k)*temp_g)*deltaUp_X(k)+(1-q_X(k)*temp_g)*deltaDown_X(k));
                U_t_L(k,i,j)=U(k,i,j)+S_L(k)/2;
                U_t_R(k,i,j)=U(k,i,j)-S_R(k)/2;
            end
        end
    end
    
    for i=1:N
        for j=1:M-1
          if(j==1)
              Uh_L(:,i,j)=U(:,i,j); 
              Uh_R(:,i,j)=U_t_L(:,i,j+1);
            elseif(j==M-1)
              Uh_L(:,i,j)=U_t_R(:,i,j);
              Uh_R(:,i,j)=U(:,i,j+1);
            else
              Uh_L(:,i,j)=U_t_R(:,i,j);
              Uh_R(:,i,j)=U_t_L(:,i,j+1);
            end   
        end 
    end
    
    
    for j=1:M
        for i=2:N-1
            for k=1:4
                deltaUp_Y(k)=U(k,i+1,j)-U(k,i,j);
                deltaDown_Y(k)=U(k,i,j)-U(k,i-1,j);
                q_Y(k)=2* deltaUp_Y(k) * deltaDown_Y(k)/(temp_e + deltaUp_Y(k)^2 + deltaDown_Y(k)^2);
                S_B(k)=0.5*q_Y(k)*((1-q_Y(k)*temp_g)*deltaUp_Y(k)+(1+q_Y(k)*temp_g)*deltaDown_Y(k));
                S_T(k)=0.5*q_Y(k)*((1+q_Y(k)*temp_g)*deltaUp_Y(k)+(1-q_Y(k)*temp_g)*deltaDown_Y(k));
                U_t_B(k,i,j)=U(k,i,j)+S_B(k)/2;
                U_t_T(k,i,j)=U(k,i,j)-S_T(k)/2;
            end
        end
    end
    
    for j=1:M
        for i=1:N-1
          if(i==1)
              Uh_B(:,i,j)=U(:,i,j); 
              Uh_T(:,i,j)=U_t_B(:,i+1,j);
            elseif(i==N-1)
              Uh_B(:,i,j)=U_t_T(:,i,j);
              Uh_T(:,i,j)=U(:,i+1,j);