C***********************************************************************
C
C
                            PROGRAM gamma
C
C
C***********************************************************************
C
C
C     Program to generate initial particle conditions.  Current 
C     options include:
C
C                1. Uniform density sphere of fixed particles.
C                2. Equilibrium Plummer sphere.
C                3. 2-D uniform surface density sheet of fixed
C                   particles.
C                4. 2-D exponential surface density disk of
C                   fixed particles.
C                5. Equilibrium LH model.
C                6. Two component equilibrium LH model.
C
C         Version 2: March 19, 1992
C
C
C=======================================================================


        IMPLICIT REAL*8(a-h,o-z)
        PARAMETER(npmax=1024000)

        CHARACTER*7 option
        CHARACTER*8 outfile

        REAL*8 msys,m,msys1,msys2

        DIMENSION x(npmax),y(npmax),z(npmax),vx(npmax),vy(npmax),
     &            vz(npmax),m(npmax)

        CALL inoption(npmax,option,np,msys,rsys,r0,np1,np2,msys1,msys2,
C            --------
     &                rsys1,rsys2,r01,r02,outfile)

        CALL initmass(option,np,msys,m,np1,np2,msys1,msys2)
C            --------

        IF(option.EQ.'plummer') THEN

           CALL initplum(np,r0,rsys,msys,x,y,z,vx,vy,vz)
C               --------
           CALL cmtrans(np,m,msys,vx,vy,vz)
C               -------
        ELSE IF(option.EQ.'sphere') THEN

           CALL initsphr(np,rsys,msys,x,y,z)
C               --------
        ELSE IF(option.EQ.'sheet') THEN

           CALL initshee(np,rsys,msys,x,y)
C               --------
        ELSE IF(option.EQ.'expdisk') THEN

           CALL initdisk(np,rsys,r0,msys,x,y)
C               --------
        ELSE IF(option.EQ.'LH') THEN

           CALL initLH(np,r0,rsys,msys,x,y,z,vx,vy,vz,m)
C               ------
c           CALL cmtrans(np,m,msys,vx,vy,vz)
C               -------

        ELSE

           CALL init2LH(np,r0,rsys,msys,x,y,z,vx,vy,vz,m,np1,np2,msys1,
C               -------
     &                  msys2,rsys1,rsys2,r01,r02)

           CALL cmtrans(np,m,msys,vx,vy,vz)
C               -------

        ENDIF

c        CALL cmtrans(np,m,msys,x,y,z)
C            -------
        CALL outdata(option,outfile,np,m,x,y,z,vx,vy,vz)
C            -------

        write(6,*) ' before outcm '

        CALL outcm(np,m,x,y,z,vx,vy,vz)
C            -----

        CALL EXIT(0)
        END

C***********************************************************************
C
C
        SUBROUTINE inoption(npm,option,np,msys,rsys,r0,np1,np2,msys1,
     &                      msys2,rsys1,rsys2,r01,r02,outfile)
C
C
C***********************************************************************
C
C
C     Subroutine to read in options from the terminal.
C
C
C=======================================================================

        IMPLICIT REAL*8(a-h,o-z)
        CHARACTER*7 option
        CHARACTER*8 outfile
      
        REAL*8 msys
	integer idum

        WRITE(6,*) ' Input option:  sphere-->uniform density sphere '
        WRITE(6,*) '                plummer-->equil. plummer model '
        WRITE(6,*) '                sheet-->2-D uniform density disk '
        WRITE(6,*) '                expdisk-->2-D exponential disk '
        WRITE(6,*) '                LH-->equil. LH model '
        WRITE(6,*) '                twoLH-->2 component LH model '
        READ(5,10) option
 10     FORMAT(1a7)

        IF(option.EQ.'twoLH') THEN

           WRITE(6,*) ' This option is not currently implemented '
           STOP

c           WRITE(6,*) 'Input number of particles in small component'
c           READ(5,*) np1
c
c           WRITE(6,*) ' Input number of particles in large component'
c           READ(5,*) np2
c        
c           np=np1+np2
c
c           IF(np.GT.npm) THEN
c              WRITE(6,*) ' np exceeds npmax '
c              CALL EXIT(0)
c           ENDIF
c
c           WRITE(6,*) ' Input mass and size of small component'
c           READ(5,*) msys1,rsys1
c
c           WRITE(6,*) ' Input mass and size of large component'
c           READ(5,*) msys2,rsys2
c
c           msys=msys1+msys2
c
c           WRITE(6,*) ' Input scale-length of small component '
c           READ(5,*) r01
c
c           WRITE(6,*) ' Input scale-length of large component '
c           READ(5,*) r02

        ELSE
           WRITE(6,*) ' Input number of particles '
           READ(5,*) np
        
           IF(np.GT.npm) THEN
              WRITE(6,*) ' np exceeds npmax '
              CALL EXIT(0)
           ENDIF

           WRITE(6,*) ' Input system mass and size '
           READ(5,*) msys,rsys

           IF(option.EQ.'plummer'.OR.option.EQ.'expdisk'.OR.option
     &        .EQ.'LH') THEN
              WRITE(6,*) ' Input scale-length '
              READ(5,*) r0
           ENDIF
        ENDIF

        WRITE(6,*) ' Input name of file '
        READ(5,20) outfile
 20     FORMAT(1a8)

        WRITE(6,*) ' input iseed '
        READ(5,*) iseed

	iseed = -ABS(iseed)
c        CALL xraninit(iseed)
	idum = ABS(iseed)
        CALL dxrand(iseed)
C            --------

        RETURN
        END

C***********************************************************************
C
C
         SUBROUTINE initmass(option,np,msys,m,np1,np2,msys1,msys2)
C
C
C***********************************************************************
C
C
C     Subroutine to initialize masses of individual particles.
C
C
C=======================================================================

        IMPLICIT REAL*8(a-h,o-z)
        CHARACTER*7 option

        REAL*8 m,msys,msys1,msys2

        DIMENSION m(np)

        IF(option.EQ.'twoLH') THEN
           DO 10 i=1,np1
              m(i)=msys1/np1
 10        CONTINUE
           DO 20 i=np1+1,np1+np2
              m(i)=msys2/np2
 20        CONTINUE
        ELSE
           DO 30 i=1,np
              m(i)=msys/np
 30        CONTINUE
        ENDIF

        RETURN
        END

C***********************************************************************
C
C
          SUBROUTINE initplum(np,r0,rsys,msys,x,y,z,vx,vy,vz)
C
C
C***********************************************************************
C
C
C     Subroutine to initialize particle positions and velocities
C     according to an equilibrium Plummer model.
C
C
C=======================================================================

        IMPLICIT REAL*8(a-h,o-z)
        REAL*8 msys
        REAL*8 dxrand
	integer idum

        DIMENSION x(np),y(np),z(np),vx(np),vy(np),vz(np)

        ntot=0
        sigma0=SQRT(msys/(6.*r0))
        twopi=2.*3.141592654
        const=(45./32.)/SQRT(5.)*(9./(14.*SQRT(5.)))**3.5

 10     xr=dxrand(idum)*rsys
        xv=dxrand(idum)*SQRT(12.)*sigma0
        e=.5*xv*xv-6.*sigma0*sigma0/SQRT(1.+xr*xr/(r0*r0))
        esig02=e/sigma0**2
        IF(esig02.GE.0..OR.esig02.LE.-6.) GO TO 10

        p=dxrand(idum)

        pn=const*(xr*xr/(r0*r0))*(xv*xv/(sigma0*sigma0))*(-esig02)**3.5

        IF(p.LE.pn) THEN
           cth=2.*(dxrand(idum)-.5)
           sth=SQRT(1.-cth*cth)
           signs=2.*(dxrand(idum)-.5)
           sth=signs*sth/ABS(signs)
           phi=twopi*dxrand(idum)
           cthv=2.*(dxrand(idum)-.5)
           sthv=SQRT(1.-cthv*cthv)
           signs=2.*(dxrand(idum)-.5)
           sthv=signs*sthv/ABS(signs)
           phiv=twopi*dxrand(idum)
           ntot=ntot+1
           x(ntot)=xr*sth*COS(phi)
           y(ntot)=xr*sth*SIN(phi)
           z(ntot)=xr*cth
           vx(ntot)=xv*sthv*COS(phiv)
           vy(ntot)=xv*sthv*SIN(phiv)
           vz(ntot)=xv*cthv
        ELSE
           GO TO 10
        ENDIF

        IF(ntot.LT.np) GO TO 10

        RETURN
        END

C***********************************************************************
C
C
          SUBROUTINE initLH(np,r0,rsys,msys,x,y,z,vx,vy,vz,m)
C
C
C***********************************************************************
C
C
C     Subroutine to initialize particle positions and velocities
C     according to an equilibrium gamma model.  The variable const
C     represents the maximum value of r**2 * v**2 * f(q), which
C     occurs at r/r_0 = 0.63798179 (g=1) and 
c     v/v_g = v_circ = r*d\Phi/dr.
C
C	const0 = 2.88385786d-3
C	const1 = 1.20223157581242
C	const2 = 5.2573d-3
C
C=======================================================================

        IMPLICIT REAL*8(a-h,o-z)
        REAL*8 msys,m,mk
        REAL*8 dxrand
	integer idum

        DIMENSION x(np),y(np),z(np),vx(np),vy(np),vz(np),m(np)
	
	write(6,*)'which gamma model (0,1,2)'
	read(5,*)gamma

	
C         if(gamma.eq.1) mk=msys*rsys*rsys/(rsys+r0)**2/np
	  mk = msys*((rsys/(rsys+r0))**(3.-gamma))/np

        DO 5 k=1,np
           m(k)=mk
 5      CONTINUE

        ntot=0

        one=1.0
        twopi=2.*4.*ATAN(one)
        pi=4.*ATAN(one)
        const1=1.20223157581242
        const0=0.2884d-2
        const2=0.52573d-2
        c0max = 0
        c1max = 0
        c2max = 0

	if(gamma.eq.0) emin = -msys/(2.*r0)
	if(gamma.eq.1) emin = -msys/r0
	if(gamma.eq.2) emin = -1.0d1

	if(gamma.eq.0) vemax = SQRT(-2.*emin)
	if(gamma.eq.1) vemax = SQRT(-2.*emin)
	if(gamma.eq.2) vemax = SQRT(-2.*emin)

 10     xr=dxrand(idum)*rsys
C        xv=dxrand(idum)*SQRT(2.*msys/r0)
        xv=dxrand(idum)*vemax

C	gamma dependent
        if(gamma.eq.0) e=0.5*xv*xv-0.5*msys*(1.-(xr/(xr+r0))**2)/r0   
        if(gamma.eq.1) e=0.5*xv*xv-msys/(xr+r0)
        if(gamma.eq.2) e=0.5*xv*xv+msys*LOG(xr/(xr+r0))/r0

        IF(e.GT.0.0.OR.e.LT.emin) GO TO 10
        q=SQRT(-r0*e/msys)

        p=dxrand(idum)

c
	if (gamma.eq.0) then
         q = q*q
         dum1 = SQRT(2.D0*q/(1.D0-2.D0*q))
         dum2 = SQRT(1.D0/(1.D0-2.D0*q))
         fq= 3.D0/(2.D0*pi**3)*(
     &    (3.D0-4.D0*q)*dum1*dum2 - 3.D0*LOG(dum1+dum2)   )
         pn=(xr*xr/(r0*r0))*(xv*xv/(msys/r0))*fq/const0
         c0tmp=(xr*xr/(r0*r0))*(xv*xv/(msys/r0))*fq
         c0max = MAX(c0max,c0tmp)
	endif
c
	if (gamma.eq.1) then
C=======================================================================
        fq=(3.*ASIN(q)+q*SQRT(1.-q*q)*(1.-2.*q*q)*(8.*q**4-8.*q**2-3.))
     &     /(1.-q**2)**2.5
         pn=(xr*xr/(r0*r0))*(xv*xv/(msys/r0))*fq/const1
         c1tmp=(xr*xr/(r0*r0))*(xv*xv/(msys/r0))*fq
         c1max = MAX(c1max,c1tmp)
	endif
c
	if (gamma.eq.2) then
         q = q*q
         dum1 = SQRT(q)
         dum2 = SQRT(2.D0)*dum1
         fq=1.D0/SQRT(8.D0*pi**5)*(
     &    EXP(2.D0*q)*erf(dum2)/SQRT(2.D0) - EXP(q)*erf(dum1) +
     &    SQRT(2.D0/pi)*( dawson(dum2) - SQRT(2.D0)*dawson(dum1) ) )
         pn=(xr*xr/(r0*r0))*(xv*xv/(msys/r0))*fq/const2
         c2tmp=(xr*xr/(r0*r0))*(xv*xv/(msys/r0))*fq
         c2max = MAX(c2max,c2tmp)
	endif
c

        IF(pn.GT.1.0) THEN
           WRITE(6,*) ' pn error in initLH ',pn
        ENDIF

        IF(p.LE.pn) THEN
           cth=2.*(dxrand(idum)-.5)
           sth=SQRT(1.-cth*cth)
           signs=2.*(dxrand(idum)-.5)
           sth=signs*sth/ABS(signs)
           phi=twopi*dxrand(idum)
           cthv=2.*(dxrand(idum)-.5)
           sthv=SQRT(1.-cthv*cthv)
           signs=2.*(dxrand(idum)-.5)
           sthv=signs*sthv/ABS(signs)
           phiv=twopi*dxrand(idum)
           ntot=ntot+1
           x(ntot)=xr*sth*COS(phi)
           y(ntot)=xr*sth*SIN(phi)
           z(ntot)=xr*cth
           vx(ntot)=xv*sthv*COS(phiv)
           vy(ntot)=xv*sthv*SIN(phiv)
           vz(ntot)=xv*cthv
        ELSE
           GO TO 10
        ENDIF

        IF(ntot.LT.np) GO TO 10

	write(6,*)0,c0max,const0
	write(6,*)1,c1max,const1
	write(6,*)2,c2max,const2

        RETURN
        END

C***********************************************************************
C
C
           SUBROUTINE init2LH(np,r0,rsys,msys,x,y,z,vx,vy,vz,m,np1,np2,
     &                        msys1,msys2,rsys1,rsys2,r01,r02)
C
C
C***********************************************************************
C
C
C     Subroutine to initialize particle positions and velocities
C     according to a two component equilibrium LH model.  The variable 
C     const represents the maximum value of r**2 * v**2 * f(q), which
C     occurs at r/r_0 = 0.63798179 and v/v_g = v_circ.
C
C
C=======================================================================

        IMPLICIT REAL*8(a-h,o-z)
        REAL*8 msys,m,msys1,msys2,mk1,mk2
        REAL*8 dxrand
	integer idum

        DIMENSION x(np),y(np),z(np),vx(np),vy(np),vz(np),m(np)

        mk1=msys1*rsys1*rsys1/(rsys1+r01)**2/np1
        mk2=msys2*rsys2*rsys2/(rsys2+r02)**2/np2

        DO 5 k=1,np1
           m(k)=mk1
 5      CONTINUE

        DO 7 k=np1+1,np1+np2
           m(k)=mk2
 7      CONTINUE

        one=1.0
        twopi=2.*4.*ATAN(one)
        const=1.20223157581242

C  Initilize spatial coordinates of component 1.
C  ---------------------------------------------

        write(6,*) ' initializing spatial coords of component 1 '

        ntot=0

 10     xr=dxrand(idum)*rsys1

        p=dxrand(idum)

        pn=27.0*r01*r01*xr/(4.0*(xr+r01)**3)

        IF(pn.GT.1.0) THEN
           WRITE(6,*) ' pn error in initLH '
        ENDIF

        IF(p.LE.pn) THEN
           cth=2.*(dxrand(idum)-.5)
           sth=SQRT(1.-cth*cth)
           signs=2.*(dxrand(idum)-.5)
           sth=signs*sth/ABS(signs)
           phi=twopi*dxrand(idum)
           ntot=ntot+1
           x(ntot)=xr*sth*COS(phi)
           y(ntot)=xr*sth*SIN(phi)
           z(ntot)=xr*cth
           IF(MOD(ntot,1000).EQ.0) WRITE(6,*) ntot
        ELSE
           GO TO 10
        ENDIF

        IF(ntot.LT.np1) GO TO 10

C  Initilize spatial coordinates of component 2.
C  ---------------------------------------------

        write(6,*) ' initializing spatial coords of component 2 '

        ntot=np1

 20     xr=dxrand(idum)*rsys2

        p=dxrand(idum)

        pn=27.0*r02*r02*xr/(4.0*(xr+r02)**3)

        IF(pn.GT.1.0) THEN
           WRITE(6,*) ' pn error in initLH '
        ENDIF

        IF(p.LE.pn) THEN
           cth=2.*(dxrand(idum)-.5)
           sth=SQRT(1.-cth*cth)
           signs=2.*(dxrand(idum)-.5)
           sth=signs*sth/ABS(signs)
           phi=twopi*dxrand(idum)
           ntot=ntot+1
           x(ntot)=xr*sth*COS(phi)
           y(ntot)=xr*sth*SIN(phi)
           z(ntot)=xr*cth
           IF(MOD(ntot,1000).EQ.0) WRITE(6,*) ntot
        ELSE
           GO TO 20
        ENDIF

        IF(ntot.LT.np1+np2) GO TO 20

C  Initilize velocity coordinates of component 1.
C  ----------------------------------------------

        write(6,*) ' initializing velocity coords of component 1 '

        ntot=0

 30     xr=SQRT(x(ntot+1)**2+y(ntot+1)**2+z(ntot+1)**2)
        xv=dxrand(idum)*SQRT(2.*msys1/r01)
        e=0.5*xv*xv-msys1/(xr+r01)
        IF(e.GT.0.0.OR.e.LT.-msys1/r01) GO TO 30
        q=SQRT(-r01*e/msys1)

        p=dxrand(idum)

        fq=(3.*ASIN(q)+q*SQRT(1.-q*q)*(1.-2.*q*q)*(8.*q**4-8.*q**2-3.))/
     &     (1.-q**2)**2.5
        pn=(xr*xr/(r01*r01))*(xv*xv/(msys1/r01))*fq/const

        IF(pn.GT.1.0) THEN
           WRITE(6,*) ' pn error in initLH '
        ENDIF

        IF(p.LE.pn) THEN
           cthv=2.*(dxrand(idum)-.5)
           sthv=SQRT(1.-cthv*cthv)
           signs=2.*(dxrand(idum)-.5)
           sthv=signs*sthv/ABS(signs)
           phiv=twopi*dxrand(idum)
           ntot=ntot+1
           vx(ntot)=xv*sthv*COS(phiv)
           vy(ntot)=xv*sthv*SIN(phiv)
           vz(ntot)=xv*cthv
           IF(MOD(ntot,1000).EQ.0) WRITE(6,*) ntot
        ELSE
           GO TO 30
        ENDIF

        IF(ntot.LT.np1) GO TO 30

C  Initilize velocity coordinates of component 2.
C  ----------------------------------------------

        write(6,*) ' initializing velocity coords of component 2 '

        ntot=np1

 40     xr=SQRT(x(ntot+1)**2+y(ntot+1)**2+z(ntot+1)**2)
        xv=dxrand(idum)*SQRT(2.*msys2/r02)
        e=0.5*xv*xv-msys2/(xr+r02)
        IF(e.GT.0.0.OR.e.LT.-msys2/r02) GO TO 40
        q=SQRT(-r02*e/msys2)

        p=dxrand(idum)

        fq=(3.*ASIN(q)+q*SQRT(1.-q*q)*(1.-2.*q*q)*(8.*q**4-8.*q**2-3.))/
     &     (1.-q**2)**2.5
        pn=(xr*xr/(r02*r02))*(xv*xv/(msys2/r02))*fq/const

        IF(pn.GT.1.0) THEN
           WRITE(6,*) ' pn error in initLH '
        ENDIF

        IF(p.LE.pn) THEN
           cthv=2.*(dxrand(idum)-.5)
           sthv=SQRT(1.-cthv*cthv)
           signs=2.*(dxrand(idum)-.5)
           sthv=signs*sthv/ABS(signs)
           phiv=twopi*dxrand(idum)
           ntot=ntot+1
           vx(ntot)=xv*sthv*COS(phiv)
           vy(ntot)=xv*sthv*SIN(phiv)
           vz(ntot)=xv*cthv
           IF(MOD(ntot,1000).EQ.0) WRITE(6,*) ntot
        ELSE
           GO TO 40
        ENDIF

        IF(ntot.LT.np1+np2) GO TO 40

        RETURN
        END

C***********************************************************************
C
C
                  SUBROUTINE cmtrans(np,m,msys,x,y,z)
C
C
C***********************************************************************
C
C
C     Subroutine to transform to center of mass coordinates.
C
C
C=======================================================================

        IMPLICIT REAL*8(a-h,o-z)
        REAL*8 m,msys,mtot

        DIMENSION x(np),y(np),z(np),m(np)

        xcm=0.
        ycm=0.
        zcm=0.
        mtot=0.

        DO 10 i=1,np
           mtot=mtot+m(i)
           xcm=xcm+m(i)*x(i)
           ycm=ycm+m(i)*y(i)
           zcm=zcm+m(i)*z(i)
 10     CONTINUE

        xcm=xcm/mtot
        ycm=ycm/mtot
        zcm=zcm/mtot

        DO 20 i=1,np
           x(i)=x(i)-xcm
           y(i)=y(i)-ycm
           z(i)=z(i)-zcm
 20     CONTINUE

        RETURN
        END

C***********************************************************************
C
C
                SUBROUTINE initsphr(np,rsys,msys,x,y,z)
C
C
C***********************************************************************
C
C
C     Subroutine to generate an initial distribution of particles
C     corresponding to a uniform density sphere.
C
C
C=======================================================================

        IMPLICIT REAL*8(a-h,o-z)
        REAL*8 msys
        REAL*8 dxrand
	integer idum

        DIMENSION x(np),y(np),z(np)

        ntot=0

 10     xr=(.5-dxrand(idum))*2.*rsys
        yr=(.5-dxrand(idum))*2.*rsys
        zr=(.5-dxrand(idum))*2.*rsys

        IF(xr*xr+yr*yr+zr*zr.GT.rsys*rsys) GO TO 10

        ntot=ntot+1
        x(ntot)=xr
        y(ntot)=yr
        z(ntot)=zr

        IF(ntot.LT.np) GO TO 10

        RETURN
        END

C***********************************************************************
C
C
                  SUBROUTINE initshee(np,rsys,msys,x,y)
C
C
C***********************************************************************
C
C
C     Subroutine to initial particle positions in a uniform 2-D
C     disk.
C
C
C=======================================================================

        IMPLICIT REAL*8(a-h,o-z)
        REAL*8 msys
        REAL*8 dxrand
	integer idum

        DIMENSION x(np),y(np)

        ntot=0

 10     xr=(.5-dxrand(idum))*2.*rsys
        yr=(.5-dxrand(idum))*2.*rsys

        IF(xr*xr+yr*yr.GT.rsys*rsys) GO TO 10

        ntot=ntot+1
        x(ntot)=xr
        y(ntot)=yr

        IF(ntot.LT.np) GO TO 10

        RETURN
        END

C***********************************************************************
C
C
                  SUBROUTINE initdisk(np,rsys,r0,msys,x,y)
C
C
C***********************************************************************
C
C
C     Subroutine to initialize particle positions in an exponential
C     disk with scale-length r0.
C
C
C=======================================================================

        IMPLICIT REAL*8(a-h,o-z)
        REAL*8 msys
        REAL*8 dxrand
	integer idum

        DIMENSION x(np),y(np)

        ntot=0
        sigma0=msys/(2.*3.141592654*r0*r0)

 10     xr=(.5-dxrand(idum))*2.*rsys
        yr=(.5-dxrand(idum))*2.*rsys

        IF(xr*xr+yr*yr.GT.rsys*rsys) GO TO 10

        r=SQRT(xr*xr+yr*yr)
        p=dxrand(idum)
        fn=EXP(-r/r0)
       
        IF(p.LE.fn) THEN
           ntot=ntot+1
           x(ntot)=xr
           y(ntot)=yr
        ELSE
           GO TO 10
        ENDIF

        IF(ntot.LT.np) GO TO 10

        RETURN
        END

C***********************************************************************
C
C
          SUBROUTINE outdata(option,outfile,np,m,x,y,z,vx,vy,vz)
C
C
C***********************************************************************
C
C
C     Subroutine to write particle data to an output file.
C
C
C=======================================================================

        IMPLICIT REAL*8(a-h,o-z)
        CHARACTER*7 option
        CHARACTER*8 outfile
        REAL*8 m

        DIMENSION x(np),y(np),z(np),vx(np),vy(np),vz(np),m(np)

        OPEN(UNIT=10,FILE=outfile,STATUS='NEW')

        time=0.
        ndim=2
        IF(option.EQ.'plummer'.OR.option.EQ.'sphere') ndim=3

        WRITE(10,111) np,time

        DO 100 i=1,np
           WRITE(10,112) m(i),x(i),y(i),z(i),vx(i),vy(i),vz(i)
 100    CONTINUE

 111    FORMAT(1x,1i8,1pe14.6)
 112    FORMAT(1x,7(1pe14.6))

        CLOSE(UNIT=10)

        RETURN
        END

C***********************************************************************
C
C
                  SUBROUTINE outcm(np,m,x,y,z,vx,vy,vz)
C
C
C***********************************************************************
C
C
C     Subroutine to print out center of mass coordinates.
C
C
C=======================================================================

        IMPLICIT REAL*8(a-h,o-z)
        REAL*8 m,mtot

        DIMENSION x(np),y(np),z(np),m(np),vx(np),vy(np),vz(np)

        xcm=0.
        ycm=0.
        zcm=0.
        vxcm=0.
        vycm=0.
        vzcm=0.
        mtot=0.

        DO 10 i=1,np
           mtot=mtot+m(i)
           xcm=xcm+m(i)*x(i)
           ycm=ycm+m(i)*y(i)
           zcm=zcm+m(i)*z(i)
           vxcm=vxcm+m(i)*vx(i)
           vycm=vycm+m(i)*vy(i)
           vzcm=vzcm+m(i)*vz(i)
 10     CONTINUE

        xcm=xcm/mtot
        ycm=ycm/mtot
        zcm=zcm/mtot
        vxcm=vxcm/mtot
        vycm=vycm/mtot
        vzcm=vzcm/mtot

        write(6,*) ' xcm = ',xcm,ycm,zcm
        write(6,*) ' vxcm = ',vxcm,vycm,vzcm

        write(6,*) ' '

        RETURN
        END

CXXX NR random routine
C
CXXX
    
      FUNCTION DXRAND(IDUM)
      INTEGER IDUM,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB
      REAL*8 DXRAND,AM,ATAB
      PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1,
     *    IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,
     *    IR2=3791,NTAB=32.,ATAB=NTAB-1.)
      INTEGER IDUM2,J,K,IV(NTAB)
      REAL*8 Y
      DATA IDUM2/123456789/, IV/NTAB*0/, Y/0./
      SAVE Y,IV,IDUM2
      IF (IDUM.LE.0) THEN
        IDUM=MAX(-IDUM,1)
        IDUM2=IDUM
        DO 11 J=NTAB,1,-1
          K=IDUM/IQ1
          IDUM=IA1*(IDUM-K*IQ1)-K*IR1
          IF (IDUM.LT.0) IDUM=IDUM+IM1
          IV(J)=IDUM
11      CONTINUE
        Y=AM*IV(1)
      ENDIF
1     CONTINUE
        K=IDUM/IQ1
        IDUM=IA1*(IDUM-K*IQ1)-K*IR1
        IF (IDUM.LT.0) IDUM=IDUM+IM1
        K=IDUM2/IQ2
        IDUM2=IA2*(IDUM2-K*IQ2)-K*IR2
        IF (IDUM2.LT.0) IDUM2=IDUM2+IM2
        J=1+INT(ATAB*Y)
        K=IV(J)-IDUM2
        IV(J)=IDUM
        IF(K.LT.1)K=K+IMM1
        Y=AM*K
        DXRAND=Y
      IF (DXRAND.EQ.0..OR.DXRAND.EQ.1.) GOTO 1
      RETURN
      END

C***********************************************************************
      FUNCTION dawson(x)
      DOUBLE PRECISION dawson,x
C***********************************************************************
      INTEGER Nmax
      DOUBLE PRECISION H,A1,A2,A3
C*      PARAMETER (Nmax=6,H=0.4E0,A1=2.E0/3.E0,A2=0.4E0,A3=2.E0/7.E0)
C*       Not sure how we should change Nmax for double precision.
      PARAMETER (Nmax=12,H=0.4D0,A1=2.D0/3.D0,A2=0.4D0,A3=2.D0/7.D0)
      INTEGER i,init,n0
      DOUBLE PRECISION d1,d2,e1,e2,sum,x2,xp,xx,c(Nmax)
      SAVE init,c
      DATA init/0/
      IF (init.EQ.0) THEN
        init=1
        DO 11 i=1,Nmax
          c(i)=EXP(-((2.D0*DBLE(i)-1.D0)*H)**2)
   11   CONTINUE
      ENDIF
      IF (ABS(x).LT.0.2D0) THEN
        x2=x**2
        dawson=x*(1.D0-A1*x2*(1.D0-A2*x2*(1.D0-A3*x2)))
      ELSE
        xx=ABS(x)
        n0=2*NINT(0.5D0*xx/H)
        xp=xx-DBLE(n0)*H
        e1=EXP(2.D0*xp*H)
        e2=e1**2
        d1=DBLE(n0+1)
        d2=d1-2.D0
        sum=0.D0
        DO 12 i=1,Nmax
          sum=sum+c(i)*(e1/d1+1.D0/(d2*e1))
          d1=d1+2.D0
          d2=d2-2.D0
          e1=e2*e1
   12   CONTINUE
        dawson=0.5641895835*SIGN(EXP(-xp**2),x)*sum
      ENDIF
      RETURN
      END




C***********************************************************************
      FUNCTION erf(x)
      DOUBLE PRECISION erf,x
CU      USES gammp
C***********************************************************************
      DOUBLE PRECISION gammp
      IF (x.LT.0.) THEN
        erf=-gammp(0.5D0,x**2)
      ELSE
        erf=gammp(0.5D0,x**2)
      ENDIF
      RETURN
      END
C***********************************************************************
      FUNCTION gammln(xx)
      DOUBLE PRECISION gammln,xx
C***********************************************************************
      INTEGER j
      DOUBLE PRECISION ser,stp,tmp,x,y,cof(6)
      SAVE cof,stp
      DATA cof,stp/76.18009172947146d0,-86.50532032941677d0,
     $  24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2,
     $  -.5395239384953d-5,2.5066282746310005d0/
      x=xx
      y=x
      tmp=x+5.5d0
      tmp=(x+0.5d0)*LOG(tmp)-tmp
      ser=1.000000000190015d0
      DO 11 j=1,6
        y=y+1.d0
        ser=ser+cof(j)/y
   11 CONTINUE
      gammln=tmp+LOG(stp*ser/x)
      RETURN
      END
C***********************************************************************
      FUNCTION gammp(a,x)
      DOUBLE PRECISION a,gammp,x
CU      USES gcf,gser
C***********************************************************************
      DOUBLE PRECISION gammcf,gamser,gln
      IF (x.LT.0.D0 .OR. a.LE.0.D0) THEN
c        CALL termin('gammp: bad arguments',-1)
      ENDIF
      IF (x.LT.a+1.D0) THEN
        CALL gser(gamser,a,x,gln)
        gammp=gamser
      ELSE
        CALL gcf(gammcf,a,x,gln)
        gammp=1.D0-gammcf
      ENDIF
      RETURN
      END
C***********************************************************************
      SUBROUTINE gcf(gammcf,a,x,gln)
      DOUBLE PRECISION a,gammcf,gln,x
CU      USES gammln
C***********************************************************************
      INTEGER Itmax
      DOUBLE PRECISION Eps,Fpmin
C*       PARAMETER (Itmax=100,Eps=3.E-7,Fpmin=1.E-30)
      PARAMETER (Itmax=150,Eps=1.D-15,Fpmin=1.D-44)
      INTEGER i
      DOUBLE PRECISION an,b,c,d,del,h,gammln
      gln=gammln(a)
      b=x+1.D0-a
      c=1.D0/Fpmin
      d=1.D0/b
      h=d
      DO 11 i=1,Itmax
        an=-i*(i-a)
        b=b+2.D0
        d=an*d+b
        IF (ABS(d).LT.Fpmin) d=Fpmin
        c=b+an/c
        IF (ABS(c).LT.Fpmin) c=Fpmin
        d=1.D0/d
        del=d*c
        h=h*del
        IF (ABS(del-1.D0).LT.Eps) GOTO 1
   11 CONTINUE
c      CALL termin('gcf: a too large, Itmax too small',-2)
    1 gammcf=EXP(-x+a*LOG(x)-gln)*h
      RETURN
      END
C***********************************************************************
      SUBROUTINE gser(gamser,a,x,gln)
      DOUBLE PRECISION a,gamser,gln,x
CU      USES gammln
C***********************************************************************
      INTEGER Itmax
      DOUBLE PRECISION Eps
C*       PARAMETER (Itmax=100,Eps=3.E-7)
      PARAMETER (Itmax=150,Eps=1.D-15)
      INTEGER n
      DOUBLE PRECISION ap,del,sum,gammln
      gln=gammln(a)
      IF (x.LE.0.D0) THEN
        IF (x.LT.0.D0) THEN
c          CALL termin('gser: x < 0',-1)
        ENDIF
        gamser=0.D0
        RETURN
      ENDIF
      ap=a
      sum=1.D0/a
      del=sum
      DO 11 n=1,Itmax
        ap=ap+1.D0
        del=del*x/ap
        sum=sum+del
        IF (ABS(del).LT.ABS(sum)*Eps) GOTO 1
   11 CONTINUE
c      CALL termin('gser: a too large, Itmax too small',-2)
    1 gamser=sum*EXP(-x+a*LOG(x)-gln)
      RETURN
      END



