Nbody6
 All Files Functions Variables
expel2.f
Go to the documentation of this file.
1  SUBROUTINE expel2(J1,J2,ICASE)
2 *
3 *
4 * Common envelope stage for chain system.
5 * ---------------------------------------
6 *
7  include 'common6.h'
8  parameter(nmx=10,nmx2=2*nmx,nmx3=3*nmx,nmx4=4*nmx,
9  & nmx8=8*nmx,nmxm=nmx*(nmx-1)/2)
10  REAL*8 m,mass,mc,mij,mkk,xcm(3),xrel(3),vcm(3),vrel(3)
11  common/chain1/ xch(nmx3),vch(nmx3),m(nmx),
12  & zz(nmx3),wc(nmx3),mc(nmx),
13  & xi(nmx3),pi(nmx3),mass,rinv(nmxm),rsum,mkk(nmx),
14  & mij(nmx,nmx),tkk(nmx),tk1(nmx),iname(nmx),nn
15  common/chainc/ xc(3,ncmax),uc(3,ncmax),bodyc(ncmax),ich,
16  & listc(lmax)
17  common/chreg/ timec,tmax,rmaxc,cm(10),namec(6),nstep1,kz27,kz30
18  common/ccoll2/ qk(nmx4),pk(nmx4),rik(nmx,nmx),SIZE(nmx),vstar1,
19  & ecoll1,rcoll,qperi,istar(nmx),icoll,isync,ndiss1
20  common/ebsave/ ebs
21  common/ksave/ k10,k20
22  REAL*8 lums(10),tscls(20),gb(10),m01,m1,tm,tn,lum1,lum2,aj1,r1,
23  & m02,m2,aj2,r2,sep,mi(2),mc1,mc2,rcc
24  REAL*8 jspin1,jspin2,menv,renv,k2
25  LOGICAL coals
26 *
27 *
28 * Define global indices such that body #I1 is giant-like.
29  IF(kstar(j1).GE.2.AND.kstar(j1).LE.9.AND.kstar(j1).NE.7)THEN
30  i1 = j1
31  i2 = j2
32  ELSE
33  i1 = j2
34  i2 = j1
35  ENDIF
36 *
37 * Save original total mass and reduced mass.
38  zmb0 = body(i1) + body(i2)
39  zmu0 = body(i1)*body(i2)/zmb0
40 *
41 * Accumulate c.m. variables for dominant bodies.
42  rij2 = 0.d0
43  vij2 = 0.d0
44  vcm2 = 0.d0
45  rdot = 0.d0
46  DO 5 k = 1,3
47  xrel(k) = x(k,i1) - x(k,i2)
48  vrel(k) = xdot(k,i1) - xdot(k,i2)
49  rij2 = rij2 + xrel(k)**2
50  vij2 = vij2 + vrel(k)**2
51  rdot = rdot + xrel(k)*vrel(k)
52  xcm(k) = (body(i1)*x(k,i1) + body(i2)*x(k,i2))/zmb0
53  vcm(k) = (body(i1)*xdot(k,i1) + body(i2)*xdot(k,i2))/zmb0
54  vcm2 = vcm2 + vcm(k)**2
55  5 CONTINUE
56 *
57 * Form binding energy per unit mass and original semi-major axis.
58  rij0 = sqrt(rij2)
59  dminc = rij0
60  hi = 0.5d0*vij2 - zmb0/rij0
61  semi0 = 2.d0/rij0 - vij2/zmb0
62  semi0 = 1.d0/semi0
63  semi00 = semi0
64 *
65  ecc2 = (1.d0 - rij0/semi0)**2 + rdot**2/(semi0*zmb0)
66  ecc = sqrt(ecc2)
67 *
68 * Obtain SEMI0 & HI from EBS because of some velocity bug (19/8/96).
69  semi0 = -0.5d0*m(k10)*m(k20)/ebs
70  hi = -0.5d0*(m(k10) + m(k20))/semi0
71 *
72 * Define non-dominant members as perturbers (JLIST(4) = 0 if NCH = 3).
73  i3 = jlist(3)
74  i4 = jlist(4)
75  jpert(1) = i3
76  np = 1
77  IF (nch.EQ.4) THEN
78  jpert(2) = i4
79  np = 2
80  END IF
81 * Evaluate potential energy of #I3/I4 before mass loss & displacement.
82  CALL nbpot(2,np,pot1)
83 *
84 * Specify basic parameters for both stars.
85  tev1 = max(tev0(i1),tev0(i2))
86  m01 = body0(i1)*zmbar
87  m1 = body(i1)*zmbar
88  mc1 = 0.d0
89  aj1 = tev1*tstar - epoch(i1)
90  jspin1 = spin(i1)*spnfac
91  kw1 = kstar(i1)
92  m02 = body0(i2)*zmbar
93  m2 = body(i2)*zmbar
94  mc2 = 0.d0
95  aj2 = tev1*tstar - epoch(i2)
96  jspin2 = spin(i2)*spnfac
97  kw2 = kstar(i2)
98  sep = semi0*su
99  ecc0 = ecc
100  iter = 0
101 *
102 * Perform common envelope evolution (note: SEP in SU).
103  1 CALL comenv(m01,m1,mc1,aj1,jspin1,kw1,
104  & m02,m2,mc2,aj2,jspin2,kw2,ecc,sep,coals)
105 *
106 * Obtain consistent radii for the stars (skip #I2 on coalescence).
107  CALL star(kw1,m01,m1,tm,tn,tscls,lums,gb,zpars)
108  CALL hrdiag(m01,aj1,m1,tm,tn,tscls,lums,gb,zpars,
109  & r1,lum1,kw1,mc1,rcc,menv,renv,k2)
110  if(kw1.ne.kstar(i1))then
111  write(38,*)' EXPEL2 TYPE CHANGE1 ',kstar(i1),kw1
112  write(38,*)' EXPEL2 TYPE CHANGE1 ',i1,name(i1),time
113  endif
114  IF(coals)THEN
115  kw2 = kstar(i2)
116  r2 = 0.d0
117  lum2 = 0.d0
118  jspin2 = 0.d0
119  ELSE
120  CALL star(kw2,m02,m2,tm,tn,tscls,lums,gb,zpars)
121  CALL hrdiag(m02,aj2,m2,tm,tn,tscls,lums,gb,zpars,
122  & r2,lum2,kw2,mc2,rcc,menv,renv,k2)
123  if(kw2.ne.kstar(i2))then
124  write(38,*)' EXPEL2 TYPE CHANGE2 ',kstar(i2),kw2
125  endif
126  ENDIF
127 *
128  dmsun = m1 + m2 - zmb0*smu
129  IF (abs(dmsun).LT.1.0d-12.AND..NOT.coals) THEN
130 * WRITE (78,7) TIME, ECC, DMSUN, SEMI0*SU-SEP
131 * 7 FORMAT (' FAIL - CHAIN T E DMS DSEP ',F10.4,F8.4,1P,2E10.2)
132 * CALL FLUSH(78)
133  iphase = 0
134  go to 100
135  END IF
136 *
137 * Check common envelope condition again after circularization (10/08).
138  IF (.NOT.coals) THEN
139  IF (ecc0.GT.0.001d0.AND.ecc.LE.0.001d0) THEN
140  q1 = m1/m2
141  q2 = 1.d0/q1
142  rl1 = rl(q1)
143  rl2 = rl(q2)
144  IF (r1.GT.rl1*sep.OR.r2.GT.rl2*sep) THEN
145  iter = iter + 1
146  IF (iter.EQ.1) THEN
147  WRITE (6,8) rl1*sep, rl2*sep
148  8 FORMAT (' ENFORCED CHAIN CE RL*SEP ',1p,2e10.2)
149  ecc0 = ecc
150  go to 1
151  END IF
152  END IF
153  END IF
154  END IF
155 *
156 * Copy new values of radius and luminosity.
157  radius(i1) = r1/su
158  radius(i2) = r2/su
159  zlmsty(i1) = lum1
160  zlmsty(i2) = lum2
161  spin(i1) = jspin1/spnfac
162  spin(i2) = jspin2/spnfac
163 *
164 * Update initial & chain masses (M2 = 0 for coalescence).
165  body0(i1) = m01/zmbar
166  body0(i2) = m02/zmbar
167  m(k10) = m1/zmbar
168  m(k20) = m2/zmbar
169  ecoll = ecoll + zmu0*hi
170  chcoll = chcoll + zmu0*hi
171  egrav = egrav + zmu0*hi
172 *
173 * Distinguish between coalescence and surviving binary.
174  epoch(i1) = tev1*tstar - aj1
175  IF(coals)THEN
176  mi(1) = m1
177  mi(2) = m2
178  body0(i2) = 0.d0
179 * Check for TZ object formed by CE evolution.
180  IF(kstar(i2).GE.13.AND.kw1.GE.13)THEN
181  ntz = ntz + 1
182  WRITE (6,10) m1, m2
183  10 FORMAT (' NEW TZ M1 M2 ',2f7.2)
184  ENDIF
185  ipair = -1
186  CALL coal(ipair,kw1,kw2,mi)
187 * Reverse case indicator to denote exit from collision routine.
188  icase = -icase
189  ELSE
190 *
191 * Update evolution times and TMDOT.
192  epoch(i2) = tev1*tstar - aj2
193  tev(i1) = tev1
194  tev0(i1) = tev(i1)
195  tev(i2) = tev1
196  tev0(i2) = tev(i2)
197  tmdot = min(tmdot,tev1)
198 *
199 * Copy new semi-major axis & masses and specify mass loss.
200  semi = sep/su
201  body(i1) = m1/zmbar
202  body(i2) = m2/zmbar
203  zmb = body(i1) + body(i2)
204  zmu = body(i1)*body(i2)/zmb
205  dm = zmb0 - zmb
206  zmass = zmass - dm
207 *
208 * Copy chain perturbers and obtain potential energy w.r.t. I1 & I2.
209 * NP = 0
210 * DO 15 L = 3,NCH
211 * NP = NP + 1
212 * JPERT(NP) = JLIST(L)
213 * 15 CONTINUE
214 * CALL NBPOT(2,NP,POT1)
215 *
216 * Distinguish between bound and hyperbolic orbit.
217  IF (ecc.LT.1.0) THEN
218 *
219 * Specify eccentricity factor for transformation to apocentre.
220  efac = sqrt((1.0 - ecc)/(1.0 + ecc))
221 *
222 * Set new coordinates and velocities for relative & absolute motion.
223  DO 20 k = 1,3
224  xrel(k) = xrel(k)*semi/rij0
225  xrel(k) = (1.0 + ecc)*xrel(k)
226  x(k,i1) = xcm(k) + body(i2)*xrel(k)/zmb
227  x(k,i2) = xcm(k) - body(i1)*xrel(k)/zmb
228  vrel(k) = sqrt(zmb/(semi*vij2))*vrel(k)
229  vrel(k) = efac*vrel(k)
230  xdot(k,i1) = vcm(k) + body(i2)*vrel(k)/zmb
231  xdot(k,i2) = vcm(k) - body(i1)*vrel(k)/zmb
232  20 CONTINUE
233 *
234  ELSE
235 * Form relative velocity by ratio of new and old pericentre value.
236  v2 = zmb*(2.0/rij0 - 1.0/semi)
237  v20 = vrel(1)**2 + vrel(2)**2 + vrel(3)**2
238  DO 30 k = 1,3
239  vrel(k) = sqrt(v2/v20)*vrel(k)
240  xdot(k,i1) = vcm(k) + body(i2)*vrel(k)/zmb
241  xdot(k,i2) = vcm(k) - body(i1)*vrel(k)/zmb
242  x(k,i1) = xcm(k) + body(i2)*xrel(k)/zmb
243  x(k,i2) = xcm(k) - body(i1)*xrel(k)/zmb
244  30 CONTINUE
245  END IF
246 *
247 * Evaluate new interaction energy and perform differential correction.
248 * CALL NBPOT(2,NP,POT2)
249 * ECOLL = ECOLL + (POT2 - POT1)
250 *
251 * Specify new energies (note BODY(I2) & ZMU = 0 for coalescence).
252  hf = -0.5d0*zmb/semi
253 *
254 * Update energy corrections (change in H and mass loss kinetic energy).
255  ecoll = ecoll - zmu*hf
256  chcoll = chcoll - zmu*hf
257  ecdot = ecdot + 0.5*dm*vcm2
258  egrav = egrav - zmu*hf
259 *
260 * Include potential energy terms due to all non-chain members.
261  potj = 0.d0
262  rx = 10.
263 * Note no initialization of neighbours done in the chain case.
264  DO 50 j = ifirst,ntot
265  IF (j.NE.i1.AND.j.NE.i2.AND.j.NE.i3.AND.j.NE.i4) THEN
266  rij2 = (x(1,j) - xcm(1))**2 +
267  & (x(2,j) - xcm(2))**2 +
268  & (x(3,j) - xcm(3))**2
269  potj = potj - body(j)/sqrt(rij2)
270  rx = min(rx,rij2)
271  END IF
272  50 CONTINUE
273 *
274 * Obtain non-dominant chain contributions after mass loss and new XREL.
275  CALL nbpot(2,np,pot2)
276 *
277 * See whether tidal terms should be included.
278  IF (kz(14).EQ.1) THEN
279  ecdot = ecdot - 0.5d0*dm*(tidal(1)*xcm(1)**2 +
280  & tidal(3)*xcm(3)**2)
281  END IF
282 *
283 * Add both potential energy contributions to yield final correction.
284  ecdot = ecdot + dm*potj + (pot2 - pot1)
285 *
286  WRITE (6,55) dm*potj, sqrt(rx), semi00*su, semi0*su, ebs
287  55 FORMAT (' EXPEL2 DM*POTJ RX A00 A0 EBS ',1p,6e10.2)
288  WRITE (6,60) (name(jlist(k)),k=1,4), kstar(i1), kstar(i2),
289  & kw1, kw2, m1, m2, dm*zmbar, ecc, r1, r2,
290  & semi0*su, semi*su, ebs, pot2-pot1
291  60 FORMAT (' CHAIN CE NAM K0* K* M1 M2 DM E R1 R2 A0 A EB DP',
292  & 4i6,4i3,3f5.1,f8.4,2f6.1,2f8.1,1p,2e9.1)
293 * Flag hyperbolic chain CE.
294  IF (semi0.LT.0.0) THEN
295  WRITE (6,61) time+toff
296  61 FORMAT (' HYPERB CHAIN CE T ',f7.1)
297  END IF
298 *
299  kstar(i1) = kw1
300  kstar(i2) = kw2
301  iphase = 0
302  IF (ecc.LE.0.001d0) nce = nce + 1
303 * Include possible NS or BH kick.
304  IF (kw1.GE.10.AND.kw1.LE.14) THEN
305 * Save chain members to prevent possible over-writing in HIVEL.
306  DO 62 l = 1,nch
307  jpert(l) = jlist(l)
308  62 CONTINUE
309  CALL kick(i1,1,kw1,dm)
310  DO 65 l = 1,nch
311  jlist(l) = jpert(l)
312  65 CONTINUE
313  END IF
314  END IF
315 *
316 * Save reference body in temporary variables (original chain masses).
317  DO 70 k = 1,3
318  xrel(k) = x(k,ich)
319  vrel(k) = xdot(k,ich)
320  xcm(k) = 0.0
321  vcm(k) = 0.0
322  70 CONTINUE
323 *
324 * Form new c.m. body for new FPOLY1/2 with continuation of chain.
325  zm = 0.0
326  DO 80 l = 1,nch
327  j = jlist(l)
328  IF (j.EQ.0) go to 80
329  zm = zm + body(j)
330  DO 75 k = 1,3
331  xcm(k) = xcm(k) + body(j)*x(k,j)
332  vcm(k) = vcm(k) + body(j)*xdot(k,j)
333  75 CONTINUE
334  80 CONTINUE
335 * Place the c.m. system temporarily in location #ICH.
336  DO 85 k = 1,3
337  x(k,ich) = xcm(k)/zm
338  xdot(k,ich) = vcm(k)/zm
339  85 CONTINUE
340 *
341 * Set individual chain masses to zero and copy sum to #ICH.
342  b1 = body(i1)
343  b2 = body(i2)
344  b3 = body(i3)
345  IF (i4.GT.0) b4 = body(i4)
346  body(i1) = 0.0
347  body(i2) = 0.0
348  body(i3) = 0.0
349  IF (i4.GT.0) body(i4) = b4
350  body(ich) = zm
351 *
352 * Loop over all perturbers (choice of perturber or neighbour list).
353  IF (dmsun.LT.0.2) THEN
354  nnb1 = listc(1) + 1
355  ELSE
356  nnb1 = list(1,ich) + 1
357  END IF
358  DO 95 l = 2,nnb1
359  IF (dmsun.LT.0.2) THEN
360  j = listc(l)
361  ELSE
362  j = list(l,ich)
363  END IF
364  CALL dtchck(time,step(j),dtk(40))
365  DO 90 kk = 1,3
366  x0dot(kk,j) = xdot(kk,j)
367  90 CONTINUE
368  sjr = stepr(j)
369  CALL fpoly1(j,j,0)
370  CALL fpoly2(j,j,0)
371 * Replace possible small regular step with current value above 1D-04.
372  IF (sjr.GT.1.0d-04) stepr(j) = sjr
373  95 CONTINUE
374 *
375 * Restore masses and reference body (needed for CHINIT).
376  body(i1) = b1
377  body(i2) = b2
378  body(i3) = b3
379  IF (i4.GT.0) body(i4) = b4
380  DO 99 k = 1,3
381  x(k,ich) = xrel(k)
382  xdot(k,ich) = vrel(k)
383  99 CONTINUE
384 *
385 * Re-define indices of colliding bodies with J1 as new c.m.
386  100 j1 = i1
387  j2 = i2
388 *
389  RETURN
390 *
391  END