/* test_triquad_int.c  test  triquad_int.c integration accuracy */
#include "triquad_int.h"
#include "phi_tri.h"
#include "phi_tri_cm.h"
#include "simeq.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#undef abs
#define abs(a) ((a)<0.0?(-(a)):(a))

/* test functions */

static double f12(double x, double y)
{
  return x+2.0*y+3.0;
}

static double f22(double x, double y)
{
  return x*x+2.0*y*y+3.0*x*y+4.0*x+5.0*y+6.0;
}

static double f32(double x, double y)
{
  return x*x*x+2.0*y*y*y+3.0*x*x*y+4.0*x*y*y+
         5.0*x*x+6.0*y*y+7.0*x*y+8.0*x+9.0*y+10.0;
}

static double f42(double x, double y)
{
  return x*x*x*x+2.0*y*y*y*y+3.0*x*x*x*y+4.0*x*x*y*y+
         5.0*x*x*x*y+6.0*x*x*x+7.0*y*y*y+8.0*x*x*y+
         9.0*x*y*y+10.0*x*x+11.0*y*y+12.0*x*y+13.0*x+
         14.0*y+15.0;
}

static double fe2(double x, double y)
{
  return exp(x+y);
}

static double triarea(double x1, double y1,
                      double x2, double y2,
                      double x3, double y3)
{
  double a;

  a = x1*y2 - x2*y1 + x2*y3 - x3*y2 + x3*y1 - x1*y3;
  return abs(a)/2.0;
}


int main(int argc, char * argv[])
{
  int n, nn, m;
  double x1, y1, x2, y2, x3, y3;
  double tarea;
  double xx[100];
  double ww[100];
  double yy[100];
  double cm1[100], cm2[100], cm3[100];
  double area, exact; /* area is actually volume here */
  int i, j, debug;

  printf("test_triquad_int.c running \n");
  x1 = 1.0;
  y1 = 1.0;
  x2 = 3.0;
  y2 = 2.5;
  x3 = 2.0;
  y3 = 0.5;
  /* for phi22 calls */
  debug = 3;
  printf("phi_tri_cm22 running \n");
  phi_tri_cm22(x1, y1, x2, y2, x3, y3, debug, cm1);
  phi_tri_cm22(x2, y2, x3, y3, x1, y1, debug, cm2);
  phi_tri_cm22(x3, y3, x1, y1, x2, y2, debug, cm3);
  for(i=0; i<6; i++) printf("cm3[%d}=%f \n", i, cm3[i]);

  for(n=1; n<=6; n++)
  {
    printf("calling triquad_int(%d, %f,%f  %f,%f  %f,%f  xx, yy, ww) \n",
	   n, x1, y1, x2, y2, x3, y3);
    if(n==2)
    {
      nn =triquad_int(n, x1, y1, x2, y2, x3, y3, xx, yy, ww);
      for(i=0; i<nn; i++)
        printf("123 i=%d, xx=%f, yy=%f, ww=%f \n", i, xx[i], yy[i], ww[i]);
      nn =triquad_int(n, x2, y2, x3, y3, x1, y1, xx, yy, ww);
      for(i=0; i<nn; i++)
        printf("231 i=%d, xx=%f, yy=%f, ww=%f \n", i, xx[i], yy[i], ww[i]);
      nn =triquad_int(n, x3, y3, x1, y1, x2, y2, xx, yy, ww);
      for(i=0; i<nn; i++)
        printf("312 i=%d, xx=%f, yy=%f, ww=%f \n", i, xx[i], yy[i], ww[i]);
    }
    else
    {
      nn =triquad_int(n, x1, y1, x2, y2, x3, y3, xx, yy, ww);
      for(i=0; i<nn; i++)
        printf("i=%d, xx=%f, yy=%f, ww=%f \n", i, xx[i], yy[i], ww[i]);
    }
    printf("\n");
  }

  printf("f12(x,y)=x+2.0*y+3.0 \n");
  for(n=1; n<=6; n++)
  {
    nn = triquad_int(n, x1, y1, x2, y2, x3, y3, xx, yy, ww);
    area = 0.0;
    for(i=0; i<nn; i++)
    {
      area += ww[i]*f12(xx[i],yy[i]);
    }
    printf("n=%d, nn=%d, area=%f \n", n, nn, area);
  }
  printf("\n");

  printf("f22(x,y)=x*x+2.0*y*y+3.0*x*y+4.0*x+5.0*y+6.0; \n");
  for(n=1; n<=6; n++)
  {
    nn = triquad_int(n, x1, y1, x2, y2, x3, y3, xx, yy, ww);
    area = 0.0;
    for(i=0; i<nn; i++)
    {
      area += ww[i]*f22(xx[i],yy[i]);
    }
    printf("n=%d, nn=%d, area=%f \n", n, nn, area);
  }
  printf("\n");


  printf("f32(x,y)=x*x*x+2.0*y*y*y+3.0*x*x*y+4.0*x*y*y+ \n");
  printf("         5.0*x*x+6.0*y*y+7.0*x*y+8.0*x+9.0*y+10.0 \n");
  for(n=1; n<=6; n++)
  {
    nn = triquad_int(n, x1, y1, x2, y2, x3, y3, xx, yy, ww);
    area = 0.0;
    for(i=0; i<nn; i++)
    {
      area += ww[i]*f32(xx[i],yy[i]);
    }
    printf("n=%d, nn=%d, area=%f \n", n, nn, area);
  }
  printf("\n");


  printf("f42(x,y)=x*x*x*x+2.0*y*y*y*y+3.0*x*x*x*y+4.0*x*x*y*y+ \n");
  printf("         5.0*x*x*x*y+6.0*x*x*x+7.0*y*y*y+8.0*x*x*y+ \n");
  printf("         9.0*x*y*y+10.0*x*x+11.0*y*y+12.0*x*y+13.0*x+ \n");
  printf("         14.0*y+15.0 \n");
  for(n=1; n<=6; n++)
  {
    nn = triquad_int(n, x1, y1, x2, y2, x3, y3, xx, yy, ww);
    area = 0.0;
    for(i=0; i<nn; i++)
    {
      area += ww[i]*f42(xx[i],yy[i]);
    }
    printf("n=%d, nn=%d, area=%f \n", n, nn, area);
  }
  printf("\n");

  printf("fe2(x,y)=exp(x+y) \n");
  for(n=1; n<=6; n++)
  {
    nn = triquad_int(n, x1, y1, x2, y2, x3, y3, xx, yy, ww);
    area = 0.0;
    for(i=0; i<nn; i++)
    {
      area += ww[i]*fe2(xx[i],yy[i]);
    }
    printf("n=%d, nn=%d, area=%f \n", n, nn, area);
  }
  printf("\n");

  printf("phi(cm1)*fe2() \n");
  for(n=1; n<=6; n++)
  {
    nn = triquad_int(n, x1, y1, x2, y2, x3, y3, xx, yy, ww);
    area = 0.0;
    for(i=0; i<nn; i++)
    {
      area += ww[i]*phi22(cm1,xx[i],yy[i])*fe2(xx[i],yy[i]);
    }
    printf("n=%d, nn=%d, area=%f \n", n, nn, area);
  }
  printf("\n");

  printf("phi(cm2)*fe2() \n");
  for(n=1; n<=6; n++)
  {
    nn = triquad_int(n, x2, y2, x3, y3, x1, y1, xx, yy, ww);
    area = 0.0;
    for(i=0; i<nn; i++)
    {
      area += ww[i]*phi22(cm2,xx[i],yy[i])*fe2(xx[i],yy[i]);
    }
    printf("n=%d, nn=%d, area=%f \n", n, nn, area);
  }
  printf("\n");

  printf("phi(cm3)*fe2() \n");
  for(n=1; n<=6; n++)
  {
    nn = triquad_int(n, x3, y3, x1, y1, x2, y2, xx, yy, ww);
    area = 0.0;
    for(i=0; i<nn; i++)
    {
      area += ww[i]*phi22(cm3,xx[i],yy[i])*fe2(xx[i],yy[i]);
    }
    printf("n=%d, nn=%d, area=%f \n", n, nn, area);
  }
  printf("\n");


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