/* gauleg.c     P145 Numerical Recipes in Fortran               */
/* compute x(i) and w(i)  i=1,n  Legendre ordinates and weights */
/* on interval  -1.0 to 1.0 (length is 2.0)                     */
/* use ordinates and weights for Gauss Legendre integration     */

#include <stdio.h>
#include <math.h>
#undef  abs
#define abs(x) ((x)<0.0?(-(x)):(x))

void gaulegf(double x1, double x2, double x[], double w[], int n)
{
  int i, j, m;
  double eps = 3.0E-14;
  double p1, p2, p3, pp, xl, xm, z, z1;
  
  m = (n+1)/2;
  xm = 0.5*(x2+x1);
  xl = 0.5*(x2-x1);
  for(i=1; i<=m; i++)
  {
    z = cos(3.141592654*((double)i-0.25)/((double)n+0.5));
    while(1)
    {
      p1 = 1.0;
      p2 = 0.0;
      for(j=1; j<=n; j++)
      {
        p3 = p2;
        p2 = p1;
        p1 = ((2.0*(double)j-1.0)*z*p2-((double)j-1.0)*p3)/
              (double)j;
      }
      pp = (double)n*(z*p1-p2)/(z*z-1.0);
      z1 = z;
      z = z1 - p1/pp;
      if(abs(z-z1) <= eps) break;
    }
    x[i] = xm - xl*z;
    x[n+1-i] = xm + xl*z;
    w[i] = 2.0*xl/((1.0-z*z)*pp*pp);
    w[n+1-i] = w[i];
  }
} /* end gaulegf */


int main(int argc, char *argv[])
{
  int i, j;
  double sum, a, b;
  double x[100], w[100];

  printf("test gauleg.c on interval -1.0 to 1.0  ordinates, weights\n");
  for(i=1; i<=15; i++)
  {
    gaulegf(-1.0, 1.0, x, w, i);
    sum = 0.0;
    for(j=1; j<=i; j++)
    {
      printf("x[%d]=%21.13E, w[%d]=%21.13E \n", j, x[j], j, w[j]);
      sum = sum + w[j];
    }
    printf("              integral(1.0, -1.0..1.0)=%21.13E", sum);
    printf("\n");
  }

  a = 0.5;
  b = 1.0;
  for(i=2; i<=10; i++)
  {
    gaulegf(a, b, x, w, i);
    sum = 0.0;
    for(j=1; j<=i; j++)
    {
      sum = sum + w[j]*sin(x[j]);
    }
    printf("integral (0.5,1.0) sin(x) dx = %21.13E\n", sum);
  }
  printf("-cos(1.0)+cos(0.5) = %21.13E\n", (-cos(1.0)+cos(0.5)));
  printf("Maple says 3.372802560E-001 \n");
  printf("\n");

  a = 0.5;
  b = 5.0;
  for(i=2; i<=10; i++)
  {
    gaulegf(a, b, x, w, i);
    sum = 0.0;
    for(j=1; j<=i; j++)
    {
      sum = sum + w[j]*exp(x[j]);
    }
    printf("integral (0.5,5.0) exp(x) dx = %21.13E\n", sum);
  }
  printf("exp(5.0)-exp(0.5) = %21.13E\n", (exp(5.0)-exp(0.5)));
  printf("Maple says 1.467644378E+002 \n");
  printf("\n");

  a = 0.5;
  b = 5.0;
  for(i=2; i<=30; i++)
  {
    gaulegf(0.5, 5.0, x, w, i);
    sum = 0.0;
    for(j=1; j<=i; j++)
    {
      sum = sum + w[j]*((pow(pow(x[j],x[j]),x[j])*(x[j]*(log(x[j])+1.0)+x[j]*log(x[j]))));
    }
    printf("integral (0.5,5.0) mess(x) dx = %21.13E\n", sum);
  }
  printf("((5.0**5.0)**5.0)-(0.5**0.5)**0.5 = %21.13E\n",
         (pow(pow(5.0,5.0),5.0)-pow(pow(0.5,0.5),0.5)));
  printf("Maple says 2.980232239E+017 \n");
  printf("\nDone.\n");
} /* end main of gauleg.c */