%
% solve the linearised Hyperbolic system in one space dimension 
%
%     /    \       /          \   /   \
%     | q1 |       |  0    -a |   | 0 |
% D_t |    | + D_x |          | = |   |
%     | q2 |       | -a     0 |   | 0 |
%     \    /       \          /   \   /
%
% for t > 0 and a <= x <= b
%
%  P = parameter
%
%  what
%  'setup'  
%  'flux'    compute flux
%
function [varargout] = Linear2x2Model( what, varargin )

  global Linear2x2_paramaters ;

  switch lower( what )

    case 'setup'
      expected( what, 1, nargin, 0, nargout ) ; % a
      Linear2x2_paramaters = clear_struct(Linear2x2_paramaters) ;
      Linear2x2_paramaters . a = varargin{1} ;

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

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

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

    case 'riemann'
      expected( what, 4, nargin, 1, nargout ) ; % UL, UR, T, X
      UL = varargin{1} ;
      UR = varargin{2} ;
      T  = varargin{3} ;
      X  = varargin{4} ;
      varargout{1} = RiemannExact( UL, UR, T, X ) ;

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

    case 'extrapolate'
      expected( what, 3, nargin, 1, nargout ) ; % U, left_bc, right_bc
      U        = varargin{1} ;
      left_bc  = varargin{2} ;
      right_bc = varargin{3} ;
      varargout{1} = Extrapolate( U, left_bc, right_bc ) ;

    otherwise
      error(['Linear2x2Model 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 Linear2x2 call to ', what, ' expected ', int2str( ein ), ' input argument, found ', int2str(nin-1) ]) ;
  end
  if eout ~= nout 
    error( ['in Linear2x2 call to ', what, ' expected ', int2str( eout ), ' output argument, found ', int2str(nout) ]) ;  
  end
end  

%  #######                      
%  #       #      #    # #    # 
%  #       #      #    #  #  #  
%  #####   #      #    #   ##   
%  #       #      #    #   ##   
%  #       #      #    #  #  #  
%  #       ######  ####  #    # 
% compute the flux
function F = Flux( US )
  global Linear2x2_paramaters ;
  a = Linear2x2_paramaters . a ;
  F = [ 0, -a ; -a, 0 ] * US ;
end

%  ######                                       
%  #     # # ###### #    #   ##   #    # #    # 
%  #     # # #      ##  ##  #  #  ##   # ##   # 
%  ######  # #####  # ## # #    # # #  # # #  # 
%  #   #   # #      #    # ###### #  # # #  # # 
%  #    #  # #      #    # #    # #   ## #   ## 
%  #     # # ###### #    # #    # #    # #    # 
%
function QS = RiemannStar( QL, QR )
  global Linear2x2_paramaters ;
  a  = Linear2x2_paramaters . a ;
  A  = (QR + QL) ./ 2 ;
  B  = (QR - QL) ./ 2 ;
  QS(1,:) = A(1,:) + B(2,:) ;
  QS(2,:) = A(2,:) + B(1,:) ;
end

function US = RiemannExact( UL, UR, T, X )
  global Linear2x2_paramaters ;
  a = Linear2x2_paramaters . a ;

  l1 = -abs(a) ;
  l2 = +abs(a) ;
  
  IDX = find( X - l1 * T >= 0 & X - l2 * T <= 0 ) ;
  US(:,IDX) = RiemannStar( UL, UR ) * ones( size(IDX) ) ;

  IDX = find( X - l1 * T < 0 ) ;
  US(:,IDX) = UL * ones( size(IDX) ) ;

  IDX = find( X - l2 * T > 0 ) ;
  US(:,IDX) = UR * ones( size(IDX) ) ;
end

%                        #####                              
%  #    #   ##   #    # #     # #####  ###### ###### #####  
%  ##  ##  #  #   #  #  #       #    # #      #      #    # 
%  # ## # #    #   ##    #####  #    # #####  #####  #    # 
%  #    # ######   ##         # #####  #      #      #    # 
%  #    # #    #  #  #  #     # #      #      #      #    # 
%  #    # #    # #    #  #####  #      ###### ###### ##### 
% 
function MS = maxSpeed( U )
  global Linear2x2_paramaters ;
  a  = Linear2x2_paramaters . a ;
  MS = abs(a) ;
end

%  #######                                                                     
%  #       #    # ##### #####    ##   #####   ####  #        ##   ##### ###### 
%  #        #  #    #   #    #  #  #  #    # #    # #       #  #    #   #      
%  #####     ##     #   #    # #    # #    # #    # #      #    #   #   #####  
%  #         ##     #   #####  ###### #####  #    # #      ######   #   #      
%  #        #  #    #   #   #  #    # #      #    # #      #    #   #   #      
%  ####### #    #   #   #    # #    # #       ####  ###### #    #   #   ###### 
%
function Uextra = Extrapolate( U, L, R )

  global Linear2x2_paramaters ;
  a = Linear2x2_paramaters . a ;

  % Left boundary
  switch lower(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
      error( ['in Linear2x2 call to Extrapolate L=', L, ' unknown boundary condition']) ;  
  end
  % Right boundary
  switch lower(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
      error( ['in Linear2x2 call to Extrapolate R=', R, ' unknown boundary condition']) ;  
  end

  % build vector of extrapolated values
  Uextra = [ ULL, UL, U, UR, URR ] ;
end
