function xint = gauss24pts(f, a, b)
% 1D Gauss-Legendre quadrature using 48 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 = 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=[...
3.238017096286936D-02
.0970046992094627
.1612223560688917
.2247637903946891
.2873624873554556
.3487558862921608
.4086864819907167
.4669029047509584
.523160974722233
.5772247260839727
.6288673967765136
.6778723796326639
.7240341309238146
.7671590325157404
.8070662040294426
.8435882616243935
.8765720202742479
.9058791367155696
.9313866907065543
.9529877031604309
.9705915925462473
.9841245837228269
.9935301722663508
.9987710072524261];
w=[...
6.473769681268392D-02
6.446616443595009D-02
6.392423858464819D-02
6.311419228625402D-02
6.203942315989266D-02
6.070443916589388D-02
5.911483969839564D-02
5.727729210040321D-02
5.519950369998416D-02
5.289018948519367D-02
5.035903555385447D-02
4.761665849249048D-02
4.467456085669428D-02
4.154508294346475D-02
3.824135106583071D-02
3.477722256477044D-02
3.116722783279809D-02
2.742650970835695D-02
2.357076083932438D-02
1.961616045735553D-02
1.557931572294385D-02
1.147723457923454D-02
7.327553901276262D-03
3.153346052305838D-03];
       
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)));
 
%fprintf('\nIntegral = %g\n', xint);
