%
%  PURPOSE:  Call a RBC table dowloaded from 
%    http://www.dam.brown.edu/people/sfield/KernelsRWZ/
%    http://www.math.unm.edu/~lau/KernelsRWZ/
%
%   Supported for Rho = 15, 30, 60, 120, 240 
%                 ell = 2 to 64
%
%  Compressed kernel is a sum of poles in Laplace frequency or
%  sum of exponentials in time domain.
%
%  Kernel = \sum_{k=1}^d \gamma_k exp(t \beta_k)
%  this routine calls tables which return \gamma_k and \beta_k
%
%  INPUT:    Rho         (areal value)/(2M) of outer boundary
%  
%            ell         spherical harmonic mode ell
%
%            M           black hole mass
%
%            type        potential: 'rw' for RW/CPM and 'z' for Zerilli
%
%            RBC_Dir     absolute or relative path for RBC kernel tables
%
%
%  OUTPUT:   kap = poles \beta_k which are negative and real
%
%            mu  = streangths \gamma_k corresponding to kap
%
%            qR and qI = real/imaginary parts of conjugate pair poles
%
%            mR and mI = streangth of these conjugate pair poles
%
%
%  AUTHOR:  Scott Field
%           Department of Physics
%           University of Maryland, Collge Park
%           August 31, 2013
% 
% function [kap,mu,qR,qI,mR,mI,nmu,nmR] = GetBoundaryTables_RWZ(Rho,ell,M,type,RBC_Dir)
  function [kap,mu,qR,qI,mR,mI,nmu,nmR] = GetBoundaryTables_RWZ(Rho,ell,M,type,RBC_Dir)

%%% complete path to folder containing downloaded table %%%
switch lower(type)
    
    case{'rw','cpm'}
        PathToRBCFolder = strcat(RBC_Dir,'/TablesHeunRBC/');
    case{'z'}
        PathToRBCFolder = strcat(RBC_Dir,'/TablesZerilliRBC/');
        
end

%%% get corresponding outer boundary in rho = r/(2M) from xb  %%%
switch(Rho)
    case(15)
       RBCFolder = 'Rho015/';
    case(30)
       RBCFolder = 'Rho030/';
    case(60)
       RBCFolder = 'Rho060/';
    case(120)
       RBCFolder = 'Rho120/';
    case(240)
       RBCFolder = 'Rho240/';
    otherwise
       error('Boundary kernel not found')
end

PathToRBCFolder = strcat(PathToRBCFolder,RBCFolder)
addpath(PathToRBCFolder)


%%% load the table %%%
if (ell < 2 || ell > 64 )
    
    error('Boundary kernel not found')
    
else
    
    EllString = num2str(ell);
    Ell = '000';
    if( numel(EllString) == 1 )
        Ell(3) = EllString(1);
    elseif( numel(EllString) == 2 )
        Ell(2) = EllString(1);
        Ell(3) = EllString(2);
    end
    
    UnixCommand = ['ls ' PathToRBCFolder '*' Ell 'D*.m  | xargs -n 1 basename']
    [status,table] = unix(UnixCommand);
    table = table(1:end-3)
    pause(4)
    [kap, mu, qR, qI, mR, mI] = eval(table);
    mu  = mu/(2*M);
    kap = kap/(2*M);
    mR  = mR/(2*M);
    mI  = mI/(2*M);
    qR  = qR/(2*M);
    qI  = qI/(2*M);
    nmu = length(mu); % number of pure real poles
    nmR = length(mR); % number of conjugate poles
    
end