hipop.f

Go to the documentation of this file.

      SUBROUTINE hipop
*
*
*       Initial hierarchical population.
*       --------------------------------
*
      include 'common6.h'
      REAL*8  xorb(2),vorb(2),xrel(3),vrel(3),px(3),qx(3),
     &        bs(mmax),xs(3,mmax),vs(3,mmax)
      REAL*8  ran2
*
*
*       Read input parameters (same usage as routine BINPOP).
      READ (5,*)  semi0, ecc0, ratio, range, icirc
      nhi = nhi0
      WRITE (6,1)  nhi, semi0, ecc0, ratio, range, icirc
    1 FORMAT (/,12x,'HIERARCHIES:    NHI =',i4,'  A =',f9.6,
     &              '  E =',f6.2,'  RATIO =',f4.1,'  RANGE =',f6.1,
     &              '  ICIRC =',i2,/)
*
      IF (nhi.GT.nbin0) THEN
          WRITE (6,2)  nbin0, nhi
    2     FORMAT (5x,'FATAL ERROR!   NBIN0 =',i4,'  NHI =',i4)
          stop
      END IF
*
*       Introduce binary components by splitting the primary.
      DO 40 i = 1,nhi
*
*       Randomize perihelion, node & inclination.
          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) 
*
*       Determine two-body elements for original binary.
          i1 = 2*i - 1
          i2 = 2*i
          rij2 = 0.0
          vij2 = 0.0
          rdot = 0.0
          DO 5 k = 1,3
              xrel(k) = x(k,i1) - x(k,i2)
              vrel(k) = xdot(k,i1) - xdot(k,i2)
              rij2 = rij2 + xrel(k)**2
              vij2 = vij2 + vrel(k)**2
              rdot = rdot + xrel(k)*vrel(k)
    5     CONTINUE
          rij = sqrt(rij2)
          zmb1 = body(i1) + body(i2)
          a1 = 2.0/rij - vij2/zmb1
          semi1 = 1.0/a1
          ecc2 = (1.0 - rij/semi1)**2 + rdot**2/(semi1*zmb1)
          ecc1 = sqrt(ecc2)
          pmin = semi1*(1.0 - ecc1)
*
*       Specify component masses (primary fraction range 0.5 - 0.9).
          q0 = 0.5 + 0.4*ran2(idum1)
          IF (ratio.EQ.1.0) q0 = 0.5
          bs(i) = body(i2)
          body(i1) = q0*body(i1)
          body(i2) = body(i1)*(1.0 - q0)/q0
          zmb = body(i1) + body(i2)
*
*       Choose random (thermalized) or fixed eccentricity.
          IF (ecc0.LT.0.0) THEN
              ecc2 = ran2(idum1)
              ecc = sqrt(ecc2)
          ELSE
              ecc = ecc0
          END IF
*
*       Select semi-major axis from uniform distribution in log(A) or SEMI0.
          iter = 0
   10     IF (range.GT.0.0) THEN
              IF (iter.LE.5) THEN
                  exp = ran2(idum1)*log10(range)
                  semi = semi1/10.0**exp
              ELSE
*       Shrink by factor 2 if no success after 5 iterations.
                  semi = 0.5*semi
              END IF
          ELSE
              semi = semi0
          END IF
*
*       Check stability criterion (maximum 12 tries with inclination effect).
          pcrit = stability(body(i1),body(i2),bs(i),ecc,ecc1,zi)
          pcrit = pcrit*semi
          iter = iter + 1
          IF (pmin.LT.pcrit.AND.iter.LE.12.AND.semi.LT.semi0) go to 10
*
          p0 = days*semi*sqrt(semi/zmb)
          p1 = days*semi1*sqrt(semi1/zmb1)
          WRITE (6,20)  iter, ecc, ecc1, pmin, pcrit, p0, p1
   20     FORMAT (' HIERARCHY:    IT E E1 PMIN PCRIT P0 P1 ',
     &                            i4,2f7.3,1p,4e9.1)
*
*       Specify relative motion at apocentre and sum binding energy.
          xorb(1) = semi*(1.0 + ecc)
          xorb(2) = 0.0
          vorb(1) = 0.0
          vorb(2) = sqrt(zmb*(1.0d0 - ecc)/(semi*(1.0d0 + ecc)))
          ebin0 = ebin0 - 0.5*body(i1)*body(i2)/semi
*
*       Transform to relative variables.
          DO 25 k = 1,3
              xrel(k) = px(k)*xorb(1) + qx(k)*xorb(2)
              vrel(k) = px(k)*vorb(1) + qx(k)*vorb(2)
   25     CONTINUE
*
*       Save old secondary and set global variables for each component.
          DO 30 k = 1,3
              xs(k,i) = x(k,i2)
              vs(k,i) = xdot(k,i2)
              x(k,i1) = x(k,i1) + body(i2)*xrel(k)/zmb
              x(k,i2) = x(k,i1) - xrel(k)
              xdot(k,i1) = xdot(k,i1) + body(i2)*vrel(k)/zmb
              xdot(k,i2) = xdot(k,i1) - vrel(k)
   30     CONTINUE
   40 CONTINUE
*
*       Move single particles down by NHI to make room for outer components.
      l = n
      DO 50 j = 2*nbin0+1,n
          i = nhi + l
          body(i) = body(l)
          DO 45 k = 1,3
              x(k,i) = x(k,l)
              xdot(k,i) = xdot(k,l)
   45     CONTINUE
          l = l - 1
   50 CONTINUE
*
*       Place hierarchical components immediately after the binaries.
      DO 60 l = 1,nhi
          i = 2*nbin0 + l
          body(i) = bs(l)
          DO 55 k = 1,3
              x(k,i) = xs(k,l)
              xdot(k,i) = vs(k,l)
   55     CONTINUE
   60 CONTINUE
*
*       Update the particle number and reset NHI (might be used elsewhere).
      n = n + nhi
      nzero = n
      ntot = n
      nhi = 0
*
      RETURN
*
      END