% Estimate advective flux difference div(a*u)^{n+1/2} at midpoint in time to second order
% This uses Fromm's method including a constant-in time source term
function w_t = Fromm(w, a, dt, h, source, n, t, DBC)
  
   global periodic;
   global exponent;
   
   n=length(w);
        
   cfl = a*dt/h;
      
   if(0 && periodic) % Try limiters

      %phi = @(theta) (1+theta)/2; % Fromm method rewritten with a "limiter"
      phi = @(theta) max(0, min(min((1+theta)/2,2),2*theta));
      my_sign = @(x) sign(sign(x)+10*eps); % Returns either +1 or -1 but never zero
   
      % Limited version:
      %theta = (w-circshift(w, 1)) ./ (circshift(w, -1)-w) ; % Division by zero!
      % Make sure we never divide by zero:
      slopes = circshift(w, -1)-w;
      slopes = max(abs(slopes),10*eps) .* my_sign(slopes); % Get away from zero by 1e-15
      theta = (w-circshift(w, 1)) ./ slopes;
      
      fluxes = a*(w + 0.5*(1-cfl)*phi(theta).*slopes) + 0.5*a*dt*source;
      w_t = (circshift(fluxes, 1) - fluxes)/h ;
          
   elseif(periodic) % Here we number the fluxes so that cell i is upwind from face i
   
      if(1) % Fromm's method

         fluxes = a*(w + 0.25*(1-cfl)*(circshift(w, -1)-circshift(w, 1)) + 0.5*dt*source);
         w_t = (circshift(fluxes, 1) - fluxes)/h ;
      
      else % Lax-Wendroff
      
         fluxes = a*(w + 0.5*(1-cfl)*(circshift(w, -1)-w) + ...
                     0.25*dt*(source+circshift(source,-1)));
         w_t = (circshift(fluxes, 1) - fluxes)/h ;
            
      end   
         
   else % Here we number faces so that the first face is #1
        % This means cell i is now downwind from face i

      fluxes = zeros(n+1,1);

      % Interior (Fromm's method = linear reconstruction with centered slopes)  
      %----------------
      fluxes(3:n) = a*(w(2:n-1)  + 0.25*(1-cfl)*(w(3:n)-w(1:n-2)) + 0.5*dt*source(2:n-1));
            
      % Left boundary (inflow):
      %----------------      
      w0_beg=DBC(t); % Dirichlet BC at the beginning of the time step
      w0_mid=DBC(t+dt/2); % Dirichlet BC at left at midpoint
      
      fluxes(1) =   a*w0_mid; % Make sure the BC is read at the *midpoint* in time
      
      % For the second face:
      if(0) % Estimate the slope by a centered difference
         % This is second order but with boundary layer in error
         fluxes(2) =   a*(w(1) + 0.50*(1-cfl)*(w(2)-w(1))     + 0.5*dt*source(1));
      else % estimate the slope by a one-sided difference to the boundary (works great)
         % A good way to think of this is as a linear reconstruction that enforces the BC at x=0
         % Remember the BC here is evaluated is at *beginning* of time step
         fluxes(2) =   a*(w(1) + 1.00*(1-cfl)*(w(1)-w0_beg)   + 0.5*dt*source(1)); 
      end   
      
      % Right boundary (outflow):
      %----------------
      if(0) % For the last face estimate slope upwind      
         % This is second-order but has a boundary layer in the error so not perfect
         fluxes(n+1) = a*(w(n)      + 0.50*(1-cfl)*(w(n)-w(n-1))     + 0.5*dt*source(n));
      else  % Use beam-warming for last cell (works great)
         % Beam warming is linear reconstruction with upwinded slopes:
         w_t(n) = -a/(2*h)*(3*w(n)-4*w(n-1)+w(n-2)) + 0.5*dt*(a/h)^2*(w(n)-2*w(n-1)+w(n-2));
         % You can think of this as setting fluxes(n+1) to whatever it needs to be
         fluxes(n+1) = fluxes(n) - w_t(n)*h; % Fix the flux      
      end   
      
      % Flux difference
      %----------------
      w_t = -(fluxes(2:n+1)-fluxes(1:n))/h;
                  
   end

end
