binpop.f

Go to the documentation of this file.

      SUBROUTINE binpop
*
*
*       Initial binary population.
*       --------------------------
*
      include 'common6.h'
      REAL*8  xorb(2),vorb(2),xrel(3),vrel(3),px(3),qx(3),bs(nmax)
      REAL*8  ran2
      DATA  eta1,eta2 /2.5,45.0/
*
*
      READ (5,*)  semi0, ecc0, ratio, range, nskip, idorm, icirc
      nbin = nbin0
      nbin1 = nbin + 1
      WRITE (6,1)  nbin, semi0, ecc0, ratio, range, nskip, idorm, icirc
    1 FORMAT (/,12x,'BINARIES:   NBIN =',i5,'  A =',f9.6,'  E =',f6.2,
     &              '  RATIO =',f4.1,'  RANGE =',f6.1,'  NSKIP =',i3,
     &              '  IDORM =',i2,'  ICIRC =',i2,/)
*
*       Check type of binary mass distribution (NSKIP, IMF2 or split c.m.).
      IF (nskip.EQ.0.OR.kz(20).GE.2) go to 10
      IF (ratio.EQ.1.0) go to 20
*
*       Select binaries from the most massive bodies (frequency NSKIP).
      ilast = (1 + nskip)*nbin
      jskip = 0
      js = 0
      jb = 1
*
*       Transfer binary masses to first NBIN locations.
      DO 6 i = 2,ilast
          jskip = jskip + 1
*       Copy binary mass of body #I to new global location.
          IF (jskip.GT.nskip) THEN
              jskip = 0
              jb = jb + 1
              body(jb) = body(i)
          ELSE
*       Save next NSKIP masses of single bodies.
              js = js + 1
              bs(js) = body(i)
          END IF
    6 CONTINUE
*
*       Restore the single bodies in subsequent locations.
      js = 0
      DO 8 i = nbin1,ilast
          js = js + 1
          body(i) = bs(js)
    8 CONTINUE
*
*       Move main variables of all single bodies.
   10 DO 15 i = n,nbin1,-1
          j = i + nbin
          body(j) = body(i)
          DO 12 k = 1,3
              x(k,j) = x(k,i)
              xdot(k,j) = xdot(k,i)
   12     CONTINUE
   15 CONTINUE
*
*       Create space for each binary component next to primary.
   20 DO 30 i = nbin,2,-1
          j = 2*i - 1
          body(j) = body(i)
          DO 25 k = 1,3
              x(k,j) = x(k,i)
              xdot(k,j) = xdot(k,i)
   25     CONTINUE
   30 CONTINUE
*
*       Introduce binary components from relative motion.
      DO 60 i = 1,nbin
*
*       Randomize perihelion, node & inclination (ZI = 0.25 before 3/99).
          pi = twopi*ran2(idum1)
          omega = twopi*ran2(idum1)
          zi = 0.5*twopi*ran2(idum1)
*
*       Set transformation elements (Brouwer & Clemence p. 35).
          px(1) = cos(pi)*cos(omega) - sin(pi)*sin(omega)*cos(zi)
          qx(1) =-sin(pi)*cos(omega) - cos(pi)*sin(omega)*cos(zi)
          px(2) = cos(pi)*sin(omega) + sin(pi)*cos(omega)*cos(zi)
          qx(2) =-sin(pi)*sin(omega) + cos(pi)*cos(omega)*cos(zi)
          px(3) = sin(pi)*sin(zi)
          qx(3) = cos(pi)*sin(zi) 
*
*       Specify component masses (copy BODY0 from IMF2 or use RATIO).
          i1 = 2*i - 1
          i2 = 2*i
          IF (kz(20).GE.2) THEN
              body(i1) = body0(i1)
              body(i2) = body0(i2)
          ELSE IF (ratio.EQ.1.0) THEN
              body(i2) = body(i1) 
          ELSE
              body(i1) = ratio*body(i1)
              body(i2) = body(i1)*(1.0 - ratio)/ratio
          END IF
*
*       Choose random (thermalized; < 0.9 for triples) or fixed eccentricity.
          IF (ecc0.LT.0.0) THEN
              ecc2 = ran2(idum1)
              ecc = sqrt(ecc2)
              IF (kz(18).GT.1) ecc = min(ecc,0.9d0)
          ELSE
              ecc = ecc0
          END IF
*
*       Generate semi-major axis (new option added 4/8/05, modified 27/7/10).
          IF (kz(8).EQ.4) THEN
*       Adopt Pavel Kroupa (1995, Eq.11b) distribution for semi-major axis.
             exp1 = exp(2.d0*ran2(idum1)/2.5d0) - 1.d0
             exp1 = sqrt(exp1*45.d0) + 1.d0
*       Specify period in yrs.
             exp1 = 10**exp1 /365.25d0
*       Transform to semi-major axis in model units by Kepler's Law.
             semi = (body(i1)+body(i2))*zmbar*exp1*exp1
             semi = semi**(1.d0/3.d0)
*       Set semi in pc and then in model units.
             semi = semi/206259.591d0
             semi = semi/rbar
*       Select semi-major axis from uniform distribution in log(A) or SEMI0.
          ELSE IF (range.GT.0.0) THEN
              exp1 = ran2(idum1)*log10(range)
              semi = semi0/10.0**exp1
          ELSE
              semi = semi0
          END IF
*
*       Check for eigen-evolution (Pavel Kroupa & Rosemary Mardling).
          IF (icirc.EQ.1) THEN
              icoll = 0
              ic0 = 0
              ic1 = 0
              ic2 = 0
              ic3 = 0
              zmb = (body(i1) + body(i2))*zmbar
*       Include minimum period (copy RANGE; at least 1 day).
              pmin = max(range,1.0d0)
              it = 0
   35         xr = ran2(idum1)
*       Generate period distribution (Pavel Kroupa: MN 277, 1491, eq.11b).
              p0 = log10(pmin) + sqrt(eta2*(exp(2.0*xr/eta1) - 1.0))
              tk = 10.0**p0
*       Invert eccentricity from thermal distribution (XR = E**2).
              xr = ran2(idum1)
              es0 = sqrt(xr)
*       Set pericentre distance in AU with period in days & mass in SU.
              rp0 = (1.0 - es0)*((tk/365.0)**2*zmb)**0.3333
*       Convert to N-body units.
              rp0 = rp0/rau
              a0 = rp0/(1.0 - es0)
              e0 = es0
*       Limit the maximum semi-major axis to 1000 AU (Oct 2004).
              IF (a0*rau.GT.1000.0) go to 35
*       Define K* = 0/1 and enhanced radii for pre-main sequence.
              kstar(i1) = 1
              kstar(i2) = 1
              IF (body(i1)*zmbar.LT.0.7) kstar(i1) = 0
              IF (body(i2)*zmbar.LT.0.7) kstar(i2) = 0
              radius(i1) = 5.0*sqrt(body(i1)*zmbar)/su
              radius(i2) = 5.0*sqrt(body(i2)*zmbar)/su
              IF (rp0.LT.max(radius(i1),radius(i2))) THEN
                  WRITE (6,38)  i1, zmb, es0, a0, rp0
   38             FORMAT (12x,'COLLISION:    I1 MB E A RP ',
     &                                       i6,f6.1,f7.3,1p,2e10.2)
                  icoll = icoll + 1
                  go to 35
              END IF
*       Perform eigen-evolution of pericentre & eccentricity for 10^6 yrs.
              tc = -1.0/tscale
              CALL tcirc(rp0,es0,i1,i2,icirc,tc)
*       Copy modified eccentricity and re-evaluate the semi-major axis.
              ecc = es0
              semi = rp0/(1.0 - ecc)
              it = it + 1
              IF (semi.GT.semi0.AND.it.LT.25) go to 35
              tk = 365.0*sqrt((semi*rau)**3/zmb)
              IF (ecc.LE.0.002) ic0 = ic0 + 1
              IF (tk.LT.pmin) ic1 = ic1 + 1
              IF (tk.LT.2.0*pmin) ic2 = ic2 + 1
              IF (tk.LT.5.0*pmin) ic3 = ic3 + 1
              WRITE (23,40)  it, i1, zmb, e0, ecc, a0, semi, tk
   40         FORMAT (12x,'BINARY:   IT I1 MB E0 E A0 A P ',
     &                               i2,i5,f5.1,2f7.3,1p,3e10.2)
              CALL flush(23)
          ELSE IF (kz(27).EQ.1.OR.kz(19).GE.3) THEN
*       Obtain tidal encounter distance (4*RADIUS) from square root relation.
              rsun = 1.0/su
              zm = max(body(i1),body(i2))*zmbar
              rt = 4.0*rsun*sqrt(zm)
              IF (kstar(i1).GT.1.OR.kstar(i2).GT.1) rt = 0.d0
*       Modify orbital elements to avoid tidal interaction or collision.
   42         IF (semi*(1.0 - ecc).LT.2.0*rt) THEN
*       Increase semi-major axis or reduce eccentricity until peri > 2*RT.
   44             IF (semi.LT.2.0*rt) THEN
                      semi = 2.0*semi
                      go to 44
                  ELSE
                      ecc = 0.9*ecc
                  END IF
                  WRITE (51,46)  i1, i2, ecc, semi, semi*(1.0-ecc), rt
   46             FORMAT (12x,'REDUCE ECC:    I1 I2 E A PM RT ',
     &                                        2i5,f7.3,1p,3e10.2)
                  CALL flush(51)
                  go to 42
              END IF
          END IF
*
*       Evolve ECC and SEMI via circularization and feeding (Kroupa 1995).
          IF (icirc.EQ.2) THEN
              CALL proto_star(zmbar,rbar,body(i1),body(i2),ecc,semi)
          END IF
*
*       Specify relative motion at apocentre.
          xorb(1) = semi*(1.0 + ecc)
*       Specify relative motion at apocentre and sum binding energy.
          xorb(1) = semi*(1.0 + ecc)
          xorb(2) = 0.0
          vorb(1) = 0.0
          zmbin = body(i1) + body(i2)
          vorb(2) = sqrt(zmbin*(1.0d0 - ecc)/(semi*(1.0d0 + ecc)))
          ebin0 = ebin0 - 0.5*body(i1)*body(i2)/semi
*
*       Transform to relative variables.
          DO 50 k = 1,3
              xrel(k) = px(k)*xorb(1) + qx(k)*xorb(2)
              vrel(k) = px(k)*vorb(1) + qx(k)*vorb(2)
   50     CONTINUE
*
*       Set global variables for each component.
          DO 55 k = 1,3
              x(k,i1) = x(k,i1) + body(i2)*xrel(k)/zmbin
              x(k,i2) = x(k,i1) - xrel(k)
              xdot(k,i1) = xdot(k,i1) + body(i2)*vrel(k)/zmbin
              xdot(k,i2) = xdot(k,i1) - vrel(k)
   55     CONTINUE
   60 CONTINUE
*
*       Update the total particle number after primary splitting or IMF2.
      IF (ratio.LT.1.0.OR.kz(20).GE.2) THEN
          n = n + nbin
          nzero = n
          ntot = n
          bodym = zmass/float(n)
          IF (nskip.GT.0) THEN
              WRITE (6,62)  (body(j),j=1,10)
              WRITE (6,64)  (body(j),j=2*nbin+1,2*nbin+10)
   62         FORMAT (/,12x,'BINARY MASSES (1-10):  ',10f9.5)
   64         FORMAT (/,12x,'SINGLE MASSES (1-10):  ',10f9.5,/) 
          END IF
      END IF
*
*       Include procedure for introducing dormant binaries.
      IF (idorm.GT.0) THEN
          DO 66 i = 1,nbin
              i1 = 2*i - 1
              i2 = i1 + 1
              zmbin = body(i1) + body(i2)
              DO 65 k = 1,3
                  x(k,i) = (body(i1)*x(k,i1) + body(i2)*x(k,i2))/zmbin
                  xdot(k,i) = (body(i1)*xdot(k,i1) +
     &                         body(i2)*xdot(k,i2))/zmbin
   65         CONTINUE
              body(i) = zmbin
   66     CONTINUE
*
*       Move the original single particles up to form compact array.
          i1 = 2*nbin + 1
          i2 = nbin
          DO 68 i = i1,n
              i2 = i2 + 1
              body(i2) = body(i)
              DO 67 k = 1,3
                  x(k,i2) = x(k,i)
                  xdot(k,i2) = xdot(k,i)
   67         CONTINUE
   68     CONTINUE
*
*       Reset particle membership and turn off binary output option (if = 1).
          n = n - nbin
          nzero = n
          ntot = n
          nbin0 = 0
          ebin0 = 0.0
          IF (kz(8).EQ.1) kz(8) = 0
      END IF
*
*       Set coordinates & velocities in c.m. rest frame.
      DO 70 k = 1,3
          cmr(k) = 0.0d0
          cmrdot(k) = 0.0d0
   70 CONTINUE
*
      DO 80 i = 1,n
          DO 75 k = 1,3
              cmr(k) = cmr(k) + body(i)*x(k,i)
              cmrdot(k) = cmrdot(k) + body(i)*xdot(k,i)
   75     CONTINUE
   80 CONTINUE
*
      DO 90 i = 1,n
          DO 85 k = 1,3
              x(k,i) = x(k,i) - cmr(k)/zmass
              xdot(k,i) = xdot(k,i) - cmrdot(k)/zmass
   85     CONTINUE
   90 CONTINUE
*
      RETURN
*
      END