/* lagrange_2d.c  test, second order  h/2 on a triangle 
 *
 * (( x-x1)^2+( y-y1)^2)   (( x-x2)^2+( y-y2)^2)   (( x-x3)^2+( y-y3)^2)  
 * --------------------- * --------------------- * ---------------------  
 * ((x0-x1)^2+(y0-y1)^2)   ((x0-x2)^2+(y0-y2)^2)   ((x0-x3)^2+(y0-y3)^2)  
 *
 * zero at x1,y1  x2,y2  x3,y3    one at  x0,y0
 *
 */
#include <stdio.h>
#include <stdlib.h>

#define lterm(x,y,xi,yi,x0,y0) (((x-xi)*(x-xi)+(y-yi)*(y-yi))/((x0-xi)*(x0-xi)+(y0-yi)*(y0-yi)))

int main(int argc, char *argv[])
{
  int i, j;
  double x0=1.0, y0=1.0;
  double x1=2.0, y1=2.0;
  double x2=2.0, y2=1.0;
  double dx=0.1, dy=0.1;
  double x01, y01, x02, y02, x12, y12; /* mid points */
  double v, x, y;
  double z[11][11], xp[11][11], yp[11][11];

  printf("lagrange_2d running \n");
  x01=(x0+x1)/2.0;
  y01=(y0+y1)/2.0;
  x02=(x0+x2)/2.0;
  y02=(y0+y2)/2.0;
  x12=(x1+x2)/2.0;
  y12=(y1+y2)/2.0;

  for(i=0; i<11; i++)
  {
    x=1.0+i*dx;
    for(j=0; j<11; j++)
    {
      y=1.0+j*dy;
      v =lterm(x,y,x1,y1,x0,y0);
      v*=lterm(x,y,x2,y2,x0,y0);
      v*=lterm(x,y,x01,y01,x0,y0);
      v*=lterm(x,y,x02,y02,x0,y0);
      v*=lterm(x,y,x12,y12,x0,y0);
      z[i][j]=v;
      xp[i][j]=x;
      yp[i][j]=y;
      if(y<=x)
      {
        printf("x=%9.5f, y=%9.5f, z=%9.5f \n", x, y, v);
      }
    }
  }
  printf("lagrange_2d finished \n");

  return 0;
} /* end lagrange_2d */