/* lagrange_3d.c  test, second order  h/2 on a tetrahedron 
 *
 * (( x-x1)^2+( y-y1)^2+( z-z1)^2)   (( x-x2)^2+( y-y2)^2+( z-z2)^2)     
 * ------------------------------- * --------------------- * ...  
 * ((x0-x1)^2+(y0-y1)^2+(z0-z1)^2)   ((x0-x2)^2+(y0-y2)^2+(z0-z2)^2)     
 *
 * zero at x1,y1,z1  x2,y2,z2      1.000  at  x0,y0,z0
 *
 */
#include <stdio.h>
#include <stdlib.h>

#define lterm3(x,y,z,xi,yi,zi,x0,y0,z0) (((x-xi)*(x-xi)+(y-yi)*(y-yi)+(z-zi)*(z-zi))/((x0-xi)*(x0-xi)+(y0-yi)*(y0-yi)+(z0-zi)*(z0-zi)))

int main(int argc, char *argv[])
{
  int i, j, k;
  double x0=1.0, y0=1.0, z0=1.0;
  double x1=2.0, y1=2.0, z1=1.0;
  double x2=2.0, y2=1.0, z2=1.0;
  double x3=2.0, y3=2.0, z3=2.0;
  double dx=0.1, dy=0.1, dz=0.1;
  double x01, y01, z01, x02, y02, z02, x03, y03, z03; /* mid points */
  double x12, y12, z12, x13, y13, z13, x23, y23, z23; /* mid points */
  double v, x, y, z;
  double s[11][11][11], xp[11][11][11], yp[11][11][11], zp[11][11][11];

  printf("lagrange_3d running \n");
  x01=(x0+x1)/2.0;
  y01=(y0+y1)/2.0;
  z01=(z0+z1)/2.0;
  x02=(x0+x2)/2.0;
  y02=(y0+y2)/2.0;
  z02=(z0+z2)/2.0;
  x03=(x0+x3)/2.0;
  y03=(y0+y3)/2.0;
  z03=(z0+z3)/2.0;
  x12=(x1+x2)/2.0;
  y12=(y1+y2)/2.0;
  z12=(z1+z2)/2.0;
  x13=(x1+x3)/2.0;
  y13=(y1+y3)/2.0;
  z13=(z1+z3)/2.0;
  x23=(x2+x3)/2.0;
  y23=(y2+y3)/2.0;
  z23=(z2+z3)/2.0;

  for(i=0; i<11; i++)
  {
    x=1.0+i*dx;
    for(j=0; j<11; j++)
    {
      y=1.0+j*dy;
      for(k=0; k<11; k++)
      {
        z=1.0+k*dz;
        v =lterm3(x,y,z,x1,y1,z1,x0,y0,z0);
        v*=lterm3(x,y,z,x2,y2,z2,x0,y0,z0);
        v*=lterm3(x,y,z,x3,y3,z3,x0,y0,z0);
        v*=lterm3(x,y,z,x01,y01,z01,x0,y0,z0);
        v*=lterm3(x,y,z,x02,y02,z02,x0,y0,z0);
        v*=lterm3(x,y,z,x03,y03,z03,x0,y0,z0);
        v*=lterm3(x,y,z,x12,y12,z12,x0,y0,z0);
        v*=lterm3(x,y,z,x13,y13,z13,x0,y0,z0);
        v*=lterm3(x,y,z,x23,y23,z23,x0,y0,z0);
        s[i][j][k]=v;
        xp[i][j][k]=x;
        yp[i][j][k]=y;
        zp[i][j][k]=z;
        if(y<=x && z<=x) /* inside tetrahedron */
        {
          printf("x=%9.6f, y=%9.6f, z=%9.6f, v=%9.6f \n", x, y, z, v);
        }
      }
    }
  }
  printf("lagrange_3d finished \n");

  return 0;
} /* end lagrange_2d */