#include "ModifiedVolume.h"

static void DumpAFace 
(neibour_descr * nd, atom_record * a, int ll,vertex * FaceVertices[],
 double FaceVolume)
{
  double xs = 0;
  double ys = 0;
  double zs = 0;
  int i ;

  printf("DumpAFace(): BEGIN\n");
  printf("  -- Face between atom %3d and neighbour %3d, which are separated by %7.3f A\n",
	     a->number, nd->atom->number, nd->distance);
  printf("  -- Vertex List: number, derived from atoms (4 IDs), coordinates (x,y,z)\n");
  for (i=0; i<= ll; i++) {
    vertex * v = FaceVertices[i];
    printf("%3d   %5d %5d %5d %5d   %7.4f %7.4f %7.4f\n",i,a->number,
	   v->a1->number,v->a2->number,v->a3->number,v->x,v->y,v->z);
    xs += x ;
    ys += y ;
    zs += z ;
  }
  xs /= (ll+1) ;
  ys /= (ll+1) ;
  zs /= (ll+1) ;
  {
    double dx = a->x - xs;
    double dy = a->y - ys;
    double dz = a->z - zs;
    double d = sqrt( dx*dx + dy*dy + dz*dz);
    double area = 3 * FaceVolume / d ;
    printf("  -- Face-Centroid= %7.4f %7.4f %7.4f\n",xs,ys,zs,area);
    printf("  -- Distance of face to central atom: %8.4f\n",d);
    printf("  -- Face-Area= %9.4f   Pyramid-Volume=%9.4f\n",area,FaceVolume);
  }
  printf("DumpAFace(): END\n");
}


void CalcVHookFcn
     (atom_record *atom,neibour_descr *neibours,int neibour_count)
{
  int ii,jj,n_vertices_count,ll ;
  double volume = 0.0 ; 
  printf ("CalcVHookFcn: BEGIN polyhedron for following atom, which has ID %5d.\n",
	  atom->number);
  write_pdb_record(stdout,atom);

  for (ii=0 ; ii < neibour_count ; ii++) {
    neibour_descr * nd = &neibours[ii] ; 
    if (nd->flag) {
      vertex * n_vertices[MAX_VERTEX] ;
      n_vertices_count= 0;
      for (jj = 0 ; jj < vertices_count; jj++) {
	vertex *ve = &vertices[jj];
	if (ve->a1 == nd || ve->a2 == nd || ve->a3 == nd ) {
	  ve->flag = 0 ;
	  n_vertices[n_vertices_count++] =ve ;
	}
      }
      
      printf("CalcVHookFcn: BEGIN face %3d\n",ii);
      
      if (n_vertices_count < 3) {
	STDERR("CalcVHookFcn: Less than 3 vertices skipping neighbor...");
      } else  { 
	double FaceVolume = 0;
	vertex * FaceVertices[MAX_VERTEX];
	vertex *first_vertex = n_vertices[0] ;
	vertex *second_vertex, *third_vertex, *previous_second = first_vertex;
	neibour_descr *xnd = first_vertex->a1;
	
	if (xnd == nd)
	  xnd = first_vertex->a2 ;

	for(jj =1 ; jj< n_vertices_count ; jj++)
	  if (n_vertices[jj]->a1 == xnd ||
	      n_vertices[jj]->a2 == xnd ||
	      n_vertices[jj]->a3 == xnd ) {
	    second_vertex = n_vertices[jj] ;
	    break ;
	  }
	first_vertex->flag = 1;
	second_vertex->flag = 1 ;
	
	FaveVertices[0] = first_vertex;
	FaceVertices[1] = second_vertex ;

	for (ll = 2 ; ll < n_vertices_count ; ll++) {
	  neibour_descr *next_nd = second_vertex->a1 ;
	  
	  if (next_nd == nd)
	    if  (second_vertex->a2 == xnd )
	      next_nd = second_vertex->a3 ;
	    else next_nd = second_vertex->a2 ;
	  else if (next_nd == xnd)
	    if  (second_vertex->a2 == nd )
	      next_nd = second_vertex->a3 ;
	    else next_nd = second_vertex->a2 ;
	  
	  for (jj = 1  ; jj < n_vertices_count ; jj++) {
	    third_vertex = n_vertices[jj] ;
	    if (!third_vertex->flag &&
		(third_vertex->a1 == next_nd ||
		 third_vertex->a2 == next_nd ||
		 third_vertex->a3 == next_nd))
	      break ;
	  }
	  if (jj == n_vertices_count) {
	    warning(ARG_TYPE_ATOM,"Did not find third_vertex", atom) ;
	    return;
	  }
	  { 
	    double t_volume = 
	      tet_volume(atom,first_vertex, second_vertex, third_vertex);
	    FaceVolume += t_volume ;
	    
	  }
	  xnd = next_nd;
	  FaceVertices[ll] = second_vertex ;
	  third_vertex->flag = 1;
	  second_vertex = third_vertex ;
	} /* for ll */ 

	DumpAFace(atom,nd,ll,FaceVertices,FaceVolume);
	printf ("CalcVHookFcn: END face %3d\n",ii);
	volume += FaceVolume ;
      }
    }
  }
  
  printf ("CalcVHookFcn: END polyhedron for atom with ID %5d and total volume %9.4f.\n\n",
	  atom->number,volume);
  atom->volume = volume ;
  atom->dumm_solvents = n_solvent ;
  SaveAtomVertices(5,atom,NULL,NULL);
}