% Start timer
tic;

% Check if restart flag exists
if (exist('Q1D_restart') == 0),
  Q1D_restart = 0;
end

% CFL number
CFL = 0.90;

if (Q1D_restart == 0),
  
  % Gas properties
  gamma = 1.4;
  Rgas  = 287.06;
  cv    = Rgas/(gamma - 1);
  cp    = cv*gamma;
  
  % Inlet total pressure and temperature
  P0inlet = 101327;
  T0inlet = 288.15;
  poutlet = P0inlet*0.9;
  
  % Maximum number of iterations
  maxiter = 100000;

  % Number of iterations per plot
  nplot = 1;

  % Set grid size and geometry info
  Nx = 40;
  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
  
  for j = 1:Nx,
    dA(j) = A(j+1) - A(j);
    Vol(1,j) = 1.0/(0.5*(A(j+1)+A(j))*dx);
    Vol(2,j) = 1.0/(0.5*(A(j+1)+A(j))*dx);
    Vol(3,j) = 1.0/(0.5*(A(j+1)+A(j))*dx);
  end
  
  % Allocate space solution and residual vectors
  U  = zeros(3,Nx);
  R  = zeros(3,Nx);
  Fh = zeros(3,1);

  % Initialize flow to stagnation conditions
  pinit   = P0inlet;
  Tinit   = T0inlet;
  uinit   = 0;
  rhoinit = pinit/(Rgas*Tinit);
  
  U(1,:) = rhoinit;
  U(2,:) = rhoinit*uinit;
  U(3,:) = pinit/(gamma-1) + 0.5*rhoinit*uinit^2;
  
  % Initialize niter, restol, curres, and maxres
  niter = 0;
  restol = 1e-4;
  maxres = 0;
  curres = 1;
end

while ((niter < maxiter) & (curres > maxres*restol)),
  % Zero residual
  R(:,:) = 0;
  
  % 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 flux at the inlet
  % can be calculated.

  j = 2;
  P0I  = P0inlet;
  T0I  = T0inlet;
  uI   = U(2,j)/U(1,j);
  a2I  = gamma*Rgas*T0I - 0.5*(gamma-1)*uI^2;
  M2I  = uI^2/a2I;
  pI   = P0I*(1+0.5*(gamma-1)*M2I)^(-gamma/(gamma-1));
  rhoI = gamma*pI/a2I;
  
  U(1,1) = rhoI;
  U(2,1) = rhoI*uI;
  U(3,1) = rhoI*(cp*T0I - pI/rhoI);
  
  % Set outlet states as follows (we assume a subsonic outlet condition):
  % Calculate flux at outlet face and distribute to cell residual
  % To do this, we will assume that the states at the inlet are
  % given as follows:
  % 
  % P0 = P0(Nx-1) (i.e. set equal to total pressure in last cell)
  % T0 = T0(Nx-1) (i.e. set equal to total temperature in last cell)
  % p  = poutlet  (i.e. prescribed)
  %
  % Given these three pieces of information, the flux at the outlet
  % can be calculated.
  
  j   = Nx-1;
  pj  = (gamma-1)*(U(3,j) - 0.5*U(2,j)^2/U(1,j));
  a2j = gamma*pj/U(1,j);
  uj  = U(2,j)/U(1,j);
  JpF = uj + 2*sqrt(a2j)/(gamma-1);
  sF  = pj/U(1,j)^(gamma);
  
  pF  = poutlet;
  
  rhoF = (pF/sF)^(1/gamma);
  a2F  = gamma*pF/rhoF;
  uF   = JpF - 2*sqrt(a2F)/(gamma-1);
  HF   = a2F/(gamma-1) + 0.5*uF^2;

  U(1,Nx) = rhoF;
  U(2,Nx) = rhoF*uF;
  U(3,Nx) = rhoF*(HF - pF/rhoF);

  % Calculate flux at interior faces and distribute to cell residuals
  for j = 2:Nx,
    Fh = Fupwind(U(:,j-1),U(:,j),gamma);
    R(:,j-1) = R(:,j-1) + Fh*A(j);
    R(:,j  ) = R(:,j  ) - Fh*A(j);
  end
  
  % Add source term to each cell residual  
  pvec   = (gamma-1)*(U(3,:) - 0.5*U(2,:).^2./U(1,:));
  R(2,:) = R(2,:) - pvec.*dA;
  
  % Assign timestep
  cvec = sqrt(gamma*pvec./U(1,:));
  lambdamax = max(abs(U(2,:)./U(1,:)) + cvec);  
  dt = CFL*dx/lambdamax;
  niter = niter + 1;
  curres = sqrt(sum((Vol(1,2:Nx-1).*R(1,2:Nx-1)).^2)*dx);
  Res(niter) = curres;  
  if (curres > maxres),
    maxres = curres;
  end
  
  % Update U using Forward Euler (ignoring update to first and last cells)
  U(:,2:Nx-1) = U(:,2:Nx-1) - dt*R(:,2:Nx-1).*Vol(:,2:Nx-1);

  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;
  P0    = p.*(1 + 0.5*(gamma-1)*M.^2).^(gamma/(gamma-1));
  Fout  = Fupwind(U(:,Nx-1),U(:,Nx),gamma);

  P0out(niter) = P0(Nx);
  mdot(niter)  = Fout(1)*A(Nx);

  if (mod(niter,nplot) == 0),
    subplot(221);
    plot(xc,M);
    xlabel('x');ylabel('M');
    
    subplot(222);
    plot(log10(Res(min(niter,max(2,niter-199)):niter)));
    xlabel('iters');ylabel('log10 Res');

    subplot(223);
    plot(P0out(min(niter,max(2,niter-199)):niter));
    xlabel('iters');ylabel('P0outlet');

    subplot(224);
    plot(mdot(min(niter,max(2,niter-199)):niter));
    xlabel('iters');ylabel('mdot outlet');
    
    drawnow;
  end
  
end

% Stop timer 
fprintf('Elapsed time = %f\n',toc);
fprintf('Iters        = %i\n',niter);

% Final plot
subplot(221);
plot(xc,M);
xlabel('x');ylabel('M');

subplot(222);
plot(log10(Res));
xlabel('iters');ylabel('log10 RMSres');

subplot(223);
plot(P0out);
xlabel('iters');ylabel('P0outlet');

subplot(224);
plot(mdot);
xlabel('iters');ylabel('mdot outlet');

fprintf('Outlet P0    = %f\n',P0out(niter));
fprintf('Outlet mdot  = %f\n',mdot(niter));