function q = gauss2 ( a, b, f )

  x1 = - sqrt ( 1.0 / 3.0 );
  x2 = + sqrt ( 1.0 / 3.0 );
  w1 = 1.0;
  w2 = 1.0;
%
%  Transform from [-1,+1] to [a,b].
%
  x1 = ( ( 1.0 - x1 ) * a + ( 1.0 + x1 ) * b ) / 2.0;
  x2 = ( ( 1.0 - x2 ) * a + ( 1.0 + x2 ) * b ) / 2.0;
  w1 = w1 * ( b - a ) / 2.0;
  w2 = w2 * ( b - a ) / 2.0;

  q = w1 * f ( x1 ) + w2 * f ( x2 );

  return
end
