extern "C" {
  void writeH(short *);
}

#include "DAGH.h"

extern "C" {
#include "hdf.h"
}

/* For debugging! */
void main() {
  GridData(3)<short> flag(0,0,0,64,64,64,1);
  GridData(2)<short> flag_2(0,0,128,128,1);
  int thedim = 3;
  flag = 0;
  flag_2 = 0;
  void Cluster3(GridData(3)<short> &flag, double, int, int, BBoxList &);
  void Cluster2(GridData(2)<short> &flag, double, int, int, BBoxList &);
  /* TEST CASES! */
  int whichcase = 1;
  (cout <<"CASE:\n" <<
    "1. Points\n" <<
    "2. Spheres\n" <<
    "3. Ring\n" <<
    "4. Orthogonal Rings\n" <<
    "5. Spherical Shell\n" <<
    "6. TwoD Ring\n" <<
    "7. TwoD Shock\n" <<
    "").flush();
  printf ("case  : ");
  scanf ("%d",&whichcase);
  /* Single point in the middle of the grid */
  double eff;
  int mw;

  printf ("eff   : ");
  scanf("%lf",&eff);
  printf ("Min Width : ");
  scanf ("%d",&mw);

  // The #&@! AIX compiler won't allow me to conditionally declare and
  // initialize these in the case statement without complaining
  int i = 0,j = 0,k = 0;
  int rad=1, i0, j0, k0;
  int radi, rado, h;

  switch (whichcase) {
  case 1:
    /* Points */
    printf ("Enter points, enter i = -1 to done:");
    while (i > -1) {
      printf ("i:");
      scanf("%d",&i);
      if (i > 0) {
	printf ("j:");
	scanf("%d",&j);
	printf ("k:");
	scanf("%d",&k);
	flag(i,j,k) = 1;
      }
    }
    break;
  case 2:
    /* Spheres */
    printf ("Enter spheres, enter rad^2 = -1 to done:\n");
    while (rad > 0) {
      printf ("rad :");
      scanf ("%d",&rad);
      if (rad > 0) {
	printf ("i0 :"); scanf("%d",&i0);
	printf ("j0 :"); scanf("%d",&j0);
	printf ("k0 :"); scanf("%d",&k0);
	for (int i=0;i<=64;i++)
	  for (int j=0;j<=64;j++)
	    for (int k=0;k<=64;k++) {
	      if ((i-i0)*(i-i0) + 
		  (j-j0)*(j-j0) +
		  (k-k0)*(k-k0) < rad*rad)
		flag(i,j,k) = 1;
	    }
      }
    }
    break;
  case 3:
    /* Ring */
    scanf("%d",&radi);
    scanf("%d",&rado);
    scanf("%d",&h);
    for ( i=0;i<=64;i++)
      for (j=0;j<=64;j++)
	for ( k=0;k<=64;k++) {
	  if (abs(k-32) < h) {
	    int tr =(i-32)*(i-32)+(j-32)*(j-32);
	    if (tr < rado*rado && tr > radi*radi)
	      flag(i,j,k) = 1;
	  }
	}
    break;
  case 4:
    /* Double Ring */
    scanf("%d",&radi);
    scanf("%d",&rado);
    scanf("%d",&h);
    for ( i=0;i<=64;i++)
      for (j=0;j<=64;j++)
	for ( k=0;k<=64;k++) {
	  if (abs(k-32) < h) {
	    int tr =(i-32)*(i-32)+(j-32)*(j-32);
	    if (tr < rado*rado && tr > radi*radi)
	      flag(i,j,k) = 1;
	  }
	  if (abs(i-32) < h) {
	    int tr =(k-32)*(k-32)+(j-32)*(j-32);
	    if (tr < rado*rado && tr > radi*radi)
	      flag(i,j,k) = 1;
	  }
	}
    break;
  case 5:
    /* Ring */
    scanf("%d",&radi);
    scanf("%d",&rado);
    scanf("%d",&h);
    for ( i=0;i<=64;i++)
      for (j=0;j<=64;j++)
	for ( k=0;k<=64;k++) {
	    int tr =(i-32)*(i-32)+(j-32)*(j-32) + (k-32)*(k-32);
	    if (tr < rado*rado && tr > radi*radi)
	      flag(i,j,k) = 1;
	}
    break;
  case 6:
    /* 2D Ring */
    thedim = 2;
    scanf("%d",&radi);
    scanf("%d",&rado);
    scanf("%d",&h);
    scanf("%d",&i0);
    scanf("%d",&j0);
    for ( i=0;i<=128;i++)
      for (j=0;j<=128;j++) {
	int tr = (i-i0)*(i-i0) + (j-j0)*(j-j0);
	if (tr < rado*rado && tr > radi*radi)
	  flag_2(i,j) = 1;
      }
    break;
  case 7:
    /* 2D Ring */
    thedim = 2;
    for ( i=3;i<=120;i++)
      for (j=-2;j<=2;j++) {
	flag_2(i,i+j) = 1;
      }
    break;

  default:
    printf ("Please choose a valid case! (you chose %d)\n",whichcase);
    exit(123);
  }
  
  /* A sphere */
  int32  *HDFdim;
  HDFdim = (int32 *)malloc(3*sizeof(int32));
  printf ("Starting Clustere\n");
  system("rm flag.hdf");
  BBoxList result;
  if (thedim == 3) {
    Cluster3(flag,eff,mw,1,result);
    float32 *data;
    data = (float32 *)malloc(flag.size()*sizeof(float32));
    for (int i=0;i<flag.size();i++)
      data[i] = (float)(flag.data())[i];
    for (i=0;i<3;i++)
      HDFdim[i] = flag.extents(i);
    DFSDsetdims(3, HDFdim);
    DFSDsetNT(DFNT_FLOAT32);
    DFSDadddata("flag.hdf",3,HDFdim,data);
    free(data);
  } else {

    Cluster2(flag_2,eff,mw,1,result);
    float32 *data;
    data = (float32 *)malloc(flag_2.size()*sizeof(float32));
    for (int i=0;i<flag_2.size();i++)
      data[i] = (float)(flag_2.data())[i];
    for (i=0;i<2;i++) {
      HDFdim[i] = flag_2.extents(i);
      printf ("Setting extents to %d\n",HDFdim[i]);
    }
    DFSDsetdims(2, HDFdim);
    DFSDsetNT(DFNT_FLOAT32);
    DFSDadddata("flag.hdf",2,HDFdim,data);
    free(data);
  }

  cout << "-------\n"<< result <<"\n\n";

  FILE *ofp = fopen("cluster.bbox","w");
  int c = 0;
  for (BBox* newB = result.first(); newB; newB = result.next()) {
    c++;			// grin
    fprintf (ofp, "%d %d %d %d %d %d %d %d %d %d %d %d %d\n",
	     0,3,0,c,
	     (*newB).extents(0),(*newB).extents(1),(*newB).extents(2),
	     (*newB).lower(0),(*newB).lower(1),(*newB).lower(2),
	     (*newB).upper(0),(*newB).upper(1),(*newB).upper(2));
  }
}