#include "dist.h"

static int VertexOK(vertex * v) 
{
  if (v == NULL) 
    return 0;
  return 1 ;
}

static int NbrOk (neibour_descr *nd) {
  if (nd == NULL) 
    return -1;
  if (nd->atom == NULL) 
    return -1 ;
  return nd->atom->number ;
}

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

  printf("DumpAFace(): BEGIN face %3d\n",ii);
  printf("  -- Face between atom %3d and neighbour %3d, which are separated by %7.3f A\n",
	     a->number, nd->atom->number, nd->distance);
  printf("  -- List of %3d vertices: number, derived from atoms (4 IDs), coord. (x,y,z)\n",
	 n_vertices_count);
  for (i=0; i< n_vertices_count; i++) {
    vertex * v = FaceVertices[i];
    if (VertexOK(v)) {
      printf("%3d   %5d %5d %5d %5d   %7.4f %7.4f %7.4f\n",i,a->number,
	     NbrOk(v->a1),NbrOk(v->a2),NbrOk(v->a3),v->x,v->y,v->z);
      xs += v->x ;
      ys += v->y ;
      zs += v->z ;
    } else {
      printf("%3d   %5d ** BAD VERTEX **\n",i,a->number);
    }
  }
  assert(n_vertices_count > 0);
  xs /= (n_vertices_count) ;
  ys /= (n_vertices_count) ;
  zs /= (n_vertices_count) ;

  {
    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   face %3d\n\n",ii);
}


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

  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 ;
	}
      }
      
      if (n_vertices_count < 3) {
	STDERR("FullDumpPoly(): Less than 3 vertices skipping neighbor...");
	exit(1);
      } 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 ;
	
	FaceVertices[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;
	  third_vertex->flag = 1;
	  FaceVertices[ll] = third_vertex ;
	  second_vertex = third_vertex ;
	} /* for ll */ 

	DumpAFace(atom,nd,ii,n_vertices_count,FaceVertices,FaceVolume);
	volume += FaceVolume ;
      } /* else */
    }
  }
  printf ("FullDumpPoly(): END polyhedron for atom with ID %5d and total volume %9.4f.\n\n\n",
	  atom->number,volume);
  atom->volume = volume ;
}

void SaveAtomVertices(int n,atom_record * a, vertex *v1, vertex *v2) {}

void LoopOverNeighbors
(atom_record *atom,neibour_descr *neibours,int neibour_count) {}


main(int argc,char **argv)
{
  SAY_HI_FIRST ;
  CITATION ; 
  default_initialization();
  {
    file_records *f = open_and_read_pdb(STRG_AFTER_OPT_W_ENV("i","in.pdb"),0);
    char *outputfile = STRG_AFTER_OPT_W_ENV("o",NULL);
    FILE *ff = outputfile == NULL ? stdout : fopen(outputfile, "w") ;

    CalcVHookFcn = FullDumpPoly ;
    NN_Calculate_Volumes(f, (int) DoubleAfterOption(argc,argv,"method",2),
			 64.0, 6.0);

  }
}