%  #    # #    # #    # ###### #####  #  ####    ##   #      
%  ##   # #    # ##  ## #      #    # # #    #  #  #  #      
%  # #  # #    # # ## # #####  #    # # #      #    # #      
%  #  # # #    # #    # #      #####  # #      ###### #      
%  #   ## #    # #    # #      #   #  # #    # #    # #      
%  #    #  ####  #    # ###### #    # #  ####  #    # ###### 
%                                      
%  ###### #      #    # #    # ######  ####  
%  #      #      #    #  #  #  #      #      
%  #####  #      #    #   ##   #####   ####  
%  #      #      #    #   ##   #           # 
%  #      #      #    #  #  #  #      #    # 
%  #      ######  ####  #    # ######  ####

function [varargout] = StandardFlux( what, varargin )

  numerical_flux_name = lower(varargin{1}) ;

  switch what
    case 'nflux'
      Flux = varargin{2} ;
      DT   = varargin{3} ;
      DX   = varargin{4} ;
      U    = varargin{5} ;
      %
      switch numerical_flux_name
      case 'laxfriedrichs'
        varargout{1} = LaxFriedrichs( Flux, DT, DX, U ) ;
      case 'force'
        varargout{1} = Force( Flux, DT, DX, U ) ;
      case 'godunovcentered'
        varargout{1} = GodunovCentered( Flux, DT, DX, U ) ;
      case 'richtmyer'
        varargout{1} = Richtmyer( Flux, DT, DX, U ) ;
      otherwise
        error( [ 'StandardFlux nflux: ', varargin{1}, ' unknown' ] ) ;
      end

    case 'get_cfl'
      switch numerical_flux_name
      case { 'laxfriedrichs', 'force', 'richtmyer' }
        varargout{1} = 1.0 ;
      case { 'slic', 'godunovcentered' }
        varargout{1} = 0.5 ;
      otherwise
        error( [ 'StandardFlux nflux: ', numerical_flux_name, ' unknown' ] ) ;
      end
  end
end
%   _                    _____     _          _      _      _         
%  | |    __ ___  __    |  ___| __(_) ___  __| |_ __(_) ___| |__  ___ 
%  | |   / _` \ \/ /____| |_ | '__| |/ _ \/ _` | '__| |/ __| '_ \/ __|
%  | |__| (_| |>  <_____|  _|| |  | |  __/ (_| | |  | | (__| | | \__ \
%  |_____\__,_/_/\_\    |_|  |_|  |_|\___|\__,_|_|  |_|\___|_| |_|___/
%
% Compute the Lax-Friedrichs flux
%
function F = LaxFriedrichs( FLUX, DT, DX, U )
  UL = U(:,2:end-2) ;
  UR = U(:,3:end-1) ;
  %%
  FL = feval( FLUX, UL) ;
  FR = feval( FLUX, UR) ;
  F  = 0.5 .* ( FL + FR + (DX/DT) .* (UL-UR) ) ;
end
%   _____                  
%  |  ___|__  _ __ ___ ___ 
%  | |_ / _ \| '__/ __/ _ \
%  |  _| (_) | | | (_|  __/
%  |_|  \___/|_|  \___\___|
%                        
% compute an intercell flux according to the FORCE scheme.
%
function F = Force( FLUX, DT, DX, U )
  UL = U(:,2:end-2) ;
  UR = U(:,3:end-1) ;
  %%
  FL = feval( FLUX, UL) ;
  FR = feval( FLUX, UR) ;
  % Richtmyer
  US  = 0.5 .* ( UL + UR + (DT/DX) .* (FL-FR) ) ;
  FRM = feval( FLUX, US ) ;
  % LaxFriedrichs
  FLF = 0.5 .* ( FL+FR + (DX/DT) .* (UL-UR) ) ;
  F   = 0.5 .* ( FRM + FLF ) ;
end
%    ____           _                         ____           _                    _ 
%   / ___| ___   __| |_   _ _ __   _____   __/ ___|___ _ __ | |_ ___ _ __ ___  __| |
%  | |  _ / _ \ / _` | | | | '_ \ / _ \ \ / / |   / _ \ '_ \| __/ _ \ '__/ _ \/ _` |
%  | |_| | (_) | (_| | |_| | | | | (_) \ V /| |__|  __/ | | | ||  __/ | |  __/ (_| |
%   \____|\___/ \__,_|\__,_|_| |_|\___/ \_/  \____\___|_| |_|\__\___|_|  \___|\__,_|
%
% compute an intercell flux according to the Godunov centred scheme (non-monotone).
% Stability: 0 < CFL Coefficient < 0.5 * sqrt(2), about 0.7
%
function F = GodunovCentered( FLUX, DT, DX, U )
  UL = U(:,2:end-2) ;
  UR = U(:,3:end-1) ;
  %%
  FL = feval( FLUX, UL ) ;
  FR = feval( FLUX, UR ) ;
  US = 0.5 .* (UL+UR) + (DT/DX) .* ( FL - FR ) ;
  % Compute the flux at the state US
  F  = feval( FLUX, US ) ;
end
%   ____  _      _     _                             
%  |  _ \(_) ___| |__ | |_ _ __ ___  _   _  ___ _ __ 
%  | |_) | |/ __| '_ \| __| '_ ` _ \| | | |/ _ \ '__|
%  |  _ <| | (__| | | | |_| | | | | | |_| |  __/ |   
%  |_| \_\_|\___|_| |_|\__|_| |_| |_|\__, |\___|_|
%                                    |___/           
% Compute intercell fluxes according to the Richtmyer method, or two step Lax-Wendroff method
% See Eq. 5.79, Chapter 5, Ref. 1
%
function F = Richtmyer( FLUX, DT, DX, U )
  UL = U(:,2:end-2) ;
  UR = U(:,3:end-1) ;
  %%
  FL = feval( FLUX, UL) ;
  FR = feval( FLUX, UR) ;
  US = 0.5 * ( UL + UR + (DT/DX) * (FL-FR) ) ;
  F  = feval( FLUX, US ) ;
end
