#include "NNvolume.h"

/* Calculate the ratio of the distance between the origin atom and its
   neighbour, which "belongs" to the origin atom.

   *********************
   Covalent Connectivity
   *********************

   Right now we are deciding whether two atoms are covalently
   connected based on the distance between them being greater than
   2.0. In contrast, a more reasonable criteria (the one Mike uses,
   for instance) for determing if two atoms are covalently bonded is
   the following:

   Two atoms X are bonded if the distance between them is less than
   1.7 Angstroms.  If one of the X is S (sulfur) use 2 A and if both
   are S use 2.8 A.  If one of the X is H, use 1.2 A.

   ***********
   Definitions
   ***********

   Distance of dividing plane from origin atom         = D1 
   Distance between origin atom and neighbor           = D = nd->distance 
   Ratio returned by routine                           = r
                                               D1      = r * D 
   Radius (either covalent or vdw) of origin atom      = R
   Radius (either covalent or vdw) of neighboring atom = N
    
   **********************
   Methods to calculate r:
   **********************

   --------
   Method 1: Richards Method A: 
   --------  r = 0.5
             The classical Voronoi method using bisection. 

   --------	     
   Method 3:
   --------	     

   eradius = ad->vdw + *probe ; 
   Eradius of origin atom = E
   Eradius of neighboring atom = e

                    D^2 + E^2 - e^2     1     e^2 - E^2
             r  =   ---------------  =  -  -  ---------
                         2 D^2          2       2 D^2

   --------	     
   Method 5: 
   --------	     

            r = 1 - R/D  

   --------
   Method 2: is Richards Method B
   --------

   --------
   Method 4: uses Richards Method B formula for covalently connected atoms (part a)
   --------  for all atoms whether or not they are covalently connected


   ******************
   Richards Method B:
   ******************

           --------------------
        a. Covalently connected: 
           --------------------

	r = R/(R+N) 

	(R/N is what Richards puts in his paper but it is wrong!)
             
C
C----   set up plane based on division by covalent radii.
C
520      F= AOCOV/(AOCOV+RCOV(JID))

              ------------------------
        ** b. Not covalently connected: 
              ------------------------


               D - (R + N)       D + R - N
      r = R  + -----------   =   ---------
                   2                2 D 
          ----------------    
                 D 

C----   set up plane based on division by vanderwaals radii.
      F=(BOX(1,J) + AOVDW - RVDW(JID))/(2.*BOX(1,J))
      CALL BXLOAD(J,F)
      GO TO 540
515      F=(BOX(1,J) + AOVDW - VDWDFT)/(2.*BOX(1,J))
      CALL BXLOAD(J,F)
      GO TO 540

   --------
   Method 6: uses Richards Method B formula for not covalently connected atoms (part b)
   --------  for all atoms whether or not they are covalently connected.
             The VDW radius is always used. This is just a test option to use
             to compare with VDW_fr_vol.f .

   --------
   Method 7: uses normal voronoi bisection (method A) for all waters,
   --------  same as method 4 otherwise 


*/

double CalcRatio7 (atom_record *atom, neibour_descr *nd)
{
  if (WATER_ATOM_P(atom)) {
    OFF_MSG_NTIMES(10,"water atom");
    return 0.5 ;
  } else if (WATER_ATOM_P(nd->atom)) {
    OFF_MSG_NTIMES(10,"water nd");
    return 0.5 ;
  } else {
    double distance = nd->distance ;
    int vdw = distance > 2.0 ;   
    double rad = vdw ? atom->definition->vdw : atom->definition->covalent;
    double nrad = vdw ? nd->atom->definition->vdw : nd->atom->definition->covalent;
    double ratio = (rad / (rad + nrad)) ;
    OFF_MSG_NTIMES(10,"ratio not 0.5");
    return ratio ;
  }
}  

double  NNcalculate_ratio (atom_record *atom, neibour_descr *nd)
{
  if (! nd->protein) {
    STDERR("I'm here ! nd->protein in calculate_ratio()");
    return atom->definition->vdw / nd->distance ;
  } else if (the_method == 7) {
    return CalcRatio7 (atom, nd) ;
  } else if (the_method == 1) {
    return 0.5 ;
  } else if (the_method == 3) {
    double dsq = nd->dsq;
    double eradsq = atom->definition->eradsq ;
    double neradsq = nd->atom->definition->eradsq;
    return (dsq + eradsq - neradsq) / (2 * dsq) ;
  } else if (the_method == 5) {
    return (1 - atom->definition->vdw / nd->distance) ;
  } else if (the_method == 6) {
    static int flag = 1 ;
    double distance = nd->distance ;
    double rad = atom->definition->vdw ;
    double nrad = nd->atom->definition->vdw ;
    double ratio =  (distance + rad - nrad) / (2 * distance) ;
    if (flag) {
      STDERR("calculate_ratio(): Using method 6! This is only for a test.");
      flag = 0;
    }
    return ratio ;
  } else {
    double distance = nd->distance ;
    int vdw = distance > 2.0 ;   
    double rad = vdw ? atom->definition->vdw : atom->definition->covalent;
    double nrad = vdw ? nd->atom->definition->vdw : nd->atom->definition->covalent;
    double ratio ;
    if (the_method == 4 || ! vdw)
      ratio = (rad / (rad + nrad)) ;
    else 
      ratio =  (distance + rad - nrad) / (2 * distance) ;
    return ratio ;

#   ifdef OFF
    if (ratio != 0.5) {
      fprintf(stderr,"ratio= %lf\n",ratio);
    }
      fprintf(stderr,"\ncalculate_ratio(): distance= %lf  vdw= %1d  the_method= %1d\n",
	      distance,vdw,the_method);
      fprintf(stderr,"                 : rad= %lf nrad= %lf  ratio= %lf\n",
	      rad,nrad,ratio);
      fprintf(stderr,"                 : cov= %lf ncov= %lf\n",
	      atom->definition->covalent,nd->atom->definition->covalent);
      fprintf(stderr,"                 : type= %s ntype= %s\n",
	      atom->definition->name,nd->atom->definition->name);
      write_pdb_record (stderr, atom,0);
      write_pdb_record (stderr, nd->atom,0);
    }
#   endif
  }
}