/*qg32 32 point Gaussian quadrature, originally given to Gary F by Jim Lattimer */
#include "qg32.h"

double qg32(
	double xl, /*lower limit to integration range*/
	double xu, /*upper limit to integration range*/
	/*following is the pointer to the function that will be evaulated*/
	double (*fct)(double) )
{
	double a, 
	  b, 
	  c, 
	  y;


	/********************************************************************************
	 *                                                                              *
	 *  32-point Gaussian quadrature                                                *
	 *  xl  : the lower limit of integration                                        *
	 *  xu  : the upper limit                                                       *
	 *  fct : the (external) function                                               *
	 *  returns the value of the integral                                           *
	 *                                                                              *
	 * simple call to integrate sine from 0 to pi                                   *
	 * double agn = qg32( 0., 3.141592654 ,  sin );                                 *
	 *                                                                              *
	 *******************************************************************************/

	a = .5*(xu + xl);
	b = xu - xl;
	c = .498631930924740780*b;
	y = .35093050047350483e-2*((*fct)(a+c) + (*fct)(a-c));
	c = b*.49280575577263417;
	y += .8137197365452835e-2*((*fct)(a+c) + (*fct)(a-c));
	c = b*.48238112779375322;
	y += .1269603265463103e-1*((*fct)(a+c) + (*fct)(a-c));
	c = b*.46745303796886984;
	y += .17136931456510717e-1*((*fct)(a+c) + (*fct)(a-c));
	c = b*.44816057788302606;
	y += .21417949011113340e-1*((*fct)(a+c) + (*fct)(a-c));
	c = b*.42468380686628499;
	y += .25499029631188088e-1*((*fct)(a+c) + (*fct)(a-c));
	c = b*.3972418979839712;
	y += .29342046739267774e-1*((*fct)(a+c) + (*fct)(a-c));
	c = b*.36609105937014484;
	y += .32911111388180923e-1*((*fct)(a+c) + (*fct)(a-c));
	c = b*.3315221334651076;
	y += .36172897054424253e-1*((*fct)(a+c) + (*fct)(a-c));
	c = b*.29385787862038116;
	y += .39096947893535153e-1*((*fct)(a+c) + (*fct)(a-c));
	c = b*.2534499544661147;
	y += .41655962113473378e-1*((*fct)(a+c) + (*fct)(a-c));
	c = b*.21067563806531767;
	y += .43826046502201906e-1*((*fct)(a+c) + (*fct)(a-c));
	c = b*.16593430114106382;
	y += .45586939347881942e-1*((*fct)(a+c) + (*fct)(a-c));
	c = b*.11964368112606854;
	y += .46922199540402283e-1*((*fct)(a+c) + (*fct)(a-c));
	c = b*.7223598079139825e-1;
	y += .47819360039637430e-1*((*fct)(a+c) + (*fct)(a-c));
	c = b*.24153832843869158e-1;
	y = b*(y + .482700442573639e-1*((*fct)(a+c) + (*fct)(a-c)));

	/* the answer */

	return( y );
}