%This script will attempt to solve the Quasi 1-d Euler equations
%
Cv=.17;                         %Specific Heat in BTU/lb/deg R
gamma=1.4;                      %Ratio of specific heats
gamma_fac=(gamma-1)/(gamma+1);  %Used coefficient
Beta=gamma-1;                   %Used coefficient
sigma=1;                        %Constant used in ROE 
finite_vol=1;                   %1=finite volume        2=finite difference
flux_method=2;                  %1=Modified Staeger     2=Roe Method
ghost_cell=2;                   %1=use ghost cell       2=no ghost cell
flow_case=2;                    %1=No flow              2=pressure dropped at exit

%The grid extends from -.5 to 1.5 and it's divided into 41 equally spaced points
IL=41;            IO=[1:IL];        IM=[2:IL-1];        %Surfaces
%The grid encloses either volumes or represents the surfaces
if finite_vol==1
    IVL=IL-1;         IVO=[1:IVL-1];      IVM=[2:IVL-2];      %IL-1 Volumes
else
    IVL=IL;         IVO=[1:IL];      IVM=[2:IL-1];            %IL Surfaces
end

xleft=-.5;        xright=1.5;       delx=(xright-xleft)/(IL-1);
x=[xleft:delx:xright]; 
xvols=x(2:IL)-.5*delx;

num_vols=IVL;

%The leading edge of the chord is at 0 and the trailing edge is at 1
I_lead_edge=(IL-1)*.25+1;         I_trail_edge=(IL-1)*.75+1;
%Get indices of before, over, and after chord
I_before_chord=[1:I_lead_edge];
I_over_chord=[I_lead_edge+1:I_trail_edge];
I_after_chord=[I_trail_edge+1:IL];

%Figure out the xvalues of the chord and the surface y value S(x)
tau=.03;                    c=1;              Surf=-.1*ones(1,IL);
[xo,yo,R]=Circular_Chord_Quas_1d(tau,c);
Surf(I_over_chord)=sqrt(R^2-(x(I_over_chord)-xo).^2)+yo-.1;
%Now get the actual width of the channel and the corresponding volumes
Surf=0-Surf;
if ghost_cell==1
  Ghost_cell=delx.*(Surf(2)+Surf(1))/2;
  Vol=[Ghost_cell,delx.*(Surf(2:IL)+Surf(1:IL-1))./2];
else
  Vol=delx.*(Surf(2:IL)+Surf(1:IL-1))./2;
end

%Set up the flux vectors.  Initial conditions are determined using Newton's
%method and isentropic relations 
%F is the flux vector---- F(1,:)=p      F(2,:)=pu   F(3,:)=e
U=zeros(3,IVL);           Vols=repmat(Vol,3,1);     Surfs=repmat(Surf,3,1);
Q=zeros(3,IVL);           Q(1,IVL)=0;               Q(3,IVL)=0;

%Constants for flow cases
if flow_case==1                   %Flow_case=1-----NO FLOW
    po=2117;                      %Entrance Pressure in lb/ft^2
    Mo=0;                         %Mach Number at Entrance=0
    To=531.2;                     %Temperature at entrance
    p=po*ones(1,IVL);              %Initial Pressure
    rho=p./((gamma-1)*Cv*To);     %Initial density
    u=zeros(1,IVL);                    %Initial Velocity
    a=sqrt(2*gamma*gamma_fac*Cv*To);
else
    po=2117;                      %Entrance Pressure in lb/ft^2
    Mo=0;                         %Mach Number at Entrance=0
    To=531.2;                     %Temperature at entrance
    p=po*ones(1,IVL);             %Initial Pressure
    p(IVL)=1400;                  %Instantaneously drop pressure
    rho=p./((gamma-1)*Cv*To);                %Initial density
    u=zeros(1,IVL);                              %Initial Velocity
    a=sqrt(2*gamma*gamma_fac*Cv*To);
end

%Set up conserved quantities
U(1,:)=rho;                       
U(2,:)=rho.*u;
U(3,:)=rho.*(Cv.*To+.5.*(u.^2));

figure;
set(gcf,'doublebuffer','on');

num_steps=1000;

for t=0:num_steps
    %Update Pressure first since it will be needed throughout
    p=(gamma-1)*U(1,:).*(U(3,:)./U(1,:)-.5*(U(2,:)./U(1,:)).^2);    %pressure
    
    %-------------------------------------------------------------STEGER-WARMING METHOD    
    if flux_method==1
      %Exit Boundary Condition
      
      %State and Flux vectors and Q vector
      U_b=(U(:,1:IVL-1)+U(:,2:IVL))./2;                           %Averages       
      F_msw=zeros(3,IVL-1);                                       %Backward Flux   Modified SW
      F_sw=zeros(3,IVL-1);                                        %Forward Flux    SW
      F=zeros(3,IVL-1);                                           %Modified Flux
      
      %Jabcobians-----MSW   
      A_b=zeros(3,3,IVL-1);
      
      A_b11=zeros(1,IVL-1);                                                        
      A_b12=ones(1,IVL-1);                                                      
      A_b13=zeros(1,IVL-1);                                                       
      
      A_b21=.5*(gamma-3).*U_b(2,:).^2./U_b(1,:).^2;                                
      A_b22=-(gamma-3).*U_b(2,:)./U_b(1,:);                            
      A_b23=Beta.*ones(1,IVL-1);                                                        
      
      A_b31=-gamma.*U_b(3,:).*U_b(2,:)./U_b(1,:).^2+Beta.*U_b(2,:).^3./U_b(1,:).^3;       
      A_b32=gamma.*U_b(3,:)./U_b(1,:)-3*Beta.*U_b(2,:).^2./(2.*U_b(1,:).^2);            
      A_b33=gamma.*U_b(2,:)./U_b(1,:);                                             
      
      A_b(1,1,:)=A_b11;             A_b(1,2,:)=A_b12;             A_b(1,3,:)=A_b13;
      A_b(2,1,:)=A_b21;             A_b(2,2,:)=A_b22;             A_b(2,3,:)=A_b23;
      A_b(3,1,:)=A_b31;             A_b(3,2,:)=A_b32;             A_b(3,3,:)=A_b33;
      
      %Jacobians------SW
      A=zeros(3,3,IVL);
      
      A11=zeros(1,IVL);                                                        
      A12=ones(1,IVL);                                                      
      A13=zeros(1,IVL);                                                       
      
      A21=.5*(gamma-3).*U(2,:).^2./U(1,:).^2;                                
      A22=-(gamma-3).*U(2,:)./U(1,:);                            
      A23=Beta.*ones(1,IVL);                                                        
      
      A31=-gamma.*U(3,:).*U(2,:)./U(1,:).^2+Beta.*U(2,:).^3./U(1,:).^3;       
      A32=gamma.*U(3,:)./U(1,:)-3*Beta.*U(2,:).^2./(2.*U(1,:).^2);            
      A33=gamma.*U(2,:)./U(1,:);                                             
      
      A(1,1,:)=A11;             A(1,2,:)=A12;             A(1,3,:)=A13;
      A(2,1,:)=A21;             A(2,2,:)=A22;             A(2,3,:)=A23;
      A(3,1,:)=A31;             A(3,2,:)=A32;             A(3,3,:)=A33;
    
      
      %Flux Split the Jacobians
      for vol_ind=2:num_vols
        %Get A+ A- from data from i and i-1 for MSW
        [Lam_v_b,Lam_e_b]=eig(A_b(:,:,vol_ind-1));
        Lam_e_p_b=(Lam_e_b+abs(Lam_e_b))./2;
        Lam_e_m_b=(Lam_e_b-abs(Lam_e_b))./2;
        %Get A+ A- from data from i and i+1 for SW
        [Lam_v_i_1,Lam_e_i_1]=eig(A(:,:,vol_ind));
        Lam_e_m=(Lam_e_i_1-abs(Lam_e_i_1))./2;                                       
        [Lam_v_i,Lam_e_i]=eig(A(:,:,vol_ind-1));
        Lam_e_p=(Lam_e_i+abs(Lam_e_i))./2;
        
        %Now get Generic Fluxes   SW=Steger-warming and MSW=Modified SW
        F_sw(:,vol_ind-1)= Lam_v_i_1*Lam_e_m*Lam_v_i_1^-1*U(:,vol_ind)+...
                           Lam_v_i*Lam_e_p*Lam_v_i^-1*U(:,vol_ind-1);
        F_msw(:,vol_ind-1)=Lam_v_b*Lam_e_p_b*Lam_v_b^-1*U(:,vol_ind-1)+...
                           Lam_v_b*Lam_e_m_b*Lam_v_b^-1*U(:,vol_ind);
        %Use pressure gradient waiting to avoid going negative in pressure
        pg=(p(vol_ind)-p(vol_ind-1))/min(p(vol_ind),p(vol_ind-1));
        wt=1/(1+pg^2);
        F(:,vol_ind-1)=wt*F_msw(:,vol_ind-1)+(1-wt)*F_sw(:,vol_ind-1);
                
        %Find the minimum delx/u+c (i.e. the greatest eigenvalue for the CFL)
        if vol_ind==2
            lam_max=max(Lam_e_b(:));
            uc_min=delx/lam_max; 
        else
           if max(Lam_e_b(:)) > lam_max
            lam_max=max(Lam_e_b(:));
            uc_min=delx/lam_max;

           end
        end
       
   end
    
      %----------------------------------ROE METHOD
   else 
    
      %State and Flux vectors and Q vector
      U_b=(U(:,1:IVL-1)+U(:,2:IVL))./2;                           %Averages       
      F_msw=zeros(3,IVL-1);                                       %Backward Flux   Modified SW
      F_sw=zeros(3,IVL-1);                                        %Forward Flux    SW
      F=zeros(3,IVL-1);                                           %Modified Flux
      
      %Get the Roe Averaged Variables
      rho=U(1,:);             u=U(2,:)./U(1,:);         e=U(3,:);
      p=(gamma-1)*U(1,:).*(U(3,:)./U(1,:)-.5*(U(2,:)./U(1,:)).^2);    %pressure

      rho_hat=sqrt(rho(1:IVL-1)).*sqrt(rho(2:IVL));                         %RA density
      u_hat=(sqrt(rho(1:IVL-1)).*u(1:IVL-1)+sqrt(rho(2:IVL)).*u(2:IVL))./...  %RA velocity
            (sqrt(rho(1:IVL-1))+sqrt(rho(2:IVL)));
      h_hat=(sqrt(rho(1:IVL-1)).*(e(1:IVL-1)+p(1:IVL-1))./rho(1:IVL-1)+...
             sqrt(rho(2:IVL)).*(e(2:IVL)+p(2:IVL))./rho(2:IVL))./...
            (sqrt(rho(1:IVL-1))+sqrt(rho(2:IVL)));
      c_hat=sqrt((gamma-1).*(h_hat-.5*u_hat.^2));
      
      %Jabcobians-----MSW   
      A_h=zeros(3,3,IVL-1);
      
      A_h11=zeros(1,IVL-1);                                                        
      A_h12=ones(1,IVL-1);                                                      
      A_h13=zeros(1,IVL-1);                                                       
      
      A_h21=.5*(gamma-3).*u_hat.^2;                                
      A_h22=-(gamma-3).*u_hat;                            
      A_h23=Beta.*ones(1,IVL-1);                                                        
      
      A_h31=(gamma/2-1).*u_hat.^3-c_hat.^2.*u_hat./(gamma-1);       
      A_h32=(3/2-gamma).*u_hat.^2+c_hat.^2./(gamma-1);            
      A_h33=gamma.*u_hat;                                             
      
      A_h(1,1,:)=A_h11;             A_h(1,2,:)=A_h12;             A_h(1,3,:)=A_h13;
      A_h(2,1,:)=A_h21;             A_h(2,2,:)=A_h22;             A_h(2,3,:)=A_h23;
      A_h(3,1,:)=A_h31;             A_h(3,2,:)=A_h32;             A_h(3,3,:)=A_h33;
      
      %Jacobians------SW
      A=zeros(3,3,IVL);
      
      A11=zeros(1,IVL);                                                        
      A12=ones(1,IVL);                                                      
      A13=zeros(1,IVL);                                                       
      
      A21=.5*(gamma-3).*U(2,:).^2./U(1,:).^2;                                
      A22=-(gamma-3).*U(2,:)./U(1,:);                            
      A23=Beta.*ones(1,IVL);                                                        
      
      A31=-gamma.*U(3,:).*U(2,:)./U(1,:).^2+Beta.*U(2,:).^3./U(1,:).^3;       
      A32=gamma.*U(3,:)./U(1,:)-3*Beta.*U(2,:).^2./(2.*U(1,:).^2);            
      A33=gamma.*U(2,:)./U(1,:);                                             
      
      A(1,1,:)=A11;             A(1,2,:)=A12;             A(1,3,:)=A13;
      A(2,1,:)=A21;             A(2,2,:)=A22;             A(2,3,:)=A23;
      A(3,1,:)=A31;             A(3,2,:)=A32;             A(3,3,:)=A33;
    
      
      %Split the Flux Difference vectors
      for vol_ind=2:num_vols
        %Get A_h+ A_h- from data from i and i-1 for ROE
        [Lam_v_h,Lam_e_h]=eig(A_h(:,:,vol_ind-1));
        %Get eigenvalues of A to avoid shock-happiness
        Lam_e_i=eig(A(:,:,vol_ind));
        Lam_e_i_1=eig(A(:,:,vol_ind-1));
        
        %Check eigenvalues 
        for e_ind=1:3
            ei_half=sigma*max([0,...
                              Lam_e_h(e_ind,e_ind)-Lam_e_i_1(e_ind),...
                              Lam_e_i(e_ind)-Lam_e_h(e_ind,e_ind)]);
            if abs(Lam_e_h(e_ind,e_ind)) < ei_half
                Lam_e_h(e_ind,e_ind)=.5*((Lam_e_h(e_ind,e_ind))^2/ei_half+ei_half);
            end

        end
        max_e_val=max(max([max(Lam_e_h)',Lam_e_i(:),Lam_e_i_1(:)]))
        %Get the absolute value ROE averaged hat matrix
        Lam_e_p_h=(Lam_e_h+abs(Lam_e_h))./2;
        Lam_e_m_h=(Lam_e_h-abs(Lam_e_h))./2;
        abs_A_h=Lam_v_h*Lam_e_p_h*Lam_v_h^-1-...
                Lam_v_h*Lam_e_m_h*Lam_v_h^-1;
        %Now get the normal fluxes
        F_i=A(:,:,vol_ind)*U(:,vol_ind);
        F_i_1=A(:,:,vol_ind-1)*U(:,vol_ind-1);
        %And finally, get the ROE flux
        F(:,vol_ind-1)=(F_i+F_i_1)./2-.5*abs_A_h*(U(:,vol_ind)-U(:,vol_ind-1));
                
        %Find the minimum delx/u+c (i.e. the greatest eigenvalue for the CFL)
        if vol_ind==2
            lam_max=max_e_val;
            uc_min=delx/lam_max; 
        else
           if max_e_val > lam_max
            lam_max=max_e_val;
            uc_min=delx/lam_max;

           end
        end
       
      end
    
  end    
      %Calculate the time step using the 
      delt=0.9*uc_min;
    
      Q(2,2:IVL-1)=(1./Vol(:,2:IVL-1)).*(Surf(:,3:IL-1)-Surf(:,2:IL-2)).*...
                    p(2:IVL-1);
      U(:,2:IVL-1)=U(:,2:IVL-1)-...
                 (delt./Vols(:,2:IVL-1)).*...
                 (F(:,2:IVL-1).*Surfs(:,3:IL-1)-F(:,1:IVL-2).*Surfs(:,2:IL-2))+...
                 delt.*Q(:,2:IVL-1);
              
      %Exit Boundary Conditions--------Works
      [U_IVL]=Quasi_1D_exit_SI_bc(delt,delx,Cv,a,To,po,gamma,...
                         p(IVL),p(IVL-1),U(:,IVL),U(:,IVL-1),Surf(IVL),Surf(IVL-1));
      U(:,IVL)=U_IVL;
      %Entrance Boundary Condition-----
      p(1) 
      [delU]=Quasi_1D_entrance_SI_bc(delt,delx,Cv,a,To,po,gamma,...
                         p(1),p(2),U(:,1),U(:,2),Surf(1),Surf(2)); 
      U(:,1)=U(:,1)+delU;
      
      %Plot the conservative variables
      c=sqrt(gamma*p./U(1,:));
      subplot(2,1,1);
      plot(xvols,(U(2,:)./U(1,:))./c);
      set(gca,'ylim',[0,2]);
      %plot(xvols,U(2,:),xvols,U(1,:),xvols,U(3,:),'*')
      subplot(2,1,2);
      plot(xvols,(p/po))
      set(gca,'ylim',[0 2]);
      pause(0.03)
      if t==100
        pause;
      end
    
      
end
    
                 
                 
end