function f = HermiteILT(F,t,n)

% HermiteILT: function for inverting the Laplace transform F(z),
% which should be supplied as an anonymous function.  The target function 
% f(t) is computed at the single input value t.  n-point Gauss-Hermite
% quadrature is used, with n less than or equal to 20. A function for 
% computing the Hermite weights and nodes (as column and row vectors, 
% respectively) should be available.   Recommendation:  hermpts as
% available in the Chebfun system.  Contour parameters mu and L have
% been pre-computed as described in the paper "Gauss-Hermite Quadrature
% for the Bromwich integral" by JAC Weideman.   
% Example of usage: with values of t and n defined, 
% >>F = @(z) exp(-1./z)./sqrt(z);
% >>f = HermiteILT(F,t,n)
% With n = 20, values close to machine precision should be obtained
% for transforms with singularities restricted to the negative real
% axis.                                   JACW, March 2019
 
 
 
% Use pre-computed parameters (for the case of singularities on or
% near the negative real axis, n <= 20), or define your own. 
[mu,L] = params(n);

z =  @(phi) (mu/t)*(1+1i*phi).^2;   % Define contour
dz = @(phi) 2i*(mu/t)*(1+1i*phi);   % and its derivative

[r,w] = hermpts(n);            % Compute Hermite points

m = floor(n/2);                % Adjust for n even/odd       
r = r(m+1:n); w = w(m+1:n);
w(1) = w(1)/2^rem(n,2);  
 
                                % Compute integrand
g = L/(2i*pi)*exp(r.^2+z(L*r)*t).*F(z(L*r)).*dz(L*r);
f = 2*real(w*g);                % and integrate

end

function [mu,L] = params(n)

if n > 20; disp('Pre-computed contour parameters not available'); return; end

mul = [
 2 0.906983 0.922150 ;
 3 1.182682 0.819084 ;
 4 1.454542 0.744987 ;
 5 1.723091 0.688037 ;
 6 1.989685 0.642548 ;
 7 2.255579 0.605112 ;
 8 2.521662 0.573598 ;
 9 2.786037 0.546588 ;
10 3.050443 0.523100 ;
11 3.314528 0.502426 ;
12 3.577201 0.484044 ;
13 3.841881 0.467559 ;
14 4.104752 0.452663 ;
15 4.367202 0.439116 ;
16 4.629882 0.426727 ;
17 4.891091 0.415335 ;
18 5.153948 0.404817 ;
19 5.416676 0.395064 ;
20 5.680141 0.385991 ;
];

mu = mul(n-1,2); L = mul(n-1,3);

end

