#include "NNvolume.h"

double volume_max_distance = 64.0 ; 
int the_method = 2 ;
static double the_edge = 2.8 ;
vertex vertices [MAX_VERTEX] ;
int vertices_count  = 0;
FILE * plot_file = NULL ;

/* This plots a line in frodo vector list format */
void plot_a_line(vertex *v1, vertex *v2)
{
  fprintf(plot_file, "%10.5Lf%10.5Lf%10.5Lf    0\n", v1->x,v1->y,v1->z);
  fprintf(plot_file, "%10.5Lf%10.5Lf%10.5Lf    1\n", v2->x,v2->y,v2->z);
}

/* The generation of a vertex from the three neighbours.  The vertex
   is at the point where :

 a1->dx * x + a1->dy * y + a1->dz * z + a1->plane_constant = 0 
 a2->dx * x + a2->dy * y + a2->dz * z + a2->plane_constant = 0 
 a3->dx * x + a3->dy * y + a3->dz * z + a3->plane_constant = 0 

The solution to this equation is :

       |a1->dx a1->dy a1->dz|
B = det|a2->dx a2->dy a2->dz|
       |a3->dx a3->dy a3->dz|

        |a1->plane_constant a1->dy a1->dz|    / 
x = -det|a2->plane_constant a2->dy a2->dz|   /  B
        |a3->plane_constant a3->dy a3->dz|  /   

        |a1->dx a1->plane_constant a1->dz|    / 
y = -det|a2->dx a2->plane_constant a2->dz|   /  B
        |a3->dx a3->plane_constant a3->dz|  /   

        |a1->dx a1->dy a1->plane_constant|    / 
z = -det|a2->dx a2->dy a2->plane_constant|   /  B
        |a3->dx a3->dy a3->plane_constant|  / 


*/


#define DET(x1,x2,x3,y1,y2,y3,z1,z2,z3) ( x1 * (y2 * z3 - y3 * z2) - \
                                          x2 * (y1 * z3 - y3 * z1) + \
					  x3 * (y1 * z2 - y2 * z1))


void generate_vertex (neibour_descr *a1, neibour_descr *a2,
		      neibour_descr *a3)
{
  vertex *ve = & vertices[vertices_count ++ ] ;
  
  if (debugging & DEBUG_VOLUME) {
    fprintf(stderr,LF_STRING_14,a1->dx,a1->dy,a1->dz,a1->plane_constant);
    fprintf(stderr,LF_STRING_14,a2->dx,a2->dy,a2->dz,a2->plane_constant);
    fprintf(stderr,LF_STRING_14,a3->dx,a3->dy,a3->dz,a3->plane_constant);
  }
  
  ve->a1 = a1 ;
  ve->a2 = a2 ;
  ve->a3 = a3 ;
  
  
  {
    double B  = DET (a1->dx, a1->dy, a1->dz,
		     a2->dx, a2->dy, a2->dz,
		     a3->dx, a3->dy, a3->dz) ;
    
    if (fabs(B) < 0.000001) {
      fprintf(stderr,
	      "generate_vertex(): fabs( B= %lf ) < 0.00000001, so skip vol. calc.\n",
	      B);
      write_pdb_record(stderr,a1->atom,0);
      write_pdb_record(stderr,a2->atom,0);
      write_pdb_record(stderr,a3->atom,0);
#     ifdef OFF
      print_an_atom(a1->atom,stderr);
      print_an_atom(a2->atom,stderr);
      print_an_atom(a3->atom,stderr);
      fprintf(stderr," %Lf \n%Lf  %Lf  %Lf  \n %Lf  %Lf  %Lf  \n%Lf  %Lf  %Lf  \n", 
	      B,a1->dx, a1->dy, a1->dz,
	      a2->dx, a2->dy, a2->dz,
	      a3->dx, a3->dy, a3->dz);

      STDERR("generate_vertex(): skipping vertex...");
      vertices_count-- ;
#     endif      
      vertices_count = -1 ;
      return ;
    }
    B = 0.0 - B ;
    
    ve->x = DET (a1->plane_constant, a1->dy, a1->dz,
		 a2->plane_constant, a2->dy, a2->dz,
		 a3->plane_constant, a3->dy, a3->dz) / B ;
    
    ve->y = DET (a1->dx, a1->plane_constant, a1->dz,
		 a2->dx, a2->plane_constant, a2->dz,
		 a3->dx, a3->plane_constant, a3->dz) / B ;

    ve->z = DET (a1->dx ,a1->dy, a1->plane_constant,
		 a2->dx, a2->dy ,a2->plane_constant,
		 a3->dx, a3->dy, a3->plane_constant) / B ;
  }
  if (debugging & DEBUG_VOLUME)
    fprintf(stderr,LF_STRING_15,ve->x , ve->y,ve->z) ;
  
}

/* To delete a vertex: if it is not the last one, copy the last one on
   top of it. */

void delete_vertex (int vert_number)
{
  vertices_count-- ;

  if ( vertices_count != vert_number) {
    vertex * ve = & vertices[vert_number] ;
    vertex * ve1= & vertices[vertices_count] ;
    
    ve->x = ve1-> x;
    ve->y = ve1-> y;
    ve->z = ve1-> z;
    ve->a1 = ve1->a1 ;
    ve->a2 = ve1->a2 ;
    ve->a3 = ve1->a3 ;
  }
}

/* Process a neighbour. First check which side of the plane the origin
   atom is, positive or negative.  Then go throgh the vertices, and
   for each one check if it is on the same side as the atom. If it is
   not, we have to get rid of the vertex, and add new vertices. These
   are added for every pair of atoms participating in forming any of
   the removed vertices. */

void process_neibour(atom_record * atom, neibour_descr * nd)
{
  double atom_off = - (nd->dx * atom->x +
		       nd->dy * atom->y +
		       nd->dz * atom->z ) ;
  int positive = (atom_off > nd->plane_constant) ;
  int add_count = 0 ;
  int ii ;
  neibour_descr *addition_vector[MAX_NBR_DESCRIPTOR] ;
  

  for(ii = 0 ; ii < vertices_count ; ii++ ) {
    vertex *ve = &vertices[ii] ;
    double vertex_off = - (nd->dx * ve->x +
			   nd->dy * ve->y +
			   nd->dz * ve->z );
    
    if (positive ? vertex_off < nd->plane_constant : 
	vertex_off > nd->plane_constant ) {
      add_count = mark_for_addition(ve->a1,ve->a2,addition_vector,add_count) ;
      add_count = mark_for_addition(ve->a1,ve->a3,addition_vector,add_count) ;
      add_count = mark_for_addition(ve->a2,ve->a3,addition_vector,add_count) ;
      ve->a1 = NULL;
    }
  }
  
  for(ii = vertices_count - 1 ; ii >= 0 ; ii--) {
    if (vertices[ii].a1 == NULL)
      delete_vertex (ii) ;
  }
  for (ii = 0 ; ii < add_count; ii += 2) {
    generate_vertex(nd, addition_vector[ii], addition_vector[ii+1]) ;
    if (vertices_count == -1)
      return  ;
  }
}

 
/* calculate the plane_constant for a plane cutting the line between
   the origin atom and the neibour at a point determined by ratio. The
   plane equation is:

nd->dx * x + nd->dy * y + nd->dz * z  + plane_constant = 0 

point_x,point_y and point_z are the coordinates of the point where the
plane cuts the line. */

void calculate_plane_constant
(atom_record * atom, neibour_descr * nd,double ratio)
{
  double point_x = nd->dx * ratio +atom->x  ;
  double point_y = nd->dy * ratio +atom->y  ;
  double point_z = nd->dz * ratio +atom->z  ;
	
  nd->plane_constant = - (point_x * nd->dx +
			  point_y * nd->dy +
			  point_z * nd->dz ) ;
  if (DEBUG_VOLUME & debugging && nd->atom != NULL) {
    if (nd->atom == NULL)
      fprintf(stderr,"Dummy ")  ;
    else 
      print_an_atom(nd->atom, stderr);
    fprintf(stderr,"  %9.5Lf  %9.5Lf  %9.5Lf  d %9.5Lf dr %9.5Lf rat %9.5Lf pc %9.5Lf\n",
	    point_x ,  point_y,point_z,nd->distance,nd->distance * ratio, ratio,nd->plane_constant) ;
  }
}

/* Generate 4 dummy atoms, and from them 4 dummy vertices. This forms
   an initial polyhedron, which is then cut down to the final
   polyhedron. */

void generate_initial_dummy_vertices 
(atom_record * atom, neibour_descr *n_atoms,int n_count)
{
  neibour_descr *dumm_nd1 = &n_atoms[n_count] ;
  neibour_descr *dumm_nd2 = &n_atoms[n_count + 1] ;
  neibour_descr *dumm_nd3 = &n_atoms[n_count + 2] ;
  neibour_descr *dumm_nd4 = &n_atoms[n_count + 3] ;

  dumm_nd1->atom = NULL ;
  dumm_nd2->atom = NULL ;
  dumm_nd3->atom = NULL ;
  dumm_nd4->atom = NULL ;
  dumm_nd1->dx =  10.0 ; dumm_nd1->dy =  10.0 ; dumm_nd1->dz =  10.0;
  dumm_nd2->dx =  10.0 ; dumm_nd2->dy = -10.0 ; dumm_nd2->dz = -10.0;
  dumm_nd3->dx = -10.0 ; dumm_nd3->dy =  10.0 ; dumm_nd3->dz = -10.0;
  dumm_nd4->dx =  -10.0 ; dumm_nd4->dy = -10.0 ; dumm_nd4->dz =  10.0;

  calculate_plane_constant(atom,dumm_nd1,0.5);
  calculate_plane_constant(atom,dumm_nd2,0.5);
  calculate_plane_constant(atom,dumm_nd3,0.5);
  calculate_plane_constant(atom,dumm_nd4,0.5);

  generate_vertex (dumm_nd1, dumm_nd2, dumm_nd3);
  generate_vertex (dumm_nd2, dumm_nd3, dumm_nd4);
  generate_vertex (dumm_nd3, dumm_nd4, dumm_nd1);
  generate_vertex (dumm_nd4, dumm_nd1, dumm_nd2);
}



/* 
   Check that the four dummies create in the above routine 
   have been eliminated from the calculation .
*/

int check_dummy_were_eliminated ( neibour_descr *n_atoms,int n_count)
{
  if (n_atoms[n_count].flag || n_atoms[n_count + 1 ].flag || 
      n_atoms[n_count + 2].flag || n_atoms[n_count + 3 ].flag )
    return 0;
  else 
    return 1 ;
}

/*  
   Calculate all the vertices for the polyhedron around the atom.
*/

int 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 ;
  }
  
  return vertices_count;

# ifdef OFF
  if (check_dummy_were_eliminated(n_atoms,n_count))
      return vertices_count ;
  else 
    return -1;
# endif
}


/* Check if the pair a1 and a2 is already in the addition_vector. If
   it isn't, add it. Else if it is already there, remove the other
   pair. */

int mark_for_addition (neibour_descr *a1, neibour_descr *a2,
			neibour_descr **addition_vector, int add_count)
{
  int ii ;

  if (a1 > a2) {
    neibour_descr * te = a1 ; 
    a1 = a2; 
    a2 = te; 
  }
  for (ii = 0 ; ii< add_count ; ii += 2)
    if (addition_vector[ii] == a1 && addition_vector[ii+1] == a2) {
      add_count -= 2 ;
      addition_vector[ii] =  addition_vector[add_count];
      addition_vector[ii+1] =  addition_vector[add_count+1];
      return add_count ;
    }
  
  addition_vector[add_count] = a1;
  addition_vector[add_count +1] = a2;
  return add_count + 2;
}

/* Calculate the volume of the tetrahedron defined by a central atom
   and 3 vertices. The volume is :

                                            | b - a |
     (b - a) . (c -a ) x (d - a) / 6  =  det| c - a | / 6
                                            | d - a |

where b,c and d are the vertices vectors, and a is the atom vector */

double tet_volume (atom_record *atom, vertex *v1, vertex *v2, vertex *v3)
{
  double ax = atom->x;
  double ay = atom->y;
  double az = atom->z;
	
  double v1x = v1->x - ax ;
  double v1y = v1->y - ay ;
  double v1z = v1->z - az ; 
  double v2x = v2->x - ax ;
  double v2y = v2->y - ay ;
  double v2z = v2->z - az ;
  double v3x = v3->x - ax ;
  double v3y = v3->y - ay ;
  double v3z = v3->z - az ;

  return fabs((DET(v1x,v1y,v1z,v2x,v2y,v2z,v3x,v3y,v3z)) / 6) ;
}


/* 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_vertices, 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 fficient way is to calculate only the area, and
   multiply by the distance of the plane to to the origin atom (and
   divide by 3).  */

void 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 ; 
  if (plot_file != NULL) {
    print_an_atom(atom,plot_file) ;
    fputc('\n', plot_file) ;
  }

  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);
	
      } 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 ;
	
	if (plot_file != NULL) 
	  plot_a_line(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 ;
	    
	    if (plot_file != NULL)
	      plot_a_line(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 ;
	}
	if (plot_file != NULL)
	  plot_a_line(first_vertex, third_vertex) ;
      }
    }
  }
  atom->volume = volume ;
  atom->dumm_solvents = n_solvent ;
}

/* calculate the volume of one particular atom */

void 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 = calculate_vertices(atom,n_atoms, neibour_count);
  
  if (vertex_count < 0) {
    atom->volume = -1.0;
    atom->dumm_solvents = -1 ;
    return ;
  }
  
  calculate_volume(atom,n_atoms,neibour_count) ;
}

/* 
   initially get rid of all the non-zero dummy waters, because that
   was what Richards, the great, does. */

int water_model_mode = 1 ;

void check_shell_resnum (atom_record *atom, void *dummy, int index)
{
   int resnum = atom->resnum ;

   switch(water_model_mode)  {
   case 1 : if (resnum != 0)
		SWITCH_ON(atom,IGNORE_ATOM) ; break; 
   case 2 : if ( resnum == -2)
		SWITCH_ON(atom,IGNORE_ATOM) ; break; 
  }
}

/* Calculate the volume of all the atoms (these that are marked so) */

void calculate_volumes (file_records *file,file_records *shell, double probe)
{
  file_records *xshell = shell ;
  if (distance_check < 0.0)
    distance_check = 8.0 ;
  setup_file_records (file , 0.0,1) ;
  setup_minmax_xyzr(file,1) ;
  the_edge = probe + probe ;
  if (shell != NULL) {
    if (shell == file)
      shell = generate_shell(file,probe ) ;
    traverse_file_atom_records (shell, check_shell_resnum, NULL);	
  }
  setup_neibouring_cube(file, distance_check) ;
  if (shell != NULL)
    add_file_to_neibouring_cube(shell) ;
  traverse_file_atom_records (file, call_with_neibours,  calculate_atom_volume);
  if (shell != xshell)
    free_file_records(shell) ;
}

/* The following function is defined for qsort. */

# ifdef OFF
int compare_neibour_descriptors(neibour_descr *a, neibour_descr *b)
# endif

int compare_neibour_descriptors(const void *aa, const void *bb)
{
  neibour_descr *a = (neibour_descr *) aa;
  neibour_descr *b = (neibour_descr *) bb;
  
  if (a->distance > b->distance)
    return 1;
  return -1;
}

int fill_atoms_for_volume(atom_record * atom, atom_record **neibours, 
			     int neibour_count, neibour_descr *descr)
{
  int ii,jj = 0 ;
  double x = atom->x ;
  double y = atom->y ;
  double z = atom->z ;
	
  double xt,xb,yt,yb,zb,zt;     /* This junk for compatibility with
				   Richards, the Great. */
  
  ii = rint((x - minimum_x) / the_edge);
  xb = minimum_x + (ii  - 0.5) * the_edge - 0.1 ;
  xt = xb + 0.2 + 2 * the_edge ;
  ii = rint((y - minimum_y) / the_edge);
  yb = minimum_y + (ii  - 0.5) * the_edge - 0.1 ;
  yt = yb + 0.2 + 2 * the_edge ;
  ii = rint((z - minimum_z) / the_edge);
  zb = minimum_z + (ii  - 0.5) * the_edge - 0.1 ;
  zt = zb + 0.2 + 2 * the_edge ;

  if (debugging & DEBUG_VOLUME)
    fprintf(stderr,
	    "limits x : %9.5Lf   %9.5Lf  \n       y : %9.5Lf   %9.5Lf  \n       z : %9.5Lf   %9.5Lf  \n", xb, xt,yb,yt,zb,zt);
  for (ii =0 ; ii < neibour_count ; ii++) {
    neibour_descr * nd = &descr[jj] ; 
    atom_record *neibour =  neibours[ii] ;
    
    if (neibour->atomtype != RECORD_DUMMY)
      nd->protein = 1;	
    else {
      if (neibour->x > xt || neibour->x < xb ||
	  neibour->y > yt || neibour->y < yb ||
	  neibour->z > zt || neibour->z < zb )
	continue ;
      
      nd->protein = 0 ;
    }
    jj++ ;
    
    nd->atom = neibour ;
    nd->flag = 0 ;
    
    { 
      double dx = neibour->x - x ;
      double dy = neibour->y - y ;
      double dz = neibour->z - z ;
      double dsq = dx * dx + dy * dy + dz * dz ;
      
      nd->dx = dx ; 
      nd->dy = dy;
      nd->dz = dz;
      nd->dsq = dsq;
      nd->distance = sqrt(dsq) ;
      calculate_plane_constant(atom, nd, NNcalculate_ratio(atom,nd)) ;
    }
  }
  qsort(descr, jj , sizeof(neibour_descr),compare_neibour_descriptors) ;
  if (DEBUG_VOLUME & debugging)
    for (ii = 0 ; ii< jj ; ii++) {
      neibour_descr *nd = &descr[ii] ;
      print_an_atom(nd->atom,stderr);
      fprintf(stderr,"%3d %9.5Lf %9.5Lf %9.5Lf %9.5Lf %9.5Lf\n", ii,
	      nd->distance,nd->dx,nd->dy,nd->dz,nd->plane_constant) ;
    }
  return jj ;
}