hicirc.f

Go to the documentation of this file.

      SUBROUTINE hicirc(RP,ES0,I1,I2,BODYI,TG,TC,ECC1,EDOT,W)
*
*
*       Eccentricity for given circularization time.
*       --------------------------------------------
*
*       Theory of Rosemary Mardling, Ap. J. XX, YYY, 1995.
*       @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
*
      include 'common6.h'
      REAL*8  ww(3),qq(3),w(2),q(2),at0(2),m21,wg(2),qg(2),
     &        bodyi(2),bod(2)
      REAL*8  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/
*
*
*       Specify index J1 as biggest radius to be used with AT0(1).
      IF (radius(i1).GE.radius(i2)) THEN
          j1 = i1
          j2 = i2
          bod(1) = bodyi(1)
          bod(2) = bodyi(2)
      ELSE
          j1 = i2
          j2 = i1
          bod(1) = bodyi(2)
          bod(2) = bodyi(1)
      END IF
*
*       Define oscillation period (dimensionless time) and damping constants.
      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
              bodi = bod(k)
              CALL giant3(ik,bodi,wg,qg,xn,ql)
              w(k) = wg(1)
              q(k) = qg(1)
              CONTINUE
          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)/bod(k)/w(k))
          at0(k) = 1.0/(ql*tl)
    5 CONTINUE
*
*       Form mass, radius & pericentre ratio (note BODYI(2) is a ghost).
      IF (radius(i1).GE.radius(i2)) THEN
          m21 = bodyi(2)/bodyi(1)
          r21 = radius(i2)/radius(i1)
      rp1 = rp/radius(i1)
      ELSE
      m21 = bodyi(1)/bodyi(2)
          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
*
*   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)
*
*       Determine eccentricity corresponding to t_{circ} = t_{grow}.
      z = tg*const/tstar + ff
      ecc1 = (-1.0 + c(1)*z - sqrt(c(2)*z**2 + c(3)*z + c(4)))
     &                                        /(c(5) + c(6)*z)
*
*       Evaluate circularization time (in units of 10**6 yrs).
      tc = tstar*(1.0 - ff)/const
*
*       Obtain (de/dt) due to tidal circularization.
      fe = 1.0 + 3.75*es0**2 + 1.875*es0**4 + (5.0/64.0)*es0**6
      fe = (9.0*twopi/10.0)*es0*(1.0 - es0**2)**1.5*fe
      edot = -const*fe
*
      RETURN
*
      END