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