qtides.f

Go to the documentation of this file.

      SUBROUTINE qtides(I1,I2,IM,SEMI,ES0)
*
*
*       Quadrupole and tidal terms for hierarchical binary.
*       ---------------------------------------------------
*
*       Theory of Rosemary Mardling, Ap. J. XX, YYY, 1995.
*       @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
*
      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/tidal/  cq(2),ct(2),cgr,dedt
      REAL*8  ww(3),qq(3),w(2),q(2),at0(2),m21,wg(2),qg(2),bod(2),
     &        a(2),b(2),c(6),qd(2),k1,k2
      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/
      SAVE itime
      DATA itime /0/
*
*
*       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) = cm(1,im)
          bod(2) = cm(2,im)
      ELSE
          j1 = i2
          j2 = i1
          bod(1) = cm(2,im)
          bod(2) = cm(1,im)
      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) THEN
              bodi = bod(k)
              CALL giant3(ik,bodi,wg,qg,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.6) 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)
          qd(k) = ql
    5 CONTINUE
*
*       Form mass, radius & pericentre ratio.
      rp = semi*(1.0 - es0)
      IF (radius(i1).GE.radius(i2)) THEN
          m21 = cm(2,im)/cm(1,im)
          r21 = radius(i2)/radius(i1)
      rp1 = rp/radius(i1)
      ELSE
      m21 = cm(1,im)/cm(2,im)
          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
*
*       Define mass ratio and apsidal motion constants.
      qin = bod(2)/bod(1)
      k1 = 0.2*twopi*q(1)**2/w(1)
      k2 = 0.2*twopi*q(2)**2/w(2)
*
*       Absorb constants to simplify derivatives (note: n^2=(1+qin)/a^3).
      cq(1) = k1*qin*sqrt(bod(1)*(1.0+qin))*radius(j1)**5
      cq(2) = k2*sqrt(bod(1)*(1.0+qin))*radius(j2)**5/qin
*       Note that BOD(1)*(1+qin) is used instead of BODY(I1) for clarity.
      ct(1) = -(1.0/twopi)*k1*qin*(1.0+qin)/(qd(1)*sqrt(w(1)))
      ct(2) = -(1.0/twopi)*k2*(1.0+qin)/qin**2/(qd(2)*sqrt(w(2)))
      ct(1) = ct(1)*sqrt(bod(1))*radius(j1)**6.5
      ct(2) = ct(2)*sqrt(bod(2))*radius(j2)**6.5
*
*       Form gravitational precession constant (Holman et al, Nature 1998).
      tgr = 3.4d+07*yrs*rau/(body(i1)*smu*sqrt(body(i1)))
*       Note: P_{gr} = 3.4E+07*(1-e^2)*P_{yr}*a_{AU}/m_{sun}.
      tgr = tgr/(1.0d+06*tstar)
*       Convert to angular velocity from N-body time units.
      cgr = twopi/tgr
*
      itime = itime + 1
      IF (itime.LE.1) THEN
      tf = time + toff + (1.0 - ff)/const
      WRITE (6,20)  cq, ct, edot, tf, cgr
   20 FORMAT (' QTIDES    CQ CT EDT TF CGR  ',1p,7e10.2)
      WRITE (6,30)  k1,k2,qd,w
   30 FORMAT (' k1 k2 QD W   ',2f10.5,1p,4e10.2)
      tb = yrs*semi*sqrt(semi/body(i1))
      tgr = tgr*(1.0 - es0**2)*semi**2.5
      tgr = 1.0d+06*tstar*tgr
      WRITE (6,40)  tb, tgr, semi*rau, es0
   40 FORMAT (' TB(yr) TGR A(AU) E  ',1p,3e10.2,0p,f8.4)
      CALL flush(6)
      END IF
*
      RETURN
*
      END