induce.f

Go to the documentation of this file.

      SUBROUTINE induce(IPAIR,EMAX,EMIN,ICIRC,TC2,ANGLE,TG,EDAV)
*
*
*       Induced eccentricity of hierarchical binary.
*       --------------------------------------------
*
      include 'common6.h'
      REAL*8  a1(3),a2(3),xrel(3),vrel(3),ei(3),hi(3),ho(3),bhat(3)
*
*
*       Define c.m. & KS indices and evaluate semi-major axis & eccentricity.
      i = n + ipair
      i1 = 2*ipair - 1
      i2 = i1 + 1
      semi = -0.5d0*body(i)/h(ipair)
      ecc2 = (1.d0 - r(ipair)/semi)**2 + tdot2(ipair)**2
     &                                   /(body(i)*semi)
      ecc = sqrt(ecc2)
*
*       Ensure that KS components are resolved (perturbed not always done).
      CALL resolv(ipair,1)
*
*       Obtain inner & outer angular momentum by cyclic notation.
      rij2 = 0.d0
      vij2 = 0.d0
      rdot = 0.d0
      a12 = 0.d0
      a22 = 0.d0
      a1a2 = 0.d0
      ri2 = 0.d0
      vi2 = 0.d0
      rvi = 0.d0
      DO 5 k = 1,3
          rij2 = rij2 + (x(k,i) - x(k,jcomp))**2
          vij2 = vij2 + (xdot(k,i) - xdot(k,jcomp))**2
          rdot = rdot + (x(k,i) - x(k,jcomp))*(xdot(k,i) -xdot(k,jcomp))
          k1 = k + 1
          IF (k1.GT.3) k1 = 1
          k2 = k1 + 1
          IF (k2.GT.3) k2 = 1
          a1(k) = (x(k1,i1) - x(k1,i2))*(xdot(k2,i1) - xdot(k2,i2))
     &          - (x(k2,i1) - x(k2,i2))*(xdot(k1,i1) - xdot(k1,i2))
          a2(k) = (x(k1,jcomp) - x(k1,i))*(xdot(k2,jcomp) - xdot(k2,i))
     &          - (x(k2,jcomp) - x(k2,i))*(xdot(k1,jcomp) - xdot(k1,i))
          a12 = a12 + a1(k)**2
          a22 = a22 + a2(k)**2
          a1a2 = a1a2 + a1(k)*a2(k)
*       Form relative vectors and scalars for inner binary.
          xrel(k) = x(k,i1) - x(k,i2)
          vrel(k) = xdot(k,i1) - xdot(k,i2)
          ri2 = ri2 + xrel(k)**2
          vi2 = vi2 + vrel(k)**2
          rvi = rvi + xrel(k)*vrel(k)
    5 CONTINUE
*
*       Evaluate orbital parameters for outer orbit.
      rij = sqrt(rij2)
      zmb = body(i) + body(jcomp)
      semi1 = 2.d0/rij - vij2/zmb
      semi1 = 1.d0/semi1
      ecc1 = sqrt((1.d0 - rij/semi1)**2 + rdot**2/(semi1*zmb))
*     PMIN1 = SEMI1*(1.d0 - ECC1)
*       Determine inclination (8 bins of 22.5 degrees).
      fac = a1a2/sqrt(a12*a22)
      fac = acos(fac)
      angle = fac
***   ANGLE = FAC*360.d0/TWOPI
      in = 1 + fac*360.0/(twopi*22.5)
*
*       Exit on hyperbolic orbit or outer eccentricity near 1.
      IF (semi1.LT.0.0.OR.ecc1.GT.0.9999) THEN
          emax = 0.d0
          tc2 = 999.d0
          tg = 1.0d+04
          edav = 1.0d+04
          go to 30
      END IF
*
*       Construct the Runge-Lenz vector (Heggie & Rasio 1995, Eq.(5)).
      ei2 = 0.d0
      DO 10 k = 1,3
          ei(k) = (vi2*xrel(k) - rvi*vrel(k))/body(i) -
     &                                                 xrel(k)/sqrt(ri2)
          ei2 = ei2 + ei(k)**2
   10 CONTINUE
*
*       Define unit vectors for inner eccentricity and angular momenta.
      cosj = 0.d0
      sjsg = 0.d0
      DO 15 k = 1,3
          ei(k) = ei(k)/sqrt(ei2)
          hi(k) = a1(k)/sqrt(a12)
          ho(k) = a2(k)/sqrt(a22)
          cosj = cosj + hi(k)*ho(k)
          sjsg = sjsg + ei(k)*ho(k)
   15 CONTINUE
*
*       Form unit vector BHAT and scalars AH & BH (Douglas Heggie, 10/9/96).
      ah = 0.d0
      bh = 0.d0
      DO 20 k = 1,3
          k1 = k + 1
          IF (k1.GT.3) k1 = 1
          k2 = k1 + 1
          IF (k2.GT.3) k2 = 1
          bhat(k) = hi(k1)*ei(k2) - hi(k2)*ei(k1)
          ah = ah + ei(k)*ho(k)
          bh = bh + bhat(k)*ho(k)
   20 CONTINUE
*
*       Evaluate the expressions A & Z.
      a = cosj*sqrt(1.d0 - ei2)
      z = (1.d0 - ei2)*(2.d0 - cosj**2) + 5.d0*ei2*sjsg**2
*
*       Obtain maximum inner eccentricity (Douglas Heggie, Sept. 1995).
      z2 = z**2 + 25.d0 + 16.d0*a**4 - 10.d0*z - 
     &     20.d0*a**2 - 8.d0*a**2*z
      emax = one6*(z + 1.d0 - 4.d0*a**2 + sqrt(z2))
      emax = max(emax,0.0001d0)
      emax = sqrt(emax)
      zi = fac
      em = sqrt(sin(zi)**2 + ei2*cos(zi)**2)
*
*       Form minimum eccentricity (Douglas Heggie, Sept. 1996).
      az = a**2 + z - 2.d0
      IF (az.GE.0.0) THEN
          az1 = 1.d0 + z - 4.d0*a**2
          emin2 = one6*(az1 - sqrt(az1**2 - 12.d0*az))
      ELSE
          emin2 = 1.d0 - 0.5d0*(a**2 + z)
      END IF
      emin2 = max(emin2,0.0001d0)
      emin = sqrt(emin2)
*
*       Set impact parameters and estimate eccentricity growth time-scale.
      pmin = semi*(1.d0 - ecc)
      pmin2 = semi*(1.d0 - emax)
      tk = twopi*semi*sqrt(semi/body(i))
      tk1 = twopi*semi1*sqrt(semi1/zmb)
      tg = tk1**2*zmb*(1.d0 - ecc1**2)**1.5/(body(jcomp)*tk)
*
      const = pfac(a,z)
      const = const*4.d0/(1.5d0*twopi*sqrt(6.d0))
*       Convert growth time to units of 10**6 yrs.
      tg = const*tg*tstar
*
*       Form doubly averaged eccentricity derivative (Douglas Heggie 9/96).
      yfac = 15.d0*body(jcomp)/(4.d0*zmb)*twopi*tk/tk1**2
      yfac = yfac*ecc*sqrt(1.d0 - ecc**2)/
     &                             (1.d0 - ecc1**2)**(1.5)
      edav = yfac*ah*bh
*
*       Determine circularization time for current & smallest pericentre.
      tc2 = 2000.d0
      tc = 2000.d0
      IF (pmin2.LT.10.0*max(radius(i1),radius(i2)).AND.
     &    kz(27).EQ.2.AND.ecc.GT.0.002) THEN
          CALL tcirc(pmin,ecc,i1,i2,icirc,tc)
          IF (icirc.GE.0) icirc = -1
*       Note: only use PMIN call because ICIRC may become >0 after first!
          CALL tcirc(pmin2,emax,i1,i2,icirc,tc2)
      ELSE IF (ecc.LT.0.002) THEN
          tc = 0.0
          tc2 = 0.0
      END IF
*
*       Print diagnostics for high inclinations and TC2 < 10**7 yrs.
      IF (in.GE.4.AND.in.LE.5.AND.tc2.LT.10.0.AND.kstar(i).NE.-2.AND.
     &    ei2.GT.4.0d-06) THEN
          WRITE (44,25)  sqrt(ei2), em, emax, kstar(i1), kstar(i2),
     &                   kstar(i), semi, pmin2, in, tg, tc, tc2, tphys
   25     FORMAT (' INDUCE:    E EMX EMAX K* SEMI PM2 IN TG TC TC2 TP ',
     &                         3f8.4,3i3,1p,2e9.1,0p,i3,f7.3,3f7.1)
          CALL flush(44)
      END IF
*
*       Use first value for now.
      tc2 = tc
*
   30 RETURN
*
      END