[Invertf,sigmaP,bP,RelTotalError,AbsTotalError,AbsTruncateError,AbsRoundoffError,ToleranceMetFlag,LaguerreCoef] =WeeksMethod(FLaplace,Timevec,RelErrorTol,SearchSwitch)

function [Invertf,sigmaP,bP,RelTotalError,AbsTotalError,AbsTruncateError,AbsRoundoffError,ToleranceMetFlag,LaguerreCoef] =WeeksMethod(FLaplace,Timevec,RelErrorTol,SearchSwitch)
%WEEKSMETHOD A Weeks' method numerical inverse Laplace transform of F(s) to f(t) 
%
%  Use:
%  [Invertf,sigmaP,bP,RelTotalError,AbsTotalError,AbsTruncateError,AbsRoundoffError,ToleranceMetFlag,LaguerreCoef] =...
%  WeeksMethod(FLaplace,Timevec,RelErrorTol,SearchSwitch)
%  
%  Input:
%  FLaplace = a character string expression for the Laplace transform space
%  function F(s) in terms of 's'
%  Timevec = times where f(t) are to be computed (vector)
%  RelErrorTol (optional) = an input total relative error tolerance for each time [default 5 percent]
%  SearchSwitch (optional) = the search method [default 0]
%   0-CPU local fminbnd, 1-CPU global, 2-GPU/JACKET global
%
%  Output: 
%  (all outputs have a length equal to that of the time input Timevec)
%  Invertf = numerically computed f(t)
%  sigmaP = automatically determined sigma or the user defined value if an input 
%  bP = automatically determined beta or the user defined value if an input
%  RelTotalError = the estimated total relative error
%  AbsTotalError = the estimated total absolute error
%  AbsTruncateError = the estimated absolute trunction error
%  AbsRoundoffError = the estimated absolute round-off error
%  ToleranceMetFlag = this flag indicates if the relative error tolerance was met: 0-fail, 1-success
%  LaguerreCoef = the Laguerre polynomial expansion coefficients
%
%  Comment:
%  The method uses the following Laguerre polynomial expansion expresion for F(s):
%  s = \sigma - b\frac{w+1}{w-1}
%  \sum_{n=0}^{\infty} a_{n}w^{n} = \frac{2b}{1-w}F(s)
%
%  The time domain function f(t) is then:
%  f(t) = \exp(\sigma t)\sum_{n=0}^{\infty}a_{n} exp{(-bt)} L_{n}(2bt)
%
%  where the Laguerre polynomials are defined by:
%  L_{n}(x) = \frac{e^{x}}{n!} \frac{d^{n}}{dx^{n}}(x^{n} e^{-x})
%
%  References: 
%  Weeks' method:
%  1) Algorithms for Parameter Selection in the Weeks Method for Inverting the Laplace Transform, 
%     J. A. C. Weideman, SIAM Journal on Scientific Computing, vol. 21, pp. 111-128, 1999. 
%     http://dip.sun.ac.za/~weideman/research/weeks.html
%  2) Application of Weeks method for the numerical inversion of
%     the Laplace transform to the matrix exponential, P. Kano, M. Brio, and J. Moloney,
%     Communications in mathematical sciences", vol. 3, no. 3, pp. 335-372, 2005.
%     http://math.arizona.edu/~brio/WEEKS_METHOD_PAGE/pkanoWeeksMethod.html
%  3) Software for an Implementation of Weeks' Method for the Inverse Laplace Transform Problem
%     S. Garbow, G. Giunta, and J. N. Lyness, A. Murli, ACM, 1988.
%     http://www.utd.edu/~cantrell/ee6481/lectures/Laplace/Garbow-88.pdf
% 
%  Author: 
%  Patrick Kano, Moysey Brio - 2011
%
%  Modification Date [M/D/Y]:
%  03/31/2011 - Version 1.0

%%%Initialize output vectors
Ntimes = length(Timevec);
Invertf = zeros(1,Ntimes);
sigmaP = zeros(1,Ntimes);
bP = zeros(1,Ntimes);
RelTotalError = zeros(1,Ntimes);
AbsTotalError = zeros(1,Ntimes);
AbsTruncateError = zeros(1,Ntimes);
AbsRoundoffError = zeros(1,Ntimes);
ToleranceMetFlag = zeros(1,Ntimes); %0-failed to meet tolerance, 1-succeeded in meet tolerance
LaguerreCoef = cell(1,Ntimes); %this can be a large amount of space 

%Define default (sigma, b) space search parameters
sigmin = 0;
sigmax = 20;
tols = [1.0 0.5 0.25]; %increase grid resolution

bmax = 40;
tolb = [1.0 0.5 0.25]; %increase grid resolution

%Define the default number of Laguerre expansion coefficients
NLag = [32 64 128 256 512];
  
switch nargin
 case 2
  SearchSwitch = 0; %default local fminbnd search
  RelErrorTol = 0.05*ones(1,Ntimes);
 case 3
  SearchSwitch = 0; %default local fminbnd search
 case 4
  %do nothing, the user has specified all input parameters to the wrapper
 otherwise
  error('WeeksMethod:nargincheck','Incorrect number of inputs.  At least 2 are required.');
end %switch 

%Check that a Jacket is installed and available.
if(SearchSwitch==2) %GPU global 
 JacketFlag = wfnCheckGPUsoft();
 if(JacketFlag==0)
  error('WeeksMethod:JACKET','JACKET is unavailable.  The user may not select a GPU computation.');
 end
end

%%%Main Search
for tidx=1:Ntimes
 for ntolidx=1:length(tols)
  for nLagidx=1:length(NLag)
   %tempindices = [tidx,ntolidx,nLagidx]
        
   [testInvertf,testsigmaP,testbP,testRelTotal,testAbsTotal,testAbsTruncate,testAbsRoundoff,testLaguerreCoef]=...  
   wfnWeeksCore(FLaplace,Timevec(tidx),NLag(nLagidx),sigmin,sigmax,bmax,SearchSwitch,tols(ntolidx),tolb(ntolidx));
          
   if testRelTotal<RelErrorTol(tidx)
    Invertf(tidx) = testInvertf;
    sigmaP(tidx) = testsigmaP;
    bP(tidx) = testbP;
    RelTotalError(tidx) = testRelTotal;
    AbsTotalError(tidx) = testAbsTotal;
    AbsTruncateError(tidx) = testAbsTruncate;
    AbsRoundoffError(tidx) = testAbsRoundoff;
    ToleranceMetFlag(tidx) = 1;
    LaguerreCoef{tidx} = testLaguerreCoef;
    break;
    
   elseif(RelErrorTol(tidx)<testRelTotal) %update
    Invertf(tidx) = testInvertf;
    sigmaP(tidx) = testsigmaP;
    bP(tidx) = testbP;
    RelTotalError(tidx) = testRelTotal;
    AbsTotalError(tidx) = testAbsTotal;
    AbsTruncateError(tidx) = testAbsTruncate;
    AbsRoundoffError(tidx) = testAbsRoundoff;
    ToleranceMetFlag(tidx) = 0;
    LaguerreCoef{tidx} = testLaguerreCoef;
   end    
  end %nLagidx
  
  if(ToleranceMetFlag(tidx)==1), break; end %out to next time index
 end %ntolidx
end %tidx    

end %function definition
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%