/* test_gaulegf1D.c  test  gaulegf.c integration accuracy */
#include "gaulegf.h"
#include "laphi.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/* test functions and their indefinate integrals */

static double f11(double x)
{
  return x+2.0;
}

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

static double f21(double x)
{
  return x*x+2.0*x+3.0;
}

static double int_f21(double x)
{
  return x*x*x/3.0+2.0*x*x/2.0+3.0*x;
}

static double f31(double x)
{
  return x*x*x+2.0*x*x+3.0*x+4.0;
}

static double int_f31(double x)
{
  return x*x*x*x/4.0+2.0*x*x*x/3.0+3.0*x*x/2.0+4.0*x;
}

static double f41(double x)
{
  return x*x*x*x+2.0*x*x*x+3.0*x*x+4.0*x+5.0;
}

static double int_f41(double x)
{
  return x*x*x*x*x/5.0+2.0*x*x*x*x/4.0+3.0*x*x*x/3.0+4.0*x*x/2.0+5.0*x;
}

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

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


int main(int argc, char * argv[])
{
  int n;
  double xmin, xmax;
  double x[100];
  double w[100];
  double area, exact;
  double dxg, xg[100];
  int i, j, m, mj;

  printf("test_gaulegf1D.c running \n");
  for(n=1; n<6; n++)
  {
    xmin = -1.0;
    xmax =  2.0;
    printf("calling gaulegf(%f, %f, x, w, n) \n", xmin, xmax, n);
    printf("n=%d, equivalent h=%f \n", n, (xmax-xmin)/n);
    gaulegf(xmin, xmax, x, w, n);
    for(i=1; i<=n; i++)
      printf("x[%d]=%f, w[%d]=%f \n", i, x[i], i, w[i]);
    printf("\n");
  }

  printf("f11(x)=x+2.0 \n");
  for(n=1; n<6; n++)
  {
    xmin = -1.0;
    xmax =  2.0;
    gaulegf(xmin, xmax, x, w, n);
    area = 0.0;
    for(i=1; i<=n; i++)
    {
      area += w[i]*f11(x[i]);
    }
    exact=int_f11(xmax)-int_f11(xmin);
    printf("n=%d, area=%f, exact=%f, error=%g \n", n, area, exact, area-exact);
  }
  printf("\n");

  printf("f21(x)=x*x+2.0*x+3.0 \n");
  for(n=1; n<6; n++)
  {
    xmin = -1.0;
    xmax =  2.0;
    gaulegf(xmin, xmax, x, w, n);
    area = 0.0;
    for(i=1; i<=n; i++)
    {
      area += w[i]*f21(x[i]);
    }
    exact=int_f21(xmax)-int_f21(xmin);
    printf("n=%d, area=%f, exact=%f, error=%g \n", n, area, exact, area-exact);
  }
  printf("\n");


  printf("f31(x)=x*x*x+2.0*x*x+3.0*x+4.0 \n");
  for(n=1; n<6; n++)
  {
    xmin = -1.0;
    xmax =  2.0;
    gaulegf(xmin, xmax, x, w, n);
    area = 0.0;
    for(i=1; i<=n; i++)
    {
      area += w[i]*f31(x[i]);
    }
    exact=int_f31(xmax)-int_f31(xmin);
    printf("n=%d, area=%f, exact=%f, error=%g \n", n, area, exact, area-exact);
  }
  printf("\n");

  printf("f41(x)=x*x*x*x+2.0*x*x*x+3.0*x*x+4.0*x+5.0 \n");
  for(n=1; n<6; n++)
  {
    xmin = -1.0;
    xmax =  2.0;
    gaulegf(xmin, xmax, x, w, n);
    area = 0.0;
    for(i=1; i<=n; i++)
    {
      area += w[i]*f41(x[i]);
    }
    exact=int_f41(xmax)-int_f41(xmin);
    printf("n=%d, area=%f, exact=%f, error=%g \n", n, area, exact, area-exact);
  }
  printf("\n");

  printf("fe1(x)=exp(x) \n");
  for(n=1; n<6; n++)
  {
    xmin = -1.0;
    xmax =  2.0;
    gaulegf(xmin, xmax, x, w, n);
    area = 0.0;
    for(i=1; i<=n; i++)
    {
      area += w[i]*fe1(x[i]);
    }
    exact=int_fe1(xmax)-int_fe1(xmin);
    printf("n=%d, area=%f, exact=%f, error=%g \n", n, area, exact, area-exact);
  }
  printf("\n");

  for(m=3; m<6; m=m+2)
  {
    mj=(m+1)/2;
    printf("phi(x,%d,%d,xg) \n", mj, m);
    xmin = -1.0;
    xmax =  2.0;
    dxg = (xmax-xmin)/m;
    if(n==1) printf("phi order=%d, equivh=%d \n", m, dxg);
    for(j=0; j<=m; j++)
    {
      xg[j] = xmin+j*dxg;
      if(n==1) printf("j=%d, xg[j]=%f \n",j, xg[j]);
    }
    for(n=1; n<8; n++)
    {
      gaulegf(xmin, xmax, x, w, n);
      area = 0.0;
      for(i=1; i<=n; i++)
      {
        area += w[i]*phi(x[i],mj,m,xg);
      }
      exact=1.125;
      printf("n=%d, area=%f, exact=%f, error=%g \n",
              n, area, exact, area-exact);
    }
    printf("\n");

    printf("phi(x,%d,%d,xg)*f11=x+2.0 \n", mj, m);
    for(n=1; n<8; n++)
    {
      gaulegf(xmin, xmax, x, w, n);
      area = 0.0;
      for(i=1; i<=n; i++)
      {
        area += w[i]*phi(x[i],mj,m,xg)*f11(x[i]);
      }
      exact=3.825;
      printf("n=%d, area=%f, exact=%f, error=%g \n",
              n, area, exact, area-exact);
    }
    printf("\n");

    printf("phi(x,%d,%d,xg)*f21=x*x+2.0*x+3.0 \n", mj, m);
    for(n=1; n<8; n++)
    {
      gaulegf(xmin, xmax, x, w, n);
      area = 0.0;
      for(i=1; i<=n; i++)
      {
        area += w[i]*phi(x[i],mj,m,xg)*f21(x[i]);
      }
      exact=8.325;
      printf("n=%d, area=%f, exact=%f, error=%g \n",
              n, area, exact, area-exact);
    }
    printf("\n");

    printf("phi(x,%d,%d,xg)*f31=x*x*x+2.0*x*x+3.0*x+4.0 \n", mj, m);
    for(n=1; n<8; n++)
    {
      gaulegf(xmin, xmax, x, w, n);
      area = 0.0;
      for(i=1; i<=n; i++)
      {
        area += w[i]*phi(x[i],mj,m,xg)*f31(x[i]);
      }
      exact=15.460714;
      printf("n=%d, area=%f, exact=%f, error=%g \n",
              n, area, exact, area-exact);
    }
    printf("\n");

    printf("phi(x,%d,%d,xg)*f41=x*x*x*x+2.0*x*x*x+3.0*x*x+4.0*x+5.0 \n",
           mj, m);
    for(n=1; n<8; n++)
    {
      gaulegf(xmin, xmax, x, w, n);
      area = 0.0;
      for(i=1; i<=n; i++)
      {
        area += w[i]*phi(x[i],mj,m,xg)*f41(x[i]);
      }
      exact=26.373214;
      printf("n=%d, area=%f, exact=%f, error=%g \n",
              n, area, exact, area-exact);
    }
    printf("\n");

    printf("phi(x,%d,%d,xg)*fe1=exp(x) \n", mj, m);
    for(n=1; n<8; n++)
    {
      gaulegf(xmin, xmax, x, w, n);
      area = 0.0;
      for(i=1; i<=n; i++)
      {
        area += w[i]*phi(x[i],mj,m,xg)*fe1(x[i]);
      }
      exact=4.262081;
      printf("n=%d, area=%f, exact=%f, error=%g \n",
              n, area, exact, area-exact);
    }
    printf("\n");
    if(m==3) printf("end of 'exact' being correct \n\n"); 
  }


  m=5;
  mj=0;
  printf("phi(x,%d,%d,xg)*fe1=exp(x) \n", mj, m);
  for(n=1; n<8; n++)
  {
    gaulegf(xmin, xmax, x, w, n);
    area = 0.0;
    for(i=1; i<=n; i++)
    {
      area += w[i]*phi(x[i],mj,m,xg)*fe1(x[i]);
    }
    exact=0.0;
    printf("n=%d, area=%f, exact=%f, error=%g \n",
            n, area, exact, area-exact);
  }
  printf("\n");

  m=5;
  mj=5;
  printf("phi(x,%d,%d,xg)*fe1=exp(x) \n", mj, m);
  for(n=1; n<8; n++)
  {
    gaulegf(xmin, xmax, x, w, n);
    area = 0.0;
    for(i=1; i<=n; i++)
    {
      area += w[i]*phi(x[i],mj,m,xg)*fe1(x[i]);
    }
    exact=0.0;
    printf("n=%d, area=%f, exact=%f, error=%g \n",
            n, area, exact, area-exact);
  }
  printf("\n");

  m=5;
  mj=0;
  printf("phip(x,%d,%d,xg)*fe1=exp(x) \n", mj, m);
  for(n=1; n<8; n++)
  {
    gaulegf(xmin, xmax, x, w, n);
    area = 0.0;
    for(i=1; i<=n; i++)
    {
      area += w[i]*phip(x[i],mj,m,xg)*fe1(x[i]);
    }
    exact=0.0;
    printf("n=%d, area=%f, exact=%f, error=%g \n",
            n, area, exact, area-exact);
  }
  printf("\n");

  m=5;
  mj=5;
  printf("phip(x,%d,%d,xg)*fe1=exp(x) \n", mj, m);
  for(n=1; n<8; n++)
  {
    gaulegf(xmin, xmax, x, w, n);
    area = 0.0;
    for(i=1; i<=n; i++)
    {
      area += w[i]*phip(x[i],mj,m,xg)*fe1(x[i]);
    }
    exact=0.0;
    printf("n=%d, area=%f, exact=%f, error=%g \n",
            n, area, exact, area-exact);
  }
  printf("\n");

  return 0;
} /* end test_gaulegf.c */