%  Model for the Burgers equations
%
%   u(t,x) + u(t,x) * u (t,x) = 0
%   t                  x
%
%   for t > 0 and -1 <= x <= 1
%
function [varargout] = BurgerModel( what, varargin )

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

  global BurgerModel_paramaters ;

  switch lower(what)

    case 'setup'
      expected( what, 0, nargin, 0, nargout ) ;
       BurgerModel_paramaters = clear_struct(BurgerModel_paramaters) ;

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

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

    case 'set_numerical_flux'
      expected( what, 1, nargin, 0, nargout ) ;
      BurgerModel_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} = BurgerModel_paramaters . DX ;    

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

    case 'get_final_time'
      expected( what, 0, nargin, 1, nargout ) ;
      varargout{1} = BurgerModel_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 '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(['BurgerModel 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 BurgerModel call to ', what, ' expected ', int2str( ein ), ' input argument, found ', int2str(nin-1) ]) ;
  end
  if eout ~= nout 
    error( ['in BurgerModel call to ', what, ' expected ', int2str( eout ), ' output argument, found ', int2str(nout) ]) ;  
  end
end

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

% compute the flux
function F = Flux( U )
  F = 0.5 * U .^ 2 ;
end

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

function U = Init(N)
  global BurgerModel_paramaters ;
  % Calculate mesh size DX
  DX = 2.0/N ;
  X  = [-1+DX/2:DX:1-DX/2] ;
  
  if isfield(BurgerModel_paramaters, 'test')
    switch BurgerModel_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(['BurgerModel::Init unknown test ', BurgerModel_paramaters . test]) ;
    end
  else
    error('BurgerModel::Init unknown test you must set it by set_test!') ;      
  end
  BurgerModel_paramaters . DX = DX ;
  BurgerModel_paramaters . X  = X ;
end


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

  global BurgerModel_paramaters ;

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

  BurgerModel_paramaters . test     = test ;
  BurgerModel_paramaters . bc_left  = bc_left ;
  BurgerModel_paramaters . bc_right = bc_right ;
  BurgerModel_paramaters . tf       = tf ;

end

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

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

% Extrapolate values for boundary condition
function [ULL, UL, UR, URR] = Extrapolate_for_BC( U )
  global BurgerModel_paramaters ;
  L = BurgerModel_paramaters . bc_left ;
  R = BurgerModel_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 BurgerModel_paramaters ;
  nflux = BurgerModel_paramaters . nflux ;
  %
  switch nflux
  case 'godunov'
    CFL = 1.0 ;
  otherwise
    CFL = StandardFlux( 'get_cfl', nflux ) ;
  end
end


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

function FLUX = NFlux( DT, U )

  global BurgerModel_paramaters ;

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

  nflux = BurgerModel_paramaters . nflux ;
  DX    = BurgerModel_paramaters . DX ;

  switch nflux
  % monolitic numerical fluxes
  case 'godunov'
    FLUX = Godunov(UL,UR) ;
  otherwise
    FLUX = StandardFlux( 'nflux', nflux, @Flux, DT, DX, U ) ;
  end
end
                                       
%  _ ___  _   __ ___    _   _   _   _  _  
% |_  |  |_) (_   |    / \ |_) | \ |_ |_) 
% |  _|_ | \ __)  |    \_/ | \ |_/ |_ | \                                       
%

%   #####                                            
%  #     #  ####  #####  #    # #    #  ####  #    # 
%  #       #    # #    # #    # ##   # #    # #    # 
%  #  #### #    # #    # #    # # #  # #    # #    # 
%  #     # #    # #    # #    # #  # # #    # #    # 
%  #     # #    # #    # #    # #   ## #    #  #  #  
%   #####   ####  #####   ####  #    #  ####    ##   %
%
%  Compute the solution of the Riemann problem
%
function F = Godunov( UL, UR )

  %%USTAR = zeros(size(UL)) ;

  IDX = find( UL > UR ) ;
    % Solution is a shock wave Compute shock speed S
    S(IDX) = 0.5*(UL(IDX) + UR(IDX)) ;
    % Sample the state along the t-axis
    II = IDX( find( S(IDX) >= 0 )) ;
    USTAR(II) = UL(II) ;
    II = IDX( find( S(IDX) <  0 )) ;
    USTAR(II) = UR(II) ;

  IDX = find( UL <= UR ) ;
    % Solution is a rarefaction wave. There are 3 cases:
    II = IDX(find( UL(IDX) >= 0 )) ;
    % Right supersonic rarefaction
    USTAR(II) = UL(II) ;

    II = IDX(find( UR(IDX) <= 0 )) ;
    % Left supersonic rarefaction
    USTAR(II) = UR(II) ;

    II = IDX(find( UL(IDX) <= 0 & UR(IDX) >= 0 )) ;
    % Transonic rarefaction
    USTAR(II) = 0 ;
  
  F = 0.5 * USTAR .^ 2 ;

end
