#include "NCube.h" 

# ifdef OFF
/* returns pointers to Unit_Cube in three dimensions x, y, & z */
Cube * SetupCube (double *bxsidel, double nb_cutoff)
{
  int 
    xx, yy, zz;

  NCube * CubedBox = (NCube *)malloc(sizeof(Cube));

  AtomList **** C;


  /* defining the number of bins */
  CubedBox->rcut = nb_cutoff;
  CubedBox->dimx = floor(bxsidel[0]/nb_cutoff);
  CubedBox->dimy = floor(bxsidel[1]/nb_cutoff);
  CubedBox->dimz = floor(bxsidel[2]/nb_cutoff);

  /* defining x, y, z size of each unit box           */
  /* using definition defined for minimum coordinates */
  CubedBox->MinX = bxsidel[0]/CubedBox->dimx;
  CubedBox->MinY = bxsidel[0]/CubedBox->dimy;
  CubedBox->MinZ = bxsidel[0]/CubedBox->dimz;

  /* print out cube divisions and # of divisions */
  printf("non bonded cutoff for unit boxes is %d\n", nb_cutoff);
  printf("for x dimension: %10.5lf gives %d bins of %10.5lf A.\n",
	 bxsidel[0], CubedBox->dimx, CubedBox->MinX);
  printf("for y dimension: %10.5lf gives %d bins of %10.5lf A.\n",
	 bxsidel[2], CubedBox->dimy, CubedBox->MinY);
  printf("for z dimension: %10.5lf gives %d bins of %10.5lf A.\n", 
	 bxsidel[3], CubedBox->dimz, CubedBox->MinZ);
  
  C = (AtomList ****)calloc(xbins, sizeof(AtomList ***));
  for(xx=0;xx<xbins;xx++)  {
    C[xx] = (AtomList ***)calloc(ybins, sizeof(AtomList **));
    for(yy=0;yy<ybins;yy++)  {
      C[xx][yy] = (AtomList **)calloc(zbins, sizeof(AtomList*));
      for(zz=0;zz<zbins;zz++)  {
	C[xx][yy][zz] = (AtomList *)malloc(sizeof(AtomList));
	C[xx][yy][zz]->NextAtom = NULL;
      }
    }
  }
  CubedBox->GridData = C;

  return CubedBox;
}
#endif

NCube * AllocateDefaultCube(void)
{
  /* allocate a 3tensor of AtomLists 
     with range t[nrl..nrh][ncl..nch][ndl..ndh] */
  long nrl=0;
  long nrh=MAXDIM;
  long ncl=0;
  long nch=MAXDIM;
  long ndl=0;
  long ndh=MAXDIM;
  NCube * C = (NCube *)calloc(1,sizeof(NCube));
  long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
  AtomListPtr *** t;
  
  /* allocate pointers to pointers to rows */
  t=(AtomListPtr ***) calloc(1,(size_t)((nrow+NR_END)*sizeof(AtomListPtr**)));
  if (!t) nrerror("allocation failure 1 in AllocateDefaultCube()");
  t += NR_END;
  t -= nrl;
  
  /* allocate pointers to rows and set pointers to them */
  t[nrl]=(AtomListPtr **) 
    calloc(1,(size_t)((nrow*ncol+NR_END)*sizeof(AtomListPtr*)));
  if (!t[nrl]) nrerror("allocation failure 2 in AllocateDefaultCube()");
  t[nrl] += NR_END;
  t[nrl] -= ncl;

  /* allocate rows and set pointers to them */
  t[nrl][ncl]=(AtomListPtr *) 
    calloc(1,(size_t)((nrow*ncol*ndep+NR_END)*sizeof(AtomListPtr)));
  if (!t[nrl][ncl]) nrerror("allocation failure 3 in AllocateDefaultCube()");
  t[nrl][ncl] += NR_END;
  t[nrl][ncl] -= ndl;

  for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
  for(i=nrl+1;i<=nrh;i++) {
    t[i]=t[i-1]+ncol;
    t[i][ncl]=t[i-1][ncl]+ncol*ndep;
    for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep;
  }

  /* return pointer to array of pointers to rows */
  C->GridData = t;
  return C;
}