/* lagrange_4d.c  test, second order  h/2 on a 4D element
 *
 * (( x-x1)^2+( y-y1)^2+( z-z1)^2+( w-w1)^2)   
 * ----------------------------------------- * ...
 * ((x0-x1)^2+(y0-y1)^2+(z0-z1)^2+(w0-w1)^2)   
 *
 * zero at x1,y1,z1,w1  x2,y2,z2,w2      1.000  at  x0,y0,z0,w0
 *
 */
#include <stdio.h>
#include <stdlib.h>

#define lterm4(x,y,z,w,xi,yi,zi,wi,x0,y0,z0,w0) (((x-xi)*(x-xi)+(y-yi)*(y-yi)+(z-zi)*(z-zi)+(w-wi)*(w-wi))/((x0-xi)*(x0-xi)+(y0-yi)*(y0-yi)+(z0-zi)*(z0-zi)+(w0-wi)*(w0-wi)))

int main(int argc, char *argv[])
{
  int i, j, k, m;
  double x0=1.0, y0=1.0, z0=1.0, w0=1.0;
  double x1=2.0, y1=2.0, z1=1.0, w1=1.0;
  double x2=2.0, y2=1.0, z2=1.0, w2=1.0;
  double x3=2.0, y3=2.0, z3=2.0, w3=1.0;
  double x4=2.0, y4=2.0, z4=2.0, w4=2.0;
  double dx=0.25, dy=0.25, dz=0.25, dw=0.25;
  double x01, y01, z01, w01, x02, y02, z02, w02; /* mid points */
  double x03, y03, z03, w03, x04, y04, z04, w04; /* mid points */
  double x12, y12, z12, w12, x13, y13, z13, w13; /* mid points */
  double x14, y14, z14, w14, x23, y23, z23, w23; /* mid points */
  double x24, y24, z24, w24, x34, y34, z34, w34; /* mid points */
  double v, x, y, z, w;

  printf("lagrange_4d running \n");
  x01=(x0+x1)/2.0;
  y01=(y0+y1)/2.0;
  z01=(z0+z1)/2.0;
  w01=(w0+w1)/2.0;
  x02=(x0+x2)/2.0;
  y02=(y0+y2)/2.0;
  z02=(z0+z2)/2.0;
  w02=(w0+w2)/2.0;
  x03=(x0+x3)/2.0;
  y03=(y0+y3)/2.0;
  z03=(z0+z3)/2.0;
  w03=(w0+w3)/2.0;
  x04=(x0+x4)/2.0;
  y04=(y0+y4)/2.0;
  z04=(z0+z4)/2.0;
  w04=(w0+w4)/2.0;
  x12=(x1+x2)/2.0;
  y12=(y1+y2)/2.0;
  z12=(z1+z2)/2.0;
  w12=(w1+w2)/2.0;
  x13=(x1+x3)/2.0;
  y13=(y1+y3)/2.0;
  z13=(z1+z3)/2.0;
  w13=(w1+w3)/2.0;
  x14=(x1+x4)/2.0;
  y14=(y1+y4)/2.0;
  z14=(z1+z4)/2.0;
  w14=(w1+w4)/2.0;
  x23=(x2+x3)/2.0;
  y23=(y2+y3)/2.0;
  z23=(z2+z3)/2.0;
  w23=(w2+w3)/2.0;
  x24=(x2+x4)/2.0;
  y24=(y2+y4)/2.0;
  z24=(z2+z4)/2.0;
  w24=(w2+w4)/2.0;
  x34=(x3+x4)/2.0;
  y34=(y3+y4)/2.0;
  z34=(z3+z4)/2.0;
  w34=(w3+w4)/2.0;

  for(i=0; i<5; i++)
  {
    x=1.0+i*dx;
    for(j=0; j<5; j++)
    {
      y=1.0+j*dy;
      for(k=0; k<5; k++)
      {
        z=1.0+k*dz;
        for(m=0; m<5; m++)
        {
          w=1.0+m*dw;
          v =lterm4(x,y,z,w,x1,y1,z1,w1,x0,y0,z0,w0);
          v*=lterm4(x,y,z,w,x2,y2,z2,w2,x0,y0,z0,w0);
          v*=lterm4(x,y,z,w,x3,y3,z3,w3,x0,y0,z0,w0);
          v*=lterm4(x,y,z,w,x4,y4,z4,w4,x0,y0,z0,w0);
          v*=lterm4(x,y,z,w,x01,y01,z01,w01,x0,y0,z0,w0);
          v*=lterm4(x,y,z,w,x02,y02,z02,w02,x0,y0,z0,w0);
          v*=lterm4(x,y,z,w,x03,y03,z03,w03,x0,y0,z0,w0);
          v*=lterm4(x,y,z,w,x04,y04,z04,w04,x0,y0,z0,w0);
          v*=lterm4(x,y,z,w,x12,y12,z12,w12,x0,y0,z0,w0);
          v*=lterm4(x,y,z,w,x13,y13,z13,w13,x0,y0,z0,w0);
          v*=lterm4(x,y,z,w,x14,y14,z14,w14,x0,y0,z0,w0);
          v*=lterm4(x,y,z,w,x23,y23,z23,w23,x0,y0,z0,w0);
          v*=lterm4(x,y,z,w,x24,y24,z24,w24,x0,y0,z0,w0);
          v*=lterm4(x,y,z,w,x34,y34,z34,w34,x0,y0,z0,w0);
          if(y<=x && z<=x && w<=x) /* inside 4D element */
          {
            printf("x=%11.9f, y=%11.9f, z=%11.9f, w=%11.9f, v=%11.9f \n",
		   x, y, z, w, v);
          }
	}
      }
    }
  }
  printf("lagrange_4d finished \n");

  return 0;
} /* end lagrange_4d */