function Aw_plus_g = Advection1(t, w, h, a, d, method) % Evaluate rhs of ODE w'=A*w+g(t)
% method=0 upwind, method=1 centered, method=2 reduced-order third-order upwind biased
% method=3 third-order upwind biased (most accurate)
% Finite volume spatial discretization to solve:
% u_t+a*u_x=d*u_xx
% with Dirichlet BCs on left and Neumann on right:
% u(0,t)=sin(-pi*t)^exponent;
% u_x(1,t)=0 if d>0
% If there were only advection the solution would be u(x,t)=sin(pi*(x-t))^exponent

global exponent;
w0=sin(pi*(0-t))^exponent; % Dirichlet BC at left

n=length(w); % Unfortunately MATLAB does not allow zero or negative indexing
% So here I treat the ghost cells separately and not with array syntax

% Fluxes start on left-most and end on right-most face
fluxes=zeros(n+1,1); % MATLAB requires column vectors for ode solvers

% Compute advective fluxes first, ignoring BCs at outflow boundary
fluxes(1) = a*w0; % Left-most face has Dirichlet BC so flux is known!
switch abs(method)
   case 0 % upwind

      fluxes(2:n+1) = a*w(1:n);

   case 1 % centered in interior, upwind at right boundary

      fluxes(2:n)=0.5*a*(w(1:n-1)+w(2:n));
      
      % This actually corresponds to upwind on the right boundary:
      wnp1=2*w(n)-w(n-1); % One-sided linear extrapolation
      fluxes(n+1) = 0.5*a*(w(n)+wnp1);

   case 2 % Third-order upwind biased in the interior

      fluxes(3:n) = a*(-w(1:n-2)+5*w(2:n-1)+2*w(3:n))/6;

      % Left boundary:
      % Linear extrapolation at the Dirichlet boundary
      % This is now second-order at the boundary (as finite-difference) and gives
      % fluxes(2)=-(1/3)*w0+w(1)+(1/3)*w(2)
      % dw(1)/dt=-(-4/3*w0+w(1)+w(2)/3)/h
      % which is now a centered discretization of advection with one staggered value (w0)
      wm1=2*w0-w(1); 
      fluxes(2) = a*(-wm1+5*w(1)+2*w(2))/6;      
      
      % Testing:
      % My Maple analysis shows the above to be second-order only for finite-difference
      % but not for finite-volume. What I get for finite-volume is
      % fluxes(2) = a*(-4*w0+5*w(1)+w(2))/2;
      % but this does not seem to work at all ?!?

      % Right boundary:      
      % Extrapolate quadratically without using the BC
      % Observe that this is now actually doing second-order upwind at the right boundary
      % dw(n)/dt=a*(w(n-2)-4*w(n-1)+3*w(n))/(2*h)
      % So this is second-order accurate at the boundary as we want it to be
      wnp1=3*w(n)-3*w(n-1)+w(n-2); % One-sided quadratic extrapolation
      fluxes(n+1) = a*(-w(n-1)+5*w(n)+2*wnp1)/6;

   otherwise
      error('method=%d not implemented yet',method);  
end

if(d>0) % Add diffusive fluxes
   % In the interior we use the usual three-point centered difference:
   fluxes(2:n)=fluxes(2:n) - d*(w(2:n)-w(1:n-1))/h;

   fluxes(1)=fluxes(1) - d*(w(1)-w0)/(0.5*h); % Dirichlet BC at left
   % Neumann BC at right means no diffusive contribution to fluxes(n+1)   
end

Aw_plus_g =  (fluxes(1:n)-fluxes(2:n+1))/h;

end
