%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
finite_vol=1;                   %1=finite volume        2=finite difference
flux_method=1;                  %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
    press=po*ones(1,IVL);              %Initial Pressure
    rho=press./((gamma-1)*Cv*To);     %Initial density
    u=zeros(1,IVL);                    %Initial Velocity
    a=sqrt(2*gamma*(gamma-1)/(gamma-1)*Cv*To);
else
    po=2117;                      %Entrance Pressure in lb/ft^2
    Mo=0;                         %Mach Number at Entrance=0
    To=531.2;                     %Temperature at entrance
    press=po*ones(1,IVL);         %Initial Pressure
    press(IVL)=1400;              %Instantaneously drop pressure
    rho=press./((gamma-1)*Cv*To);                %Initial density
    u=zeros(1,IVL);                              %Initial Velocity
    a=sqrt(2*gamma*(gamma-1)/(gamma-1)*Cv*To);
end

%Set up conserved quantities
U(1,:)=rho;                       
U(2,:)=rho.*u;
U(3,:)=rho.*(Cv.*To+.5.*(u.^2));
plot(press)
num_steps=300;
for t=0:num_steps
    
    %Boundary Conditions
    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_p=zeros(3,IVL-2);                                      %Forward Flux
      F_m=zeros(3,IVL-2);                                      %Backward Flux
   
      %Jabcobiabn      
      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;
      
    end
      
    %Flux Split the Jacobians
    for vol_ind=2:num_vols-1
        %Get A+ A- from data from i and i-1
        [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
        [Lam_v_f,Lam_e_f]=eig(A_b(:,:,vol_ind));
        Lam_e_p_f=(Lam_e_f+abs(Lam_e_f))./2;
        Lam_e_m_f=(Lam_e_f-abs(Lam_e_f))./2;
        %Now get Generic Fluxes   F_p=F_i+1/2       F_m=F_i-1/2
        F_p(:,vol_ind-1)=Lam_v_f*Lam_e_p_f*Lam_v_f^-1*U(:,vol_ind)+...
                         Lam_v_f*Lam_e_m_f*Lam_v_f^-1*U(:,vol_ind+1);
        F_m(:,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);
                   
        
        %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_f(:)) > lam_max
            lam_max=max(Lam_e_f(:));
            uc_min=delx/lam_max;

           end
        end
       
    end
    
    %Calculate the time step using the 
    delt=0.9*uc_min;
    
    U
    p=(gamma-1)*U(1,:).*(U(3,:)./U(1,:)-.5*U(2,:).^2);     %pressure
    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_p.*Surfs(:,3:IL-1)-F_m.*Surfs(:,2:IL-2))+...
                  delt.*Q(:,2:IVL-1);
    p(IVL)
    pause;
   
    %Exit Boundary Conditions--------Works
    [U_IVL]=Quasi_1D_exit_SI_bc(delt,delx,Cv,a,To,po,gamma,...
                         U(2,IVL),U(2,IVL-1),U(1,IVL),U(1,IVL-1),Surf(IVL),Surf(IVL-1));
    U(:,IVL)=U_IVL;
    %Entrance Boundary Condition-----
    [delU]=Quasi_1D_entrance_SI_bc(delt,delx,Cv,a,To,po,gamma,...
                         U(2,1),U(2,2),U(1,1),U(1,2),Surf(1),Surf(2)); 
    U(:,1)=U(:,1)+delU;
  
    %Plot the conservative variables
    %plot(xvols,U(2,:),xvols,U(1,:),xvols,U(3,:),'*')
    plot(p)
  
    
end
    
                 
                 
