%  solve the linearised equations of Gas Dynamics in one space dimension 
%
%   u(t,x) + c u (t,x) = 0
%   t           x
%
%   for t > 0 and -1 <= x <= 1
%
function [varargout] = LinearAdvectionModel( what, varargin )

  % LinearAdvectionModel_paramaters
  %   c        = velocity
  %   test     = which test is used
  %   nflux    = which numerical flux is used
  %   bc_left  = left bounday condition
  %   bc_right = right boundaruy condition
  %

  global LinearAdvectionModel_paramaters ;

  switch lower(what)

    case 'setup'
      expected( what, 1, nargin, 0, nargout ) ;
       LinearAdvectionModel_paramaters = clear_struct(LinearAdvectionModel_paramaters) ;
       LinearAdvectionModel_paramaters . speed = varargin{1} ;

    case 'set_test'
      expected( what, 1, nargin, 0, nargout ) ; % TEST
      ChooseTest( varargin{1} ) ;

    case 'set_limiter'
      expected( what, 1, nargin, 0, nargout ) ;
      LinearAdvectionModel_paramaters . limiter  = lower( varargin{1} ) ;

    case 'set_numerical_flux'
      expected( what, 1, nargin, 0, nargout ) ;
      LinearAdvectionModel_paramaters . nflux = lower( varargin{1} ) ;

    case 'get_dim'
      expected( what, 0, nargin, 1, nargout ) ;
      varargout{1} = 1 ;    

    case 'get_dx'
      expected( what, 0, nargin, 1, nargout ) ;
      varargout{1} = LinearAdvectionModel_paramaters . DX ;    

    case 'get_x'
      expected( what, 0, nargin, 1, nargout ) ;
      varargout{1} = LinearAdvectionModel_paramaters . X ;    

    case 'get_final_time'
      expected( what, 0, nargin, 1, nargout ) ;
      varargout{1} = LinearAdvectionModel_paramaters . tf ;    

    case 'get_cfl'
      expected( what, 0, nargin, 1, nargout ) ;
      varargout{1} = MaxCFL() ;    

    case 'names'
      expected( what, 0, nargin, 1, nargout ) ;
      varargout{1} = 'u' ;    

    case 'flux'
      expected( what, 1, nargin, 1, nargout ) ; % U
      varargout{1} = Flux( varargin{1} ) ;

    case 'init'
      expected( what, 1, nargin, 1, nargout ) ; % N
      varargout{1} = Init( varargin{1} ) ;

    case 'exact'
      expected( what, 2, nargin, 2, nargout ) ; % TIME, N 
      [ varargout{1} varargout{2} ] = Exact( varargin{1}, varargin{2}  ) ;

    case 'muscle_nflux'
      expected( what, 2, nargin, 1, nargout ) ; % UL, UR
      varargout{1} = Godunov( varargin{1}, varargin{2} ) ;    

    case 'conservative_to_physical'
      varargout{1} = varargin{1} ;

    case 'physical_to_conservative'
      varargout{1} = varargin{1} ;
      
    case 'conservative_to_characteristic'
      varargout{1} = varargin{1} ;

    case 'characteristic_to_conservative'
      varargout{1} = varargin{1} ;

    case 'dt_from_cfl'
      expected( what, 2, nargin, 1, nargout ) ; % U,  CFLcoeff 
      varargout{1} = DT_from_CFL( varargin{1}, varargin{2} ) ;

    case 'extrapolate_for_bc'
      expected( what, 1, nargin, 4, nargout ) ; % U
      [ varargout{1}, varargout{2}, varargout{3}, varargout{4} ] = Extrapolate_for_BC( varargin{1} ) ;

    case 'nflux'
      expected( what, 2, nargin, 1, nargout ) ; % DT, U
      varargout{1} = NFlux( varargin{1}, varargin{2} ) ;

    otherwise
      error(['LinearAdvectionModel bad argument: ', what]) ;
  end
end

function SE = clear_struct(S)
  if isstruct(S)
    SE = rmfield( S, fieldnames(S) ) ;
  else
    SE = S ;
  end
end

function expected( what, ein, nin, eout, nout )
  if ein ~= nin-1 
    error( ['in LinearAdvectionModel call to ', what, ' expected ', int2str( ein ), ' input argument, found ', int2str(nin-1) ]) ;
  end
  if eout ~= nout 
    error( ['in LinearAdvectionModel call to ', what, ' expected ', int2str( eout ), ' output argument, found ', int2str(nout) ]) ;  
  end
end

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

% compute the flux
function F = Flux( U )
  global LinearAdvectionModel_paramaters ;
  speed  = LinearAdvectionModel_paramaters . speed ;
  F      = speed .* U ;
end

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

function U = Init(N)
  global LinearAdvectionModel_paramaters ;
  % Calculate mesh size DX
  DX = 2.0/N ;
  X  = [-1+DX/2:DX:1-DX/2] ;
  
  if isfield(LinearAdvectionModel_paramaters, 'test')
    switch LinearAdvectionModel_paramaters . test
    case 'classroomproblem'
      DX = 8.0/N ;
      X  = [0+DX/2:DX:8-DX/2] ; % center of the cell
      N1 = floor(3*N/8) ;
      N2 = N - N1 ;
      U  = [ones(1,N1) zeros(1,N2) ] ;
    case 'square'
      N3 = floor(N/3) ;
      M  = N - N3 - N3 ;
      U  = [zeros(1,N3) ones(1,M) zeros(1,N3) ] ;
    case 'smooth'
      U  = exp( -8 .* X .^ 2 ) ;
        otherwise
      error(['LinearAdvectionModel::Init unknown test ', LinearAdvectionModel_paramaters . test]) ;
    end
  else
    error('LinearAdvectionModel::Init unknown test you must set it by set_test!') ;      
  end
  LinearAdvectionModel_paramaters . DX = DX ;
  LinearAdvectionModel_paramaters . X  = X ;
end


%  #####                                     #######                     
% #     # #    #  ####   ####   ####  ######    #    ######  ####  ##### 
% #       #    # #    # #    # #      #         #    #      #        #   
% #       ###### #    # #    #  ####  #####     #    #####   ####    #   
% #       #    # #    # #    #      # #         #    #           #   #   
% #     # #    # #    # #    # #    # #         #    #      #    #   #   
%  #####  #    #  ####   ####   ####  ######    #    ######  ####    #
%
function ChooseTest(test)

  global LinearAdvectionModel_paramaters ;

  switch test 
    case 'classroomproblem'
      tf       = 3.5 ;
      bc_left  = 'free' ;
      bc_right = 'free' ;
    case 'square'
      tf       = 4 ;
      bc_left  = 'cyclic' ;
      bc_right = 'cyclic' ;
    case 'smooth'
      tf       = 4 ;
      bc_left  = 'cyclic' ;
      bc_right = 'cyclic' ;
   otherwise
      error(['LinearAdvectionModel:ChooseTest test name',test] ) ;
  end

  LinearAdvectionModel_paramaters . test     = test ;
  LinearAdvectionModel_paramaters . bc_left  = bc_left ;
  LinearAdvectionModel_paramaters . bc_right = bc_right ;
  LinearAdvectionModel_paramaters . tf       = tf ;

end



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

function [X,U] = Exact(TIME,N)
  global LinearAdvectionModel_paramaters ;
  speed = LinearAdvectionModel_paramaters . speed ;

  % Calculate mesh size DX
  DX = 2.0/N ;
  X  = [-1+DX/2:DX:1-DX/2] ;
  X  = X - speed*TIME ; % translate the mesh

  % put the solution in [-1..1]
  while min(X) < -1
    IDX = find( X < -1 ) ;
    X(IDX) = X(IDX)+2 ;
  end

  while max(X) > 1
    IDX = find( X > 1 ) ;
    X(IDX) = X(IDX)-2 ;
  end
  
  switch LinearAdvectionModel_paramaters . test
  case 'square'
    U      = zeros(1,N) ;
    IDX    = find( X >= -1.0/3.0 & X <= 1.0/3.0 ) ;
    U(IDX) = 1 ;
  case 'smooth'
    % use gauss quadrature to estimate status
    XL = X + DX/sqrt(3)/2;
    XR = X - DX/sqrt(3)/2 ;
    U  = (exp( -8 .* XL .^ 2 )+exp( -8 .* XR .^ 2 ))/2 ;
  end
  % ripristina la mesh
  X  = [-1+DX/2:DX:1-DX/2] ;
end

%  ######  #######                                              #####  ####### #       
%  #     #    #            ###### #####   ####  #    #         #     # #       #       
%  #     #    #            #      #    # #    # ##  ##         #       #       #       
%  #     #    #            #####  #    # #    # # ## #         #       #####   #       
%  #     #    #            #      #####  #    # #    #         #       #       #       
%  #     #    #            #      #   #  #    # #    #         #     # #       #       
%  ######     #            #      #    #  ####  #    #          #####  #       ####### 
%                  #######                             #######                         
% Compute DT from CFL
function DT = DT_from_CFL( U, CFLcoeff )
  global LinearAdvectionModel_paramaters ;
  speed  = LinearAdvectionModel_paramaters . speed ;
  DX     = LinearAdvectionModel_paramaters . DX ;
  DT     = CFLcoeff * DX / abs(speed) ;
end

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

% Extrapolate values for boundary condition
function [ULL, UL, UR, URR] = Extrapolate_for_BC( U )
  global LinearAdvectionModel_paramaters ;
  L = LinearAdvectionModel_paramaters . bc_left ;
  R = LinearAdvectionModel_paramaters . bc_right ;
  % Left boundary
  switch L
  case 'free'   % Apply transmissive boundary conditions
    UL  = U(:,1) ;
    ULL = U(:,2) ;
  case 'cyclic' % Apply cyclic boundary conditions
    UL  = U(:,end) ;
    ULL = U(:,end-1) ;
  otherwise     % Apply reflective boundary conditions
    UL  = [ U(1,1) ; -U(2,1) ] ;
    ULL = [ U(1,2) ; -U(2,2) ] ;
  end
  % Right boundary
  switch R
  case 'free'   % Apply transmissive boundary conditions
    UR  = U(:,end) ;
    URR = U(:,end-1) ;
  case 'cyclic' % Apply cyclic boundary conditions
    UR  = U(:,1) ;
    URR = U(:,2) ;
  otherwise    % Apply reflective boundary conditions
    UR  = [ U(1,end)   ; -U(2,end) ] ;
    URR = [ U(1,end-1) ; -U(2,end-1) ] ;
  end
end

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

function CFL = MaxCFL
  global LinearAdvectionModel_paramaters ;
  nflux = LinearAdvectionModel_paramaters . nflux ;
  %
  switch nflux
  case 'torobillett'
    CFL = 2.0 ;
  case { 'fromm', 'warmingbeam', 'waf', 'godunov' }
    CFL = 1.0 ;
  otherwise
    CFL = StandardFlux( 'get_cfl', nflux ) ;
  end
end


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

function FLUX = NFlux( DT, U )

  global LinearAdvectionModel_paramaters ;

  ULL = U(:,1:end-3) ;
  UL  = U(:,2:end-2) ;
  UR  = U(:,3:end-1) ;
  URR = U(:,4:end-0) ;

  nflux   = LinearAdvectionModel_paramaters . nflux ;
  DX      = LinearAdvectionModel_paramaters . DX ;
  limiter = LinearAdvectionModel_paramaters . limiter ;

  switch nflux
  % monolitic numerical fluxes
  case 'godunov'
    FLUX = Godunov(UL,UR) ;
  case 'torobillett'
    FLUX = ToroBillett(DT,ULL,UL,UR,URR) ;
  case 'fromm'
    FLUX = Fromm(DT,ULL,UL,UR,URR) ;
  case 'warmingbeam'
    FLUX = WarmingBeam(DT,ULL,UL,UR,URR) ;
  case 'waf'
    FLUX = Waf(DT,U) ;
  otherwise
    FLUX = StandardFlux( 'nflux', nflux, @Flux, DT, DX, U, limiter ) ;
  end
end
                                       
%  _ ___  _   __ ___    _   _   _   _  _  
% |_  |  |_) (_   |    / \ |_) | \ |_ |_) 
% |  _|_ | \ __)  |    \_/ | \ |_/ |_ | \                                       
%

%   #####                                            
%  #     #  ####  #####  #    # #    #  ####  #    # 
%  #       #    # #    # #    # ##   # #    # #    # 
%  #  #### #    # #    # #    # # #  # #    # #    # 
%  #     # #    # #    # #    # #  # # #    # #    # 
%  #     # #    # #    # #    # #   ## #    #  #  #  
%   #####   ####  #####   ####  #    #  ####    ##   %
%
%  Compute the solution of the Riemann problem
%
function F = Godunov( UL, UR )
  global LinearAdvectionModel_paramaters ;
  speed = LinearAdvectionModel_paramaters . speed ;
  % the flux at interface i is computed
  if speed > 0
    F = speed .* UL ; % using the left state
  else
    F = speed .* UR ; % using the right state
  end
end
%  #######                      ######                               
%     #     ####  #####   ####  #     # # #      #      ###### #####  ##### 
%     #    #    # #    # #    # #     # # #      #      #        #      #
%     #    #    # #    # #    # ######  # #      #      #####    #      #
%     #    #    # #####  #    # #     # # #      #      #        #      #
%     #    #    # #   #  #    # #     # # #      #      #        #      #
%     #     ####  #    #  ####  ######  # ###### ###### ######   #      #
%
% Compute intercell fluxes according to the Toro-Billett first order CFL-2 scheme
% See Eq. 13.22, Chapter 13, Ref. 1
function F = ToroBillett(DT,ULL,UL,UR,URR)
  global LinearAdvectionModel_paramaters ;
  speed = LinearAdvectionModel_paramaters . speed ;
  DX    = LinearAdvectionModel_paramaters . DX ;
  CN    = abs(speed*DT/DX) ;
  
  if CN <= 1
    % Conventional Godunov upwind flux
    F = Godunov(UL,UR) ;
  else
    a = (CN-1)/CN ;
    b = 1-a ;
    if speed > 0
      % Information travels from the left
      F = a .* Flux(ULL) + b .* Flux(UL) ;
    else
      % Information travels from the right
      F = a .* Flux(URR) + b .* Flux(UR) ;
    end
  end
end

%  __  _  _  _        _     _   _   _   _  _  
% (_  |_ /  / \ |\ | | \   / \ |_) | \ |_ |_) 
% __) |_ \_ \_/ | \| |_/   \_/ | \ |_/ |_ | \ 

%  #######                             
%  #       #####   ####  #    # #    # 
%  #       #    # #    # ##  ## ##  ## 
%  #####   #    # #    # # ## # # ## # 
%  #       #####  #    # #    # #    # 
%  #       #   #  #    # #    # #    # 
%  #       #    #  ####  #    # #    # 
%
% Compute intercell fluxes according to the Fromm first order upwind method
% See Eq. 13.35, Chapter 13, Ref. 1
%
function F = Fromm(DT,ULL,UL,UR,URR)
  global LinearAdvectionModel_paramaters ;
  speed = LinearAdvectionModel_paramaters . speed ;
  DX    = LinearAdvectionModel_paramaters . DX ;
  CN    = speed*DT/DX ;
  if CN > 0 
    F = speed .* (UL + 0.25 .* (1-CN) .* (UR - ULL)) ;
  else
    F = speed .* (UR - 0.25 .* (1+CN) .* (URR - UR)) ;
  end ;
end

%  #     #                                      ######                       
%  #  #  #   ##   #####  #    # # #    #  ####  #     # ######   ##   #    # 
%  #  #  #  #  #  #    # ##  ## # ##   # #    # #     # #       #  #  ##  ## 
%  #  #  # #    # #    # # ## # # # #  # #      ######  #####  #    # # ## # 
%  #  #  # ###### #####  #    # # #  # # #  ### #     # #      ###### #    # 
%  #  #  # #    # #   #  #    # # #   ## #    # #     # #      #    # #    # 
%   ## ##  #    # #    # #    # # #    #  ####  ######  ###### #    # #    # 
%                                                                        
% Compute intercell fluxes according to the Warming-Beam method
% See Eq. 13.21, Chapter 13, Ref. 1
%
function FLUX = WarmingBeam(DT,ULL,UL,UR,URR)
  global LinearAdvectionModel_paramaters ;
  speed = LinearAdvectionModel_paramaters . speed ;
  DX    = LinearAdvectionModel_paramaters . DX ;
  CN    = speed*DT/DX ;
  if speed > 0 
    FLUX = (0.5*speed) .* ( (3-CN) .* UL - (1-CN) .* ULL ) ;
  else
    FLUX = (0.5*speed) .* ( (3+CN) .* UR - (1+CN) .* URR ) ;
  end
end

%  #     #    #    ####### 
%  #  #  #   # #   #       
%  #  #  #  #   #  #       
%  #  #  # #     # #####   
%  #  #  # ####### #       
%  #  #  # #     # #       
%   ## ##  #     # #       
%
%  Compute the WAF flux with the exact Riemann solver,
%  to be used in explicit conservative formula, subroutine UPDATE

function FLUX = Waf( DT, U )
  global LinearAdvectionModel_paramaters ;
  speed   = LinearAdvectionModel_paramaters . speed ;
  DX      = LinearAdvectionModel_paramaters . DX ;
  limiter = LinearAdvectionModel_paramaters . limiter ;
  CN      = speed*DT/DX ;

  % Compute local jump DLOC. This is the change across the wave of speed c.
  DLOC = U(:,3:end-1) - U(:,2:end-2) ;

  % Identify the upwind direction
  if speed > 0 % Information travels from the left
    IUPW = -1 ;
  else         % Information travels from the right
    IUPW = +1 ;
  end

  % Compute the upwind jump DUPW. This depends on the direction of the wave
  % in the solution of the Riemann problem RP(UL, UR), namely the sign of c
  DUPW = U(:,3+IUPW:end-1+IUPW) - U(:,2+IUPW:end-2+IUPW) ;
  ONES = ones (size(DLOC)) ;

  % Apply TVD condition to find WAF limiter function WAFLIM
  WAFLIM = WafLimiter( DLOC, DUPW, CN, limiter ) ;
  WL     = 0.5 .* (ONES + WAFLIM) ;
  WR     = 0.5 .* (ONES - WAFLIM) ;

  % Compute fluxes on data
  FD = Flux( U(:,2:end-1) ) ;

  % Compute WAF Flux FLUX(I) corresponding to RP(I,I+1)
  FLUX = WL .* FD(1:end-1) + WR .* FD(2:end) ; 

end

function WAFLIM = WafLimiter( DLOC, DUPW, CN, limiter )

  TOLLIM = 1E-6 ;

  % Select limiter function WAFLIM
  ZEROS = zeros(size(DLOC)) ;
  ONES  = ones (size(DLOC)) ;

  % Compute RATIO=DUPW/DLOC of upwind to local changes 
  % some workaround to avoid division by zero are added
  RATIO = sign(DLOC).*DUPW./max(abs(DLOC),TOLLIM.*ONES) ;

  switch ( limiter )
    case 'godunov'
      WAFLIM = sign(CN) .* ONES ;
      return ;
    case 'secondorder'
      WAFLIM = CN .* ONES ;
      return ;
    case 'superbee'
      B = max( min(2.*RATIO, ONES), min(RATIO, 2.*ONES) ) ;
    case 'vanleer'
      B = max(RATIO,ZEROS) ;
      B = 2.*B./(ONES+B) ;
    case 'vanalbada'
      B = (RATIO + RATIO.^2 ) ./ (ONES + RATIO.^2) ;
    case 'minmod'
      B = min(RATIO, ONES) ;
  end

  % Transform to WAF limiter
  WAFLIM = sign(CN) .* ( ONES - (ONES - abs(CN)) .* max(ZEROS, B) ) ;
end

