tcirc.f

Go to the documentation of this file.

      SUBROUTINE tcirc(RP,ES0,I1,I2,ICIRC,TC)
*
*
*       Circularization time.
*       ---------------------
*
*       Theory of Piet Hut, A & A 99, 126 (1981).
*       Developed by Rosemary Mardling (31/1/97).
*       @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
*
      include 'common6.h'
      REAL*8  ww(3),qq(3),w(2),q(2),at0(2),m21,wg(2),qg(2),
     &        wscale(2),qscale(2),a(2),b(2),c(6)
      DATA  ww  /2.119,3.113,8.175/
      DATA  qq  /0.4909,0.4219,0.2372/
      DATA  a  /6.306505,-7.297806/
      DATA  b  /32.17211,13.01598/
      DATA  c  /5.101417,24.71539,-9.627739,1.733964,
     &                            -2.314374,-4.127795/
      DATA  eccm,eccm1 /0.002,0.002001/
      SAVE itime
      DATA itime /0/
*
*
*       Set large circularization time for merged binary or BHs.
      IF (radius(i1).EQ.0.0d0.OR.radius(i2).EQ.0.0d0.OR.
     &   (kstar(i1).EQ.14.AND.kstar(i2).EQ.14)) THEN
          tc = 1.0d+10
          go to 30
      END IF
*
*       Specify index J1 as biggest radius to be used with AT0(1).
      IF (radius(i1).GE.radius(i2)) THEN
          j1 = i1
          j2 = i2
      ELSE
          j1 = i2
          j2 = i1
      END IF
*
*       Include an appropriate numerical factor (Chris Tout 07/08).
      ar = rp/((1.0 - es0)*radius(j1))
      ff = 1.0
      IF (kstar(j1).EQ.1.AND.body(j1)*smu.GT.1.25) THEN
          ff = 100.0
      ELSE IF (kstar(j1).EQ.4.OR.kstar(j1).EQ.7) THEN
          ff = 100.0
      ELSE IF (kstar(j1).GE.10) THEN
          ff = 1.0d+17/ar**2
      END IF
      q12 = body(j1)/body(j2)
      tauc = ff*2.0*q12**2/(1.0 + q12)*ar**8
*       Exit on TC > 1.0D+08 yr (Chris Tout, Cambody Book 2008).
      IF (tauc.GT.1.0d+08) THEN
          tc = 1.0d+10
          go to 30
      END IF
*
*       Define oscillation period (dimensionless time) and damping constants.
      xn = 0.0
      DO 5 k = 1,2
          IF (k.EQ.1) THEN
              ik = j1
          ELSE
              ik = j2
          END IF
*       Specify polytropic index for each star (n = 3, 2 or 3/2).
          IF (kstar(ik).EQ.3.OR.kstar(ik).EQ.5.OR.
     &        kstar(ik).EQ.6.OR.kstar(ik).EQ.9) THEN
              ipair = kvec(i1)
              CALL giant(ipair,ik,wg,qg,wscale,qscale,xn,ql)
              w(k) = wg(1)
              q(k) = qg(1)
          ELSE
              ql = 1.0d+04
              ip = 3
              IF (kstar(ik).GE.3) ip = 2
              IF (kstar(ik).EQ.4.OR.kstar(ik).EQ.7) ip = 3
              IF (kstar(ik).EQ.8) ip = 3
              IF (kstar(ik).EQ.0) ip = 1
              w(k) = ww(ip)
              q(k) = qq(ip)
          END IF
          tl = twopi*radius(ik)*sqrt(radius(ik)/body(ik)/w(k))
          at0(k) = 1.0/(ql*tl)
    5 CONTINUE
*
*       Form mass, radius & pericentre ratio.
      IF (radius(i1).GE.radius(i2)) THEN
          m21 = body(i2)/body(i1)
          r21 = radius(i2)/radius(i1)
      rp1 = rp/radius(i1)
      ELSE
      m21 = body(i1)/body(i2)
          r21 = radius(i1)/radius(i2)
      rp1 = rp/radius(i2)
      END IF
*
*   Evaluate damping coefficient.
      rr = rp1*(1.0 + es0)
      const = 2.0*(at0(1)*(q(1)/w(1))**2*(1.0 + m21)*m21 + 
     &             at0(2)*(q(2)/w(2))**2*((1.0 + m21)/m21**2)*r21**8)/
     &                                                         rr**8
*
*       Adopt WD scaling for any NS/BH to avoid numerical problem.
      IF (kstar(i1).GE.13.OR.kstar(i2).GE.13) THEN
          const = 1.0d-04*const
      END IF
*
*   Form rational function approximation to Hut solution.
      ff = (( a(2)*es0 + a(1))*es0 + 1.0 )/
     &     (( b(2)*es0 + b(1))*es0 + 1.0 )
*     FF = MIN(FF,0.999D0)
*
*       See whether we only want the modified eccentricity (routine BINPOP).
      IF (icirc.LE.0) go to 10
*
      time0 = 0.0
*       Obtain the new eccentricity.
      z = (time - time0)*const + ff
      ecc = (-1.0 + c(1)*z - sqrt(c(2)*z**2 + c(3)*z + c(4)))
     &                                   /(c(5) + c(6)*z)
*
      ecc = max(ecc,eccm)
      rp = rp*(1.0 + es0)/(1.0 + ecc)
      es0 = ecc
      go to 30
*
*       Evaluate circularization time (in units of 10**6 yrs).
   10 tc = tstar*(1.0 - ff)/const
*
*       Activate tidal indicator if TC < 2x10**9 yrs or hyperbolic orbit.
      IF (tc.LT.2000.0.OR.es0.GT.1.0) THEN
          ip = kvec(i1)
          itime = itime + 1
*       Note possibility of counter exceeding the limit.
          IF (itime.GT.2000000000) itime = 0
          IF (icirc.EQ.0.AND.kz(27).EQ.2.AND.itime.LT.100) THEN
              semi = -0.5*body(n+ip)/h(ip)
              WRITE (6,20)  i1, nchaos, es0, rp1, m21, tc, semi, xn
   20         FORMAT (' TCIRC:    I1 NCH E RP M21 TC A n ',
     &                            2i5,f8.4,f8.1,f6.2,1p,2e10.2,0p,f5.1)
          END IF
          icirc = 1
*       Define Roche search indicator for circularized orbit (ECCM1 > 0.002).
          IF (es0.LE.eccm1.AND.kstar(n+ip).EQ.0) THEN
              kstar(n+ip) = 10
          END IF
      ELSE
*       Note ICIRC = -1 for some calls.
          icirc = 0
      END IF
*
   30 RETURN
*
      END