SUBROUTINE getaphpot (rho, theta, xchi, pot)
USE Numeric_Kinds_Module
USE CommonInfo_Module
USE Jacobi_Module
Use Masses_Module
!            P U R P O S E   O F    S U B R O U T I N E
! This routine calculates interparticle distances in Bohr and THEN
! calls surface to determine the interaction potential in Hartree's.
!            I N P U T    A R G U M E N T S
!      rho    hyperradius in Bohr.          (0 <  rho  < infinity)
!      theta  hyperangle theta in radians.  (0 < theta < pi/2)
!      xchi   hyperangle chi in radians.    (0 <  chi  < 2pi)
!            O U T P U T    A R G U M E N T S
!      pot    potential energy V(rho,theta,chi) in Hartree atomic units.
! <><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>

IMPLICIT NONE

!            I N T E G E R S

INTEGER :: k,i,ip,im
!            R E A L S
REAL (dp) :: pot, rho, theta, xchi
REAL (dp) :: r(3), s(3), sbig(3), xm(3), dscale(3), chij(3,3), chi3(3)
REAL (dp) :: mtot, mu
!-----------------------------------------------------------------------
!  Calculate the interpartical distances and store them in the arrary r.
!-----------------------------------------------------------------------
   mtot = amass + bmass + cmass
   mu = sqrt (amass * bmass * cmass / mtot)
   dscale (1) = sqrt ( (amass / mu) * (1.d0 - (amass / mtot) ) )
   dscale (2) = sqrt ( (bmass / mu) * (1.d0 - (bmass / mtot) ) )
   dscale (3) = sqrt ( (cmass / mu) * (1.d0 - (cmass / mtot) ) )

   xm(1)=amass
   xm(2)=bmass
   xm(3)=cmass
   DO i=1, 3
     chij(i,i)=0.0d0
   ENDDO
   DO i=1, 3
     ip=i+1
     im=i-1
     IF(i .eq. 3)ip=1
     IF(i .eq. 1)im=3
     chij(i,ip)=acos(-mu/(dscale(i)*dscale(ip)*xm(im)))
     chij(ip,i)=-chij(i,ip)
   ENDDO

   chi3 (1) = xchi
   chi3 (2) = chi3 (1) - chij (2, 1)
   chi3 (3) = chi3 (1) - chij (3, 1)

 DO k = 1, 3
   s (k) = rho * sqrt (0.5d0 - 0.5d0 * sin (theta) * cos (2.d0 * chi3 (k) ) )
   sbig (k) = rho * sqrt (0.5d0 + 0.5d0 * sin (theta) * cos (2.d0  * chi3 (k) ) )

   r (k) = dscale (k) * s (k)

!   fac (k) = sin (theta) * cos (2.d0 * chi3 (k) )
!   den (k) = sqrt (1.d0 - (fac (k) * fac (k) ) )
!   cthetabig (k) = (sin (theta) * sin (2.d0 * chi3 (k) ) )   / den (k)
!   sthetabig (k) = sqrt (1.d0 - (cthetabig (k) * cthetabig (k) ) )
!   sb (k) = (sthetabig (k) * s (k) * sin (2.d0 * chi3 (k) ) ) / qbig
!   carrier = ( (s (k) * sin (2.d0 * chi3 (k) ) ) * cthetabig ( k) )
!   cb (k) = (carrier + (sbig (k) * cos (2.d0 * chi3 (k) ) ) )  / qbig

 END DO

!-----------------------------------------------------------------------
! This CALL returns the interaction potential in Hartree atomic units.
!-----------------------------------------------------------------------

CALL surface (pot, r)
RETURN
END SUBROUTINE getaphpot