//+++ v. 2016-05-02

#ifndef integrals_newest_h
#define integrals_newest_h

#include <gsl/gsl_integration.h>
#include <gsl/gsl_sf_bessel.h>


//++++++++++++++++++++++++++++++++++++++++++++++++++++
//+++  QNG non-adaptive Gauss-Kronrod integration  +++++++ [1]
//++++++++++++++++++++++++++++++++++++++++++++++++++++
/*
This function applies the Gauss-Kronrod 10-point, 21-point, 43-point and 87-point integration rules in succession until an estimate of the integral of f over (a,b) is achieved within the desired absolute and relative error limits, epsabs and epsrel. The function returns the final approximation, result, an estimate of the absolute error, abserr and the number of function evaluations used, neval. The Gauss-Kronrod rules are designed in such a way that each rule uses all the results of its predecessors, in order to minimize the total number of function evaluations.
*/

int integral_QNG(const gsl_function * f, double a, double b, double epsabs, double epsrel, double * result, double *abserr, size_t * neval)
{
	return gsl_integration_qng (f,a,b,epsabs,epsrel,result,abserr,neval);
}

int integral_QNG(double (*fun)(double, void*), void* param, double a, double b, double epsabs, double epsrel, double * result, double *abserr, size_t * neval)
{
    gsl_function f;
    f.params=param;
    f.function=fun;
    
    int res = integral_QNG(&f,a,b,epsabs,epsrel,result,abserr,neval);
    return res;
}

//++++++++++++++++++++++++++++++++++++++++++++++++++++++
//+++ QAGS adaptive integration with !!! singularities !!! +++++ [2]
//++++++++++++++++++++++++++++++++++++++++++++++++++++++
/*
The presence of an integrable singularity in the integration region causes an adaptive routine to concentrate new subintervals around the singularity. As the subintervals decrease in size the successive approximations to the integral converge in a limiting fashion. This approach to the limit can be accelerated using an extrapolation procedure. The QAGS algorithm combines adaptive bisection with the Wynn epsilon-algorithm to speed up the integration of many types of integrable singularities. 
This function applies the Gauss-Kronrod 21-point integration rule adaptively until an estimate of the integral of f over (a,b) is achieved within the desired absolute and relative error limits, epsabs and epsrel. The results are extrapolated using the epsilon-algorithm, which accelerates the convergence of the integral in the presence of discontinuities and integrable singularities. The function returns the final approximation from the extrapolation, result, and an estimate of the absolute error, abserr. The subintervals and their results are stored in the memory provided by workspace. The maximum number of subintervals is given by limit, which may not exceed the allocated size of the workspace.
*/

//++++++++++++++++++++++++++++++++++++++++++++
int integral_QAGS (const gsl_function * f, double a, double b, double epsabs, double epsrel, size_t limit, double *result, double * abserr)
{
	gsl_integration_workspace * w = gsl_integration_workspace_alloc (limit);
	int res=gsl_integration_qags (f,a,b, epsabs,epsrel,limit, w, result, abserr);
	gsl_integration_workspace_free(w);
	return res;
}

int integral_QAGS (double (*fun)(double, void*), void* param, double a, double b, double epsabs, double epsrel, size_t limit, double *result, double * abserr)
{
    gsl_function f;
    f.params=param;
    f.function=fun;
    
    int res= integral_QAGS (&f,a,b,epsabs,epsrel,limit,result,abserr);
    return res;
}

//++++++++++++++++++++++++++++++++++++++++++++++++++++++
//+++ QAG adaptive integration with !!! singularities !!! +++++ [3]
//++++++++++++++++++++++++++++++++++++++++++++++++++++++
/*
This function applies an integration rule adaptively until an estimate of the integral of f over (a,b) is achieved within the desired absolute and relative error limits, epsabs and epsrel. The function returns the final approximation, result, and an estimate of the absolute error, abserr. The integration rule is determined by the value of key, which should be chosen from the following symbolic names,

              GSL_INTEG_GAUSS15  (key = 1)
              GSL_INTEG_GAUSS21  (key = 2)
              GSL_INTEG_GAUSS31  (key = 3)
              GSL_INTEG_GAUSS41  (key = 4)
              GSL_INTEG_GAUSS51  (key = 5)
              GSL_INTEG_GAUSS61  (key = 6)

corresponding to the 15, 21, 31, 41, 51 and 61 point Gauss-Kronrod rules. The higher-order rules give better accuracy for smooth functions, while lower-order rules save time when the function contains local difficulties, such as discontinuities.

On each iteration the adaptive integration strategy bisects the interval with the largest error estimate. The subintervals and their results are stored in the memory provided by workspace. The maximum number of subintervals is given by limit, which may not exceed the allocated size of the workspace. 

*/

//++++++++++++++++++++++++++++++++++++++++++++
int integral_QAG (const gsl_function * f, double a, double b, double epsabs, double epsrel, size_t limit, double *result, double * abserr, int key)
{
	static gsl_integration_workspace * w;
	static int wStatus;

               if (wStatus!=111)
	{
	gsl_integration_workspace_free(w);	
	w= gsl_integration_workspace_alloc (limit);	
	wStatus=111;
	}

	gsl_integration_qag (f,a,b, epsabs,epsrel,limit,key, w, result, abserr); int res=w->size;
//	gsl_integration_workspace_free(w);
	return res;
}

int integral_QAG (double (*fun)(double, void*), void* param, double a, double b, double epsabs, double epsrel, size_t limit, double *result, double * abserr, int key)
{
    gsl_function f;
    f.params=param;
    f.function=fun;
    
    return integral_QAG (&f,a,b, epsabs,epsrel,limit,result,abserr,key);
}


int integral_QAG_fast (const gsl_function * f, double a, double b, double epsabs, double epsrel, size_t limit, double *result, double * abserr, int key, gsl_integration_workspace * w)
{
	int res=gsl_integration_qag (f,a,b, epsabs,epsrel,limit,key, w, result, abserr);
	return w->size;
}

int integral_QAG_fast (double (*fun)(double, void*), void* param, double a, double b, double epsabs, double epsrel, size_t limit, double *result, double * abserr, int key, gsl_integration_workspace * w)
{
    gsl_function f;
    f.params=param;
    f.function=fun;
    
    return integral_QAG_fast (&f,a,b, epsabs,epsrel,limit,result,abserr,key, w);	
}
//+++++++++++++++++++++++++++++++++++++++++++++
//+++ QAGP adaptive integration with known singular points ++ [4]
//+++++++++++++++++++++++++++++++++++++++++++++
/*

This function applies the adaptive integration algorithm QAGS taking account of the user-supplied locations of singular points. The array pts of length npts should contain the endpoints of the integration ranges defined by the integration region and locations of the singularities. For example, to integrate over the region (a,b) with break-points at x_1, x_2, x_3 (where a < x_1 < x_2 < x_3 < b) the following pts array should be used

              pts[0] = a
              pts[1] = x_1
              pts[2] = x_2
              pts[3] = x_3
              pts[4] = b

with npts = 5.

If you know the locations of the singular points in the integration region then this routine will be faster than QAGS. 
*/
//++++++++++++++++++++++++++++++++++++++++++++
int integral_QAGP (const gsl_function * f, double * pts, size_t npts, double epsabs, double epsrel, size_t limit, double *result, double * abserr)
{
	gsl_integration_workspace * w = gsl_integration_workspace_alloc (limit);
	int res=gsl_integration_qagp(f, pts, npts, epsabs, epsrel,limit, w, result, abserr);
	gsl_integration_workspace_free(w);
	return res;
}


int integral_QAGP (double (*fun)(double, void*), void* param, double * pts, size_t npts, double epsabs, double epsrel, size_t limit, double *result, double * abserr)

{
    gsl_function f;
    f.params=param;
    f.function=fun;
    
	int res= integral_QAGP (&f, pts,npts,epsabs,epsrel,limit,result,abserr);
    return res;
}

//+++++++++++++++++++++++++++++++++++++++++++++
//+++ QAGI adaptive integration on infinite intervals  [-inf,+inf] + [5]
//+++++++++++++++++++++++++++++++++++++++++++++
/*

This function computes the integral of the function f over the infinite interval (-\infty,+\infty). The integral is mapped onto the semi-open interval (0,1] using the transformation x = (1-t)/t.

It is then integrated using the QAGS algorithm. The normal 21-point Gauss-Kronrod rule of QAGS is replaced by a 15-point rule, because the transformation can generate an integrable singularity at the origin. In this case a lower-order rule is more efficient. 

*/

//++++++++++++++++++++++++++++++++++++++++++++

int integral_QAGI (gsl_function * f, double epsabs, double epsrel, size_t limit, double *result, double * abserr)
{
	gsl_integration_workspace * w = gsl_integration_workspace_alloc (limit);
	int res=gsl_integration_qagi(f,epsabs,epsrel,limit, w, result, abserr);
	gsl_integration_workspace_free(w);
	return res;
}

int integral_QAGI (double (*fun)(double, void*), void* param, double epsabs, double epsrel, size_t limit, double *result, double * abserr)
{
    gsl_function f;
    f.params=param;
    f.function=fun;
    
	int res= integral_QAGI (&f, epsabs,epsrel,limit,result,abserr);
    return res;
}
//+++++++++++++++++++++++++++++++++++++++++++++
//++ QAGIU adaptive integration on infinite intervals [a,+inf] ++ [6]
//+++++++++++++++++++++++++++++++++++++++++++++
/*

This function computes the integral of the function f over the semi-infinite interval (a,+\infty). The integral is mapped onto the semi-open interval (0,1] using the transformation x = a + (1-t)/t, and then integrated using the QAGS algorithm. 

*/

//++++++++++++++++++++++++++++++++++++++++++++

int integral_QAGIU (gsl_function * f, double a, double epsabs, double epsrel, size_t limit, double *result, double * abserr)
{
	gsl_integration_workspace * w = gsl_integration_workspace_alloc (limit);
	int res=gsl_integration_qagiu(f,a,epsabs,epsrel,limit, w, result, abserr);
	gsl_integration_workspace_free(w);
	return res;
}

int integral_QAGIU (double (*fun)(double, void*), void* param, double a, double epsabs, double epsrel, size_t limit, double *result, double * abserr)
{
    gsl_function f;
    f.params=param;
    f.function=fun;
    
	int res= integral_QAGIU (&f,a,epsabs,epsrel,limit,result,abserr);
    return res;
}
//+++++++++++++++++++++++++++++++++++++++++++++
//++ QAGIL adaptive integration on infinite intervals [-inf,b] ++ [7]
//+++++++++++++++++++++++++++++++++++++++++++++
/*

This function computes the integral of the function f over the semi-infinite interval (-\infty,b). The integral is mapped onto the semi-open interval (0,1] using the transformation x = b - (1-t)/t, and then integrated using the QAGS algorithm. 

*/

//++++++++++++++++++++++++++++++++++++++++++++

int integral_QAGIL (gsl_function * f, double b, double epsabs, double epsrel, size_t limit, double *result, double * abserr)
{
	gsl_integration_workspace * w = gsl_integration_workspace_alloc (limit);
	int res=gsl_integration_qagil(f,b,epsabs,epsrel,limit, w, result, abserr);
	gsl_integration_workspace_free(w);
	return res;
}

int integral_QAGIL (double (*fun)(double, void*), void* param, double b, double epsabs, double epsrel, size_t limit, double *result, double * abserr)
{
    gsl_function f;
    f.params=param;
    f.function=fun;
    
	int res= integral_QAGIL (&f,b,epsabs,epsrel,limit,result,abserr);
    return res;
}
//+++++++++++++++++++++++++++++++++++++++++++++
//++ QAGIL adaptive integration on infinite intervals [-inf,b] ++ [8]
//+++++++++++++++++++++++++++++++++++++++++++++
/*

This function computes the Cauchy principal value of the integral of f over (a,b), with a singularity at c,

I = \int_a^b dx f(x) / (x - c)

The adaptive bisection algorithm of QAG is used, with modifications to ensure that subdivisions do not occur at the singular point x = c. When a subinterval contains the point x = c or is close to it then a special 25-point modified Clenshaw-Curtis rule is used to control the singularity. Further away from the singularity the algorithm uses an ordinary 15-point Gauss-Kronrod integration rule. 
*/

//++++++++++++++++++++++++++++++++++++++++++++

int integral_QAGWC (gsl_function * f, double a, double b, double c, double epsabs, double epsrel, size_t limit, double *result, double * abserr)
{
	gsl_integration_workspace * w = gsl_integration_workspace_alloc (limit);
	int res=gsl_integration_qawc(f,a,b,c,epsabs,epsrel,limit, w, result, abserr);
	gsl_integration_workspace_free(w);
	return res;
}

int integral_QAGWC (double (*fun)(double, void*), void* param, double a, double b, double c, double epsabs, double epsrel, size_t limit, double *result, double * abserr)
{
    gsl_function f;
    f.params=param;
    f.function=fun;
    
	int res= integral_QAGWC (&f,a,b,c,epsabs,epsrel,limit,result,abserr);
    return res;
}
//+++++++++++++++++++++++++++++++++++++++++++++
//++ QAWO adaptive integration for oscillatory functions [sine] + [9]
//+++++++++++++++++++++++++++++++++++++++++++++
/*

The QAWO algorithm is designed for integrands with an oscillatory factor  \sin(\omega x). In order to work efficiently the algorithm requires a table of Chebyshev moments which must be pre-computed with calls to the functions below. 

This function allocates space for a gsl_integration_qawo_table struct and its associated workspace describing a sine weight function W(x) with the parameters (\omega, L),

          W(x) = sin(omega x)

The parameter L must be the length of the interval over which the function will be integrated L = b - a. The choice of sine or cosine is made with the parameter sine.
The gsl_integration_qawo_table is a table of the trigonometric coefficients required in the integration process. The parameter n determines the number of levels of coefficients that are computed. Each level corresponds to one bisection of the interval L, so that n levels are sufficient for subintervals down to the length L/2^n. The integration routine gsl_integration_qawo returns the error GSL_ETABLE if the number of levels is insufficient for the requested accuracy. 

This function uses an adaptive algorithm to compute the integral of f over (a,b) with the weight function \sin(\omega x) defined by the table wf,

              I = \int_a^b dx f(x) sin(omega x)

The results are extrapolated using the epsilon-algorithm to accelerate the convergence of the integral. The function returns the final approximation from the extrapolation, result, and an estimate of the absolute error, abserr. The subintervals and their results are stored in the memory provided by workspace. The maximum number of subintervals is given by limit, which may not exceed the allocated size of the workspace.

Those subintervals with “large” widths d where d\omega > 4 are computed using a 25-point Clenshaw-Curtis integration rule, which handles the oscillatory behavior. Subintervals with a “small” widths where d\omega < 4 are computed using a 15-point Gauss-Kronrod integration. 
*/

//++++++++++++++++++++++++++++++++++++++++++++

int integral_QAWO_sine (gsl_function * f, const double a, const double b, double omega, const double epsabs,const double epsrel, const  size_t limit, double *result, double * abserr, size_t n)
{
	double L=b-a;

	gsl_integration_workspace * w = gsl_integration_workspace_alloc (limit);

	gsl_integration_qawo_table * wf=gsl_integration_qawo_table_alloc (omega, L, GSL_INTEG_SINE, n);

	int res= gsl_integration_qawo (f,a,epsabs,epsrel,limit, w, wf, result, abserr);

	gsl_integration_qawo_table_free (wf);
	gsl_integration_workspace_free(w);
	return res;
}

int integral_QAWO_sine (double (*fun)(double, void*), void* param, double a, double b, double omega, const double epsabs,const double epsrel, const  size_t limit, double *result, double * abserr, size_t n)
{
    gsl_function f;
    f.params=param;
    f.function=fun;
    
	int res= integral_QAWO_sine (&f,a,b,omega,epsabs,epsrel,limit,result,abserr,n);
    return res;
}

//+++++++++++++++++++++++++++++++++++++++++++++++
//++ QAWO adaptive integration for oscillatory functions [cosine] + [11]
//+++++++++++++++++++++++++++++++++++++++++++++++
/*

The QAWO algorithm is designed for integrands with an oscillatory factor  \cos(\omega x). In order to work efficiently the algorithm requires a table of Chebyshev moments which must be pre-computed with calls to the functions below. 

This function allocates space for a gsl_integration_qawo_table struct and its associated workspace describing a sine weight function W(x) with the parameters (\omega, L),

          W(x) = cos(omega x)

The parameter L must be the length of the interval over which the function will be integrated L = b - a. The choice of sine or cosine is made with the parameter sine.
The gsl_integration_qawo_table is a table of the trigonometric coefficients required in the integration process. The parameter n determines the number of levels of coefficients that are computed. Each level corresponds to one bisection of the interval L, so that n levels are sufficient for subintervals down to the length L/2^n. The integration routine gsl_integration_qawo returns the error GSL_ETABLE if the number of levels is insufficient for the requested accuracy. 

This function uses an adaptive algorithm to compute the integral of f over (a,b) with the weight function \sin(\omega x) defined by the table wf,

              I = \int_a^b dx f(x) cos(omega x)

The results are extrapolated using the epsilon-algorithm to accelerate the convergence of the integral. The function returns the final approximation from the extrapolation, result, and an estimate of the absolute error, abserr. The subintervals and their results are stored in the memory provided by workspace. The maximum number of subintervals is given by limit, which may not exceed the allocated size of the workspace.

Those subintervals with “large” widths d where d\omega > 4 are computed using a 25-point Clenshaw-Curtis integration rule, which handles the oscillatory behavior. Subintervals with a “small” widths where d\omega < 4 are computed using a 15-point Gauss-Kronrod integration. 
*/

//++++++++++++++++++++++++++++++++++++++++++++

int integral_QAWO_cosine (gsl_function * f, const double a, const double b, double omega, const double epsabs,const double epsrel, const  size_t limit, double *result, double * abserr, size_t n)
{
	double L=b-a;

	gsl_integration_workspace * w = gsl_integration_workspace_alloc (limit);

	gsl_integration_qawo_table * wf=gsl_integration_qawo_table_alloc (omega, L, GSL_INTEG_COSINE, n);

	int res= gsl_integration_qawo (f,a,epsabs,epsrel,limit, w, wf, result, abserr);

	gsl_integration_qawo_table_free (wf);
	gsl_integration_workspace_free(w);
	return res;
}

int integral_QAWO_cosine (double (*fun)(double, void*), void* param, const double a, const double b, double omega, const double epsabs,const double epsrel, const  size_t limit, double *result, double * abserr, size_t n)
{
    gsl_function f;
    f.params=param;
    f.function=fun;
    
	int res= integral_QAWO_cosine (&f,a,b,omega,epsabs,epsrel,limit,result,abserr,n);
    return res;
}


//+++++++++++++++++++++++++++++++++++++++++++++
//++ QAWF adaptive integration for Fourier integrals [sine] ++++ [12]
//+++++++++++++++++++++++++++++++++++++++++++++
/*

This function attempts to compute a Fourier integral of the function f over the semi-infinite interval [a,+\infty).

              I = \int_a^{+\infty} dx f(x) sin(omega x)

The parameter \omega is taken from the table wf (the length L can take any value, since it is overridden by this function to a value appropriate for the fourier integration). The integral is computed using the QAWO algorithm over each of the subintervals,

              C_1 = [a, a + c]
              C_2 = [a + c, a + 2 c]
              ... = ...
              C_k = [a + (k-1) c, a + k c]

where c = (2 floor(|\omega|) + 1) \pi/|\omega|. The width c is chosen to cover an odd number of periods so that the contributions from the intervals alternate in sign and are monotonically decreasing when f is positive and monotonically decreasing. The sum of this sequence of contributions is accelerated using the epsilon-algorithm.

This function works to an overall absolute tolerance of abserr. The following strategy is used: on each interval C_k the algorithm tries to achieve the tolerance

              TOL_k = u_k abserr

where u_k = (1 - p)p^{k-1} and p = 9/10. The sum of the geometric series of contributions from each interval gives an overall tolerance of abserr.

If the integration of a subinterval leads to difficulties then the accuracy requirement for subsequent intervals is relaxed,

              TOL_k = u_k max(abserr, max_{i<k}{E_i})

where E_k is the estimated error on the interval C_k.

The subintervals and their results are stored in the memory provided by workspace. The maximum number of subintervals is given by limit, which may not exceed the allocated size of the workspace. The integration over each subinterval uses the memory provided by cycle_workspace as workspace for the QAWO algorithm. 
*/

//++++++++++++++++++++++++++++++++++++++++++++

int integral_QAWF_sine (gsl_function * f, const double a, double omega, const double epsabs, const  size_t limit, double *result, double * abserr, size_t n)
{
	gsl_integration_workspace * w = gsl_integration_workspace_alloc (limit);
	gsl_integration_workspace * w_cycle = gsl_integration_workspace_alloc (limit);

	gsl_integration_qawo_table * wf=gsl_integration_qawo_table_alloc (omega, 1.0, GSL_INTEG_SINE, n);

	int res= gsl_integration_qawf (f,a,epsabs,limit, w,w_cycle, wf, result, abserr);

	gsl_integration_qawo_table_free (wf);
	gsl_integration_workspace_free(w);
	gsl_integration_workspace_free(w_cycle);
	return res;
}

int integral_QAWF_sine (double (*fun)(double, void*), void* param, const double a, double omega, const double epsabs, const  size_t limit, double *result, double * abserr, size_t n)
{
    gsl_function f;
    f.params=param;
    f.function=fun;
    
	int res= integral_QAWF_sine (&f,a,omega,epsabs,limit,result,abserr,n);
    return res;
}

//+++++++++++++++++++++++++++++++++++++++++++++
//++ QAWF adaptive integration for Fourier integrals [cosine] ++ [13]
//+++++++++++++++++++++++++++++++++++++++++++++
/*

This function attempts to compute a Fourier integral of the function f over the semi-infinite interval [a,+\infty).

              I = \int_a^{+\infty} dx f(x) cos(omega x)

The parameter \omega is taken from the table wf (the length L can take any value, since it is overridden by this function to a value appropriate for the fourier integration). The integral is computed using the QAWO algorithm over each of the subintervals,

              C_1 = [a, a + c]
              C_2 = [a + c, a + 2 c]
              ... = ...
              C_k = [a + (k-1) c, a + k c]

where c = (2 floor(|\omega|) + 1) \pi/|\omega|. The width c is chosen to cover an odd number of periods so that the contributions from the intervals alternate in sign and are monotonically decreasing when f is positive and monotonically decreasing. The sum of this sequence of contributions is accelerated using the epsilon-algorithm.

This function works to an overall absolute tolerance of abserr. The following strategy is used: on each interval C_k the algorithm tries to achieve the tolerance

              TOL_k = u_k abserr

where u_k = (1 - p)p^{k-1} and p = 9/10. The sum of the geometric series of contributions from each interval gives an overall tolerance of abserr.

If the integration of a subinterval leads to difficulties then the accuracy requirement for subsequent intervals is relaxed,

              TOL_k = u_k max(abserr, max_{i<k}{E_i})

where E_k is the estimated error on the interval C_k.

The subintervals and their results are stored in the memory provided by workspace. The maximum number of subintervals is given by limit, which may not exceed the allocated size of the workspace. The integration over each subinterval uses the memory provided by cycle_workspace as workspace for the QAWO algorithm. 
*/

//++++++++++++++++++++++++++++++++++++++++++++

int integral_QAWF_cosine (gsl_function * f, const double a, double omega, const double epsabs, const  size_t limit, double *result, double * abserr, size_t n)
{
	gsl_integration_workspace * w = gsl_integration_workspace_alloc (limit);
	gsl_integration_workspace * w_cycle = gsl_integration_workspace_alloc (limit);

	gsl_integration_qawo_table * wf=gsl_integration_qawo_table_alloc (omega, 1.0, GSL_INTEG_COSINE, n);

	int res= gsl_integration_qawf (f,a,epsabs,limit, w,w_cycle, wf, result, abserr);

	gsl_integration_qawo_table_free (wf);
	gsl_integration_workspace_free(w);
	gsl_integration_workspace_free(w_cycle);
	return res;
}

int integral_QAWF_cosine (double (*fun)(double, void*), void* param, const double a, double omega, const double epsabs, const  size_t limit, double *result, double * abserr, size_t n)
{
    gsl_function f;
    f.params=param;
    f.function=fun;
    
	int res= integral_QAWF_cosine (&f,a,omega,epsabs,limit,result,abserr,n);
    return res;
}

/*  GSL >=1.14
//+++++++++++++++++++++++++++++++++++++++++++++
//++ Gauss-Legendre integration ++++++++++++++++++++++ [14]
//+++++++++++++++++++++++++++++++++++++++++++++


The fixed-order Gauss-Legendre integration routines are provided for fast integration of smooth functions with known polynomial order. The n-point Gauss-Legendre rule is exact for polynomials of order 2*n-1 or less. For example, these rules are useful when integrating basis functions to form mass matrices for the Galerkin method. Unlike other numerical integration routines within the library, these routines do not accept absolute or relative error bounds. 



//++++++++++++++++++++++++++++++++++++++++++++

int integral_GLfixed (const gsl_function * f, const double a, double b, double *result, size_t n)
{
	gsl_integration_glfixed_table * t;
	t=gsl_integration_gl_table_alloc (n);
	result= gsl_integration_gl(f,a,b, t);
	gsl_integration_gl_table_alloc (t);
	return res;
}
*/

#endif