#include <assert.h>
#include "ModifiedVolume.h"

static file_records * global_file ;

/* Added by graham@molmovdb */
void SaveAtomVertices (int i, atom_record *a, vertex *v1, vertex *v2) {}

void LoopOverNeighbors (atom_record *atom, neibour_descr *n_atoms, int n_count)
{}

CalcVHookFcnType CalcVHookFcn = NULL;

void PrintVolumeParameters(void)
{
  STDERR("PrintVolumeParameters():");
  fprintf(stderr,"CalcVHookFcn= %5d  the_method= %1d  \nvolume_max_distance= %7.4f distance_check= %7.4f\n",
	 (int) CalcVHookFcn,the_method,volume_max_distance,distance_check);
  fprintf(stderr,"vertices_count= %4d plot_file= %d\n",
	  vertices_count,(int) plot_file);
}
  

/* [ 4 ] called from [ 2 ]

   Calculate the volume of an atom.  First mark all neibours that are
   still used. Then for each neighbour that is still used (nd),
   collect all the relevant vertices into atom_vertWINDOWID=46137348
ices, and also set
   their flag.  Take one of the vertices, first_vertex. Find one of
   its neighbours; call it second_vertex (they will share another
   neighbour except nd, xnd). Then repeaet: find the next neighbour of
   second_vertex; call it third_vertex (it shares with second_vertex a
   neighbour, next_nd, which is the third neighbour in second_vertex
   which is not xnd or nd). Calculate the volume of the tetrahedron
   (first_vertex,second_vertex, third_vertex, origin_atom) and collect
   it. Set second_vertex to third_vertex.

   We do it this way because Richards does, so it is easier to compare
   the code. A more efficient way is to calculate only the area, and
   multiply by the distance of the plane to to the origin atom (and
   divide by 3).  

   Modified version that calls SaveAtomVertices() : */

void N_calculate_volume 
(atom_record *atom,neibour_descr *neibours,int neibour_count)
{
  int ii,jj,n_vertices_count,ll ;
  int n_solvent = 0 ;
  double volume = 0.0 ; 
  SaveAtomVertices(1,atom, NULL, NULL) ;
  
  for (ii=0 ; ii < neibour_count ; ii++) {
    neibour_descr * nd = &neibours[ii] ; 
    if (nd->flag) {
      vertex * n_vertices[MAX_VERTEX] ;
      if (!nd->protein) n_solvent++;
      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  (DEBUG_VOLUME & debugging)
	fprintf(stderr, "Neighbour       %d   \n",ii);
      
      if (n_vertices_count < 3) {
	serious_error(2, "Less than 3 vertices for a neighbour ", 
		      n_vertices_count);
#       ifdef OFF
	STDERR("Less than 3 vertices skipping neighbor...");
#       endif
      } else  { 
	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 ;
	
	SaveAtomVertices(2,NULL,first_vertex, 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);
	    volume += t_volume ;
	    
	    SaveAtomVertices(3,NULL,second_vertex, third_vertex) ;
	    
	    if (debugging & DEBUG_VOLUME)
	      fprintf(stderr,"       %d   %d   %d    %9.5Lf\n", 
		      first_vertex -vertices, second_vertex -vertices ,
		      third_vertex - vertices , t_volume) ;
	  }
	  xnd = next_nd;
	  third_vertex->flag = 1;
	  second_vertex = third_vertex ;
	}
	SaveAtomVertices(4,NULL,first_vertex,third_vertex);
      }
    }
  }
  
  atom->volume = volume ;
  atom->dumm_solvents = n_solvent ;
  SaveAtomVertices(5,atom,NULL,NULL);
}

/* [ 3 ] Called from [ 2 ] .

   Calculate all the vertices for the
   polyhedron around the atom.  This is a modified version of
   calculate vertices that stores the maximum distance to an atom in
   the atom->surface field */

int N_calculate_vertices 
(atom_record *atom, neibour_descr *n_atoms,int n_count)
{
  int ii;
  double ax = atom->x;
  double ay = atom->y;
  double az = atom->z;

  double max_distance = 0.0 ;
  vertices_count = 0 ;
  	
  
  generate_initial_dummy_vertices(atom, n_atoms,n_count) ;
  for (ii = 0 ; ii < n_count ; ii++) {
    process_neibour(atom,&n_atoms[ii] ) ;
    if (vertices_count == -1)
      return -1 ;
  }
  
  for (ii = 0 ; ii < vertices_count ; ii++) {
    vertex *ve = &vertices[ii] ;
    double deltax = ve->x - ax ;
    double deltay = ve->y - ay ;
    double deltaz = ve->z - az ;
    double dsq = deltaz * deltaz + deltay * deltay + deltax * deltax ;
    
    if (dsq > max_distance)
      max_distance = dsq ; 
    if (DEBUG_VOLUME & debugging)
      fprintf(stderr, "%2d   %2d   %2d   %2d  %9.5Lf  %9.5Lf  %9.5Lf\n",
	      ii, ve->a1 - n_atoms, ve->a2 - n_atoms,
	      ve->a3 - n_atoms,
	      ve->x,ve->y,ve->z) ;
    ve->a1->flag = 1;
    ve->a2->flag = 1;
    ve->a3->flag = 1;
  }
  if (max_distance > volume_max_distance) {
    soft_warning(ARG_TYPE_ATOM, "Furthest vertex too far ", atom);
    return -1 ;
  }
  
  atom->surface = sqrt(max_distance) ;
  return vertices_count;
}

/* [ 2 ] Called from [ 1 ]

   Calculate the volume of one particular atom.
   This is a modified version that calls my routines. 
*/

void N_calculate_atom_volume(atom_record *atom, atom_record **neibours, int neibour_count)
{
  neibour_descr n_atoms[MAX_NBR_DESCRIPTOR] ; 
  vertex  vertices[MAX_VERTEX] ;
  int vertex_count ;
  
  neibour_count 
    = fill_atoms_for_volume(atom , neibours, neibour_count,n_atoms) ;
  
  vertex_count = N_calculate_vertices(atom,n_atoms, neibour_count);
  
  if (vertex_count < 0) {
    atom->volume = -1.0;
    atom->dumm_solvents = -1 ;
    return ;
  }
  
  N_calculate_volume(atom,n_atoms,neibour_count) ;
}

/* [ 1 ] called 1st .

   Sets up the neighbor calculation and calculates the volume of a
   whole file . */

void N_Calculate_Volumes 
(file_records *file, int TheMethod, double VolMaxDist, double DistCheck)
{
# ifdef OFF
  volume_max_dist = 200 ;   /* This setting will very liberally calculate */
# endif
  if (distance_check < 0.0)
    distance_check = 8.0 ;  /* This variable is defined in neibour-cube.c */
  the_method = TheMethod ;
  volume_max_distance = VolMaxDist;
  distance_check = DistCheck ;
  setup_file_records (file , 0.0,1) ; /* Used to set up the cube. */
  setup_minmax_xyzr(file,1) ; /* Used to set up the cube. */
  setup_neibouring_cube(file, distance_check) ;
  PrintVolumeParameters();
  traverse_file_atom_records (file, 
			      call_with_neibours, N_calculate_atom_volume);
}


/*
-------------------------------------------------------------------------------
*/

void compact_print_file(file_records * f )
{
  int ii;
  for (ii=0; ii<f->atomnum; ii++) {
    atom_record * a = f->atoms[ii] ;
    atom_definition * ad = a->definition ;
    fprintf (stderr,"[%4.4d %s%3.3d:%4.4s ",
	    a->number,a->resname,a->resnum,a->name) ;
    fprintf (stderr,"%2.2d %4s %4.1f]",
	    ad->type_index,atom_type_name(ad->type_index),a->x);

    if (ii % 3 == 0) fprintf(stderr,"\n");
  }
}





/* [ 2b ] Called from [ 1b ]

   Calculate the volume of one particular atom.
   This is a modified version that calls my routines. 
   Plus it calls LoopOverNeighbors().
*/

void NN_calculate_atom_volume(atom_record *atom, atom_record **neibours, int neibour_count)
{
  neibour_descr n_atoms[MAX_NBR_DESCRIPTOR] ; 
  vertex  vertices[2000] ;
  int vertex_count ;
  
  if (debugging & DEBUG_VOLUME) {
    STDERR("NN_calculate_atom_volume():");
    fprintf(stderr,"Central atom is %s%3.3d:%4.4s\n",
	    atom->resname,atom->resnum,atom->name,neibour_count) ;
  }

  assert(neibour_count < MAX_NBR_DESCRIPTOR);
  neibour_count 
    = fill_atoms_for_volume(atom , neibours, neibour_count,n_atoms) ;

  
  if (debugging & DEBUG_VOLUME) {
    STDERR("About to call N_calculate_vertices()");
    compact_print_file(global_file);
    fprintf(stderr,"with atom= [%4.4d]%s%3.3d:%4.4s\n",
	    atom->number,atom->resname,atom->resnum,atom->name) ;
    ITOE(neibour_count);
  }
  vertex_count = N_calculate_vertices(atom,n_atoms, neibour_count);
  
  if (vertex_count < 0) {
    atom->volume = -1.0;
    atom->dumm_solvents = -1 ;
    return ;
  }
  if (CalcVHookFcn) {
    (*CalcVHookFcn)(atom,n_atoms,neibour_count) ;
  } else {
    N_calculate_volume(atom,n_atoms,neibour_count) ;
    LoopOverNeighbors(atom,n_atoms,neibour_count);
  }
}

/* [ 1b ] called 1st :
   Sets up the neighbor calculation and calculates the volume of a whole file . 
   Calls LoopOverNeighbors() & SaveAtomVertices() .
*/
void NN_Calculate_Volumes 
(file_records *file, int TheMethod, double VolMaxDist, double DistCheck)
{
# ifdef OFF
  volume_max_dist = 200 ;   /* This setting will very liberally calculate */
# endif
  if (distance_check < 0.0)
    distance_check = 8.0 ;  /* This variable is defined in neibour-cube.c */
  the_method = TheMethod ;
  volume_max_distance = VolMaxDist;
  distance_check = DistCheck ;
  setup_file_records (file , 0.0,1) ; /* Used to set up the cube. */
  setup_minmax_xyzr(file,1) ; /* Used to set up the cube. */

  setup_neibouring_cube(file, distance_check) ;
  PrintVolumeParameters();
  global_file = file ;
  traverse_file_atom_records (file, 
			      call_with_neibours, NN_calculate_atom_volume);
}