# include "NCube.h"

int compar(AtomList ** e1, AtomList ** e2)   /* Compare funct for qsort */
{
  if ((*e1)->Dist <  (*e2)->Dist) return(-1);
  if ((*e1)->Dist == (*e2)->Dist) return(0);
  if ((*e1)->Dist >  (*e2)->Dist) return(1);
}

typedef int (*comp)(const void *, const void *);

AtomList * SortList(AtomList * head)
{
# define MAX_NBR 1000
  AtomList * Myarray[MAX_NBR];
  int j,n;
  AtomList * thehead;
  
  n=0;

# ifdef OFF
  fprintf(stderr,"Sorting. \n");
# endif
  thehead = head;
  
  while (head != NULL) {
    Myarray[n]=head;
    n++;
    head=head->NextAtom; 
  } 
  
  qsort(Myarray,n,sizeof(AtomList *),(comp)compar);        /*  !!!!!! */
  
  for (j=0; j<n; j++) {
    Myarray[j]->NextAtom=Myarray[j+1];
#   ifdef OFF
    printf(" %1.3f \n", Myarray[j]->Dist);
#   endif
  }
  
  head=Myarray[0];
  Myarray[n-1]->NextAtom=NULL;
  return head ;
# undef MAX_NBR
}

FindProteinDimensions(NCube * Pcube, file_records * f)
{
   double maxX,maxY,maxZ,dx,dy,dz;
   double minX,minY,minZ,cubelen;
   atom_record *atom;
   AtomList *temp;
   int i;

   maxX=0;   maxY=0;   maxZ=0;     /* Arbitrary initial settings */
   minX=100; minY=100; minZ=100;
    
   fprintf(stderr,"Calculating Protein Size. \n");
 
   for (i=0; i<f->atomnum; i++) {
      atom = f->atoms[i];
        
        if (maxX < atom->x) maxX= atom->x;
        if (maxY < atom->y) maxY= atom->y;
        if (maxZ < atom->z) maxZ= atom->z;

        if (minX > atom->x) minX= atom->x;
        if (minY > atom->y) minY= atom->y;
        if (minZ > atom->z) minZ= atom->z;


    }  /* for i */

  Pcube->dimx=maxX - minX;     
  Pcube->dimy=maxY - minY;
  Pcube->dimz=maxZ - minZ;

  Pcube->MinX = minX;
  Pcube->MinY = minY;
  Pcube->MinZ = minZ;
}

void PrintAtomList(AtomList * head)
{
  char buf[BUFSIZ];
  while(head != NULL) {
    sprintfAtomShortForm(buf,head->atom) ;
    printf("%30s %8.3f\n",buf,head->Dist);
    head = head->NextAtom;
  }
}

void PrintNeighbors(AtomList * head, atom_record * a)
{
  char buf[BUFSIZ];
  sprintfAtomShortForm(buf,a);
  printf("PrintNeigbors(): Neighbors of atom %20s \n",buf);
  PrintAtomList(head);
}

AssignAtomsToCubes(NCube * Pcube, file_records * f)
{

  int i,cx,cy,cz;
  atom_record *atom;
  double cubelen;

  cubelen=5+Pcube->rcut;
  
  for (i=0; i<f->atomnum; i++) {
    atom = f->atoms[i];
    
    cx=(int)ceil( ((atom->x)-Pcube->MinX)/cubelen);  
    cy=(int)ceil( ((atom->y)-Pcube->MinY)/cubelen);
    cz=(int)ceil( ((atom->z)-Pcube->MinZ)/cubelen);
    
    if (cx==0) cx=1;
    if (cy==0) cy=1;
    if (cz==0) cz=1;
    
    AddAtomToList(Pcube,atom,&cx,&cy,&cz);
  }       /* for i */
}


AddAtomToList(NCube * Pcube, atom_record * atom, int * i, int * j, int * k)
{
  AtomList *temp;
  
  temp=Pcube->GridData[*i][*j][*k];
  
  if ( temp == NULL) {
    temp=(AtomList *)malloc(sizeof(AtomList));
    temp->NextAtom = NULL;
    temp->atom = atom;
    Pcube->GridData[*i][*j][*k] = temp;
  } else { 
    temp=(AtomList *)malloc(sizeof(AtomList));
    temp->NextAtom = Pcube->GridData[*i][*j][*k];
    temp->atom = atom;
    Pcube->GridData[*i][*j][*k] = temp;
  }  
}

AtomList * N_GetMyNeighbors (atom_record * atom, NCube * Pcube, double rcut)
{
  double dx,dy,dz,ex,ey,ez;
  double adist, vdw1, vdw2, thecut;   
  atom_record *twoatom;
  AtomList *temp, *temp2, *head;
  int ix,iy,iz,cx,cy,cz,j;
  double cubelen;
  
  cubelen=5+Pcube->rcut;
  
  j=0;
  head = NULL;
  
  if (atom != NULL) {
    
    ex=atom->x;
    ey=atom->y;
    ez=atom->z;

    cx=(int)ceil( (ex-Pcube->MinX)/cubelen);  
    cy=(int)ceil( (ey-Pcube->MinY)/cubelen);
    cz=(int)ceil( (ez-Pcube->MinZ)/cubelen);
    
    if (cx==0) cx=1;
    if (cy==0) cy=1;
    if (cz==0) cz=1;
    
       for (ix=cx-1; ix<cx+2; ix++) {
	 for (iy=cy-1; iy<cy+2; iy++) {
	   for (iz=cz-1; iz<cz+2; iz++) {
	     
	     if ( (ix>0) && (iy>0) && (iz>0) && (ix<=Pcube->dimx) 
		 && (iy<=Pcube->dimy) && (iz<=Pcube->dimz) )  {
	       
	       if (Pcube->GridData[ix][iy][iz] != NULL) 
                 temp2 = Pcube->GridData[ix][iy][iz];
	       else temp2 = NULL;
	       
	       while (temp2 != NULL)  {
		 
		 twoatom=temp2->atom;
		 assert(twoatom != NULL);
		 
		 vdw2=twoatom->definition->vdw;
		 
		 dx=twoatom->x;
		 dy=twoatom->y;
		 dz=twoatom->z;
		 
		 adist=sqrt( ((ex-dx)*(ex-dx)) +
			    ((ey-dy)*(ey-dy)) +
			    ((ez-dz)*(ez-dz)) );
		 
		 thecut=vdw1+vdw2+rcut;           /* Cutoff distance */
		 
		 if ((adist < thecut) && (atom->number != twoatom->number)) {
		   j++;     
		   if (head == NULL) {
		     head=(AtomList *)malloc(sizeof(AtomList));
		     head->NextAtom = NULL;
		     head->atom = twoatom; 
		     head->Dist = adist;
		   } else  {
		     temp=(AtomList *)malloc(sizeof(AtomList));
		     temp->atom = twoatom;
		     temp->NextAtom = head;
		     temp->Dist = adist;
		     head = temp;
		   }
		 } 
		 
		 temp2 = temp2->NextAtom;
		 
	       } /* while  */
	       
	     } /* if ix  */
	   } /* for iz */
	 } /* for iy */
       } /* for ix */
    
  } /* if atom */
  return(head);  
}

void N_do_Cubeinit(NCube * Pcube)
{
  int i,j,k;
  fprintf(stderr,"N_do_Cubeinit(): Cube initialization.\n");   
  for (i=0; i<MAXDIM; i++)
    for (j=0; j<MAXDIM; j++)
      for (k=0; k<MAXDIM; k++)  {
	Pcube->GridData[i][j][k]=NULL;
	Pcube->GridData[i][j][k]=NULL;
      }
}

void N_SetupCube(NCube * Pcube, file_records * f) 
{
  double dx,dy,dz,cubelen;
  
  FindProteinDimensions(Pcube, f);
  
  dx=Pcube->dimx;     
  dy=Pcube->dimy;
  dz=Pcube->dimz;
  
  cubelen=5+Pcube->rcut;       /* Set the size of the cubes. */
                               /* vdw of Carbon =2.0 A.      */
  Pcube->dimx = (int)ceil(dx/cubelen);
  Pcube->dimy = (int)ceil(dy/cubelen);
  Pcube->dimz = (int)ceil(dz/cubelen);

  if (Pcube->dimx==0) Pcube->dimx=1;
  if (Pcube->dimy==0) Pcube->dimy=1;
  if (Pcube->dimz==0) Pcube->dimz=1;
 
  fprintf(stderr,"N_SetupCube(): Creating Cube of Size: ");
  fprintf(stderr," %d  %d  %d \n",Pcube->dimx,Pcube->dimy,Pcube->dimz);

  AssignAtomsToCubes(Pcube, f);
} /* Setup cube */


void CalculateAllNeighbors(file_records * f, double rcut)
{
  int ii;
  NCube *Pcube = AllocateDefaultCube();

  STDERR("N_CalculateAllNeighbors() begins");
  Pcube->rcut = rcut;             /* ****** Cut off value ****** */
  N_do_Cubeinit(Pcube);
  N_SetupCube(Pcube, f);

  for (ii=0; ii<f->atomnum; ii++) {
    atom_record * a = f->atoms[ii] ;
    AtomList * b = N_GetMyNeighbors(a,Pcube,rcut);
    a->Contacts = SortList(b);
  }
  STDERR("N_CalculateAllNeighbors() ends");
}