/* With the functions in this file  1-dimensional definite integrals 
   for a complex function can be calculated using  a Gauss-Legendre algorithm. */


#include <int_complex.h>


#define CFUNC(x) ((*cfunc)(x))  /* func is a pointer to the function 
                                 to be integrated */
#define CFUNC2(x) ((*cfunc2)(x))  /* func2 is a pointer to the function 
                                 to be integrated */
#define CFUNC3(x) ((*cfunc3)(x))  /* func3 is a pointer to the function 
                                 to be integrated */
#define CFUNC4(x) ((*cfunc4)(x))  /* func4 is a pointer to the function 
                                 to be integrated */




/* The function cint_gauss() calculates a one-dimensional definite integral 
   for a complex function using a gaussian quadrature formula 
   (in the present case we use the GAUSS-LEGENDRE quadrature formula.     
   - a and b are the lower and upper integration limits respectively,
   - n stands for the n-point quadrature formula (all polynomials
     up to degree (2 * n - 1) are integrated exactly),
   - x[1...n] and w[1...n] contain the abscissa and weights of the
     n-point Gauss-Legendre quadrature formula and are calculated
     by the function gauss_leg(). */

dcomplex cint_gauss(dcomplex (*cfunc)(double), double a, double b, int n)
{
     int j;
     double re_sum = 0., im_sum = 0.;
     dcomplex function;
     static double *x, *w;  
     
     x = dvector(1, n);
     w = dvector(1, n);
     
     gauss_leg(a, b, x, w, n);
     
     for (j = 1; j <= n; j++)
     { 
       function = CFUNC(x[j]);
       re_sum += w[j] * function.re;
       im_sum += w[j] * function.im;
     }  
     
     free_dvector(x, 1, n);
     free_dvector(w, 1, n);
     
     return Complex(re_sum, im_sum); 
} 


dcomplex cint_gauss2(dcomplex (*cfunc2)(double), double a, double b, int n)
{
     int j;
     double re_sum = 0., im_sum = 0.;
     dcomplex function;
     static double *x, *w;  
     
     x = dvector(1, n);
     w = dvector(1, n);
     
     gauss_leg(a, b, x, w, n);
     
     for (j = 1; j <= n; j++)
     { 
       function = CFUNC2(x[j]);
       re_sum += w[j] * function.re;
       im_sum += w[j] * function.im;
     }  
     
     free_dvector(x, 1, n);
     free_dvector(w, 1, n);
     
     return Complex(re_sum, im_sum); 
} 


dcomplex cint_gauss3(dcomplex (*cfunc3)(double), double a, double b, int n)
{
     int j;
     double re_sum = 0., im_sum = 0.;
     dcomplex function;
     static double *x, *w;  
     
     x = dvector(1, n);
     w = dvector(1, n);
     
     gauss_leg(a, b, x, w, n);
     
     for (j = 1; j <= n; j++)
     { 
       function = CFUNC3(x[j]);
       re_sum += w[j] * function.re;
       im_sum += w[j] * function.im;
     }  
     
     free_dvector(x, 1, n);
     free_dvector(w, 1, n);
     
     return Complex(re_sum, im_sum); 
} 


dcomplex cint_gauss4(dcomplex (*cfunc4)(double), double a, double b, int n)
{
     int j;
     double re_sum = 0., im_sum = 0.;
     dcomplex function;
     static double *x, *w;  
     
     x = dvector(1, n);
     w = dvector(1, n);
     
     gauss_leg(a, b, x, w, n);
     
     for (j = 1; j <= n; j++)
     { 
       function = CFUNC4(x[j]);
       re_sum += w[j] * function.re;
       im_sum += w[j] * function.im;
     }  
     
     free_dvector(x, 1, n);
     free_dvector(w, 1, n);
     
     return Complex(re_sum, im_sum); 
}