function C=ttdConvolveWithG1d(time,Tdata,Cbdata,Gamma,Delta,lambdaDecay,isSteady)

% USAGE: C = ttdConvolveWithG1d(time,Tdata,Cbdata,Gamma,Delta,[lambdaDecay],[isSteady])
% Returns the convolution of G1d with the tracer boundary condition
% at time t with flow parameters Gamma and Delta.
% INPUTS:
%  time = time (same units as Tdata and Delta)
%  Tdata = vector of times at which boundary data are available
%  Cbdata = vector of boundary data (of length Tdata)
%  Gamma = mean of 1-d TTD
%  Delta = width of 1-d TTD
%  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[G1d(t')*Cb(t-t'),{t',0,t}]
% If lambdaDecay is specified, the integral is:
% C(t) = Integral[G1d(t')*exp(-lambdaDecay*t')*Cb(t-t'),{t',0,t}]
% If isSteady=1, the integral is:
% C(t) = Cb*Integral[G1d(t')*exp(-lambdaDecay*t'),{t',0,Inf}]

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

% Process arguments
if nargin<6
  lambdaDecay=0;
end
if isempty(lambdaDecay)
  lambdaDecay=0;
end

if nargin<7
  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 length(Gamma)~=1 || length(Delta)~=1
    error('ERROR: Gamma and Delta must be scalars')
  end  

% Assume that Cb = 0 for t < Tdata(1) so shift origin
  time=time-Tdata(1);
  Tdata=Tdata-Tdata(1);

  C=repmat(0,size(time));
  for it=1:length(time)
    t=time(it);
    if Delta>=2 || Gamma>t
	  if ~verLessThan('matlab','7.5')
		C(it)=quadgk(@(tp)convolutionIntegrand(tp),0,t,'AbsTol',quadTol);
      else	  
		C(it)=quadva(@(tp)convolutionIntegrand(tp),[0 t],[],quadTol);
	  end
    else % hack to deal with sharp peak
      tp1=max([Gamma-30*Delta 0]);
      tp2=min([Gamma+30*Delta t]);
	  if ~verLessThan('matlab','7.5')
        C(it)=quadgk(@(tp)convolutionIntegrand(tp),0,t,'AbsTol',quadTol,'WayPoints',[tp1 tp2]);      
      else
        tt=[0 tp1 tp2 t]';      
        CC=0;        
		for i=1:length(tt)-1
		  CC=CC+quadva(@(tp)convolutionIntegrand(tp),[tt(i) tt(i+1)],[],quadTol);
		end
        C(it)=CC;              
      end
    end
  end
else % steady state tracer
  if Delta>=1
    if ~verLessThan('matlab','7.5')
      C = Cbdata*quadgk(@(tp)steadyIntegrand(tp),0,Inf,'AbsTol',quadTol);
    else
      C = Cbdata*quadva(@(tp)steadyIntegrand(tp),[0 Inf],[],quadTol);
    end    
  else % hack to deal with sharp peak
	tp1=max([Gamma-30*Delta 0]);
	tp2=Gamma+30*Delta;
	if ~verLessThan('matlab','7.5')
	  C = Cbdata*quadgk(@(tp)steadyIntegrand(tp),0,Inf,'AbsTol',quadTol,'WayPoints',[tp1 tp2]);
	else
	  tt=[0 tp1 tp2 Inf]';      
	  C=0;        
	  for i=1:length(tt)-1
        C=C+Cbdata*quadva(@(tp)steadyIntegrand(tp),[tt(i) tt(i+1)],[],quadTol);
	  end
	end
  end  
end

% private (nested) functions
  function GC=convolutionIntegrand(tp)
    GC=ttdG1d(tp,Gamma,Delta).*exp(-lambdaDecay*tp).*interp1(Tdata,Cbdata,t-tp);
  end
  
  function GC=steadyIntegrand(tp)
    GC=ttdG1d(tp,Gamma,Delta).*exp(-lambdaDecay*tp);
%     [tp' GC']
  end  
end