function q = gauss4 ( a, b, f )

  n = 4;

  x = zeros ( n, 1 );
  x(1) = -0.861136311594052575223946488893;
  x(2) = -0.339981043584856264802665759103;
  x(3) = 0.339981043584856264802665759103;
  x(4) = 0.861136311594052575223946488893;

  w = zeros ( n, 1 );
  w(1) = 0.347854845137453857373063949222;
  w(2) = 0.652145154862546142626936050778;
  w(3) = 0.652145154862546142626936050778;
  w(4) = 0.347854845137453857373063949222;
%
%  Transform from [-1,+1] to [a,b].
%
  x = ( ( 1.0 - x ) * a + ( 1.0 + x ) * b ) / 2.0;
  w = w * ( b - a ) / 2.0;

  q = w' * f ( x );

  return
end
