function v=Convolution(F,g,dt,N,B,mu0,nu,sigma)
% CONVOLUTION fast convolution algorithm by Lubich and Schaedle
%  v=Convolution(F,g,N,B,mu0,nu,sigma) computes approximations of the
%  convolution integrals int(f(j*dt-tau)*g(t),tau=0..j*dt),
%  j=1,2,...,length(g)-1, where f is defined by its Laplace transform
%  F, g is a vector, and dt is a time step, using the fast algorithm
%  of Lubich and Schaedle, with 2*N+1 quadrature points, base B and
%  paramters mu, nu and sigma. The result is the vector v of the same
%  length as g. 

Nt=length(g)-1;
[la0,w0]=NodesAndWeights(N,mu0/dt,nu,sigma);
Lg=ceil(log(1+Nt*(B-1))/log(B)-10*eps)-1; % number of Talbot contours
tau=zeros(Lg,1);                          % 10*eps to compensate for
for l=1:Lg                                % roundoff errors
  Tl=(2*B^l-1)*dt;
  [la(:,l),w(:,l)]=NodesAndWeights(N,mu0/Tl,nu,sigma);
  Bs(l)=2+sum(B.^(1:l-1));
end;
phi1=sum(w0.*feval(F,la0)./la0.*exp(la0*dt));
phi2=sum(w0.*feval(F,la0)./la0.^2.*exp(la0*dt));

v(1)=phi1*g(1)+phi2*(g(2)-g(1))/dt;
y=zeros(2*N+1,Lg);                        % to store ODE solutions
for n=2:Nt
  t=n*dt;                                 % do the last dt
  for l=1:Lg                              % advance all ODEs
    y(:,l)=y(:,l)+(exp(dt*la(:,l))-1)./(dt*la(:,l))...
           .*(dt*la(:,l).*y(:,l)+dt*g(n-1)+dt*(g(n)-g(n-1))/dt./la(:,l))...
           -dt*(g(n)-g(n-1))/dt./la(:,l);
    if l==1                               % always activate first
      ya(:,l)=y(:,l); tau(l)=tau(l)+dt;
    else
      if mod(n,B^(l-1))==1                % save for future use
        ys(:,l)=y(:,l); y(:,l-1)=zeros(2*N+1,1);
      elseif n>=Bs(l) & mod(n-Bs(l),B^(l-1))==0     % activate
        ya(:,l)=ys(:,l); tau(l)=tau(l)+B^(l-1)*dt;
      end;
    end;
  end;
  v(n)=phi1*g(n)+phi2*(g(n+1)-g(n))/dt;
  L=ceil(log(1+n*(B-1))/log(B)-10*eps)-1; % number of active Talbot contours
  for l=1:L                               % sum remaining contributions
    v(n)=v(n)+sum(w(:,l).*feval(F,la(:,l)).*exp((t-tau(l))*la(:,l)).*ya(:,l));
  end;
end;

function [la,w]=NodesAndWeights(N,mu,nu,sigma);
% NODESANDWEIGHTS computes nodes and weights of Talbot contour
                      
th=(-N:N)*pi/(N+1);              
for j=1:length(th)
  if th(j)~=0
    la(j,1)=sigma+mu*(th(j)*cot(th(j))+i*nu*th(j));
    w(j,1)=1/(2*i*(N+1))*mu*((cos(th(j))*sin(th(j))-th(j))/sin(th(j))^2+i*nu);
  else                                    % need to handle th=0 separately
    la(j,1)=sigma+mu;
    w(j,1)=mu*nu/(2*(N+1));
  end;
end;
       
        

       
F=inline('(sqrt(pi./s))','s');            
dt=1/100; N=12; B=2;
g=ones(2+sum(B.^(1:3)),1);
mu0=8; nu=0.6; sigma=0;
v=ConvolutionOld(F,g,dt,N,B,mu0,nu,sigma);