kozai.f

Go to the documentation of this file.

      SUBROUTINE kozai(IPAIR,IM,ECC1,SEMI1,ITERM)
*
*
*       Kozai relations & decision-making.
*       ----------------------------------
*
      include 'common6.h'
      common/binary/  cm(4,mmax),xrel(3,mmax),vrel(3,mmax),
     &                hm(mmax),um(4,mmax),umdot(4,mmax),tmdis(mmax),
     &                namem(mmax),nameg(mmax),kstarm(mmax),iflagm(mmax)
      common/koztmp/  emax
      REAL*8  ei(3),hi(3),ho(3),a1(3),a2(3)
      SAVE it
      DATA it /0/
*
*
*       Define KS & c.m. indices for outer orbit with ECC1 & SEMI1.
      i1 = 2*ipair - 1
      i2 = i1 + 1
      i = n + ipair
      iterm = 0
*
*       Evaluate inner semi-major axis and eccentricity.
      ri = 0.0
      td2 = 0.0
      DO 1 k = 1,4
          ri = ri + um(k,im)**2
          td2 = td2 + 2.0*um(k,im)*umdot(k,im)
    1 CONTINUE
      zmb = cm(1,im) + cm(2,im)
      semi = -0.5*zmb/hm(im)
      ecc2 = (1.0 - ri/semi)**2 + td2**2/(zmb*semi)
      ecc = sqrt(ecc2)
*
*       Skip circularized binaries.
      IF (ecc.LT.0.1) THEN
          IF (tev(i).GE.1.0d+10) tev(i) = time + 1.0
          go to 50
      END IF
*
*       Form useful scalars (XREL & VREL consistent with UM & UMDOT).
      a12 = 0.0
      a22 = 0.0
      a1a2 = 0.0
      ri2 = 0.0
      vi2 = 0.0
      rv = 0.0
      DO 5 k = 1,3
          k1 = k + 1
          IF (k1.GT.3) k1 = 1
          k2 = k1 + 1
          IF (k2.GT.3) k2 = 1
          a1(k) = xrel(k1,im)*vrel(k2,im) - xrel(k2,im)*vrel(k1,im)
          a2(k) = (x(k1,i2) - x(k1,i1))*(xdot(k2,i2) - xdot(k2,i1))
     &          - (x(k2,i2) - x(k2,i1))*(xdot(k1,i2) - xdot(k1,i1))
          a12 = a12 + a1(k)**2
          a22 = a22 + a2(k)**2
          a1a2 = a1a2 + a1(k)*a2(k)
          ri2 = ri2 + xrel(k,im)**2
          vi2 = vi2 + vrel(k,im)**2
          rv = rv + xrel(k,im)*vrel(k,im)
    5 CONTINUE
*
*       Construct the Runge-Lenz vector (Heggie & Rasio 1995, Eq.(5)).
      ei2 = 0.0
      DO 10 k = 1,3
          ei(k) = (vi2*xrel(k,im) - rv*vrel(k,im))/body(i1) -
     &                              xrel(k,im)/sqrt(ri2)
          ei2 = ei2 + ei(k)**2
   10 CONTINUE
      ei2 = min(ei2,0.999999d0)
*
*       Define unit vectors for inner eccentricity and angular momenta.
      cosj = 0.0
      sjsg = 0.0
      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
*
*       Evaluate the expressions A & Z.
      a = cosj*sqrt(1.0 - ei2)
      z = (1.0 - ei2)*(2.0 - cosj**2) + 5.0*ei2*sjsg**2
*
*       Obtain maximum inner eccentricity (Douglas Heggie, Sept. 1995).
      z2 = z**2 + 25.0 + 16.0*a**4 - 10.0*z - 20.0*a**2 - 8.0*a**2*z
      emax = one6*(z + 1.0 - 4.0*a**2 + sqrt(z2))
      emax = max(emax,0.0001d0)
      emax = sqrt(emax)
*       Skip max eccentricity below 0.90.
      IF (emax.LT.0.0) go to 50
*
*       Form minimum eccentricity (Douglas Heggie, Sept. 1996).
      az = a**2 + z - 2.0
      IF (az.GE.0.0) THEN
          az1 = 1.0 + z - 4.0*a**2
          emin2 = one6*(az1 - sqrt(az1**2 - 12.0*az))
      ELSE
          emin2 = 1.0 - 0.5*(a**2 + z)
      END IF
      emin2 = max(emin2,0.0001d0)
      emin = sqrt(emin2)
*
*       Estimate Kozai time-scale (N-body units).
      tk = twopi*semi*sqrt(semi/body(i1))
      tk1 = twopi*abs(semi1)*sqrt(abs(semi1)/body(i))
      tkoz = tk1**2*body(i)*(1.0 - ecc1**2)**1.5/(body(i2)*tk)
*
*       Evaluate numerical precession factor in Myr (elliptic integral).
      const = pfac(a,z)
*       Note published factor 2/(3*pi) (Kiseleva, Eggleton & Mikkola MN/98).
      const = const*4.0/(1.5*twopi*sqrt(6.0))
      tkoz = const*tkoz*tstar
*
*       Obtain inclination (see routine INCLIN).
      fac = min(cosj,1.0d0)
      angle = acos(fac)
      zi = angle*360.0/twopi
*
*       Check optional output (every 100 orbits or each time for #42 > 1).
      it = it + 1
      IF ((mod(it,100).EQ.0.AND.emax.GT.0.99).OR.kz(42).GT.1) THEN
          WRITE (42,30)  time+toff, name(i1), zi, emin, emax, ecc,
     &                   semi1/semi, tkoz, tmdis(im)
   30     FORMAT (' KOZAI    T NAM1 INC EM EX E A1/A TKOZ TMD ',
     &                       f9.3,i7,f7.1,f6.2,2f8.4,f6.1,1p,2e9.1)
          CALL flush(42)
          IF (it.GT.2000000000) it = 0
      END IF
*
*       Search for circularization candidates with large EMAX & inclination.
      IF (emax.GT.0.90.AND.zi.GT.40.0) THEN
          CALL findj(i1,jg,im2)
*       Include safety test to avoid negative JG or JG > N.
          IF (jg.LT.0.OR.jg.GT.n) go to 50
*       Consider termination condition for tidal circularization (#27 = 2).
          IF (semi*(1.0 - emax).LT.2.0*(radius(i1) + radius(jg)).AND.
     &        kz(27).EQ.2) THEN
*
*       Estimate circularization time in Myr (C. Tout, Cambody lecture 2006).
              IF (radius(jg).GT.radius(i1)) THEN
                  j1 = jg
                  j2 = i1
                  q1 = cm(2,im)/cm(1,im)
              ELSE
                  j1 = i1
                  j2 = jg
                  q1 = cm(1,im)/cm(2,im)
              END IF
*             SEMI0 = SEMI*(1.0 - ECC2)
              semi0 = semi*(1.0 - emax**2)
              tcirc = 2.0*q1**2/(1.0 + q1)*(semi0/radius(j1))**8
              tcirc = 1.0d-06*tcirc
*             WRITE (45,40)  TIME+TOFF, NAME(I1), ECC, SEMI, SEMI0,
*    &                       TCIRC
*  40         FORMAT (' CIRCULARIZED    T NM E A0 A TC ',
*    &                                  F9.2,I6,F8.4,1P,4E10.2)
*       Note possibility of switching if nothing happens.
              rp = semi*(1.0 - emax)
              tt = min(tkoz,tcirc)
              WRITE (6,45)  emax, semi, rp, radius(i1)+radius(jg), tt
   45         FORMAT (' KOZAI TIDE    EX A RP R1+R2 TT ',
     &                                f10.6,1p,e12.4,3e10.2)
*             ITERM = 1
*       Perform circularization on small Kozai period or circ time.
              IF (tkoz.LT.1.0.OR.tcirc.LT.1) iterm = 2
              tev0(i) = time
              tev(i) = time + 1.0
          END IF
      END IF
*
   50 RETURN
*
      END