Nbody6
 All Files Functions Variables
spiral.f
Go to the documentation of this file.
1  SUBROUTINE spiral(IPAIR)
2 *
3 *
4 * Tidal circularization of binary orbit.
5 * --------------------------------------
6 *
7 * Rational function approximations of solution to the Hut
8 * evolution equation with spin. Ref: A & A 99, 126, eqn (A15).
9 * Developed by Rosemary Mardling (31/1/97 & 25/5/00). See 12/red/36.
10 *
11  include 'common6.h'
12  common/modes/ eb0(ntmax),zj0(ntmax),ecrit(ntmax),ar(ntmax),
13  & br(ntmax),eosc(4,ntmax),edec(ntmax),tosc(ntmax),
14  & rp(ntmax),es(ntmax),cm(2,ntmax),iosc(ntmax),
15  & namec(ntmax)
16  common/slow0/ range,islow(10)
17  common/spins/angmom0,rg2(2),m21,r21,semi0,c1,c2,c3,c4,c5,semi
18  REAL*8 ww(3),qq(3),w(2),q(2),at0(2),m21,wg(2),qg(2),wscale(2),
19  & qscale(2),a(2),b(2),c(6),meanmotion,rj(2),rol(2)
20  REAL*8 m0,mc1,mc2,corerd
21  EXTERNAL corerd
22  DATA ww /2.119,3.113,8.175/
23  DATA qq /0.4909,0.4219,0.2372/
24  DATA a /6.306505,-7.297806/
25  DATA b /32.17211,13.01598/
26  DATA c /5.101417,24.71539,-9.627739,1.733964,
27  & -2.314374,-4.127795/
28  DATA eccm /0.002/
29  SAVE ione,iwarn
30  DATA ione,iwarn /0,0/
31 *
32 *
33 * Check for just table updating at KS termination.
34  irem = 0
35  IF (ipair.LT.0) THEN
36  i = -ipair
37  DO 1 k = 1,nchaos
38  IF (namec(k).EQ.name(i)) irem = k
39  1 CONTINUE
40  ione = 0
41  IF (irem.GT.0) go to 30
42 * Include case of removing a specified chaos index (I <= NCHAOS).
43  IF (i.LE.nchaos) THEN
44  irem = i
45  go to 30
46  END IF
47  go to 100
48  END IF
49 *
50 * Define c.m. & KS indices and search current names for chaos index.
51  i = n + ipair
52  i1 = 2*ipair - 1
53  i2 = i1 + 1
54  nspir = nspir + 1
55  ic = 0
56  DO 2 k = 1,nchaos
57  IF (namec(k).EQ.name(i)) ic = k
58  2 CONTINUE
59 *
60 * Try repair procedure if NCHAOS reduced to zero or NAMEC missing.
61  IF (nchaos.EQ.0.OR.ic.EQ.0) THEN
62  ic = 1
63  nchaos = nchaos + 1
64  semi = -0.5*body(i)/h(ipair)
65  ecc2 = (1.0-r(ipair)/semi)**2 + tdot2(ipair)**2/(body(i)*semi)
66  ecc = sqrt(ecc2)
67 * Re-initialize basic elements and NAMEC.
68  es(ic) = ecc
69  rp(ic) = semi*(1.0 - ecc)
70  namec(ic) = name(i)
71  WRITE (6,4) name(i1), name(i), list(1,i1), ecc, r(ipair)
72  4 FORMAT (' SPIRAL REPAIR NM NMI NP E RP ',
73  & 2i7,i4,f8.4,1p,e10.2)
74  END IF
75 *
76 * Set non-zero index on call from CHRECT (KSTAR = -KSTAR) for CE skip.
77  istar = 0
78  IF (kstar(i).GT.0) THEN
79  kstar(i) = -kstar(i)
80  istar = 1
81  END IF
82 *
83 * Include case of chain chaos without identified NAMEC.
84  IF (ic.EQ.0.AND.nchaos.GT.0) THEN
85  WRITE (6,3) nchaos, kstar(i), ipair, name(i1), name(i2),
86  & step(2*ipair-1), step(i)
87  3 FORMAT (' WARNING! SPIRAL NCH K* KS NAME STEP1 STEPI',
88  & 3i4,2i6,1p,2e10.2)
89 * Search for component names (set in CHAOS2).
90  DO 5 k = 1,nchaos
91  IF (namec(k).EQ.name(i1).OR.namec(k).EQ.name(i2)) ic = k
92  5 CONTINUE
93 *
94 * Include safety check just in case.
95  IF (ic.EQ.0) THEN
96  WRITE (6,6) name(i), kstar(i), (namec(k),k=1,nchaos)
97  6 FORMAT (' DANGER! SPIRAL NAMI K* NAMEC ',15i6)
98  stop
99  END IF
100  namec(ic) = name(i)
101  END IF
102 *
103 * Copy spiral parameters and set semi-major axis & period.
104  time0 = tosc(ic)
105  rp0 = rp(ic)
106  es0 = es(ic)
107  hi = h(ipair)
108  semi0 = -0.5*body(i)/h(ipair)
109  tk = semi0*sqrt(semi0/body(i))
110 *
111 * Use actual elements for perturbed binary.
112  IF (list(1,i1).GT.0) THEN
113  ecc2 = (1.0 - r(ipair)/semi0)**2 +
114  & tdot2(ipair)**2/(body(i)*semi0)
115  ecc0 = sqrt(ecc2)
116  es0 = ecc0
117  rp0 = semi0*(1.0 - ecc0)
118  rp(ic) = rp0
119  END IF
120 *
121 * Specify index J1 as biggest radius to be used with AT0(1).
122  IF (radius(i1).GE.radius(i2)) THEN
123  j1 = i1
124  j2 = i2
125  ELSE
126  j1 = i2
127  j2 = i1
128  END IF
129 *
130 * Define oscillation period (dimensionless time) and damping constants.
131  xn = 0.0
132  DO 10 k = 1,2
133  ik = i1 + k - 1
134  IF (k.EQ.1) THEN
135  ik = j1
136  ELSE
137  ik = j2
138  END IF
139 * Specify polytropic index for each star (n = 3, 2 or 3/2).
140  IF (kstar(ik).EQ.3.OR.kstar(ik).EQ.5.OR.
141  & kstar(ik).EQ.6.OR.kstar(ik).EQ.9) THEN
142  CALL giant(ipair,ik,wg,qg,wscale,qscale,xn,ql)
143  w(k) = wg(1)
144  q(k) = qg(1)
145 * Note: rg2 should really be given by k2 in HRDIAG.
146  rg2(k)= 0.1*(1.0 - cm(k,ic)/body(ik))
147  ELSE
148  ql = 1.0e+04
149  ip = 3
150  IF (kstar(ik).GE.3) ip = 2
151  IF (kstar(ik).EQ.4.OR.kstar(ik).EQ.7) ip = 3
152  IF (kstar(ik).EQ.8) ip = 3
153  IF (kstar(ik).EQ.0) ip = 1
154  w(k) = ww(ip)
155  q(k) = qq(ip)
156  rg2(k)= 0.21
157  IF (kstar(ik).LE.2.OR.kstar(ik).EQ.7) rg2(k) = 0.1
158  IF (kstar(ik).EQ.4) rg2(k)= 0.1*(1.0 - cm(k,ic)/body(ik))
159  END IF
160  tl = twopi*radius(ik)*sqrt(radius(ik)/body(ik)/w(k))
161  at0(k) = 1.0/(ql*tl)
162  10 CONTINUE
163 *
164 * Form mass, radius & pericentre ratio.
165  IF (radius(i1).GE.radius(i2)) THEN
166  m21 = body(i2)/body(i1)
167  r21 = radius(i2)/radius(i1)
168  rp1 = rp(ic)/radius(i1)
169  rad = radius(i1)
170  ELSE
171  m21 = body(i1)/body(i2)
172  r21 = radius(i1)/radius(i2)
173  rp1 = rp(ic)/radius(i2)
174  rad = radius(i2)
175  END IF
176 *
177 * Define initial angular momentum from the scaled semi-major axis.
178  semi0 = rp1/(1.0 - es0)
179 *
180 * Form the initial mean motion in N-body units.
181  meanmotion = sqrt((body(i1)+body(i2))/(rad*semi0)**3)
182 *
183 * Convert from angular momentum to omega (denoted spin1 & spin2).
184  IF (kstar(j1).LE.2.OR.(kstar(j1).GE.7.AND.kstar(j1).NE.9)) THEN
185  spin1 = spin(j1)/(rg2(1)*body(j1)*radius(j1)**2)
186  ELSE
187  kw = kstar(j1)
188  m0 = body0(j1)*smu
189  mc1 = cm(1,ic)*smu
190  IF (mc1.LE.0.0d0.OR.mc1.GT.m0) THEN
191  mc1 = 0.3 + 0.1*float(kw - 3)
192  IF(kw.EQ.9) mc1 = min(0.3d0,0.95*m0)
193  cm(1,ic) = mc1/zmbar
194  END IF
195  zdum = 2.0d0
196  rc1 = corerd(kw,mc1,m0,zdum)/su
197  spin1 = spin(j1)/(rg2(1)*body(j1)*radius(j1)**2 +
198  & 0.21*mc1/smu*rc1**2)
199  END IF
200  IF (kstar(j2).LE.2.OR.(kstar(j2).GE.7.AND.kstar(j2).NE.9)) THEN
201  spin2 = spin(j2)/(rg2(2)*body(j2)*radius(j2)**2)
202 * Produce diagnostic information for recent WD.
203  IF (kstar(j2).GE.10.AND.kstar(j2).LE.12.AND.
204  & time.LT.tev0(j2) + 0.01) THEN
205  WRITE (95,780) ttot, name(j2), name(j1), kstar(j1),
206  & spin1, spin2, meanmotion, spin(j2)
207  780 FORMAT (' WD SPIN T NM K*1 ROT1 ROT2 <n> S2 ',
208  & f8.1,2i6,i4,1p,4e10.2)
209  CALL flush(95)
210  END IF
211  ELSE
212  kw = kstar(j2)
213  m0 = body0(j2)*smu
214  mc2 = cm(2,ic)*smu
215  IF (mc2.LE.1.0d-10.OR.mc2.GT.m0) THEN
216  mc2 = 0.3 + 0.1*float(kw - 3)
217  IF(kw.EQ.9) mc2 = min(0.3d0,0.95*m0)
218  cm(2,ic) = mc2/zmbar
219  END IF
220  zdum = 2.0d0
221  rc2 = corerd(kw,mc2,m0,zdum)/su
222  spin2 = spin(j2)/(rg2(2)*body(j2)*radius(j2)**2 +
223  & 0.21*mc2/smu*rc2**2)
224  END IF
225 *
226  IF (ione.EQ.0) THEN
227  ione = ione + 1
228  IF (abs(spin1).GT.0.d0.AND.abs(spin2).GT.0.d0) THEN
229  ts = 1.0d+06*365.0*twopi*tstar
230  WRITE (6,790) ts/spin1, ts/spin2, ts/meanmotion,
231  & days*tk, radius(j1)*su, radius(j2)*su
232  790 FORMAT (' BEGIN TROT P TK R* ',1p,2e9.1,0p,4f9.2)
233  END IF
234  END IF
235 *
236 * Scale the spins by mean motion and define angular momentum.
237  spin10=spin1/meanmotion
238  spin20=spin2/meanmotion
239  angmom0=(m21/(1+m21))*semi0**2*sqrt(1-es0**2)+rg2(1)*spin10+
240  & m21*r21**2*rg2(2)*spin20
241 *
242 * Evaluate damping coefficients (Mardling & SJA, M.N. 321, 398, 2001).
243  cf = 54.0*twopi/5.0
244  c1 = cf*(at0(1)*(q(1)/w(1))**2*(1.0 + m21)*m21)/semi0**8
245  c2 = cf*(at0(2)*(q(2)/w(2))**2*((1.0 + m21)/m21**2)*r21**8)/
246  & semi0**8
247  c3 = (cf/9.0)*(at0(1)*(q(1)/w(1))**2*m21**2)/rg2(1)/semi0**6
248  c4 = (cf/9.0)*(at0(2)*(q(2)/w(2))**2/m21**2)*r21**6/rg2(2)/
249  & semi0**6
250 *
251 * Adopt WD scaling for any NS or BH to avoid numerical problem.
252  IF (kstar(i1).GE.13) THEN
253  c1 = 1.0d-04*c1
254  END IF
255  IF (kstar(i2).GE.13) THEN
256  c2 = 1.0d-04*c2
257  END IF
258 *
259  IF (time - time0.LE.0.0d0) go to 100
260 *
261 * Obtain dominant terms of the derivatives (cf. routine HUT/DERIV2).
262  fac=1-es0**2
263  udot1=-(es0/fac**6.5)*(c1 + c2)
264  udot2=(1.0/fac)**6*c3
265  udot3=(1.0/fac)**6*c4
266 *
267 * Choose the step from smallest time-scale (! time0 < time possible).
268  taux = min(abs(1.0/udot1),abs(1.0/udot2),abs(1.0/udot3))
269  nstep = 1 + 100.0*sqrt(abs(time - time0)/taux)
270  nstep = min(nstep,100)
271 * Include extra steps for long intervals and/or AGB type 5 or 6.
272  IF (time - time0.GT.0.1) nstep = nstep + 20
273  IF (kstar(j1).EQ.5.OR.kstar(j1).EQ.6) nstep = nstep + 10
274 * Evaluate the equilibrium angular velocity.
275  e2 = es0**2
276  f2 = ((0.3125*e2 + 5.625)*e2 + 7.5)*e2 + 1.0
277  f5 = (0.375*e2 + 3.0)*e2 + 1.0
278  omeq = f2*meanmotion/(f5*fac**1.5)
279 * Increase number of steps on slow primary rotation.
280  IF (omeq - spin1.GT.0.2*omeq) nstep = nstep + 10
281  dtau1 = abs(time - time0)/float(nstep)
282 *
283 * Check circularization time after merger or large interval.
284  IF (iphase.EQ.7.OR.time - time0.GT.1.0) THEN
285  icirc = -1
286  a0 = -0.5*body(i)/h(ipair)
287  qperi = a0*(1.0 - es0)
288  CALL tcirc(qperi,es0,i1,i2,icirc,tc)
289  dt = min(0.1d0*tc,1.0d0)
290 * Ensure careful integration with DT=0.1*TC for rapid evolution.
291  IF (tc.LT.10.0.AND.(kstar(j1).GT.1.OR.ecc.GT.0.1)) THEN
292  nstep = 500*(1.0 + 100.0/tc)
293  nstep = min(nstep,5000)
294  WRITE (6,12) name(j1), nstep, es0, qperi, time-time0, tc
295  12 FORMAT (' SPIRAL RESET NM # E QP T-T0 TC ',
296  & 2i6,f7.3,1p,3e10.2)
297  END IF
298  dtau1 = abs(dt)/float(nstep)
299  END IF
300 *
301  ione = ione + 1
302  IF (ione.LE.2.OR.mod(ione,10000).EQ.0) THEN
303  WRITE (6,13) list(1,i1), name(i1),es0, omeq, meanmotion,
304  & spin1, spin2
305  13 FORMAT (' SPIRAL NP NM es0 omeq <n> spin ',
306  & i4,i6,f9.5,1p,4e10.2)
307  END IF
308 *
309 * Integrate equations for eccentricity and angular velocities.
310  call hut(es0,spin10,spin20,ecc,spin1,spin2,nstep,dtau1)
311 *
312 * Ensure that no overshooting takes place.
313  IF (spin10.LT.1.0.and.spin1.gt.1.0) THEN
314  IF (mod(iwarn,100).EQ.0) THEN
315  WRITE (66,14) ttot, name(j1), kstar(j1), nstep, ecc,
316  & rp0*su/(1.0-es0), spin1-1.0
317  14 FORMAT (' SPIN WARNING T NM K* # E A S-1 ',
318  & f8.1,i7,i4,i6,f7.3,1p,2e10.1)
319  CALL flush(66)
320  END IF
321  spin1 = 1.0
322  iwarn = iwarn + 1
323  END IF
324  IF (spin20.LT.1.0.and.spin2.gt.1.0) THEN
325  IF (mod(iwarn,100).EQ.0) THEN
326  WRITE (66,14) ttot, name(j2), kstar(j2), nstep, ecc,
327  & rp0*su/(1.0-es0), spin2-1.0
328  CALL flush(66)
329  END IF
330  spin2 = 1.0
331  iwarn = iwarn + 1
332  END IF
333 *
334 * Re-scale the semi-major axis and angular velocities to N-body units.
335  semi = rad*semi0*semi
336  spin1 = meanmotion*spin1
337  spin2 = meanmotion*spin2
338  ecc = max(ecc,0.001d0)
339  ecc = min(ecc,0.999d0)
340 *
341 * Convert back to angular momenta.
342  IF (kstar(j1).LE.2.OR.(kstar(j1).GE.7.AND.kstar(j1).NE.9)) THEN
343  spin(j1) = rg2(1)*body(j1)*radius(j1)**2*spin1
344  ELSE
345  spin(j1) = (rg2(1)*body(j1)*radius(j1)**2 +
346  & 0.21*mc1/smu*rc1**2)*spin1
347  END IF
348  IF (kstar(j2).LE.2.OR.(kstar(j2).GE.7.AND.kstar(j2).NE.9)) THEN
349  spin(j2) = rg2(2)*body(j2)*radius(j2)**2*spin2
350  ELSE
351  spin(j2) = (rg2(2)*body(j2)*radius(j2)**2 +
352  & 0.21*mc2/smu*rc2**2)*spin2
353  END IF
354 *
355 * Obtain the tidal contributions from integration.
356  a0 = rp0/(1.0 - es0)
357  dh = -0.5*body(i)*(1.0/semi - 1.0/a0)
358 *
359 * Update energy and semi-major axis.
360  h(ipair) = h(ipair) + dh
361  semi = -0.5*body(i)/h(ipair)
362 *
363 * Set new pericentre from final elements and update reference values.
364  IF (list(1,i1).EQ.0) THEN
365  qperi = semi*(1.0 - ecc)
366  ELSE
367  qperi = rp0*(1.0 + es0)/(1.0 + ecc)
368  END IF
369  rp(ic) = qperi
370  es(ic) = ecc
371  tosc(ic) = time
372 *
373 * Define synchronous state at the end and adopt ECC = ECCM.
374  IF (ecc.LT.eccm) THEN
375  ione = 0
376  kstar(i) = 10
377  ncirc = ncirc + 1
378  tb = days*semi*sqrt(semi/body(i))
379 * Ensure the dominant star has synchronous rotation.
380  IF(kstar(j1).LE.2.OR.(kstar(j1).GE.7.AND.kstar(j1).NE.9))THEN
381  spin(j1) = rg2(1)*body(j1)*radius(j1)**2*meanmotion
382  ELSE
383  spin(j1) = (rg2(1)*body(j1)*radius(j1)**2 +
384  & 0.21*mc1/smu*rc1**2)*meanmotion
385  END IF
386  dom = (spin1 - meanmotion)/meanmotion
387  WRITE (6,15) ttot, ipair, es0, ecc, rp0, r(ipair), h(ipair),
388  & tb, xn, dom
389  15 FORMAT (' END SPIRAL T KS E0 E RP0 R H P n DOM/OM ',
390  & f9.2,i4,2f8.4,1p,4e10.2,0p,f5.1,f7.3)
391  ecc = eccm
392  END IF
393 *
394 * Ensure apocentre or arbitrary phase is replaced by pericentre.
395  err = abs(rp0 - r(ipair))/rp0
396  IF (r(ipair).GT.semi.OR.err.GT.1.0e-05) THEN
397 * Reduce eccentric anomaly by pi for inward motion.
398  IF (tdot2(ipair).LT.0.0d0) THEN
399  CALL ksapo(ipair)
400  END IF
401 * Transform from outward motion (anomaly < pi) to exact pericentre.
402  IF (gamma(ipair).LT.1.0d-04) THEN
403 * Note consistency problem for large GAMMA (HDOT effect not included).
404  tt0 = time
405  CALL ksperi(ipair)
406 * Restore TIME just in case (possible bug in scheduling of #I1).
407  time = tt0
408  END IF
409  END IF
410 *
411 * Form KS coordinate scaling factor from pericentre ratio.
412  c1 = sqrt(qperi/rp0)
413 *
414 * Set KS velocity scaling from energy relation with RP0 instead of R.
415  c2 = sqrt((body(i) + h(ipair)*qperi)/(body(i) + hi*rp0))
416 *
417 * Scale KS variables to yield the prescribed elements (set TDOT2 >= 0).
418  r(ipair) = 0.0d0
419  tdot2(ipair) = 0.0
420  DO 18 k = 1,4
421  u(k,ipair) = c1*u(k,ipair)
422  udot(k,ipair) = c2*udot(k,ipair)
423  u0(k,ipair) = u(k,ipair)
424  r(ipair) = r(ipair) + u(k,ipair)**2
425  tdot2(ipair) = tdot2(ipair) + 2.0*u(k,ipair)*udot(k,ipair)
426  18 CONTINUE
427  tdot2(ipair) = max(tdot2(ipair),0.0d0)
428 *
429 * Perform energy correction to maintain conservation.
430  zmu = body(i1)*body(i2)/body(i)
431  ecoll = ecoll + zmu*(hi - h(ipair))
432  egrav = egrav + zmu*(hi - h(ipair))
433 *
434 * Rectify large eccentricity deviation from integrated value.
435  ecc2 = (1.0 - r(ipair)/semi)**2 + tdot2(ipair)**2/(body(i)*semi)
436  eccf = sqrt(ecc2)
437  IF (abs(ecc - eccf).GT.0.01) THEN
438  CALL ksrect(ipair)
439  END IF
440 *
441 * Check Roche time and update RADIUS of primary frequently.
442  hi = h(ipair)
443 * Specify semi-major axis and binding energy for evaluating Roche time.
444  ac = qperi*(1.0 + ecc)
445  h(ipair) = -0.5*body(i)/ac
446  CALL trflow(ipair,dtr)
447  h(ipair) = hi
448 *
449 * Obtain Roche radius and stellar radius for each component.
450  j = i1
451  DO 20 k = 1,2
452  q1 = body(j)/(body(i) - body(j))
453  rol(k) = rl(q1)
454  rj(k) = radius(j)
455  j = i2
456  20 CONTINUE
457 *
458 * Determine indices for primary & secondary star (donor & accretor).
459  IF (rj(1)/rol(1).GE.rj(2)/rol(2)) THEN
460  j1 = i1
461  j2 = i2
462  ELSE
463  j1 = i2
464  j2 = i1
465  END IF
466 *
467 * Form critical mass ratio using approximate core mass.
468  IF(kstar(j1).EQ.2)THEN
469  qc = 4.d0
470  ELSE IF(kstar(j1).GE.3.AND.kstar(j1).LE.6)THEN
471  zpars7 = 0.3
472  zmc = cm(1,ic)
473  IF (j1.EQ.i2) zmc = cm(2,ic)
474  qc = (1.67d0-zpars7+2.d0*(zmc/body(j1))**5)/2.13d0
475 *
476 * Consider using condition of Hjellming & Webbink, 1987, ApJ, 318, 794.
477 * QC = 0.362D0 + 1.D0/(3.D0*(1.D0 - MASSC(1)/MASS(1)))
478  ELSE IF(kstar(j1).EQ.8.OR.kstar(j1).EQ.9)THEN
479  qc = 0.784d0
480  ELSE
481  qc = 1000.0
482  ENDIF
483 *
484 * Adopt common envelope evolution for eccentric binaries with Q1 > QC.
485  q1 = body(j1)/body(j2)
486  IF (n.GT.0) go to 99 ! Fudge 10/12 to avoid COAL without ISTAT(KCASE).
487  IF (dtr.LT.0.1/tstar.AND.q1.GT.qc.AND.istar.EQ.0) THEN
488  kspair = ipair
489  time = tblock
490  iqcoll = 1
491  CALL expel(j1,j2,icase)
492 * Exit before updating TEV in the unlikely case of coalescence.
493  IF (iphase.EQ.-1) go to 100
494  tev(i1) = tev(i) + 2.0*stepx
495  tev(i2) = tev(i1)
496  go to 100
497  END IF
498  99 CONTINUE ! Fudge leads to SPIRAL TERM and ENFORCED ROCHE.
499 *
500 * Enforce termination at constant angular momentum on Roche condition.
501  IF (dtr.LE.0.1/tstar.OR.
502  & (step(i1).GT.100.0*tk.AND.kstar(j1).GE.5)) THEN
503 * Include also enforcement for AGB stars in unperturbed state.
504  af = semi*(1.0 - ecc**2)/(1.0 - eccm**2)
505  h(ipair) = -0.5*body(i)/af
506  ecoll = ecoll + zmu*(hi - h(ipair))
507  egrav = egrav + zmu*(hi - h(ipair))
508  kstar(i) = 10
509  WRITE (6,21) es0, ecc, af/semi, dtr, semi, af
510  21 FORMAT (' WARNING! SPIRAL TERM E0 E A/A0 DTR A AF ',
511  & 2f8.4,f6.2,1p,3e10.2)
512 * Ensure the dominant star has synchronous rotation.
513  meanmotion = sqrt(body(i)/af**3)
514  IF(kstar(j1).LE.2.OR.(kstar(j1).GE.7.AND.kstar(j1).NE.9))THEN
515  spin(j1) = rg2(1)*body(j1)*radius(j1)**2*meanmotion
516  ELSE
517  spin(j1) = (rg2(1)*body(j1)*radius(j1)**2 +
518  & 0.21*mc1/smu*rc1**2)*meanmotion
519  END IF
520  dom = (spin1 - meanmotion)/meanmotion
521  WRITE (6,26) name(j1), q1, qc, dom
522  26 FORMAT (' ENFORCED ROCHE NM(J1) Q1 QC DOM/OM',i8,2f8.2,f7.3)
523 * Modify KS variables to yield circular velocity with energy H.
524  CALL expand(ipair,r(ipair))
525  CALL resolv(ipair,1)
526  END IF
527 *
528 * Initialize KS for perturbed motion (cf. backwards step in KSPERI).
529  IF (list(1,i1).GT.0) THEN
530 * Ensure unperturbed motion for small GAMMA.
531  IF (gamma(ipair).LT.1.0d-10.AND.es0.LT.0.95) THEN
532  list(1,i1) = 0
533  END IF
534  IF (gamma(ipair).LT.5.0d-07.AND.es0.LT.0.3) THEN
535  list(1,i1) = 0
536  END IF
537  CALL resolv(ipair,1)
538  imod = kslow(ipair)
539  CALL kspoly(ipair,imod)
540  END IF
541 *
542 * Rectify the orbit to yield consistent variables (only at end).
543  IF (kstar(i).EQ.10) THEN
544  CALL ksrect(ipair)
545 *
546 * Re-evaluate eccentricity after rectification.
547  semi = -0.5*body(i)/h(ipair)
548  ecc2 = (1.0-r(ipair)/semi)**2 + tdot2(ipair)**2/(body(i)*semi)
549  ecc = sqrt(ecc2)
550 *
551 * Deform the orbit to small eccentricity (H = const).
552  IF (ecc.GT.eccm) THEN
553  CALL deform(ipair,ecc,eccm)
554 *
555  WRITE (6,22) ecc, eccm, time-time0, semi,
556  & body(i1)*zmbar, body(i2)*zmbar
557  22 FORMAT (' DEFORM SPIRAL E EF T-TOSC SEMI M1 M2 ',
558  & 2f8.4,f9.4,1p,e10.2,0p,2f6.2)
559  END IF
560 *
561 * Determine indices for primary & secondary star (donor & accretor).
562  IF (body(i1)/radius(i1)**3.LT.body(i2)/radius(i2)**3) THEN
563  j1 = i1
564  j2 = i2
565  ELSE
566  j2 = i1
567  j1 = i2
568  END IF
569 *
570 * Define mass ratio and evaluate Roche radius for the primary.
571  q0 = body(j1)/body(j2)
572  q1 = q0**0.3333
573  q2 = q1**2
574  rl1 = 0.49*q2/(0.6*q2 + log(1.0d0 + q1))*semi
575 *
576 * Update Roche look-up time (does not depend on choice of primary).
577  CALL trflow(ipair,dtr)
578  tev(i) = time + dtr
579  tk = days*semi*sqrt(semi/body(i))
580 *
581  WRITE (6,24) ic, name(j1), name(j2), kstar(j1), kstar(j2),
582  & semi*su, rl1*su, radius(j1)*su, radius(j2)*su,
583  & tphys, tk, dtr
584  24 FORMAT (' ROCHE CHECK IC NAM K* A RL R* TP TK DTR ',
585  & i4,2i6,2i4,4f7.1,f9.2,1p,2e9.1)
586  zm1 = body(j1)*smu
587  CALL trdot(j1,dtm,zm1)
588  tev(j1) = time + dtm
589 * IF (DTR.LT.0.1/TSTAR) THEN
590 * TEV(I1) = TEV(I) + 2.0*STEPX
591 * TEV(I2) = TEV(I1)
592 * END IF
593 * Include enforcement of Roche coalescence to prevent shrinkage stop.
594  IF (dtr.EQ.0.0d0.AND.semi.LT.0.5*radius(j1)) THEN
595  kspair = ipair
596  iqcoll = 1
597  WRITE (6,23) ipair
598  23 FORMAT (' ENFORCED COAL KS = ',i4)
599 * Delay coalescence until detected by ROCHE or UNPERT.
600 * CALL CMBODY(R(IPAIR),2)
601 * IF (IPHASE.LT.0) GO TO 100
602  END IF
603  END IF
604 *
605 * Include occasional diagnostics of spiral evolution (every output).
606  IF (abs(time - tadj).LT.step(i1)) THEN
607  IF (rp1.LT.3.0.AND.(ecc.LT.0.1.OR.kstar(j1).GE.3)) THEN
608  np = list(1,i1)
609  tk = days*semi*sqrt(semi/body(i))
610  WRITE (6,25) name(i1), name(i2), ipair, nchaos,
611  & kstar(i1), kstar(i2), np, ecc, rp1, tk
612  25 FORMAT (' SPIRAL NAM KS NCH K* NP E RP P ',
613  & 2i6,i5,4i4,f8.4,f6.2,f8.2)
614  END IF
615  END IF
616 *
617 * Include enforced circularization for difficult cases.
618  IF ((ecc.LT.0.6.AND.gamma(ipair).LT.2.0d-07).OR.
619  & (ecc.LT.0.8.AND.gamma(ipair).LT.1.0d-08).OR.
620  & (ecc.LT.0.95.AND.gamma(ipair).LT.1.0d-09)) THEN
621  icirc = 0
622  CALL tcirc(qperi,ecc,i1,i2,icirc,tc)
623  IF (tc.GT.200.0.OR.(tc.GT.10.AND.ecc.LT.0.7)) THEN
624  CALL ksperi(ipair)
625  CALL ksapo(ipair)
626  CALL deform(ipair,ecc,eccm)
627  kstar(i) = 10
628  ncirc = ncirc + 1
629  CALL resolv(ipair,1)
630  CALL kspoly(ipair,1)
631  np = list(1,i1)
632  WRITE (6,28) ipair, np, name(i1), ecc, tc, gamma(ipair)
633  28 FORMAT (' ENFORCED CIRC KS NP NM E TC G ',
634  & 2i4,i6,f7.3,1p,2e9.1)
635  END IF
636  END IF
637 *
638 * Set standard binary on large ECC and TC if small GAMMA (large omeq).
639  IF (ecc.GT.0.95.AND.gamma(ipair).LT.1.0d-09) THEN
640  icirc = 0
641  CALL tcirc(qperi,ecc,i1,i2,icirc,tc)
642  IF (tc.GT.1.0d+04) THEN
643  kstar(i) = 0
644  irem = ic
645  WRITE (6,29) name(i1), list(1,i1), ecc, tc
646  29 FORMAT (' FROZEN CIRC NM NP E TC ',i7,i4,f9.5,1p,e9.1)
647  go to 30
648  END IF
649  END IF
650 *
651 * Check for terminated or escaped chaotic binaries at end of spiral.
652  30 IF (kstar(i).EQ.10.OR.irem.GT.0) THEN
653  IF (irem.GT.0) go to 50
654  j = 1
655 * See whether case #J indicates current or escaped/disrupted KS binary.
656  40 DO 45 jpair = 1,npairs
657  IF (namec(j).EQ.name(n+jpair)) THEN
658 * Update #J if KSTAR > 0, otherwise consider next member.
659  IF (kstar(n+jpair).GT.0) THEN
660  go to 50
661  ELSE
662  go to 70
663  END IF
664  END IF
665  45 CONTINUE
666 *
667 * Skip during multiple regularizations (KSTAR not visible).
668  IF (nsub.GT.0) THEN
669  go to 70
670  END IF
671 *
672 * Skip removal if chaos binary is member of single/double merger.
673  IF (nmerge.GT.0) THEN
674  DO 48 jpair = 1,npairs
675  IF (name(n+jpair).LT.0) THEN
676  IF (namec(j).EQ.nzero + name(2*jpair-1).OR.
677  & namec(j).EQ.nzero + name(2*jpair).OR.
678  & namec(j).LT.-2*nzero.OR.
679  & namec(j).EQ.name(2*jpair)) THEN
680  WRITE (71,46) nchaos, namec(j), name(n+jpair),
681  & name(2*jpair)
682  46 FORMAT (' MERGED SPIRAL NCH NAM ',i4,3i7)
683  CALL flush(71)
684  go to 70
685  END IF
686  END IF
687  48 CONTINUE
688  END IF
689 *
690 * Update chaos variables for #IC and any disrupted or escaped binaries.
691  50 nchaos = nchaos - 1
692 * Copy chaos index in case of KS termination.
693  IF (irem.GT.0) THEN
694  j = irem
695  END IF
696 *
697  DO 60 l = j,nchaos
698  l1 = l + 1
699  DO 55 k = 1,4
700  eosc(k,l) = eosc(k,l1)
701  55 CONTINUE
702  eb0(l) = eb0(l1)
703  zj0(l) = zj0(l1)
704  ecrit(l) = ecrit(l1)
705  ar(l) = ar(l1)
706  br(l) = br(l1)
707  edec(l) = edec(l1)
708  tosc(l) = tosc(l1)
709  rp(l) = rp(l1)
710  es(l) = es(l1)
711  cm(1,l) = cm(1,l1)
712  cm(2,l) = cm(2,l1)
713  iosc(l) = iosc(l1)
714  namec(l) = namec(l1)
715  60 CONTINUE
716 * Ensure next chaos location contains zero core masses.
717  cm(1,nchaos+1) = 0.0
718  cm(2,nchaos+1) = 0.0
719 * Consider the same location again after each removal (J <= NCHAOS).
720  j = j - 1
721  70 j = j + 1
722  IF (j.LE.nchaos.AND.irem.EQ.0) go to 40
723  END IF
724 *
725 * Check optional diagnostics on transition to ECC = 0 or SLEEP.
726  IF (kz(8).GT.3.AND.kstar(i).NE.-2) THEN
727 * Note case IPAIR < 0 from CHRECT with I denoting c.m.
728  IF (ipair.LT.0) ipair = i - n
729  CALL binev(ipair)
730  END IF
731 *
732 
733  100 RETURN
734 *
735  END