function C=ttdConvolveWithGv(time,Tdata,Cbdata,Pe,tau,epsilon,lambdaDecay,isSteady)

% USAGE: C = ttdConvolveWithGv(time,Tdata,Cbdata,Pe,tau,[epsilon],[lambdaDecay],[isSteady])
% Returns the convolution of Gv with the tracer boundary condition
% at time t with flow parameters Pe and tau.
% INPUTS:
%  time = time (same units as Tdata and tau)
%  Tdata = vector of times at which boundary data are available
%  Cbdata = vector of boundary data (of length Tdata)
%  Pe = Peclet number
%  tau = mean age
%  epsilon = optional parameter used to split the convolution 
%            integral. If omitted or is an empty array ([]) a default value 
%            of 1 year is used.
%  lambdaDecay = optional decay constant [1/year] for a radioactive tracer.
%                If omitted or is an empty array ([]) a default value 
%                of 0 is used.
%  isSteady = optional flag indicating whether the tracer is a steady state 
%             tracer or not. If isSteady=1, the steady state concentration is 
%             computed. The arguments 'time' and 'Tdata' are then ignored, and 
%             Cbdata is assumed to be a scalar.

% Compatibility note:
% This function uses an adaptive quadrature routine to evaluate the convolution integral. 
% The builtin 'quad' function is not very robust and is also quite slow. Matlab 7.5 introduced 
% a vectorized quad ('quadgk.m') that is much much better. However, if you don't have version 
% 7.5, you can download an alternative called 'quadva.m' written by Lawrence Shampine of SMU. 
% You can download it from his website: http://faculty.smu.edu/shampine/quadva.m or email me 
% for a copy (spk@ldeo.columbia.edu).
% This routine checks whether you are using version 7.5. If you are, then it will use quadgk.m. 
% Otherwise it will use quadva.m.
                
% Notes:
% C(t) = Integral[Gv(t')*Cb(t-t'),{t',0,t}]
% C(t) = I1(eps) + I2(eps)
%   I1(eps) = Integral[Gv(t')*Cb(t-t'),{t',0,eps}]
%           ~ <Cb>*Integral[Gv(t'),{t',0,eps}]
%        where, 
%          <Cb> = mean of Cb(t) in interval [0,eps]
%   I2(eps) = Integral[Gv(t')*Cb(t-t'),{t',eps,t}]
% If lambdaDecay is specified, the integral is:
% C(t) = Integral[Gv(t')*exp(-lambdaDecay*t')*Cb(t-t'),{t',0,t}]
% If isSteady=1, the integral is:
% C(t) = Cb*Integral[Gv(t')*exp(-lambdaDecay*t'),{t',0,Inf}]

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

% Process arguments
if nargin<6
  epsilon=1;
end
if isempty(epsilon)
  epsilon=1;
end

if nargin<7
  lambdaDecay=0;
end
if isempty(lambdaDecay)
  lambdaDecay=0;
end

if nargin<8
  isSteady=0;
end
if isempty(isSteady)
  isSteady=0;
end

quadTol=1e-10;

if ~isSteady % transient tracer
  if any(time<Tdata(1)) | any(time>Tdata(end))
    error('ERROR: time must be bracketed by Tdata(1) and Tdata(end)')
  end  
  if any(time<epsilon)
    error('ERROR: time must be >= epsilon')
  end  
  if length(Pe)~=1 | length(tau)~=1
    error('ERROR: Pe and tau must be scalars')
  end  

% Assume that Cb = 0 for t < Tdata(1) so shift origin
  time=time-Tdata(1);
  Tdata=Tdata-Tdata(1);
  if ~verLessThan('matlab','7.5')  
    Gv1=1-quadgk(@(tp)ttdGv(tp,Pe,tau),epsilon,Inf,'AbsTol',quadTol);
  else
    Gv1=1-quadva(@(tp)ttdGv(tp,Pe,tau),[epsilon Inf],[],quadTol);  
  end
  C=repmat(0,size(time));
  for it=1:length(time)
    t=time(it);
    Cbbar=mean(interp1(Tdata,Cbdata,[t-epsilon t]'));
    I1=Gv1*Cbbar;    
    if ~verLessThan('matlab','7.5')      
      I2=quadgk(@(tp)convolutionIntegrand(tp),epsilon,t,'AbsTol',quadTol);    
    else
      I2=quadva(@(tp)convolutionIntegrand(tp),[epsilon t],[],quadTol);        
    end
    C(it)=I1+I2;
  end
else % steady state tracer
  if ~verLessThan('matlab','7.5')  
    Gv1=1-quadgk(@(tp)ttdGv(tp,Pe,tau),epsilon,Inf,'AbsTol',quadTol);
  else
    Gv1=1-quadva(@(tp)ttdGv(tp,Pe,tau),[epsilon Inf],[],quadTol);  
  end
  expbar=mean(exp(-lambdaDecay*linspace(0,epsilon,10)));
  I1=Gv1*expbar;
  if ~verLessThan('matlab','7.5')    
    I2=quadgk(@(tp)steadyIntegrand(tp),epsilon,Inf,'AbsTol',quadTol);  
  else
    I2=quadva(@(tp)steadyIntegrand(tp),[epsilon Inf],[],quadTol);    
  end
  C=Cbdata*(I1+I2);
end

% private (nested) functions
  function GC=convolutionIntegrand(tp)
    GC=ttdGv(tp,Pe,tau).*exp(-lambdaDecay*tp).*interp1(Tdata,Cbdata,t-tp);
  end
  
  function GC=steadyIntegrand(tp)
    GC=ttdGv(tp,Pe,tau).*exp(-lambdaDecay*tp);
  end  
end