Nbody6
 All Files Functions Variables
kstide.f
Go to the documentation of this file.
1  SUBROUTINE kstide(IPAIR,QPERI)
2 *
3 *
4 * Tidal or GR interaction of KS pair.
5 * -----------------------------------
6 *
7  include 'common6.h'
8  REAL*8 de(2)
9  CHARACTER*8 which1
10  INTEGER is(2)
11  SAVE ttide,ione
12  DATA eccm,eccm2,ttide,ione /0.002,0.00000399,0.0d0,0/
13  SAVE sum
14  DATA sum /0.0d0/
15 *
16 *
17 * Define indices of KS components.
18  i1 = 2*ipair - 1
19  i2 = i1 + 1
20 *
21 * Set c.m. index & reduced mass and increase event counter.
22  i = n + ipair
23  zmu = body(i1)*body(i2)/body(i)
24  ndiss = ndiss + 1
25  rks = r(ipair)
26 *
27 * Form current semi-major axis & eccentricity (not at peri).
28  semi = -0.5d0*body(i)/h(ipair)
29  ecc2 = (1.0 - r(ipair)/semi)**2 + tdot2(ipair)**2/(semi*body(i))
30  ecc = sqrt(ecc2)
31  am0 = semi*(1.0d0 - ecc**2)
32  peri = semi*(1.0d0 - ecc)
33  hi = h(ipair)
34 *
35 * Distinguish between sequential circularization, standard chaos & GR.
36  r1 = max(radius(i1),radius(i2))
37  IF (kz(27).EQ.1) THEN
38  zf = 4.0
39  IF (ecc.GT.0.95) zf = 50.0
40 * Skip if penetration is not significant (prevents frequent calls).
41  IF (abs(qperi - zf*r1).LT.0.01*qperi) go to 50
42  ELSE IF (kz(27).EQ.2) THEN
43  is(1) = kstar(i1)
44  is(2) = kstar(i2)
45 * Obtain kinetic energy loss due to tidal interaction (DE > 0).
46 * CALL TIDES(QPERI,BODY(I1),BODY(I2),RADIUS(I1),RADIUS(I2),IS,DE)
47 *
48 * Evaluate chaos parameters and check for disruption.
49  CALL chaos(ipair,i1,i2,qperi,ecc,is,zmu,rks,semi1,ecc1,idis)
50  IF (idis.EQ.1) go to 45
51 *
52 * Update KS variables for CHAOS or TERMINATED CHAOS (H > 0).
53  IF (kstar(i).EQ.-1.OR.h(ipair).GT.0.0) THEN
54 * Ensure rectification to determine true pericentre (bug fix 07/08).
55  CALL ksrect(ipair)
56  CALL ksperi(ipair)
57  go to 10
58  END IF
59 * Skip on WIDE CHAOS or NEW SPIRAL.
60  go to 45
61  ELSE IF (kz(27).EQ.3) THEN
62  CALL tides3(qperi,body(i1),body(i2),vstar,h(ipair),ecc,de)
63  END IF
64 *
65 * Restore circularization index if needed (exit from CHAIN).
66  IF (ecc.LE.eccm.AND.kstar(i).LT.10) THEN
67  kstar(i) = 10
68  go to 50
69  END IF
70 *
71 * Consider sequential circularization or GR evolution.
72  IF (kz(27).EQ.1) THEN
73 * Suppress the old PT procedure (DH => ECC from AM0 = const).
74 * ECC2 = ECC**2 + 2.0D0*AM0*DH/BODY(I)
75 * ECC2 = MAX(ECC2,ECCM2)
76 * Accept circularized orbit if ACIRC < 4*R1 (use maximum radius).
77  am0 = semi*(1.0 - ecc**2)
78  ecc1 = sqrt(eccm2)
79  acirc = am0/(1.0 - eccm2)
80  IF (acirc.LT.zf*r1) THEN
81  semi1 = acirc
82  ELSE
83 * Obtain E1 from A1*(1 - E1**2) = AM0 using A1*(1 - E1) = 4*R1.
84  ecc1 = am0/(zf*r1) - 1.0
85  ecc1 = max(ecc1,eccm)
86 * Set new semi-major axis from angular momentum conservation.
87  semi1 = am0/(1.0 - ecc1**2)
88  END IF
89 * Form the corresponding energy change.
90  dh = 0.5*body(i)*(1.0/semi - 1.0/semi1)
91  de(1) = -zmu*dh
92  de(2) = 0.0
93  ELSE
94 * Include safety check on energy loss to prevent new SEMI < R.
95  dh = -(de(1) + de(2))/zmu
96  IF (h(ipair) + dh.LT.-0.5*body(i)/r(ipair)) THEN
97  dh = -0.5*body(i)/r(ipair) - h(ipair)
98  de(1) = -zmu*dh
99  de(2) = 0.0
100  END IF
101  semi1 = -0.5*body(i)/(h(ipair) + dh)
102  ecc1 = 1.0 - peri/semi1
103  ecc1 = max(ecc1,eccm)
104  END IF
105 *
106 * Skip possible hyperbolic case.
107  IF (h(ipair) + dh.GT.0.0) go to 50
108 *
109 * Update total energy loss (and H after obtaining peri).
110  ecoll = ecoll + (de(1) + de(2))
111  e(10) = e(10) + (de(1) + de(2))
112  egrav = egrav + (de(1) + de(2))
113 * Determine pericentre variables U & UDOT by backwards integration.
114  CALL ksperi(ipair)
115  h(ipair) = h(ipair) + dh
116 *
117 * Print first and last energy change and update indicator.
118  IF (kz(27).EQ.1.AND.(kstar(i).EQ.0.OR.kstar(i).EQ.9)) THEN
119  p = days*semi1*sqrt(semi1/body(i))
120  IF (kstar(i).EQ.0.AND.ecc1.GT.eccm) THEN
121  which1 = 'NEW CIRC'
122  kstar(i) = 9
123  WRITE (6,8) which1, name(i1), name(i2), kstar(i1),
124  & kstar(i2), ttot, ecc, ecc1, p, semi1, r1
125  ELSE IF (ecc1.LE.eccm) THEN
126  which1 = 'END CIRC'
127  kstar(i) = 10
128  ncirc = ncirc + 1
129  WRITE (6,8) which1, name(i1), name(i2), kstar(i1),
130  & kstar(i2), ttot, ecc, ecc1, p, semi1, r1
131  END IF
132  8 FORMAT (' ',a8,' NAM K* T E0 EF P AF R* ',
133  & 2i6,2i4,f9.2,2f8.3,f7.1,1p,2e10.2)
134  END IF
135 *
136 * Set new pericentre.
137  10 peri1 = semi1*(1.0d0 - ecc1)
138 *
139 * Form KS coordinate scaling factor from pericentre ratio.
140  c1 = sqrt(peri1/peri)
141 *
142 * Specify KS velocity scaling (conserved J, chaos or GR treatment).
143  IF (kz(27).EQ.1) THEN
144  c2 = 1.0/c1
145  ELSE
146 * Note that PERI should not change in GR case (hence same C2).
147  c2 = sqrt((body(i) + h(ipair)*peri1)/(body(i) + hi*peri))
148  END IF
149 *
150 * See whether circular orbit condition applies.
151 * IF (ECC1.LE.ECCM) THEN
152 * AM = SEMI1*(1.0D0 - ECC1**2)
153 * C2 = SQRT(AM/AM0)/C1
154 * END IF
155 *
156 * Transform KS variables to yield the prescribed elements.
157  r(ipair) = 0.0d0
158  DO 15 k = 1,4
159  u(k,ipair) = c1*u(k,ipair)
160  udot(k,ipair) = c2*udot(k,ipair)
161  u0(k,ipair) = u(k,ipair)
162  r(ipair) = r(ipair) + u(k,ipair)**2
163  15 CONTINUE
164 *
165 * Form new perturber list after significant energy loss.
166  np0 = list(1,i1)
167  IF (abs(semi1/semi).LT.0.5) THEN
168  CALL kslist(ipair)
169 * Ensure that single perturber differs from #I itself.
170  IF (np0.EQ.0.AND.list(1,i1).GT.0) THEN
171  IF (list(2,i1).EQ.i) THEN
172  list(2,i1) = ifirst
173  END IF
174  END IF
175  END IF
176 *
177 * Ensure chaos treatment at every pericentre by perturbed motion.
178  IF (kz(27).EQ.2.AND.list(1,i1).EQ.0) THEN
179  list(1,i1) = 1
180  list(2,i1) = n
181  np0 = 1
182  END IF
183 *
184 * Copy perturber list to JPERT and predict perturbers.
185  nnb = list(1,i1)
186  DO 16 l = 1,nnb
187  jpert(l) = list(l+1,i1)
188  j = list(l+1,i1)
189  CALL xvpred(j,0)
190  16 CONTINUE
191 *
192 * Evaluate potential energy of perturbers.
193  jlist(1) = i1
194  jlist(2) = i2
195  CALL nbpot(2,nnb,pot1)
196 *
197 * Re-initialize KS polynomials at pericentre for perturbed case.
198  t0(i1) = time
199  IF (np0.GT.0.OR.ecc.LE.eccm) THEN
200  CALL resolv(ipair,1)
201  imod = kslow(ipair)
202  CALL kspoly(ipair,imod)
203  END IF
204 *
205 * Obtain potential energy wrt new binary and apply tidal correction.
206  CALL nbpot(2,nnb,pot2)
207  dp = pot2 - pot1
208  ecoll = ecoll + dp
209 *
210 * Produce diagnostic output for interesting new case (first time).
211  IF (ecc.GT.0.99.AND.abs(ecc - ecc1).GT.0.01.AND.ione.EQ.0) THEN
212  WRITE (6,20) name(i1), name(i2), semi1, ecc, ecc1, hi,
213  & qperi, dh, dp
214  20 FORMAT (' NEW KSTIDE NAM AF E0 EF HI QP DH DP ',
215  & 2i5,1p,e10.2,0p,2f8.3,f9.1,1p,3e10.2)
216  ttide = time + toff
217  ione = ione + 1
218  END IF
219  IF (time + toff.GT.ttide + dtadj) ione = 0
220 *
221 * Check for hierarchical configuration with eccentric inner binary.
222  IF (kz(27).EQ.2.AND.ecc.GT.0.95.AND.hi.LT.0.0) THEN
223  np1 = list(1,i1) + 1
224  DO 30 l = 2,np1
225  j = list(l,i1)
226  rij2 = 0.0
227  vij2 = 0.0
228  rdot = 0.0
229  DO 25 k = 1,3
230  rij2 = rij2 + (x(k,i) - x(k,j))**2
231  vij2 = vij2 + (xdot(k,i) - xdot(k,j))**2
232  rdot = rdot + (x(k,i) - x(k,j))*(xdot(k,i) -xdot(k,j))
233  25 CONTINUE
234  rip = sqrt(rij2)
235  a1 = 2.0/rip - vij2/(body(i) + body(j))
236  a1 = 1.0/a1
237  IF (1.0/a1.GT.0.1/semi) THEN
238  ecc2 = (1.0 - rip/a1)**2 +
239  & rdot**2/(a1*(body(i) + body(j)))
240  rp = a1*(1.0 - sqrt(ecc2))
241  ra = semi*(1.0 + ecc)
242  sr = rp/ra
243  icirc = -1
244  jcomp = j
245  CALL induce(ipair,emax,emin,icirc,tc,angle,tg,edav)
246  WRITE (6,28) name(j), h(ipair), semi, a1, rp, edav,
247  & sqrt(ecc2), emax, sr
248  28 FORMAT (' HIERARCHY NMJ H A0 A1 RP EDAV E1 EX SR ',
249  & i7,f7.0,1p,4e9.1,0p,2f8.4,f6.1)
250  END IF
251  30 CONTINUE
252  END IF
253 *
254 * Set one unperturbed period for small apocentre perturbation (#27=1).
255  ga = gamma(ipair)*(semi1*(1.0 + ecc1)/r(ipair))**3
256  IF (ga.LT.gmin.AND.kz(27).EQ.1) THEN
257  step(i1) = twopi*semi1*sqrt(semi1/body(i))
258  list(1,i1) = 0
259  END IF
260 *
261 * Ensure T'' = 0 for pericentre test in KSINT & dissipation in UNPERT.
262  IF (tdot2(ipair).LT.0.0d0) THEN
263  tdot2(ipair) = 0.0d0
264  END IF
265 *
266 * Count any hyperbolic captures.
267  IF (semi.LT.0.0.AND.semi1.GT.0.0) THEN
268  ntide = ntide + 1
269  qps = qperi/r1
270  WRITE (6,35) time+toff, name(i1), name(i2), ecc, ecc1, qps,
271  & semi1
272  35 FORMAT (' NEW CAPTURE T NM E EF QP/R* A1 ',
273  & f9.2,2i6,2f9.4,1p,2e10.2)
274  END IF
275 *
276 * Record diagnostics for new synchronous orbit and activate indicator.
277  IF (ecc.GT.eccm.AND.ecc1.LT.eccm.AND.kz(27).LE.2) THEN
278  nsync = nsync + 1
279  esync = esync + zmu*h(ipair)
280  kstar(i) = 10
281  WRITE (6,38) name(i1), name(i2), kstar(i1), kstar(i2), semi1,
282  & ecc, ecc1, hi, qperi, r1
283  38 FORMAT (' CIRCULARIZED NAM K* AF E0 EF HI QP R* ',
284  & 2i6,2i3,1p,e10.2,0p,2f8.3,f9.1,1p,2e10.2)
285 * Check time interval until Roche overflow.
286  IF (kz(34).GT.0) THEN
287  CALL trflow(ipair,dtr)
288  tev(i) = time + dtr
289 * Update TEV but leave ROCHE call for treatment in MDOT (16/08/2006).
290 * IF (DTR.LT.STEP(I1)) THEN
291 * CALL ROCHE(IPAIR)
292 * GO TO 50
293 * END IF
294  END IF
295  END IF
296 *
297 * See whether a low-eccentricity synchronous state has been reached.
298  rcoll = 0.75*(radius(i1) + radius(i2))
299  IF (abs(semi1).LT.1.5*rcoll.AND.ecc.LT.eccm.AND.
300  & kstar(i).LT.10) THEN
301  kstar(i) = 10
302  WRITE (6,40) ecc1, semi1, r(ipair), rcoll, r1
303  40 FORMAT (' INACTIVE PHASE E A R RC R* ',f7.3,1p,4e10.2)
304 * Check time interval until Roche overflow.
305  IF (kz(34).GT.0) THEN
306  CALL trflow(ipair,dtr)
307  tev(i) = time + dtr
308 * Update TEV but leave ROCHE call for treatment in MDOT (16/08/2006).
309 * IF (DTR.LT.STEP(I1)) THEN
310 * CALL ROCHE(IPAIR)
311 * GO TO 50
312 * END IF
313  END IF
314  END IF
315 *
316 * Include warning in case of eccentricity increase (PT only).
317  ecc2 = 1.0 - r(ipair)/semi1
318  IF (kz(27).EQ.1.AND.ecc2.GT.max(ecc,eccm)+1.0d-04) THEN
319  WRITE (6,42) ttot, ipair, ecc2, ecc, r(ipair), semi1
320  42 FORMAT (' WARNING! E > E0 T IP E E0 R A ',
321  & f10.4,i5,2f8.4,1p,2e10.2)
322  END IF
323 *
324 * Reduce radius by DR/R to delay dissipation for small energy loss.
325  IF (kz(27).EQ.1.AND.de(1)+de(2).LT.1.0d-07*zmu*abs(hi)) THEN
326  j1 = i1
327  IF (radius(i2).GT.radius(i1)) j1 = i2
328 *** IF (TEV(J1) - TIME.GT.0.01*TEV(J1)) TEV(J1) = TIME
329  dr = (0.99*4.0*radius(j1) - qperi)/qperi
330  yf = max(abs(dr),0.01d0)
331  radius(j1) = (1.0 - min(yf,0.1d0))*radius(j1)
332  de1 = (de(1) + de(2))/(zmu*abs(h(ipair)))
333  WRITE (6,44) ttot, kstar(j1), h(ipair), qperi, dr, de1
334  44 FORMAT (' REDUCED RADIUS T K* H QP DR/R DE/E ',
335  & f9.2,i3,f10.1,1p,3e9.1)
336  END IF
337 *
338 * Combine the two stars inelastically in case of chaos disruption.
339  45 IF (kz(27).EQ.2.AND.idis.GT.0) THEN
340  WRITE (6,48) ttot, ipair, list(1,i1), ecc, semi, qperi,
341  & radius(i1), radius(i2)
342  48 FORMAT (' DISRUPTED CHAOS T KS NP E A QP R* ',
343  & f9.2,i6,i4,f8.4,1p,4e10.2)
344  CALL ksperi(ipair)
345  CALL xvpred(i,0)
346  kspair = ipair
347  iqcoll = 1
348  CALL cmbody(r(ipair),2)
349  END IF
350 *
351  50 RETURN
352 *
353  END