% CFL number
CFL = 0.90;

% Gas properties
gamma = 1.4;
Rgas  = 287.06;
  
% Inlet total pressure and temperature
P0inlet = 101327;
T0inlet = 288.15;

% Set outlet static pressure
poutlet = 0.9*P0inlet;

% Set grid size and geometry info
Nx = 20;
x  = linspace(0,4,Nx+1);
xc = 0.5*(x(1:Nx) + x(2:Nx+1));
dx = x(2) - x(1);

A0 = 10.0;
A1 =  8.0;
A2 = 12.0;

A = zeros(size(x));

for j = 1:Nx+1,
  xl = x(j);
  if (xl <= 1),
    A(j) = A0;
  elseif (xl <= 2),
    A(j) = A0 + 0.5*(A1-A0)*(1 - cos(pi*(xl-1)));
  elseif (xl <= 3),
    A(j) = A1 + 0.5*(A2-A1)*(1 - cos(pi*(xl-2)));
  else,
    A(j) = A2;
  end
end

% Allocate space solution and residual vectors
U  = zeros(3,Nx);
R  = zeros(3,Nx);

% Initialize flow to stagnation conditions

FILL THIS IN 

% Set initial iteration value
niter = 0;

while ( FILL IN CONVERGENCE CHECK BASED ON niter AND ON residual )
  
  % Set inlet states as follows (we assume a subsonic inlet condition):
  % 
  % P0 = P0inlet (i.e. prescribed)
  % T0 = T0inlet (i.e. prescribed)
  % u  = u(2)    (i.e. set velocity equal to velocity in 2nd cell)
  %
  % Given these three pieces of information, the state vector at the inlet
  % can be calculated.

  FILL THIS IN
  
  U(1,1) = ...
  U(2,1) = ...
  U(3,1) = ...
  
      
  % Set outlet states as follows (we assume a subsonic outlet condition):
  % To do this, we will assume that the states at the outlet are
  % given as follows:
  % 
  % Jp = Jp(Nx-1)
  % S  =  S(Nx-1)
  % p  = poutlet 
  %
  % Given these three pieces of information, the state vector at the outlet
  % can be calculated.
  
  FILL THIS IN
  
  U(1,Nx) = ...
  U(2,Nx) = ...
  U(3,Nx) = ...

      
  % Calculate residuals for all interior cells

  FILL THIS IN
  
  
  % Assign timestep
  
  FILL THIS IN
  
  
  % Update U for interior cells using Forward Euler

  FILL THIS IN

  
  % Update iteration counter and residual 

  FILL THIS IN
  
  
  
  % Post-process and plot solution
  p   = (gamma-1)*(U(3,:) - 0.5*U(2,:).^2/U(1,:));
  rho = U(1,:);
  u   = U(2,:)./rho;
  a   = sqrt(gamma*p./rho);
  M   = u./a;
  
  subplot(211);
  plot(xc,M);hold on;
  plot(xc,ones(size(M)),'--');hold off;
  xlabel('x');ylabel('M');
    
  subplot(212);
  plot(log10(Res(min(niter,max(2,niter-199)):niter)));
  xlabel('iters');ylabel('log10 Res');
  drawnow;
  
end

% Final plot
subplot(211);plot(xc,M);hold on;
plot(xc,ones(size(M)),'--');hold off;
axis([min(xc),max(xc),0,2]);
subplot(212);plot(log10(Res));