# include "data.h"

int gridify (int argc, char **argv) {

  int i, Nx, Ny, Xb, Yb;
  float Xmin, Xmax, dX, Ymin, Ymax, dY;
  float *buf, *val, *x, *y, *z;
  int *Nval;
  Buffer *bf;
  Vector *vx, *vy, *vz;

  if (argc != 11) {
    fprintf (stderr, "USAGE: vectobuf x y z buffer Xmin Xmax dX Ymin Ymax dY\n");
    return (FALSE);
  }
  
  if ((vx = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
  if ((vy = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);
  if ((vz = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) return (FALSE);
  if ((bf = SelectBuffer (argv[4], ANYBUFFER, TRUE)) == NULL) return (FALSE);

  if (vx[0].Nelements != vy[0].Nelements) return (FALSE);
  if (vx[0].Nelements != vz[0].Nelements) return (FALSE);

  Xmin = atof (argv[5]);
  Xmax = atof (argv[6]);
  dX   = atof (argv[7]);

  Ymin = atof (argv[8]);
  Ymax = atof (argv[9]);
  dY   = atof (argv[10]);

  Nx = (Xmax - Xmin) / dX + 1;
  Ny = (Ymax - Ymin) / dY + 1;
  
  /* these steps are a standard 'createbuffer', define? */
  bf[0].bitpix = -32;
  bf[0].unsign = FALSE;
  bf[0].bscale = 1.0;
  bf[0].bzero  = 0.0;
  strcpy (bf[0].file, "(empty)");

  bf[0].header.bitpix = bf[0].bitpix;
  bf[0].header.bscale = bf[0].bscale;
  bf[0].header.bzero  = bf[0].bzero; 
  bf[0].header.unsign = bf[0].unsign;
  bf[0].header.Naxes = 2;
  bf[0].header.Naxis[0] = Nx;
  bf[0].header.Naxis[1] = Ny;
  free (bf[0].header.buffer);  /* ALLOCATEd in SelectBuffer */
  free (bf[0].matrix.buffer);  /* ALLOCATEd in SelectBuffer */
  fits_create_header (&bf[0].header);
  fits_create_matrix (&bf[0].header, &bf[0].matrix);

  ALLOCATE (val, float, Nx*Ny);
  bzero (val, Nx*Ny*sizeof(float));
  ALLOCATE (Nval, int, Nx*Ny);
  bzero (Nval, Nx*Ny*sizeof(int));

  x = vx[0].elements;
  y = vy[0].elements;
  z = vz[0].elements;
  for (i = 0; i < vx[0].Nelements; i++, x++, y++, z++) {
    Xb = (*x - Xmin) / dX;
    Yb = (*y - Ymin) / dY;
    if (Xb >= Nx) continue;
    if (Yb >= Ny) continue;
    val[Xb + Yb*Nx] = *z;
    Nval[Xb + Yb*Nx]++;
  }

  buf = (float *) bf[0].matrix.buffer;
  for (i = 0; i < Nx*Ny; i++) {
    if (Nval[i] == 0) {
      buf[i] = 0;
      continue;
    }
    buf[i] = val[i] / Nval[i];
  }

  free (val);
  free (Nval);

  return (TRUE);

}
