/* line_int.c  integrate a function over a line                 */
/*             this is numerical quadrature                     */
/*             basically compute an area under line function    */
/*             compare methods for accuracy                     */

/* first trapezoid rule, varying step size, int f(x) dx         */
/* then using Gauss-Legendre coordinates and weights            */

/* remember np points or nodes means np-1 intervals or segments */
/* np will be the number of times the function is computed      */

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

static void gaulegf(double x1, double x2, double x[], double w[], int n);

static double f1(double x)
{
  return x+2.0; /* linear */
}

static double f2(double x)
{
  return x*x + 2.0*x + 3.0; /* quadratic, second order */
}

static double f3(double x)
{
  return x*x*x + 2.0*x*x + 3.0*x + 4.0; /* cubic, third order */
}

static double fe(double x)
{
  return exp(x);
}

int main(int argc, char *argv[])
{
  /* simple trapezoidal integration from x=1 to x=3*/
  double area, sum;
  double x, h;
  double xmin = 1.0;
  double xmax = 3.0;
  double f1exact = 8.0;
  double f2exact = 68.0/3.0;
  double f3exact = 57.0+1.0/3.0;
  double feexact;
  double xx[200], ww[200]; /* for Gauss-Legendre */
  double val;
  int p, np;
  int i, j;

  feexact = exp(xmax)-exp(xmin);
  printf("line_int.c integrate to find area \n");
  printf("Numerical quadrature, x=%g to x=%g \n", xmin, xmax);
  printf("\n");

  printf("Using trapezoidal method \n");
  np = 4;
  h = (xmax-xmin)/(double)(np-1);
  printf("f1(x)=x+2.0  np=%d,  h=%g \n", np, h);
  area = (f1(xmin) + f1(xmax))/2.0;
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*h;
    area = area + f1(x);
  }
  area = (xmax-xmin)*area/(double)(np-1);
  printf("%d function evaluations \n", np);
  printf("computed area=%g, exact area=%g, error=%g \n",
         area, f1exact, area-f1exact);
  printf("\n");

  np = 4;
  h = (xmax-xmin)/(double)(np-1);
  printf("f2(x)=x*x+2.0*x+3.0  np=%d,  h=%g \n", np, h);
  area = (f2(xmin) + f2(xmax))/2.0;
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*h;
    area = area + f2(x);
  }
  area = (xmax-xmin)*area/(double)(np-1);
  printf("%d function evaluations \n", np);
  printf("computed area=%g, exact area=%g, error=%g \n",
         area, f2exact, area-f2exact);

  np = 8;
  h = (xmax-xmin)/(double)(np-1);
  printf("f2(x)=x*x+2.0*x+3.0  np=%d,  h=%g \n", np, h);
  area = (f2(xmin) + f2(xmax))/2.0;
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*h;
    area = area + f2(x);
  }
  area = (xmax-xmin)*area/(double)(np-1);
  printf("%d function evaluations \n", np);
  printf("computed area=%g, exact area=%g, error=%g \n",
         area, f2exact, area-f2exact);

  np = 16;
  h = (xmax-xmin)/(double)(np-1);
  printf("f2(x)=x*x+2.0*x+3.0  np=%d,  h=%g \n", np, h);
  area = (f2(xmin) + f2(xmax))/2.0;
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*h;
    area = area + f2(x);
  }
  area = (xmax-xmin)*area/(double)(np-1);
  printf("%d function evaluations \n", np);
  printf("computed area=%g, exact area=%g, error=%g \n",
         area, f2exact, area-f2exact);
  printf("\n");

  np = 4;
  h = (xmax-xmin)/(double)(np-1);
  printf("f3(x)=x*x*x+2.0*x*x+3.0*x+4.0  np=%d,  h=%g \n", np, h);
  area = (f3(xmin) + f3(xmax))/2.0;
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*h;
    area = area + f3(x);
  }
  area = (xmax-xmin)*area/(double)np;
  printf("%d function evaluations \n", np);
  printf("computed area=%g, exact area=%g, error=%g \n",
         area, f3exact, area-f3exact);

  np = 16;
  h = (xmax-xmin)/(double)(np-1);
  printf("f3(x)=x*x*x+2.0*x*x+3.0*x+4.0  np=%d,  h=%g \n", np, h);
  area = (f3(xmin) + f3(xmax))/2.0;
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*h;
    area = area + f3(x);
  }
  area = (xmax-xmin)*area/(double)(np-1);
  printf("%d function evaluations \n", np);
  printf("computed area=%g, exact area=%g, error=%g \n",
         area, f3exact, area-f3exact);

  np = 64;
  h = (xmax-xmin)/(double)(np-1);
  printf("f3(x)=x*x*x+2.0*x*x+3.0*x+4.0  np=%d,  h=%g \n", np, h);
  area = (f3(xmin) + f3(xmax))/2.0;
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*h;
    area = area + f3(x);
  }
  area = (xmax-xmin)*area/(double)(np-1);
  printf("%d function evaluations \n", np);
  printf("computed area=%g, exact area=%g, error=%g \n",
         area, f3exact, area-f3exact);
  printf("\n");

  np = 4;
  h = (xmax-xmin)/(double)(np-1);
  printf("fe(x)=exp(x)  np=%d,  h=%g \n", np, h);
  area = (fe(xmin) + fe(xmax))/2.0;
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*h;
    area = area + f3(x);
  }
  area = (xmax-xmin)*area/(double)np;
  printf("%d function evaluations \n", np);
  printf("computed area=%g, exact area=%g, error=%g \n",
         area, feexact, area-feexact);

  np = 16;
  h = (xmax-xmin)/(double)(np-1);
  printf("fe(x)=exp(x)  np=%d,  h=%g \n", np, h);
  area = (fe(xmin) + fe(xmax))/2.0;
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*h;
    area = area + fe(x);
  }
  area = (xmax-xmin)*area/(double)(np-1);
  printf("%d function evaluations \n", np);
  printf("computed area=%g, exact area=%g, error=%g \n",
         area, feexact, area-feexact);

  np = 64;
  h = (xmax-xmin)/(double)(np-1);
  printf("fe(x)=exp(x)  np=%d,  h=%g \n", np, h);
  area = (fe(xmin) + fe(xmax))/2.0;
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*h;
    area = area + fe(x);
  }
  area = (xmax-xmin)*area/(double)(np-1);
  printf("%d function evaluations \n", np);
  printf("computed area=%g, exact area=%g, error=%g \n",
         area, feexact, area-feexact);
  printf("\n");

  printf("Using Gauss Legendre method \n");

  np = 4;
  printf("f1(x)=x+2.0  np=%d \n", np);
  gaulegf(xmin, xmax, xx, ww, np);
  area = 0.0;
  for(p=1; p<=np; p++)
  {
    area = area + ww[p]*f1(xx[p]);
  }
  printf("%d function evaluations \n", np);
  printf("computed area=%g, exact area=%g, error=%g \n",
         area, f1exact, area-f1exact);
  printf("\n");

  np = 4;
  printf("f2(x)=x*x+2.0*x+3.0  np=%d \n", np);
  gaulegf(xmin, xmax, xx, ww, np);
  area = 0.0;
  for(p=1; p<=np; p++)
  {
    area = area + ww[p]*f2(xx[p]);
  }
  printf("%d function evaluations \n", np);
  printf("computed area=%g, exact area=%g, error=%g \n",
         area, f2exact, area-f2exact);
  printf("\n");

  np = 4;
  printf("f3(x)=x*x*x+2.0*x*x+3.0*x+4  np=%d \n", np);
  gaulegf(xmin, xmax, xx, ww, np);
  area = 0.0;
  for(p=1; p<=np; p++)
  {
    area = area + ww[p]*f3(xx[p]);
  }
  printf("%d function evaluations \n", np);
  printf("computed area=%g, exact area=%g, error=%g \n",
         area, f3exact, area-f3exact);
  printf("\n");

  np = 4;
  printf("fe(x)=exp(x)  np=%d \n", np);
  gaulegf(xmin, xmax, xx, ww, np);
  area = 0.0;
  for(p=1; p<=np; p++)
  {
    area = area + ww[p]*fe(xx[p]);
  }
  printf("%d function evaluations \n", np);
  printf("computed area=%g, exact area=%g, error=%g \n",
         area, feexact, area-feexact);
  printf("\n");

  np = 8;
  printf("fe(x)=exp(x)  np=%d \n", np);
  gaulegf(xmin, xmax, xx, ww, np);
  area = 0.0;
  for(p=1; p<=np; p++)
  {
    area = area + ww[p]*fe(xx[p]);
  }
  printf("%d function evaluations \n", np);
  printf("computed area=%g, exact area=%g, error=%g \n",
         area, feexact, area-feexact);
  printf("\n");


  return 0;
} /* end main  of line_int.c */

static 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 */
/* end line_int.c */