/* test_tri_basis.c   test evaluation of derivatives */
#include <stdio.h>
#include <stdlib.h>
#include "tri_basis.h"
#include "deriv.h"

int main(int argc, char * argv[])
{
  double x, y, up, upp;
  double xx[12], yy[12];
  double u[12], u2[12][12];
  double T[6];
  int i, j, k, m;
  double hx=0.2;
  double hy=0.25;
  int a1[ 5], b1; /* first derivative */
  int a2[ 7], b2; /* second derivative */

  double dx, dxx;
  double dy, dyy, dxy;

  printf("test_tri_basis.c running \n");
  deriv(1, 5,2,a1,&b1);
  deriv(2, 7,3,a2,&b2);

  x=0.1;
  for(i=0; i<12; i++)
  {
    xx[i]=x;
    x=x+hx;
  }

  y=0.2;
  for(i=0; i<12; i++)
  {
    yy[i]=y;
    y=y+hy;
  }

  T[0]=0.2;
  T[1]=0.3;
  T[2]=1.5;
  T[3]=3.2;
  T[4]=3.3;
  T[5]=0.5;

  for(i=0; i<5; i++)
  {
    for(j=0; j<5; j++)
    {
      u2[i][j]=phi2p(T,xx[i],yy[j]); /* base phi2p */
    }
  }
  dx=(1.0/(b1*hx))*(a1[0]*u2[0][2]+a1[1]*u2[1][2]+a1[2]*u2[2][2]+
                    a1[3]*u2[3][2]+a1[4]*u2[4][2]);
  up=phi2x(T,xx[2],yy[2]);
  printf("phi2p,phi2x=%f, dx=%f, err=%f \n", up, dx, up-dx);

  dy=(1.0/(b1*hy))*(a1[0]*u2[2][0]+a1[1]*u2[2][1]+a1[2]*u2[2][2]+
                    a1[3]*u2[2][3]+a1[4]*u2[2][4]);
  up=phi2y(T,xx[2],yy[2]);
  printf("phi2p,phi2y=%f, dy=%f, err=%f \n", up, dy, up-dy);


  for(i=0; i<7; i++)
  {
    for(j=0; j<7; j++)
    {
      u2[i][j]=phi2p(T,xx[i],yy[j]); /* base phi2p */
    }
  }

  upp=phi2xx(T,xx[3],yy[3]);
  dxx=(1.0/(b2*hx*hx))*(a2[0]*u2[0][3]+a2[1]*u2[1][3]+a2[2]*u2[2][3]+
                        a2[3]*u2[3][3]+a2[4]*u2[4][3]+
                        a2[5]*u2[5][3]+a2[6]*u2[6][3]);
  printf("phi2p,phi2xx=%f, dxx=%f, err=%f \n", upp, dxx, upp-dxx);

  for(i=0; i<5; i++)
  {
    u[i]=(1.0/(b1*hy))*(a1[0]*u2[i][0]+a1[1]*u2[i][1]+a1[2]*u2[i][2]+
                        a1[3]*u2[i][3]+a1[4]*u2[i][4]);
  }
  dxy=(1.0/(b1*hx))*(a1[0]*u[0]+a1[1]*u[1]+a1[2]*u[2]+a1[3]*u[3]+a1[4]*u[4]);
  up=phi2xy(T,xx[2],yy[2]);
  printf("phi2p,phi2xy=%f, dxy=%f, err=%f \n", up, dxy, up-dxy);

  upp=phi2yy(T,xx[3],yy[3]);
  dyy=(1.0/(b2*hy*hy))*(a2[0]*u2[3][0]+a2[1]*u2[3][1]+a2[2]*u2[3][2]+
                        a2[3]*u2[3][3]+a2[4]*u2[3][4]+
                        a2[5]*u2[3][5]+a2[6]*u2[3][6]);
  printf("phi2p,phi2yy=%f, dyy=%f, err=%f \n", upp, dyy, upp-dyy);

  for(i=0; i<5; i++)
  {
    for(j=0; j<5; j++)
    {
      u2[i][j]=phi2m(T,xx[i],yy[j]); /* base phi2m */
    }
  }
  dx=(1.0/(b1*hx))*(a1[0]*u2[0][2]+a1[1]*u2[1][2]+a1[2]*u2[2][2]+
                    a1[3]*u2[3][2]+a1[4]*u2[4][2]);
  up=phi2mx(T,xx[2],yy[2]);
  printf("phi2m,phi2mx=%f, dx=%f, err=%f \n", up, dx, up-dx);

  dy=(1.0/(b1*hy))*(a1[0]*u2[2][0]+a1[1]*u2[2][1]+a1[2]*u2[2][2]+
                    a1[3]*u2[2][3]+a1[4]*u2[2][4]);
  up=phi2my(T,xx[2],yy[2]);
  printf("phi2m,phi2my=%f, dy=%f, err=%f \n", up, dy, up-dy);


  for(i=0; i<7; i++)
  {
    for(j=0; j<7; j++)
    {
      u2[i][j]=phi2m(T,xx[i],yy[j]); /* base phi2m */
    }
  }

  upp=phi2mxx(T,xx[3],yy[3]);
  dxx=(1.0/(b2*hx*hx))*(a2[0]*u2[0][3]+a2[1]*u2[1][3]+a2[2]*u2[2][3]+
                        a2[3]*u2[3][3]+a2[4]*u2[4][3]+
                        a2[5]*u2[5][3]+a2[6]*u2[6][3]);
  printf("phi2m,phi2mxx=%f, dxx=%f, err=%f \n", upp, dxx, upp-dxx);

  for(i=0; i<5; i++)
  {
    u[i]=(1.0/(b1*hy))*(a1[0]*u2[i][0]+a1[1]*u2[i][1]+a1[2]*u2[i][2]+
                        a1[3]*u2[i][3]+a1[4]*u2[i][4]);
  }
  dxy=(1.0/(b1*hx))*(a1[0]*u[0]+a1[1]*u[1]+a1[2]*u[2]+a1[3]*u[3]+a1[4]*u[4]);
  up=phi2mxy(T,xx[2],yy[2]);
  printf("phi2m,phi2mxy=%f, dxy=%f, err=%f \n", up, dxy, up-dxy);

  upp=phi2myy(T,xx[3],yy[3]);
  dyy=(1.0/(b2*hy*hy))*(a2[0]*u2[3][0]+a2[1]*u2[3][1]+a2[2]*u2[3][2]+
                        a2[3]*u2[3][3]+a2[4]*u2[3][4]+
                        a2[5]*u2[3][5]+a2[6]*u2[3][6]);
  printf("phi2m,phi2myy=%f, dyy=%f, err=%f \n", upp, dyy, upp-dyy);


  printf("test_tri_basis.c finished \n");
  return 0;
} /* end test_tri_basis.c */