function q=integrate(f1,k1,f2,k2,h,a,b)
% q=integrate(f1,k1,f2,k2,h,a,b)
% compute the integral from a to b of f1(k1,h,x)*f2(k2,h,x)dx
% using 3-point Gaussian quadrature
% f1 and f2 are strings with the names of functions in them
% f1 and f2 have signature y=f(k,h,x)
% The functions f1 and f2 are assumed to be piecewise polynomials
%   with pieces [a+n*h,a+(n+1)*h]
% If the product f1*f2 is no higher than fifth degree, the
% 3-point Gauss formula used here is exact.

% M. M. Sussman
% $Id: integrate.m,v 1.1 2011/01/14 22:36:18 mike Exp $

% parameters for 3-point Gauss quadrature over [-1,1]
% quantities are halved because length of [-1,1] is 2
quadpt=0.774596669241483377035853079956/2;
wa=5/18;
wb=4/9;
wc=5/18;

if b<=a
  error('integrate: interval [a,b] must have positive length')
end
if h<=0
  error('integrate: h must be >0')
end

xleft=a;
xright=min(a+h,b);
len=xright-xleft;

q=0;
while len> 100*eps*h    % guard against roundoff

  xb=(xleft+xright)/2;
  xa=xb-len*quadpt;
  xc=xb+len*quadpt;

  q=q+ len*(wa*feval(f1,k1,h,xa)*feval(f2,k2,h,xa) + ...
            wb*feval(f1,k1,h,xb)*feval(f2,k2,h,xb) + ...
            wc*feval(f1,k1,h,xc)*feval(f2,k2,h,xc) );

  xleft=xright;
  xright=min(xleft+h,b);
  len=xright-xleft;

end
