Nbody6
 All Files Functions Variables
expel.f
Go to the documentation of this file.
1  SUBROUTINE expel(J1,J2,ICASE)
2 *
3 *
4 * Common envelope stage of interacting stars.
5 * -------------------------------------------
6 *
7  include 'common6.h'
8  common/modes/ eb0(ntmax),zj0(ntmax),ecrit(ntmax),ar(ntmax),
9  & br(ntmax),eosc(4,ntmax),edec(ntmax),tosc(ntmax),
10  & rp(ntmax),es(ntmax),cm(2,ntmax),iosc(ntmax),
11  & namec(ntmax)
12  REAL*8 m01,m1,r1,aj1,m02,m2,r2,aj2,sep,mi(2)
13  REAL*8 tscls(20),lums(10),gb(10),tm,tn,lum1,lum2,mc1,mc2,rcc
14  REAL*8 jspin1,jspin2,menv,renv,k2
15  CHARACTER*8 which1
16  LOGICAL coals
17 *
18 *
19 * Define global indices such that body #I1 is giant-like.
20  IF(kstar(j1).GE.2.AND.kstar(j1).LE.9.AND.kstar(j1).NE.7)THEN
21  i1 = j1
22  i2 = j2
23  ELSE
24  i1 = j2
25  i2 = j1
26  ENDIF
27 *
28 * Update the stars to latest previous time.
29  tev1 = max(tev0(i1),tev0(i2))
30 *
31 * Specify basic parameters for both stars.
32  m01 = body0(i1)*zmbar
33  m1 = body(i1)*zmbar
34  mc1 = 0.d0
35  aj1 = tev1*tstar - epoch(i1)
36  jspin1 = spin(i1)*spnfac
37  kw1 = kstar(i1)
38  m02 = body0(i2)*zmbar
39  m2 = body(i2)*zmbar
40  mc2 = 0.d0
41  aj2 = tev1*tstar - epoch(i2)
42  jspin2 = spin(i2)*spnfac
43  kw2 = kstar(i2)
44  iter = 0
45 *
46 * Save current semi-major axis.
47  ipair = kspair
48  i = n + ipair
49  semi0 = -0.5d0*body(i)/h(ipair)
50  sep = semi0*su
51  ecc2 = (1.d0-r(ipair)/semi0)**2+tdot2(ipair)**2/(semi0*body(i))
52  ecc = sqrt(ecc2)
53  ecc0 = ecc
54  dm2 = 0.0
55  ep2 = epoch(i2)
56  age2 = aj2
57 *
58 * Perform common envelope evolution (note: SEP in SU).
59  1 CALL comenv(m01,m1,mc1,aj1,jspin1,kw1,
60  & m02,m2,mc2,aj2,jspin2,kw2,ecc,sep,coals)
61 *
62 * Obtain consistent radii for the stars (skip #I2 on coalescence).
63  CALL star(kw1,m01,m1,tm,tn,tscls,lums,gb,zpars)
64  CALL hrdiag(m01,aj1,m1,tm,tn,tscls,lums,gb,zpars,
65  & r1,lum1,kw1,mc1,rcc,menv,renv,k2)
66  IF(kw1.NE.kstar(i1))THEN
67  WRITE(38,*)' EXPEL TYPE CHANGE1 ',kstar(i1),kw1
68  ENDIF
69  IF(coals)THEN
70  kw2 = kstar(i2)
71  r2 = 0.d0
72  lum2 = 0.d0
73  jspin2 = 0.d0
74  ELSE
75  CALL star(kw2,m02,m2,tm,tn,tscls,lums,gb,zpars)
76  CALL hrdiag(m02,aj2,m2,tm,tn,tscls,lums,gb,zpars,
77  & r2,lum2,kw2,mc2,rcc,menv,renv,k2)
78  IF(kw2.NE.kstar(i2))THEN
79  WRITE(38,*)' EXPEL TYPE CHANGE2 ',kstar(i2),kw2
80  ENDIF
81  END IF
82 *
83 * Check coalescence condition again (bug 7/1/09).
84  IF (sep.LT.r1.OR.sep.LT.r2) THEN
85  coals = .true.
86  END IF
87 *
88 * Skip small mass loss unless coalescence (May 2003).
89  dms = m1 + m2 - body(i)*smu
90  IF (abs(dms).LT.1.0d-10.AND..NOT.coals) THEN
91 * WRITE (77,7) TIME, ECC, DMS, SEMI0*SU-SEP
92 * 7 FORMAT (' FAIL! T ECC DMS DSEP ',F12.4,F8.4,1P,2E10.2)
93 * CALL FLUSH(77)
94  iphase = 0
95  go to 60
96  END IF
97 *
98 * Copy new values of radius, luminosity & semi-major axis.
99  radius(i1) = r1/su
100  radius(i2) = r2/su
101  zlmsty(i1) = lum1
102  zlmsty(i2) = lum2
103  spin(i1) = jspin1/spnfac
104  spin(i2) = jspin2/spnfac
105 * Accept negative semi-major axis for fly-by (bug fix 12/08).
106  semi = sep/su
107 *
108 * Set reduced mass & binding energy and update masses & time T0(J1).
109  zmu0 = body(i1)*body(i2)/body(i)
110  hi = h(ipair)
111  body0(i1) = m01/zmbar
112  body0(i2) = m02/zmbar
113  t0(j1) = time
114 * Add current energy (final state: coalescence or close binary).
115  ecoll = ecoll + zmu0*hi
116  egrav = egrav + zmu0*hi
117 *
118 * Distinguish between coalescence and surviving binary.
119  epoch(i1) = tev1*tstar - aj1
120  IF(coals)THEN
121  mi(1) = m1
122  mi(2) = m2
123  body0(i2) = 0.d0
124 * Check for TZ object formed by CE evolution.
125  IF(kstar(i2).GE.13.AND.kw1.GE.13)THEN
126  ntz = ntz + 1
127  WRITE (6,5) kw1, m1, m2
128  5 FORMAT (' NEW TZ K* M1 M2 ',i4,2f7.2)
129  ENDIF
130  CALL coal(ipair,kw1,kw2,mi)
131  ELSE
132 *
133 * Update evolution times and TMDOT.
134  epoch(i2) = tev1*tstar - aj2
135  tev(i1) = tev1
136  tev0(i1) = tev(i1)
137  tev(i2) = tev1
138  tev0(i2) = tev(i2)
139  tmdot = min(tmdot,tev1)
140 * Shrink the orbit in case of no coalescence.
141  body(i1) = m1/zmbar
142  body(i2) = m2/zmbar
143 * Form mass loss & new binary mass and re-define binding energy.
144  dm = body(i) - (body(i1) + body(i2))
145 * Save mass loss for potential energy correction after ENFORCED CE.
146  dm2 = dm2 + dm
147  zmass = zmass - dm
148  body(i) = body(i) - dm
149  h(ipair) = -0.5d0*body(i)/semi
150 *
151 * Subtract final binding energy of surviving binary.
152  zmu = body(i1)*body(i2)/body(i)
153  ecoll = ecoll - zmu*h(ipair)
154  egrav = egrav - zmu*h(ipair)
155 *
156 * Set argument for new pericentre velocity (need to have H*A*(1-E)).
157  rx = r(ipair)/(1.0 - ecc)
158 * Modify KS variables consistent with new eccentricity and energy H.
159  CALL expand(ipair,rx)
160  CALL resolv(ipair,1)
161 *
162  IF (ecc0.LT.1.0) THEN
163  which1 = ' BINARY '
164  ELSE
165  which1 = ' HYPERB '
166  nhypc = nhypc + 1
167  END IF
168  WRITE (6,10) which1, name(i1), name(i2), kstar(i1),kstar(i2),
169  & kw1, kw2, m1, m2, dm*zmbar, ecc0, ecc, r1, r2,
170  & semi0*su, semi*su
171  10 FORMAT (a8,'CE NAM K0* K* M1 M2 DM E0 E R1 R2 A0 A ',
172  & 2i6,4i3,3f5.1,2f8.4,2f7.1,1p,e9.1,0p,f7.1)
173 *
174 * Check common envelope condition again after circularization (09/08).
175  IF (ecc0.GT.0.001d0.AND.ecc.LE.0.001d0) THEN
176  q1 = m1/m2
177  q2 = 1.d0/q1
178  rl1 = rl(q1)
179  rl2 = rl(q2)
180  IF (r1.GT.rl1*sep.OR.r2.GT.rl2*sep) THEN
181  iter = iter + 1
182  IF (iter.EQ.1) THEN
183  WRITE (6,8) rl1*sep, rl2*sep
184  8 FORMAT (' ENFORCED CE RL*SEP ',1p,2e10.2)
185  ecc0 = ecc
186  semi0 = semi
187  sep = semi*su
188 * Correct for mass loss in case of coalescence which ignores DM2.
189  IF (dm*smu.LT.0.05) THEN
190  CALL ficorr(i,dm)
191  ELSE
192  CALL fcorr(i,dm,0)
193  END IF
194  dm2 = 0.0
195  go to 1
196  END IF
197  END IF
198  END IF
199 *
200 * Copy mass loss (one or two COMENV) and new types.
201  IF (iter.GT.0) dm = dm2
202  kstar(i1) = kw1
203  kstar(i2) = kw2
204 *
205 * Obtain neighbour list for force corrections and polynomials.
206  nnb = list(1,i)
207  DO 11 l = 1,nnb+1
208  ilist(l) = list(l,i)
209  11 CONTINUE
210 *
211 * Ensure rare case of massless remnant will escape at next output.
212  IF (kw1.EQ.15.OR.kw2.EQ.15) THEN
213  CALL ksterm
214  iphase = -1
215  i = ifirst
216  IF (body(i).GT.0.0d0) i = ifirst + 1
217  t0(i) = tadj + dtadj
218  step(i) = 1.0d+06
219  ri = sqrt(x(1,i)**2 + x(2,i)**2 + x(3,i)**2)
220  vi = sqrt(xdot(1,i)**2 + xdot(2,i)**2 + xdot(3,i)**2)
221  DO 12 k = 1,3
222  x0(k,i) = min(1.0d+04+x(k,i),1000.0*rscale*x(k,i)/ri)
223  x(k,i) = x0(k,i)
224  x0dot(k,i) = sqrt(0.004*zmass/rscale)*xdot(k,i)/vi
225  xdot(k,i) = x0dot(k,i)
226  f(k,i) = 0.0d0
227  fdot(k,i) = 0.0d0
228  d2(k,i) = 0.0d0
229  d3(k,i) = 0.0d0
230  12 CONTINUE
231  WRITE (6,13) name(i), kstar(i), body(i)*zmbar
232  13 FORMAT (' MASSLESS GHOST NAM K* M ',i6,i4,1p,e10.2)
233 * Include special treatment for velocity kick of KS binary.
234  ELSE IF (kw1.GE.10) THEN
235 *
236  IF (ecc.LE.0.001d0) kstar(n+ipair) = 10
237 * Re-initialize the KS solution with new perturber list.
238  CALL kslist(ipair)
239  CALL kspoly(ipair,1)
240 * Evaluate binary disruption velocity and introduce velocity kick.
241  vi2 = xdot(1,i)**2 + xdot(2,i)**2 + xdot(3,i)**2
242  kstar(i1) = -kstar(i1)
243  CALL kick(ipair,0,kw1,dm)
244 *
245 * Include potential and kinetic energy corrections due to mass loss.
246  CALL poti(i,potj)
247  ecdot = ecdot + dm*potj + 0.5d0*dm*vi2
248 *
249 * See whether tidal terms should be included (standard or scaled).
250  IF (kz(14).GT.0) THEN
251  IF (kz(14).LE.2) THEN
252  ecdot = ecdot - 0.5d0*dm*(tidal(1)*x(1,i)**2 +
253  & tidal(3)*x(3,i)**2)
254  ELSE
255  body(i) = body(i) + dm
256  CALL xtrnlv(i,i)
257  ecdot = ecdot + ht
258  body(i) = body(i) - dm
259  CALL xtrnlv(i,i)
260  ecdot = ecdot - ht
261  END IF
262  END IF
263 *
264 * Form new c.m. variables.
265  t0(i) = time
266  DO 20 k = 1,3
267  x(k,i) = (body(i1)*x(k,i1) + body(i2)*x(k,i2))/body(i)
268  xdot(k,i) = (body(i1)*xdot(k,i1) +
269  & body(i2)*xdot(k,i2))/body(i)
270  x0(k,i) = x(k,i)
271  x0dot(k,i) = xdot(k,i)
272  20 CONTINUE
273 *
274 * Obtain new F & FDOT and time-steps for neighbours.
275  DO 30 l = 2,nnb+1
276  j = ilist(l)
277  DO 25 k = 1,3
278  x0dot(k,j) = xdot(k,j)
279  25 CONTINUE
280  CALL fpoly1(j,j,0)
281  CALL fpoly2(j,j,0)
282  30 CONTINUE
283  tprev = time - stepx
284 *
285 * Terminate KS binary and assign kick velocity to single star #I.
286  i = i1 + 2*(npairs - ipair)
287  CALL ksterm
288  CALL kick(i,1,kw1,dm)
289 * Initialize new KS polynomials after velocity kick (H > 0 is OK).
290  icomp = ifirst
291  jcomp = ifirst + 1
292  CALL ksreg
293  iphase = -1
294 *
295 * Provide diagnostic output if binary survives the kick.
296  i = ntot
297  kstar(i) = 0
298  IF (h(npairs).LT.0.0) THEN
299  semi = -0.5d0*body(i)/h(npairs)
300  WRITE (6,35) kw1, r(npairs)/semi, semi*su,
301  & body(i)*zmbar, (xdot(k,i)*vstar,k=1,3)
302  35 FORMAT (' WD/NS BINARY KW R/A A M V ',
303  & i4,3f7.2,3f7.1)
304  END IF
305  ELSE
306 * Set new binary indicator and Roche look-up time (ECC > 0 is OK).
307  kstar(i) = 0
308  CALL trflow(ipair,dtr)
309  tev(i) = time + dtr
310  tmdot = min(tev(i),tmdot)
311 *
312 * Update KSTAR & chaos variables for SYNCH after circularization.
313  IF (ecc.LE.0.001d0) THEN
314  kstar(i) = 10
315  nce = nce + 1
316  ic = 0
317 * Note: NAMEC may be set in routine CHAOS without call to SPIRAL.
318  DO 36 l = 1,nchaos
319  IF (name(i).EQ.namec(l)) ic = l
320  36 CONTINUE
321 * Include safety test on no membership just in case (not sure!).
322  IF (ic.EQ.0) THEN
323  nchaos = nchaos + 1
324  ic = nchaos
325  namec(ic) = name(i)
326  IF (nchaos.GT.ntmax) THEN
327  WRITE (6,38) name(i1), nchaos, ecc
328  38 FORMAT (' FATAL ERROR! EXPEL NM NCH E ',
329  & i6,i4,f8.4)
330  stop
331  END IF
332  END IF
333  IF (ic.GT.0) THEN
334  rp(ic) = semi*(1.0 - ecc)
335  es(ic) = ecc
336  tosc(ic) = time
337  END IF
338  END IF
339 *
340 * Include prediction to end of current block-time before new KS.
341  time = tblock
342  IF (t0(i).LT.time.AND.dm.GT.0.0d0) THEN
343  CALL xvpred(i,-2)
344  t0(i) = time
345  CALL dtchck(time,step(i),dtk(40))
346  DO 40 k = 1,3
347  x0(k,i) = x(k,i)
348  x0dot(k,i) = xdot(k,i)
349  40 CONTINUE
350  END IF
351 *
352 * Set IPHASE < 0 and re-initialize KS solution with new perturber list.
353  iphase = -1
354  CALL kslist(ipair)
355  CALL kspoly(ipair,1)
356 * Avoid negative radial velocity at pericentre (bug fix 12/08).
357  IF (tdot2(ipair).LT.0.0d0) tdot2(ipair) = 0.0
358 *
359 * Correct neighbour forces and update system energy loss due to DM.
360  IF (dm.GT.0.0) THEN
361  IF (dm*smu.LT.0.05) THEN
362  CALL ficorr(i,dm)
363  ELSE
364  CALL fcorr(i,dm,0)
365  END IF
366  END IF
367 *
368  IF (dm*smu.GT.0.1) THEN
369 * Include body #I at the end (counting from location #2).
370  nnb2 = nnb + 2
371  ilist(nnb2) = i
372 *
373 * Obtain new F & FDOT and time-steps.
374  DO 50 l = 2,nnb2
375  j = ilist(l)
376  IF (l.EQ.nnb2) THEN
377  j = i
378  ELSE
379  CALL xvpred(j,-2)
380  CALL dtchck(time,step(j),dtk(40))
381  DO 45 k = 1,3
382  x0dot(k,j) = xdot(k,j)
383  45 CONTINUE
384  END IF
385  CALL fpoly1(j,j,0)
386  CALL fpoly2(j,j,0)
387  50 CONTINUE
388  END IF
389  tprev = time - stepx
390  END IF
391  END IF
392 *
393 * Include cleaning of CHAOS arrays (cf. KSTERM for current COAL).
394  l = 0
395  ic = 0
396  52 l = l + 1
397  kk = 0
398  DO 55 i = n+1,ntot
399  IF (name(i).EQ.namec(l)) kk = 1
400  55 CONTINUE
401  nch1 = nchaos
402 * Perform removal of NAMEC on failed search.
403  IF (kk.EQ.0.AND.kstar(n+ipair).LT.0) THEN
404  ii = -l
405  CALL spiral(ii)
406  END IF
407  ic = ic + 1
408 * Consider the same location again after array update.
409  IF (nchaos.LT.nch1) l = l - 1
410  IF (ic.LT.nchaos+1) go to 52
411 *
412 * Reverse case indicator for exit from collision routine.
413  60 icase = -icase
414 *
415  RETURN
416 *
417  END