Numerical integration with Gaussian Quadrature

#pragma rtGlobals=3		// Use modern global access method and strict wave access.

//Perform Gaussian Quadrature Integration of a given function.
//This is slightly different to the inbuilt Integrate1D in that one can pass in a wave containing wave references as extra input 
//to the function to be integrated!
//The function to be integrated can also be threadsafe to speed up calculation.

//There is also a threadsafe version of the integration code (integrate1D is not threadsafe). In this case the function to be integrated also needs to be threadsafe.

//Lookup tables are used to store commonly used values of the Gauss-Legendre nodes.
//The GaussLegendre node calculation is based on code from the NIST SANS macros, which was originally based on code from
//Numerical recipes.
//Type test() for a quick test of the code.
//Written by Andrew Nelson July 2013

Function GQIntegrate1D(fcn, min_x, max_x, w, N, [tol, rtol, adaptive, startiter])
	//a function for performing Gaussian Quadrature integration with a set order of quadrature points
	//INPUT
	//fcn 			- the function to be integrated.  Has the signature:
	//						Function fcn(w, y, x)
	//							wave/wave w
	//							wave y, x
	//							y = f(x)
	//						End
	//					If the function is declared threadsafe, then the integration is done in a parallel manner.   The function is written
	//					in much the same manner as an all-at-once fit function.  i.e. All the abscissa (xpoints) are passed in at once, and the function
	//					is supposed to populate all the y values at once.
	//min_x		- the lower limit for the integration
	//max_x		- the upper limit for the integration
	//w			- wave containing wave references to pass into the function fcn.  Use this to pass data into your function
	//N				- number of quadrature points to be used. If the adaptive variable is set then this is the maximum order examined.
	//
	//OPTIONAL (for adaptive quadrature)
	//				NOTE: adaptive quadrature will only take place if the adaptive variable is set to non-zero.
	//tol			- the Iteration stops when the error between successive iterates is less than tol OR the relative change is less than rtol. default = 1e-13
	//rtol			- the Iteration stops when the error between successive iterates is less than tol OR the relative change is less than rtol.
	//startiter		- starting order for adaptive gaussian quadrature. Default = 5
	//adaptive		- Specify this variable to non-zero if you wish to do adaptive quadrature.
	
	
	string fcn
	Variable min_x, max_x		
	Wave w						
	Variable N, tol, rtol, adaptive, startiter

	// local variables
	variable ii, jj, va, vb, isthreadsafe, NORD, err, val, maxiter
	string info = "", funcname = "", fri = ""
	
	if(paramisdefault(tol))
		tol = 1.e-13
	endif
	if(paramisdefault(rtol))
		rtol = 1.49e-08
	endif
	if(paramisdefault(startiter))
		startiter = 5
	endif

	if(adaptive == 0)
		maxiter = N
		startiter = N
		make/free/n=(maxiter - startiter + 1)/d W_gaussQuad
	else
		maxiter = N
		make/o/n=(maxiter - startiter + 1)/d W_gaussQuad
		W_gaussQuad = NaN
	endif
	
	funcname = fcn
	info = functioninfo(funcname)
	isthreadsafe = stringmatch(stringbykey("THREADSAFE", info), "YES")

	//limits of integration are input to function
	va = min_x
	vb = max_x

	for(NORD = startiter, jj = 0 ; NORD <= maxiter ; NORD += 1,  jj += 1)
		//calculate the Gauss Legendre abscissae and weights
		Wave GLpoints = gauss_legendre_tbl(NORD)	
		
		//place to put evaluated data  and x point to be evaluated
		make/n=(NORD)/d/free yy, zi
		yy = 0
	
		// calculate Gauss points on integration interval (q-value for evaluation)
		Multithread zi[] = ( GLpoints[p][0] * (vb - va) + vb + va ) / 2.0
		
		//evaluate function at all the points
		if(isthreadsafe)
			Funcref TGaussQuad_proto cfcn = $fcn
			fri = funcrefinfo(cfcn)
			if(numberbykey("ISPROTO", fri) == 0)
				cfcn(w, yy, zi)
			else
				print "Function does not satisfy the signature:"
				print " Threadsafe Function TGaussQuad_proto(w, y, x)"
				return NaN
			endif
		else
			Funcref GaussQuad_proto cfcn1 = $fcn
			fri = funcrefinfo(cfcn1)
			if(numberbykey("ISPROTO", fri) == 0)
				cfcn1(w, yy, zi)
			else
				print "Function does not satisfy the signature:"
				print  "Function GaussQuad_proto(w, y, x)"
				return NaN
			endif
		endif
		
		//multiply by weights to calculate partial sum
		Multithread yy[] *= GLpoints[p][1]
		W_gaussQuad[jj] = (vb - va) / 2.0 * sum(yy)
		val = W_gaussQuad[jj]
		if(jj)
			err = abs(val - W_gaussQuad[jj - 1])
			if(err < tol || err < rtol * abs(val))
				break
			endif
		endif
	endfor
	
	//calculate value to return
	return  val
End

Threadsafe Function TGQIntegrate1D(fcn, min_x, max_x, w, N, [tol, rtol, adaptive, startiter])
	//a threadsafe function for performing Gaussian Quadrature integration with a set order of quadrature points
	//INPUT
	//fcn 			- the function to be integrated.  Has the signature:
	//						Threadsafe Function fcn(w, y, x)
	//							wave/wave w
	//							wave y, x
	//							y = f(x)
	//						End
	//min_x			- the lower limit for the integration
	//max_x		- the upper limit for the integration
	//w			- wave containing wave references to pass into the function fcn.  Use this to pass data into your function
	//N				- number of quadrature points to be used. If the adaptive variable is set then this is the maximum order examined.
	//
	//OPTIONAL (for adaptive quadrature)
	//				NOTE: adaptive quadrature will only take place if the adaptive variable is set to non-zero.
	//tol			- the Iteration stops when the error between successive iterates is less than tol OR the relative change is less than rtol. default = 1e-13
	//rtol			- the Iteration stops when the error between successive iterates is less than tol OR the relative change is less than rtol.
	//startiter		- starting order for adaptive gaussian quadrature. Default = 5
	//adaptive		- Specify this variable to non-zero if you wish to do adaptive quadrature.
	
	string fcn
	Variable min_x, max_x		
	Wave w						
	Variable N, tol, rtol, adaptive, startiter

	// local variables
	variable ii, jj, va, vb, isthreadsafe, NORD, err, val, maxiter
	string info = "", funcname = "", fri = ""
	
	if(paramisdefault(tol))
		tol = 1.e-13
	endif
	if(paramisdefault(rtol))
		rtol = 1.49e-08
	endif
	if(paramisdefault(startiter))
		startiter = 5
	endif
	
	if(adaptive == 0)
		maxiter = N
		startiter = N
		make/free/n=(maxiter - startiter + 1)/d W_gaussQuad
	else
		maxiter = N
		make/o/n=(maxiter - startiter + 1)/d W_gaussQuad
		W_gaussQuad = NaN
	endif
	
	//limits of integration are input to function
	va = min_x
	vb = max_x

	for(NORD = startiter, jj = 0 ; NORD <= maxiter ; NORD += 1,  jj += 1)
		//calculate the Gauss Legendre abscissae and weights
		Wave GLpoints = gauss_legendre(NORD)
	
		//place to put evaluated data  and x point to be evaluated
		make/n=(NORD)/d/free yy, zi

		// calculate Gauss points on integration interval (q-value for evaluation)
		Multithread zi[] = ( GLpoints[p][0] * (vb - va) + vb + va ) / 2.0
	
		//evaluate function at each point
		Funcref TGaussQuad_proto cfcn = $fcn
		fri = funcrefinfo(cfcn)
		if(numberbykey("ISPROTO", fri) == 0)
			cfcn(w, yy, zi)
		else
			print "Function does not satisfy the signature:"
			print " Threadsafe Function TGaussQuad_proto(w, y, x)"
			return NaN
		endif
	
		//multiply by weights to calculate partial sum
		Multithread yy[] *= GLpoints[p][1]
		W_gaussQuad[jj] = (vb - va) / 2.0 * sum(yy)
		val = W_gaussQuad[jj] 
		
		if(jj)
			err = abs(val - W_gaussQuad[jj - 1])
			if(err < tol || err < rtol * abs(val))
				break
			endif
		endif
	endfor
	
	//calculate value to return
	return  val
End

Function GaussQuad_proto(w, y, x)
	Wave/wave w
	Wave y, x
	print  "in GaussQuad_proto function"
	return Inf
end

Threadsafe Function TGaussQuad_proto(w, y, x)
	Wave/wave w
	wave y, x
	print  "in TGaussQuad_proto function"
	return Inf
end

Function/wave gauss_legendre_tbl(N)
	variable N
	///routine from Numerical Recipes to calculate weights and abscissae for 
	// Gauss-Legendre quadrature, over the range (-1, 1)
	//
	//Use this function in preference to gauss_legendre. This function creates a lookup table and is faster
	//to use in the longrun, compared to gauss_legendre which calculates new values each time
	//INPUT - order for Gauss Legendre (>> 1)
	//
	//RETURNS - free wave containing Gauss abscissa (first column) and weights (second column).
	//			The wave has dimensions [N][2]
	
	//	variable timer = startmstimer
	variable newsize, original_idx_size
	
	if(!datafolderexists("root:packages:GLquad"))
		DFREF saveDFR = GetDataFolderDFR()
		newdatafolder/o root:packages
		newdatafolder/o/s root:packages:GLquad
		make/n=0/i/u/o idx
		make/n=(0, 0)/d/o wt = 0, zi = 0
		SetDataFolder saveDFR
	endif
	
	Wave wt = root:packages:GLquad:wt
	Wave zi = root:packages:GLquad:zi
	Wave idx = root:packages:GLquad:idx

	findvalue/I=(N)/z idx
	if(V_Value == -1)
		//haven't created those values yet
		original_idx_size = numpnts(idx)
		if(N > dimsize(wt, 0))
			newsize = N
		else
			newsize = dimsize(wt, 0)
		endif
		redimension/n=(original_idx_size + 1) idx
		redimension/n=(newsize, original_idx_size + 1) zi, wt
		Wave vals = gauss_legendre(N)
		idx[numpnts(idx) - 1] = N
		zi[0, N - 1][numpnts(idx) - 1] = vals[p][0]
		wt[0, N - 1][numpnts(idx) - 1] = vals[p][1]
		//		print stopmstimer(timer)/1e6	
		return vals	
	else
		make/n=(N, 2)/d/free GLpoints
		GLpoints[][0] = zi[p][V_Value]
		GLpoints[][1] = wt[p][V_Value]
		//		print stopmstimer(timer)/1e6	
		return GLpoints
	endif
End
	
Threadsafe Function/wave gauss_legendre(N)
	///routine from Numerical Recipes to calculate weights and abscissae for 
	// Gauss-Legendre quadrature.
	//Code modified from GaussUtils_v40.ipf, written by Steve Kline @NCNR NIST
	//
	//INPUT - order for Gauss Legendre (>> 1)
	//
	//RETURNS - free wave containing Gauss abscissa (first column) and weights (second column).
	//			The wave has dimensions [N][2]
	
	variable N

	Variable x1, x2
	variable m, jj, ii
	variable z1, z, xm, xl, pp, p3, p2, p1
	Variable eps = 3e-11
	
	make/free/n=(N, 2)/d GLpoints

	x1 = -1
	x2 = 1
	
	m = (N+1)/2	
	xm = 0.5 * (x2 + x1)
	xl = 0.5 * (x2 - x1)
	
	for (ii = 1; ii <= m; ii += 1) 
		z = cos(pi * (ii - 0.25)/(N + 0.5))
		do 
			p1 = 1.0
			p2 = 0.0
			for (jj = 1;jj <= N;jj += 1) 
				p3 = p2
				p2 = p1
				p1 = ((2.0 * jj - 1.0) * z * p2 - (jj - 1.0) * p3) / jj
			endfor
			pp = N * (z * p1 - p2) / (z * z -1.0)
			z1 = z
			z = z1 - p1 / pp
		while (abs(z - z1) > EPS)
		GLpoints[ii - 1][0] = xm - xl * z
		GLpoints[N  - ii][0] = xm + xl * z
		GLpoints[ii - 1][1] = 2.0 * xl / ((1.0 - z * z) * pp * pp)
		GLpoints[N  - ii][1] = GLpoints[ii - 1][1]
	Endfor
	
	return GLpoints
End

Function test_quadratic(w, y, x)
	//Lets integrate x^2
	Wave/wave w
	wave y, x
	y = x^2
End

Function test_sine(w, y, x)
	//Lets integrate x^2
	Wave/wave w
	wave y, x
	y = sin(1/x)
End

Function test_sine_IGOR(x)
variable x
return sin(1/x)
End

Function test()
	variable x1, x2, true, err
	Wave/wave w
	//quadratic
	x1 = 0
	x2 = 1
	print "Integrating x^2 between", x1, "and", x2
	print/d GQIntegrate1D("test_quadratic", x1, x2, w, 100)
	print "with adaptive quadrature"
	print/d GQIntegrate1D("test_quadratic", x1, x2, w, 100, adaptive = 1, startiter = 3)
	Wave W_gaussQuad
	print/d W_gaussQuad
	
	//sin(x)
	x1 = 0.001
	x2 = pi / 2
	print "Integrating sin(1/x) between", x1, "and", x2
	print "IGOR's version"
	true = integrate1D(test_sine_igor, 0.001, pi/2, 2, 0)
	print/d true
	err = GQIntegrate1D("test_sine", x1, x2, w, 100)
	print/d true - err
	print "with adaptive quadrature"
	err =  GQIntegrate1D("test_sine", x1, x2, w, 1000, adaptive = 1, startiter = 2)
	print true - err
	Wave W_gaussQuad
	print/d W_gaussQuad
End