function xint = gauss24pts(f, a, b)
% 1D Gauss-Legendre quadrature using 24 points.
% INPUT:
% f = integrand function (a string) must be in an m-file
% a = lower limit of integration
% b = upper limit of integration
% Local variables:
% r = positive roots of Gauss-Legendre quadrature 
% w = wts of Gauss-Legendre quadrature
% xint = value of the integral
% OUTPUT:
% xint = value of the Gauss-Legendre quadrature
% K. Ming Leung, 08/11/03

% 24 pts
r = [
    6.40568928626056260827e-02
    1.91118867473616309153e-01
    3.15042679696163374398e-01
    4.33793507626045138478e-01
    5.45421471388839535649e-01
    6.48093651936975569268e-01
    7.40124191578554364260e-01
    8.20001985973902921981e-01
    8.86415527004401034190e-01
    9.38274552002732758539e-01
    9.74728555971309498199e-01
    9.95187219997021360195e-01 ];

w = [
    1.27938195346752156976e-01
    1.25837456346828296117e-01
    1.21670472927803391202e-01
    1.15505668053725601353e-01
    1.07444270115965634785e-01
    9.76186521041138882720e-02
    8.61901615319532759152e-02
    7.33464814110803057346e-02
    5.92985849154367807461e-02
    4.42774388174198061695e-02
    2.85313886289336631809e-02
    1.23412297999871995469e-02 ];


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)+feval(f,B-Ar)));
