%5.3 ilaplacecode.m
%This code was contributed by Andre Weideman
function f = talbotmod(F,t,N)
% The function f = talbot(F,t,N) implements Talbot's
% *modified* method for the inversion of the Laplace
% Transform. F is the transform to be inverted, supplied
% as an inline function of z. (The code assumes that F is
% real when z is real (*).) N is the number of F evaluations.
% The method is based on a 2N panel midpoint rule approximation,
% but because of the symmetry (*) it can be reduced to N evaluations
% of F. The parameters sigma, mu, and nu are tuned for optimality
% in the case where all the singularities of F are on or near the
% negative imaginary axis. In this case a convergence rate of
% about exp(-2.7N) can be expected, i.e., N = 14 should usually be
% sufficient for full precision in MATLAB. The code does not yet
% make provision for multiple values of t.
%
% Example of usage:
% >> F = inline('1./(z.*sqrt(z+1))')
% >> f = talbotmod(F,1,10)
% JAC Weideman, September 2006
sigma = -1.2244*N/t; mu = 1.0035*N/t; nu = 0.5272; % define optimal params
alpha = 0.6407; % new parameter
theta = (2*[0:N-1]+1)*pi/(2*N); % equispaced theta vals
z = sigma + mu*(theta.*cot(alpha*theta) + nu*i*theta); % modified contour
dz = mu*(cot(alpha*theta)-alpha*theta.*csc(alpha*theta).^2 + nu*i); % dz/dtheta
f = 1/N*imag(sum(exp(z*t).*F(z).*dz)); % midpoint sum
y