Nbody6
qtides.f
Go to the documentation of this file.
1  SUBROUTINE qtides(I1,I2,IM,SEMI,ES0)
2 *
3 *
4 * Quadrupole and tidal terms for hierarchical binary.
5 * ---------------------------------------------------
6 *
7 * Theory of Rosemary Mardling, Ap. J. XX, YYY, 1995.
8 * @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
9 *
10  include 'common6.h'
11  common/binary/ cm(4,mmax),xrel(3,mmax),vrel(3,mmax),
12  & hm(mmax),um(4,mmax),umdot(4,mmax),tmdis(mmax),
13  & namem(mmax),nameg(mmax),kstarm(mmax),iflagm(mmax)
14  common/tidal/ cq(2),ct(2),cgr,dedt
15  REAL*8 ww(3),qq(3),w(2),q(2),at0(2),m21,wg(2),qg(2),bod(2),
16  & a(2),b(2),c(6),qd(2),k1,k2
17  DATA ww /2.119,3.113,8.175/
18  DATA qq /0.4909,0.4219,0.2372/
19  DATA a /6.306505,-7.297806/
20  DATA b /32.17211,13.01598/
21  DATA c /5.101417,24.71539,-9.627739,1.733964,
22  & -2.314374,-4.127795/
23  SAVE itime
24  DATA itime /0/
25 *
26 *
27 * Specify index J1 as biggest radius to be used with AT0(1).
28  IF (radius(i1).GE.radius(i2)) THEN
29  j1 = i1
30  j2 = i2
31  bod(1) = cm(1,im)
32  bod(2) = cm(2,im)
33  ELSE
34  j1 = i2
35  j2 = i1
36  bod(1) = cm(2,im)
37  bod(2) = cm(1,im)
38  END IF
39 *
40 * Define oscillation period (dimensionless time) and damping constants.
41  DO 5 k = 1,2
42  IF (k.EQ.1) THEN
43  ik = j1
44  ELSE
45  ik = j2
46  END IF
47 * Specify polytropic index for each star (n = 3, 2 or 3/2).
48  IF (kstar(ik).EQ.3.OR.kstar(ik).EQ.5) THEN
49  bodi = bod(k)
50  CALL giant3(ik,bodi,wg,qg,xn,ql)
51  w(k) = wg(1)
52  q(k) = qg(1)
53  ELSE
54  ql = 1.0d+04
55  ip = 3
56  IF (kstar(ik).GE.3) ip = 2
57  IF (kstar(ik).EQ.4.OR.kstar(ik).EQ.6) ip = 3
58  IF (kstar(ik).EQ.0) ip = 1
59  w(k) = ww(ip)
60  q(k) = qq(ip)
61  END IF
62  tl = twopi*radius(ik)*sqrt(radius(ik)/bod(k)/w(k))
63  at0(k) = 1.0/(ql*tl)
64  qd(k) = ql
65  5 CONTINUE
66 *
67 * Form mass, radius & pericentre ratio.
68  rp = semi*(1.0 - es0)
69  IF (radius(i1).GE.radius(i2)) THEN
70  m21 = cm(2,im)/cm(1,im)
71  r21 = radius(i2)/radius(i1)
72  rp1 = rp/radius(i1)
73  ELSE
74  m21 = cm(1,im)/cm(2,im)
75  r21 = radius(i1)/radius(i2)
76  rp1 = rp/radius(i2)
77  END IF
78 *
79 * Evaluate damping coefficient.
80  rr = rp1*(1.0 + es0)
81  const = 2.0*(at0(1)*(q(1)/w(1))**2*(1.0 + m21)*m21 +
82  & at0(2)*(q(2)/w(2))**2*((1.0 + m21)/m21**2)*r21**8)/
83  & rr**8
84 *
85 * Form rational function approximation to Hut solution.
86  ff = (( a(2)*es0 + a(1))*es0 + 1.0 )/
87  & (( b(2)*es0 + b(1))*es0 + 1.0 )
88  ff = min(ff,0.999d0)
89 *
90 * Determine eccentricity corresponding to t_{circ} = t_{grow}.
91 * Z = TG*CONST/TSTAR + FF
92 * ECC1 = (-1.0 + C(1)*Z - SQRT(C(2)*Z**2 + C(3)*Z + C(4)))
93 * & /(C(5) + C(6)*Z)
94 *
95 * Evaluate circularization time (in units of 10**6 yrs).
96  tc = tstar*(1.0 - ff)/const
97 *
98 * Obtain (de/dt) due to tidal circularization.
99  fe = 1.0 + 3.75*es0**2 + 1.875*es0**4 + (5.0/64.0)*es0**6
100  fe = (9.0*twopi/10.0)*es0*(1.0 - es0**2)**1.5*fe
101  edot = -const*fe
102 *
103 * Define mass ratio and apsidal motion constants.
104  qin = bod(2)/bod(1)
105  k1 = 0.2*twopi*q(1)**2/w(1)
106  k2 = 0.2*twopi*q(2)**2/w(2)
107 *
108 * Absorb constants to simplify derivatives (note: n^2=(1+qin)/a^3).
109  cq(1) = k1*qin*sqrt(bod(1)*(1.0+qin))*radius(j1)**5
110  cq(2) = k2*sqrt(bod(1)*(1.0+qin))*radius(j2)**5/qin
111 * Note that BOD(1)*(1+qin) is used instead of BODY(I1) for clarity.
112  ct(1) = -(1.0/twopi)*k1*qin*(1.0+qin)/(qd(1)*sqrt(w(1)))
113  ct(2) = -(1.0/twopi)*k2*(1.0+qin)/qin**2/(qd(2)*sqrt(w(2)))
114  ct(1) = ct(1)*sqrt(bod(1))*radius(j1)**6.5
115  ct(2) = ct(2)*sqrt(bod(2))*radius(j2)**6.5
116 *
117 * Form gravitational precession constant (Holman et al, Nature 1998).
118  tgr = 3.4d+07*yrs*rau/(body(i1)*smu*sqrt(body(i1)))
119 * Note: P_{gr} = 3.4E+07*(1-e^2)*P_{yr}*a_{AU}/m_{sun}.
120  tgr = tgr/(1.0d+06*tstar)
121 * Convert to angular velocity from N-body time units.
122  cgr = twopi/tgr
123 *
124  itime = itime + 1
125  IF (itime.LE.1) THEN
126  tf = time + toff + (1.0 - ff)/const
127  WRITE (6,20) cq, ct, edot, tf, cgr
128  20 FORMAT (' QTIDES CQ CT EDT TF CGR ',1p,7e10.2)
129  WRITE (6,30) k1,k2,qd,w
130  30 FORMAT (' k1 k2 QD W ',2f10.5,1p,4e10.2)
131  tb = yrs*semi*sqrt(semi/body(i1))
132  tgr = tgr*(1.0 - es0**2)*semi**2.5
133  tgr = 1.0d+06*tstar*tgr
134  WRITE (6,40) tb, tgr, semi*rau, es0
135  40 FORMAT (' TB(yr) TGR A(AU) E ',1p,3e10.2,0p,f8.4)
136  CALL flush(6)
137  END IF
138 *
139  RETURN
140 *
141  END