!****************************************************************
! This module is a component of VINE, copyright by the authors, *
! as listed below. Please see the COPYING and LICENSE files     *
! in the parent directory for additional information.           *
! This notification is present as an attempt to follow the      *
! recommendations for licensing published on the FSF page.      * 
!***********************************************************************
! Various common functions used in the tree walk/build etc             *
!                                                                      *
! hsum2       : Squared sum of the node sizes.                         *
! r2conv      : Squared multipole convergence radius                   *
! r6conv      : sixth power of multipole convergence radius            *
!               (for Gadget MAC--differently defined than r2conv)      *
! dxi,dyi,dzi : Distance vectors, appropriately redefined              *
!               to account for periodic boundaries, if                 *
!               neccessary (Easier than ghost particles)               *
! rr2neig     : squared distance between two tree points,              *
!               possibly accounting for periodic boundaries.           *
!               (i.e. rr2neig uses dxi, dyi and dzi)                   *
! rr2         : squared distance between two tree points               *
!               not accounting for periodic boundaries.                *
!                                                                      *
! Why am I using statement functions? They are archaic aren't they?    *
!   Well yes they are, but the reason is that as far as I can tell,    *
!   compilers don't always inline normal functions, even if required   *
!   to do so by using a command line option to do it. This means       *
!   there can be quite a bit of overhead for calling them. Very bad    *
!   for this set of functions, as you can probably guess.              *
!   Since a statement function is not much different than fortran's    *
!   version of a macro expansion it doesn't have the inlining issues,  *
!   and will therefore be faster. Perhaps this argument will requre    *
!   modification in f90/f95/f2000                                      *
!                                                                      *
! NB: Do NOT use dxi,dyi or dzi for anything except                    *
! the distance calculations for neighbors. They can                    *
! give wrong results for the gravity calculation                       *
! if periodic boundaries are used, because the                         *
! distances may result in particles appearing outside                  *
! of the simulation box. The Ewald correction already                  *
! accounts for this--it separates the components from                  *
! inside the box, from the duplications outside the box.               *
! Thus you would 1) be doing a double counting error if                *
! you accidentally used a node from outside the box                    *
! and 2) miscalculate the force from the node itself,                  *
! since it was supposed to be in the box by assumption                 *
! in the Ewald analysis.                                               *
!                                                                      *
! Written by A. Nelson                                                 *
!                                                                      *
!***********************************************************************

      integer iii,jjj
      double precision dxi,dyi,rr2neig,rr2
#ifndef PERIODIC
!
! Free boundary distance definitions. No difference between
!                               rr2neig and rr2 required.
!
      dxi    (iii,jjj)=treept(1,iii)-treept(1,jjj)
      dyi    (iii,jjj)=treept(2,iii)-treept(2,jjj)
#if NDIM==2
      rr2neig(iii,jjj)= dxi(iii,jjj)**2 + dyi(iii,jjj)**2
      rr2    (iii,jjj)= dxi(iii,jjj)**2 + dyi(iii,jjj)**2
#else
      double precision dzi
      dzi    (iii,jjj)=treept(3,iii)-treept(3,jjj)
      rr2neig(iii,jjj)= dxi(iii,jjj)**2 + dyi(iii,jjj)**2 
     &                                  + dzi(iii,jjj)**2
      rr2    (iii,jjj)= dxi(iii,jjj)**2 + dyi(iii,jjj)**2
     &                                  + dzi(iii,jjj)**2
#endif


#else
!***********************************************************************
! For neighbor calculations with periodic boundaries, redefine         *
! the distance to be the closest of the distances to the               *
! node itself or its reproduction in one of the bounding               *
! periodic boxes that define its duplicates.                           *
!                                                                      *
!      Four cases:                                                     *
!          -rbox   < xi - xj < -rbox/2     then add rbox               *
!          -rbox/2 < xi - xj <  0          then add nothing            *
!           0      < xi - xj <  rbox/2     then add nothing            *
!           rbox/2 < xi - xj <  rbox       then subtract rbox          *
!                                                                      *
!      anint function calculates the appropriate if case.              *
!        It rounds distance/rbox to -1,0 or 1 as appropriate,          *
!        then we use that value to multiply times rbox and             *
!        add or subract the result                                     *
!                                                                      *
! NOTE: This definition of separation distance implicitly assumes      *
! that the box in which all particles live is centered on the origin.  *
! If this is not true, you will need to change the definition by       *
! subtracting off the actual center of the box as follows:             *
!                                                                      *
! dx = xi - xj - box_origin - boxlen*anint( (xi-xj)/boxlen)            *
!                                                                      *
!***********************************************************************
      dxi(iii,jjj)=treept(1,iii) - treept(1,jjj) - 
     &         boxsize(1)*anint((treept(1,iii)-treept(1,jjj))*boxinv(1))

      dyi(iii,jjj)=treept(2,iii) - treept(2,jjj) - 
     &         boxsize(2)*anint((treept(2,iii)-treept(2,jjj))*boxinv(2))

#if NDIM==2
      rr2neig(iii,jjj)= dxi(iii,jjj)**2 + dyi(iii,jjj)**2
      rr2    (iii,jjj)= (treept(1,iii)-treept(1,jjj))**2 +
     &                  (treept(2,iii)-treept(2,jjj))**2 
#else
      double precision dzi
      dzi(iii,jjj)=treept(3,iii) - treept(3,jjj) - 
     &         boxsize(3)*anint((treept(3,iii)-treept(3,jjj))*boxinv(3))

      rr2neig(iii,jjj)= dxi(iii,jjj)**2 + dyi(iii,jjj)**2 
     &                                  + dzi(iii,jjj)**2
      rr2    (iii,jjj)= (treept(1,iii)-treept(1,jjj))**2 +
     &                  (treept(2,iii)-treept(2,jjj))**2 +
     &                  (treept(3,iii)-treept(3,jjj))**2 
#endif

#endif

      double precision hsum2,r2conv,r6conv
!
! squared smoothing length
! If you change the definition of what a neighbor is
! then you need to change this definition as well, in order
! to find all the neighbors. For example, the set up now
! defines neighbors as those within 2h of each other. If you 
! want a different definition (say within 2.5h) then multiply
! this hsum2 definition by a factor 'hfac', defined as (5/4)**2=25d0/16d0
!
      hsum2(iii,jjj)= (treept(ihnode,iii)+treept(ihnode,jjj))**2
!
! squared multipole convergence distance
!   size of prospective multipole node i added to
!   the (conservative) physical size of node j.
!
      r2conv(iii,jjj)= (treept(iconvrg,iii)+treept(ihnode,jjj))**2
!
! Sixth power of multipole convergence radius for Gadget MAC
! see Springel etal New Astr. 6 (2001), equation 18.
!
      r6conv(iii,jjj)= rr2(iii,jjj)**3
