cCode to construct initial conditions for KyotoII
cOutput, on unit 8, is in "astrophysical units" (solar mass, km/s, pc)
      implicit none
      integer i,idum,j,nb,ns
c i,j        loop counter
c idum       argument of random number generator 
c nb         number of binaries
c ns         number of single stars
      real a,alist,amax,amin,b,cosi,ea,eadot,g,gamma(4),e,e2,elist,
     &     emax,inc,kt,loga,m(2),ma,mass,mbin,msun,mm,mu(2),nbu2msun,
     &     nbu2pc,nbu2kps,node,p(3,2),pc,peri,pi,q,qmin,radius(2),ran2,
     &     rg,rking,rop,rrel(3),rsun,rsun2pc,rt,rtnewt,select,x,x1,x2,
     &     xacc,vop,vrel(3)
c a          semi-major axis
c alist      diagnostic list of semi-major axes
c amax,amin  maximum and minimum semi-major axis
c b          semi-minor axis
c cosi       cosine of inclination of binary orbit
c e          eccentricity
c e2         e**2
c ea         eccentric anomaly
c eadot      rate of change of ea
c elist      diagnostic list of eccentricities
c emax       maximum eccentricity (so that a(1-e) > 2(R1+R2))
c g          gravitation constant in astrophysical units
c gamma(4)   parameters of KTG3 mass function 
c inc        inclination of binary orbit
c kt         kT
c loga       log (base 10) of semi-major axis
c m(2)       masses of binary components
c ma         mean anomaly
c mass       mass of single star or binary or binary component
c mbin       total mass of one binary
c mm         mean motion in binary ("n")
c msun       solar mass in g
c mu         fractional mass of binary component
c nbu2kps    factor converting velocity in N-body units to km/s
c nbu2msun   factor converting mass in N-body units to solar masses
c nbu2pc     factor converting length in N-body units to pc
c node       longitude of ascending node of binary orbit
c p          direction cosines converting orbital plane to cluster frame; 
c               cf. Brouwer & Clemence, p.33, with altered notation
c pc         1pc in cm
c peri       argument of pericentre of binary orbit
c pi         pi
c q          binary mass ratio (<1)
c qmin       minimum value of q set by stellar radii
c radius(2)  radii of binary components
c ran2       name of random number generator
c rg         galactocentric radius in pc
c rking      edge radius of initial King model in N-body units (from output 
c               of mkking)
c rop        components of relative position of binary components in orbital 
c               plane
c rrel       component of relative position of binary components in cluster 
c               frame
c rsun       solar radius in cm
c rsun2pc    factor converting solar radius to pc
c rt         tidal radius of cluster in pc (as specified)
c rtnewt     name of Newton-Raphson solver
c select     name of routine for finding median, etc
c x          random variable distributed uniformly on [0,1] 
c x1, x2     range of solution of Kepler's equation
c xacc       precision for solution of Kepler's equation
c vop        components of relative velocity of binary components in orbital 
c vrel       component of relative velocity of binary components in cluster 
c               frame
      double precision r,rbar(3), totmass,v,vbar(3)
c r          component of position vector
c rbar       component of position vector of binary barycentre
c totmass    total cluster mass
c v          component of velocity
c               plane
c vbar       component of velocity of binary barycentre
      parameter (msun = 1.989e33, ns = 12288, nb = 4096, 
cThe following line is non-ansi compliant, and should perhaps be replaced by 
c     &     pc = 3.085678e18, pi = 3.1415926536, rking = 4.37094, 
cThanks to Peter Teuben for this and other comments 13/8/1
     &     pc = 3.085678e18, pi = 4.*atan(1.), rking = 4.37094, 
     &     rsun = 6.9599e10, rt = 27.9)
      dimension mass(ns+2*nb),r(ns+2*nb),rop(2),v(ns+2*nb),vop(2),
     &     alist(nb),elist(nb)
      parameter (g = 6.670e-8*msun*1e-10/pc, nbu2pc = rt/rking, 
     &     rsun2pc = rsun/pc)
      data totmass/0.d0/,gamma/0.19,1.55,0.050,0.6/,idum/1/
c For gamma see KTG, 1993, MN, 262, 575, Table 10
      external funcd
c funcd      returns value and derivative for solution of Kepler's equation
      common ma,e
cFirst, create masses for single stars (solar masses)
      do i = 1,ns
 10      continue
         x = ran2(idum)
         mass(i) = 0.08 + (gamma(1)*x**gamma(2) + gamma(3)*x**gamma(4))*
     &        (1.0-x)**(-0.58)
         if (mass(i).le.0.1.or.mass(i).ge.100) go to 10
         totmass = totmass + mass(i)
      enddo
cSecond, create masses for binary components (solar masses)
      do i = ns+2,ns+2*nb,2
 20      continue
         x = ran2(idum)
cSee HTAP, 2001, MN, 323, 633, eq.(1) with x -> 1-x
         mbin = 0.33*(1.0/(x**0.75 + 0.04*x**0.25) - x**2/1.04)
         if (mbin.le.0.2.or.mbin.ge.200) go to 20
         totmass = totmass + mbin
cSee HTAP, eq.(2), with different maximum mass
         qmin = max(0.1/(mbin - 0.1),(mbin - 100.)/100.)
         q = qmin + ran2(idum)*(1. - qmin)
         mass(i-1) = q*mbin/(1.+q)
         mass(i) = mbin - mass(i-1)
      enddo
cWrite data for single stars, in units of solar mass, km/s, pc
      nbu2kps = sqrt(g*totmass/nbu2pc)
      nbu2msun = totmass
      open (7, file='initial_phase_list')
      open (8, file='ic.dat')      
      do i = 1,ns
         read (7,*) (r(j),j=1,3),(v(j),j=1,3)
         write (8,'(1pg14.6,6g24.16)') mass(i),(nbu2pc*r(j),j=1,3),
     &        (nbu2kps*v(j),j=1,3)
      enddo
cCompute and write data for binaries; calculations in astrophysical units
      kt = (1./(6*(ns + nb)))*nbu2msun*nbu2kps**2
      do i = 1,nb
cDetermine radii
         do j = 1,2
            m(j) = mass(ns + 2*(i-1) + j)
cFormulae from EFT, 1989, ApJ, 347, 998
            if (m(j).lt.1.334) then
               radius(j) = (0.1148*m(j)**1.25 + 0.8604*m(j)**3.25)/
     &              (0.04651 + m(j)**2)
            else
               radius(j) = (1.968*m(j)**2.887 - 0.7388*m(j)**1.679)/
     &              (1.821*m(j)**2.337 - 1)
            endif
            radius(j) = radius(j)*rsun2pc
         enddo
cDetermine semi-major axis, eccentricity, and place on orbit
         amax = g*m(1)*m(2)/(20*kt)
         amin = 2*(radius(1)+radius(2))
         if (amin.ge.amax) then
            write (6,*) 'amin>=amax!'
            stop
         endif
         loga = log10(amin) + ran2(idum)*log10(amax/amin)
         a = 10**loga
         alist(i) = a
         emax = 1.0 - amin/a
         e2 = emax**2*ran2(idum)
         e = sqrt(e2)
         elist(i) = e
         b = a*sqrt(1-e2)
         ma = 2*pi*ran2(idum)
         x1 = 0.
         x2 = 2*pi
         xacc = 1.e-6
         ea = rtnewt(funcd,x1,x2,xacc)
         rop(1) = a*(cos(ea) - e)
         rop(2) = b*sin(ea)
         mm = sqrt(g*(m(1)+m(2))/a**3)
         eadot = mm/(1 - e*cos(ea))
         vop(1) = -a*sin(ea)*eadot
         vop(2) = b*cos(ea)*eadot
cConvert to cluster frame
         cosi = 2*ran2(idum)-1
         inc = acos(cosi)
         node = 2*pi*ran2(idum)
         peri = 2*pi*ran2(idum)
         p(1,1) = cos(peri)*cos(node) - sin(peri)*sin(node)*cosi
         p(2,1) = cos(peri)*sin(node) + sin(peri)*cos(node)*cosi
         p(3,1) = sin(peri)*sin(inc)
         p(1,2) = -sin(peri)*cos(node) - cos(peri)*sin(node)*cosi
         p(2,2) = -sin(peri)*sin(node) + cos(peri)*cos(node)*cosi
         p(3,2) = cos(peri)*sin(inc)
         do j = 1,3
            rrel(j) = p(j,1)*rop(1) + p(j,2)*rop(2)
            vrel(j) = p(j,1)*vop(1) + p(j,2)*vop(2)
         enddo
         read (7,*) (rbar(j),j=1,3),(vbar(j),j=1,3)
         do j = 1,3
            rbar(j) = rbar(j)*nbu2pc
            vbar(j) = vbar(j)*nbu2kps
         enddo            
         mu(1) = m(1)/(m(1)+m(2))
         mu(2) = 1 - mu(1)
         write (8,'(1pg14.6,6g24.16)') mass(ns+2*i-1),
     &        (rbar(j) + mu(2)*rrel(j),j=1,3),
     &        (vbar(j) + mu(2)*vrel(j),j=1,3)
         write (8,'(1pg14.6,6g24.16)') mass(ns+2*i),
     &        (rbar(j) - mu(1)*rrel(j),j=1,3),
     &        (vbar(j) - mu(1)*vrel(j),j=1,3)
      enddo
cCompute galactocentric radius for v = 220km/s
      rg = rt**1.5*sqrt(3*220**2/(g*totmass))
      write (6,*) totmass, rg, kt, select(nb/2,nb,alist),
     &     select(nb/2,nb,elist)
      stop
      end
c
c
      real function funcd(x,f,df)
      implicit none
      real x,f,df,ma,e
      common ma,e
      f = ma - x + e*sin(x)
      df = -1 + e*cos(x)
      return
      end
c
c
      FUNCTION select(k,n,arr)
      INTEGER k,n
      real arr(n)
      INTEGER i,ir,j,l,mid
      REAL a,temp
      l=1
      ir=n
1     if(ir-l.le.1)then
        if(ir-l.eq.1)then
          if(arr(ir).lt.arr(l))then
            temp=arr(l)
            arr(l)=arr(ir)
            arr(ir)=temp
          endif
        endif
        select=arr(k)
        return
      else
        mid=(l+ir)/2
        temp=arr(mid)
        arr(mid)=arr(l+1)
        arr(l+1)=temp
        if(arr(l+1).gt.arr(ir))then
          temp=arr(l+1)
          arr(l+1)=arr(ir)
          arr(ir)=temp
        endif
        if(arr(l).gt.arr(ir))then
          temp=arr(l)
          arr(l)=arr(ir)
          arr(ir)=temp
        endif
        if(arr(l+1).gt.arr(l))then
          temp=arr(l+1)
          arr(l+1)=arr(l)
          arr(l)=temp
        endif
        i=l+1
        j=ir
        a=arr(l)
3       continue
          i=i+1
        if(arr(i).lt.a)goto 3
4       continue
          j=j-1
        if(arr(j).gt.a)goto 4
        if(j.lt.i)goto 5
        temp=arr(i)
        arr(i)=arr(j)
        arr(j)=temp
        goto 3
5       arr(l)=arr(j)
        arr(j)=a
        if(j.ge.k)ir=j-1
        if(j.le.k)l=i
      endif
      goto 1
      END
C  (C) Copr. 1986-92 Numerical Recipes Software ]vrt43D04.
c
c
      real FUNCTION RAN2(IDUM)
      PARAMETER (M=714025,IA=1366,IC=150889,RM=1.4005112E-6)
      DIMENSION IR(97)
      DATA IFF /0/
      IF(IDUM.LT.0.OR.IFF.EQ.0)THEN
        IFF=1
        IDUM=MOD(IC-IDUM,M)
        DO 11 J=1,97
          IDUM=MOD(IA*IDUM+IC,M)
          IR(J)=IDUM
11      CONTINUE
        IDUM=MOD(IA*IDUM+IC,M)
        IY=IDUM
      ENDIF
      J=1+(97*IY)/M
      IF(J.GT.97.OR.J.LT.1)PAUSE
      IY=IR(J)
      RAN2=IY*RM
      IDUM=MOD(IA*IDUM+IC,M)
      IR(J)=IDUM
      RETURN
      END
c
c
      real FUNCTION RTNEWT(FUNCD,X1,X2,XACC)
      PARAMETER (JMAX=20)
      RTNEWT=.5*(X1+X2)
      DO 11 J=1,JMAX
        CALL FUNCD(RTNEWT,F,DF)
        DX=F/DF
        RTNEWT=RTNEWT-DX
        IF((X1-RTNEWT)*(RTNEWT-X2).LT.0.)PAUSE 'jumped out of brackets'
        IF(ABS(DX).LT.XACC) RETURN
11    CONTINUE
      PAUSE 'RTNEWT exceeding maximum iterations'
      END