C
      REAL*4 FUNCTION ECONF(AN,NA)
C
C Function ECONF: Calculate the energy per charge from a set of
C                 charges with locations specified by AN(1,2,..NA)
C                 at fixed distance from the origin. 
C
C Update history:
C    V1.0    Original version                          8-Oct-1992
C    V2.0    Internal calculations R*8                23-Oct-1992
C
      IMPLICIT NONE
C
C Arguments:
C    AN = R*4(?) = vector of angles
C    NA = I*4    = number of angles in vector
C
      REAL*4    AN(1)
      INTEGER*4 NA
C
C Parameters:
C
      INCLUDE 'PARAMETERS.INC'
C
C Local variables:
C   ESUM = R*8 = summed energy
C   I1   = I*4 = counter
C   I2   = I*4 = counter
C   IA   = I*4 = counter over pair particles
C   IB   = I*4 = counter over pair particles
C   NN   = I*4 = number of particles
C   XA   = R*8 = location of particle a
C   YA   = R*8 = location of particle a
C   ZA   = R*8 = location of particle a
C   XB   = R*8 = location of particle b
C   YB   = R*8 = location of particle b
C   ZB   = R*8 = location of particle b
C   RAB  = R*8 = separation of a and b
C   XAB  = R*8 = separation of a and b
C   YAB  = R*8 = separation of a and b
C   ZAB  = R*8 = separation of a and b
C
      INTEGER*4 I1,I2,IA,IB,NN
      REAL*8    XA,XB,YA,YB,ZA,ZB
      REAL*8    ESUM,RAB,XAB,YAB,ZAB
C
C Initialise
C
      NN  =NA/2
      ESUM=0.0D0
C
C For each pair of charges, sum their mutual energy in units of
C    q^2/4 pi epsilon_0
C
      DO IA=1,NN
C
C Location of first particle in pair (Cartesian coordinates)
C
        I1=2*IA-1
        I2=I1+1
        XA=DBLE(SIND(AN(I1))*COSD(AN(I2)))
        YA=DBLE(SIND(AN(I1))*SIND(AN(I2)))
        ZA=DBLE(COSD(AN(I1)))
        DO IB=(IA+1),NN
C
C Location of second particle in pair (Cartesian coordinates)
C
          I1=2*IB-1
          I2=I1+1
          XB=DBLE(SIND(AN(I1))*COSD(AN(I2)))
          YB=DBLE(SIND(AN(I1))*SIND(AN(I2)))
          ZB=DBLE(COSD(AN(I1)))
C 
C Particle separations
C
          XAB=XA-XB
          YAB=YA-YB
          ZAB=ZA-ZB
          RAB=DSQRT(XAB*XAB+YAB*YAB+ZAB*ZAB)
C
C Guard against zero separations
C
          IF (SNGL(RAB) .LT. RABMIN) THEN
            ESUM=ESUM+DBLE(EABMAX)
          ELSE
            ESUM=ESUM+(1.0D0/RAB)
          ENDIF
C
        ENDDO
      ENDDO
C
C Normalize by number of particles
C
      ECONF=ESUM/DBLE(FLOAT(NN))
C
C Return result
C
      RETURN
C
      END