SUBROUTINE CollinearFixedRho
USE Numeric_Kinds_Module
USE Numbers_Module
USE Parms_Module
USE FileUnits_Module
USE Masses_Module
USE MassFactor2_Module
USE Narran_Module
USE BasisRegions_Module

IMPLICIT NONE
INTEGER, PARAMETER :: nc=600
INTEGER k, ic, isector
REAL(KIND=WP_Kind) pot, rho, theta, xchi, sin, cos, dc
REAL(KIND=WP_Kind), ALLOCATABLE :: PotChi(:)
REAL(KIND=WP_Kind), PARAMETER :: rhovec(4)=(/2.d0, 3.13d0, 6.d0, 20.d0/)
INTRINSIC sqrt, sin, cos
EXTERNAL surface

DO ISector=1,4
  Basis_Dist(ISector) = rhovec(Isector)
ENDDO

OPEN(Unit=639,File=TRIM(OutDIR)//'GraphicsOut/CollinearFixedRho.csv',Form='FORMATTED')
WRITE(639,'(A,",",300(A4,1PE10.2,","))')"xchi", ("Rho=",Basis_Dist(ISector),ISector=1,4) !NSectors)
theta=pi/two
dc=(chij(1,1)-chij(2,1))/(nc-1)/2
ALLOCATE(PotChi(1:NSectors))
!-----------------------------------------------------------------------
!  Calculate the interparticle distances and store them in the array r.
!-----------------------------------------------------------------------
DO ic=0,nc
  DO ISector=1,4 !NSectors
    rho=Basis_Dist(ISector)
    xchi=ic*dc
    chi3(1)=xchi
    chi3(2)=chi3(1)-chij(2,1)
    chi3(3)=chi3(1)-chij(3,1)
    DO k=1,narran
       s(k)=rho*sqrt(Half-Half*sin(theta)*cos(Two*chi3(k)))
       r(k)=d(k)*s(k)
    ENDDO

!-----------------------------------------------------------------------
! This CALL returns the interaction potential in Hartree atomic units.
!-----------------------------------------------------------------------
    CALL surface (pot, r)
    potchi(ISector)=pot
  ENDDO
  !WRITE(639,'(300(1PE14.6,","))')xchi*180.d0/Pi, (potchi(ISector),ISector=1,4) !NSectors)  
  WRITE(639,'(300(1PE14.6,","))')xchi, (potchi(ISector),ISector=1,4) !NSectors)
ENDDO
Deallocate(PotChi)
CLOSE(Unit=639)

STOP "CollinearFixedRho"
ENDSUBROUTINE CollinearFixedRho