function xint = gauss4(f, a, b, p)
% 1D Gauss-Legendre quadrature using 4 (8) points.
% INPUT:
% f = function (a string name for an m-file) tales
% a vector argument & returns a vector of values
% a = lower limit of integration
% b = upper limit of integration
% p = is a parameter of the function f
% Local variables:
% r = 24 roots of Gauss-Legendre quadrature
% w = 24 wts of Gauss-Legendre quadrature
% xint = value of the integral
% OUTPUT:
% xint = value of the 24-pt Gauss-Legendre quadrature
% K. Ming Leung, 03/15/03

r=[0.960289856497536 0.796666477413627 0.525532409916329 0.183434642495650];
w=[0.101228536290376 0.222381034453374 0.313706645877887 0.362683783378362];

A = (b-a)/2; % parameter A for current segment
B = (b+a)/2; % parameter B for current segment
Ar = A*r;
xint = A*sum(w.*(feval(f,B+Ar,p)+feval(f,B-Ar,p)));