# include "NRVects.h"

/*  ------                                      */
int fixwat_(float waters[], long * nwat, long * ifix)
/*  ------                                      */

/* 
Takes as input an array of water positions (water) and pointers to two
integers, which describe the number of waters and the index of the
water to fix. Note,

* the first water index in waters is assumed to be 1 and the last is 
  nwat. 

* nwat and ifix relate to the number of waters, not
  the number of atoms, so if nwat=10, waters would contain
  30 atoms and 90 floats. 

* waters[] has the following format:

0     1     2     3     4     5     6     7     8     9    
O-1-x O-1-y O-1-z H-1-x H-1-y H-1-z h-1-x h-1-y h-1-z O-2-x 

10    11    12    13    14    15    16    17
O-2-y O-2-z H-2-x H-2-y H-2-z h-2-x h-2-y h-2-z 


The goal of this routine is to orient all the waters so that water
ifix has the following orientation:

   ^
   |
  -y   (z into page)
   |

   O  ----x--->
 /   \
H     h

*/

#define BAD_INDEX (fprintf(stderr,"%s [%d]: BAD_INDEX\n",__FILE__,__LINE__) ) 
#define IDX(nn,at,xyz) \
 ((nn-1)*9 + \
 (at  == 'O' ? 0 : (at  == 'H' ? 3 : (at  == 'h' ? 6 : BAD_INDEX ))) + \
 ( xyz-1))

{
  AllocateNRVect3 h, x, y, z, o;
  NRMatrix33P  R = dmatrix(1,3,1,3); 
  int i = *ifix ;
  
  o[1] = waters[IDX(i,'O',1)] ;  
  o[2] = waters[IDX(i,'O',2)] ;    
  o[3] = waters[IDX(i,'O',3)] ;    
  
  h[1] = waters[IDX(i,'h',1)] - o[1] ;
  h[2] = waters[IDX(i,'h',2)] - o[2] ;  
  h[3] = waters[IDX(i,'h',3)] - o[3] ; 

  y[1] = 0.5 * (waters[IDX(i,'h',1)] + waters[IDX(i,'H',1)] ) - o[1] ;
  y[2] = 0.5 * (waters[IDX(i,'h',2)] + waters[IDX(i,'H',2)] ) - o[2] ;
  y[3] = 0.5 * (waters[IDX(i,'h',3)] + waters[IDX(i,'H',3)] ) - o[3] ;

  NR_Normalize(h);
  NR_Normalize(y);
  NR_perp_dir_1eq2x3(z,h,y);
  NR_perp_dir_1eq2x3(x,z,y);
  
  NR_MakeMat1WithVect234(R,x,y,z);
  Transpose(R);
  
  for (i=0; i<*nwat; i++) {
    h[1] = waters[3*i]   - o[1] ;
    h[2] = waters[3*i+1] - o[2] ;
    h[3] = waters[3*i+2] - o[3] ; 
    VectRotate_1by2 (h, R) ;
    waters[3*i  ] = h[1] ;
    waters[3*i+1] = h[2] ;
    waters[3*i+2] = h[3] ;
  }
}