function [G,Tc]=ttdGvConvolutionMatrix(T,dt,pe,tau)

% USAGE: [G,Tc]=ttdGvConvolutionMatrix(T,dt,pe,tau)
% Returns a convolution matrix with kernel ttdGv(pe,tau) over the interval [0:dt:T]. 
% INPUTS:
%   T  = upper time limit of convolution (scalar)
%   dt = discretization step (same units as T) for convolution (scalar)
%        N=T/dt must be an integer
%   pe = Peclet number for ttdGv
%   tau = Tau for ttdGv
% OUTPUTS:
%   G = sparse convolution matrix of size (N+1,N), such that G*Cbc is the convolution 
%       of ttdGv and Cbc. G*Cbc gives the solution at times t=0,dt,...,T (returned in 
%       vector Tc), given Cbc at times t=dt/2,3*dt/2,...,T-dt/2. 
%       Note: the first row of G is always zero. Thus the returned 
%       solution at time t=0 (first element of G*Cb) is always zero. 
%   Tc = vector of times at which the solution G*Cbc is defined.

% Samar Khatiwala (spk@ldeo.columbia.edu)

N=T/dt;
if rem(T,dt)
  error('T must be divisible by dt')
end
t=[dt/2:dt:T-dt/2]';
Tc=[0:dt:T]';

quadTol=1e-10;

g=ttdGv(t,pe,tau);

% compute 'value' of Gv at dt/2 such that the value of the convolution integral from
% 0 to dt is correctly computed.
if ~verLessThan('matlab','7.5')
  g(1)=(1-quadgk(@(tp)ttdGv(tp,pe,tau),dt,Inf,'AbsTol',quadTol))/dt;
else
  g(1)=(1-quadva(@(tp)ttdGv(tp,pe,tau),[dt Inf],[],quadTol))/dt;  
end
  
G=sparse(N+1,N);
G(2:end,:)=dt*tril(toeplitz(g));