function [integral,est] = ModifiedTalbot( F,t,opt)

% This function computes an inverse Laplace transform using the cotangent 
% version of the Talbot contour. For small N precomputed parameter values 
% are used and for larger N parameters optimized for stability issues 
% are used to counteract the roundoff error.
% 
% BD & JACW, 2014/03/24
%  
% INPUT
% 
%   F           Function handle of the transform which shall be inverted 
%               numerically      
%   t           Single time value for which the inverse shall be evaluated
% 
% OPTIONAL 
% 
%   opt         Containing a struct with arbitrary many of the following
% 	
%       shift_transform Frequency shift of the transform
%                            0 (default)
%       re_control      Determines if roundoff error parameters should be 
%                       taken into account.
%                            re_control = 1 if so (default)
%                            re_control = 0 if not
%                       If singularities are located off the negative real 
%                       axis, this parameter should always be 0.
%       accuracy_goal   Stops the evaluation if the error estimate reaches 
%                       the given goal.
%                           1e-10 (default)				
%       max_nodes       Maximal number of nodes used
%                           100 (default)
%       fig             Sets if a figure of the convergence history is 
%                       shown.
%                           0 no plot (default)
%                           1 plot 
%       exact_result    Exact result of the inverse of F at time t for 
%                       testing purposes 
% 
% OUTPUT
% 
%   integral    Value of the approximation	
%   est         Estimate for the relative error of 'integral' to the
%               exact solution or (if 'exact_result' is assigned) 
%               relative error of 'integral' to 'exact_result'
%  
% EXAMPLES
% 
%   ModifiedTalbot(@(z) 1./(z+1),1) 
% 	ModifiedTalbot(@(z) 1./sqrt(z),2,struct('accuracy_goal',1e-6,'fig',1))
%
%   proper treating of branchcuts for 1./sqrt(z.^2+1)
%   ModifiedTalbot(@(z) 1./(sqrt(z-1i).*sqrt(z+1i)),1e-2)
%

% In case of no additional inputs
if nargin<3, opt=struct();end

% Getting values specified in options, if not specified, set default
[shift, re, acc, NMax, fig, fe]=processOptions(opt,...
    {'shift_transform','re_control','accuracy_goal','max_nodes','fig','exact_result'},...
    {  0,1,1e-10,100,0,[] });

% Precomputed optimal parameters for Talbot's contour
a= .6407; m = .5017; s = .6122; n = .2645;

% Fix the correction term if a transform shift is involved
shiftI = exp(shift*t);

% Compute the inverse using Talbot's method and the precomputed parameters
% until the accuracy goal is reached or an increase of the error due to 
% roundoff errors is observed
NN1 = 4:2:NMax;
NN2Start = NMax+2;
E1 = [];

% First value for comparing the approximation
integrals=1;
for N = NN1
    theta = ([0:N/2-1]+0.5)*2*pi/N;
    z  = N/t*(-s + m*theta.*cot(a.*theta)+n*1i*theta);
    Dz = N/t*(n*1i+m*cot(a*theta)-m*a*theta.*(csc(a*theta)).^2);
    fa = shiftI*2/N*imag(sum(exp(t*z).*Dz.*F(z+shift)));
    
    % Check if fe is assigned
    % TRUE      compare with exact result
    % FALSE     compare with previous computed value of the integral
    if isempty(fe)
        E1 = [E1; abs((fa-integrals(end))/fa);];
        integrals = [integrals;fa];
    else
        E1 = [E1; abs((fa-fe)/fe);];
        integrals = [integrals;fa];
    end  
    
    % If given accuracy is reached stop the evaluation of the loop
    % The second RE-parameter loop will not be evaluated due to the
    % previous setting of NN2Start
    if E1(end)<acc
        NN1 = 4:2:N;
        break;
    end
    
    % If error starts to increase two times in a row, at least 24 nodes are
    % used and the errors are of a reasonable size at which we can speak of
    % convergence (i.e. < 1e-2 ), switch on RE-control
    if (re && N >=24)
        if (E1(end) > E1(end-1) && E1(end)>E1(end-2) && E1(end-2)<1e-2)
            NN1 = 4:2:N;
            NN2Start = N+2;
            break;
        end
    end
end

% If we consider roundoff error parameters, the following section will
% begin right after an increase of the error due to roundoff issues is
% recognized. 
if re
    E2 = []; 
    NN2 = NN2Start:2:NMax;

    % Calculate an approximation to the constant in the O expressions when
    % comparing the estimated rate of convergence with the estimated
    % roundoff error impact
    K = eps*exp(1.529*(NN2Start-4));
    for N = NN2
        theta = ([0:N/2-1]+0.5)*2*pi/N;
    
        % Calculate the optimal value for c for roundoff parameters
        fz = @(c) c + (log(eps/K))/N - 2*a*c^2*(sin (a*pi))^2/(2*a*c^2*...
            (sin (a*pi))^2 - pi*sin (2*a*pi)*sinh (a*c)^2) + (2*c*...
            (sin (a*pi))^2*(sinh (a*c))^2)/(a*(2*a*c^2*(sin (a*pi))^2 ...
            -pi*sin (2*a*pi)*(sinh (a*c))^2));
    	c = fzero(fz,1);    
       
        % Determine the roundoff parameters and calculate the approximation
        B = c*(sin(a*pi))^2/(2*a*c^2*(sin(a*pi))^2-pi*sin(2*pi*a)*(sinh(a*c))^2);
        s = 2*a*c^2*B;
        m = 2*(sinh(a*c))^2*B;
        n = (sinh(2*a*c)-2*a*c)*B;
        z  = N/t*(-s + m*theta.*cot(a.*theta)+n*1i*theta);
    	Dz = N/t*(n*1i+m*cot(a*theta)-m*a*theta.*(csc(a*theta)).^2);
    	fa = shiftI*2/N*imag(sum(exp(t*z).*Dz.*F(z+shift)));
    
        % Check if fe is assigned
    	% TRUE      compare with exact result
        % FALSE     compare with previous computed value of the integral
      	if isempty(fe)
            E2 = [E2; abs((fa-integrals(end))/fa);];
            integrals = [integrals;fa];
        else
            E2 = [E2; abs((fa-fe)/fe);];
            integrals = [integrals;fa];
        end
    
    	% If given accuracy is reached stop the evaluation of the loop
      	if E2(end)<acc
            NN2 = (NN2Start):2:N;
            break;
        end
    end  
end

% Plot results
if fig
    if re
        plotResults(NN1,E1,NMax,t,NN2,E2)
    else
        plotResults(NN1,E1,NMax,t)
    end
end

% Return solution
if re
    E = [E1;E2];
else
    E = E1;
end
[est, pos] = min(E);
integral = integrals(pos+1);
end







function [varargout]=processOptions (opt,names,stdvalues)
% Reads options or default values out of struct.
% See http://www-m3.ma.tum.de/Software/ComplexVisualToolbox for further
% information

varargout=cell(1,length(names));

for i=1:length(names)
   if (isfield(opt,names{i}))
      varargout{i}=opt.(names{i}); 
   else
      varargout{i}=stdvalues{i};
   end
end
end



function plotResults(NN1,E1,NMax,t,NN2,E2)
    figure;    
    semilogy(NN1,E1,'k-'); hold on;
    semilogy(NN1,E1,'k.','MarkerSize',16);

    % Theoretically predicted rate of convergence
    semilogy(NN1,exp(-1.358*NN1),'b-.','LineWidth',2);
    if nargin == 6
    	semilogy(NN2,E2,'k-'); shg; hold on;
    	semilogy(NN2,E2,'k.','MarkerSize',14);   
    end
    axis([0 NMax+2 1e-16 1e0]); grid on;
    title(['t = '  num2str(t)],'FontSize',14,'interp','LaTeX')
end
