function xint = gauss48pts(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 = 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=[...
1.627674484960297D-02
4.881298513604973D-02
8.129749546442556D-02
.1136958501106659
.1459737146548969
.1780968823676186
.2100313104605672
.24174315616384
.2731988125910492
.3043649443544963
.3352085228926254
.3656968614723136
.3957976498289086
.4254789884073005
.454709422167743
.4834579739205964
.5116941771546677
.5393881083243575
.5665104185613972
.593032364777572
.6189258401254686
.6441634037849671
.6687183100439161
.6925645366421715
.7156768123489676
.7380306437444001
.7596023411766475
.7803690438674332
.8003087441391408
.8194003107379316
.8376235112281871
.8549590334346014
.8713885059092965
.8868945174024204
.9014606353158523
.9150714231208981
.9277124567223087
.9393703397527552
.9500327177844377
.9596882914487426
.9683268284632642
.9759391745851365
.9825172635630147
.9880541263296237
.9925439003237626
.9959818429872093
.9983643758631817
.9996895038832307 ];
w=[...
3.255061449236317D-02
3.251611871386884D-02
3.244716371406427D-02
3.234382256857593D-02
3.220620479403025D-02
3.203445623199266D-02
3.182875889441101D-02
3.158933077072717D-02
3.131642559686135D-02
3.101033258631384D-02
3.067137612366915D-02
.0302999154208276
2.989634413632838D-02
2.946108995816791D-02
2.899461415055524D-02
2.849741106508539D-02
2.797000761684833D-02
2.741296272602924D-02
2.682686672559176D-02
2.621234073567241D-02
2.557003600534936D-02
2.490063322248361D-02
2.420484179236469D-02
2.348339908592622D-02
2.273706965832937D-02
2.196664443874435D-02
.0211729398921913
2.035679715433332D-02
1.951908114014502D-02
1.866067962741147D-02
1.778250231604526D-02
1.688547986424517D-02
1.597056290256229D-02
1.503872102699494D-02
1.409094177231486D-02
1.312822956696157D-02
1.215160467108832D-02
.0111621020998385
1.016077053500842D-02
9.148671230783386D-03
8.126876925698759D-03
7.096470791153865D-03
6.058545504235961D-03
5.014202742927518D-03
3.964554338444687D-03
2.910731817934947D-03
1.853960788946922D-03
7.967920655520124D-04 ];

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);