Nbody6
 All Files Functions Variables
trflow.f
Go to the documentation of this file.
1  SUBROUTINE trflow(IPAIR,DTR)
2 *
3 *
4 * Time until Roche overflow.
5 * --------------------------
6 *
7 *
8  include 'common6.h'
9  REAL*8 tscls(20),lums(10),gb(10),tm,tn
10  REAL*8 m0,m1,age,lum,mc,rcc,q,rl1,rl2
11  REAL*8 menv,renv,k2
12  REAL*8 aj,rhigh,thigh,t,tt,rr,rx,tlow,rlow,del,der,eps,tol
13  parameter(eps = 1.0d-06, tol = 1.0d-04)
14  REAL*8 rl,rtmsf,rgbf,rgbdf,ragbf,ragbdf,rzhef
15  EXTERNAL rl,rtmsf,rgbf,rgbdf,ragbf,ragbdf,rzhef
16  LOGICAL itev
17 *
18 *
19 * Set index of first component and semi-major axis.
20  itev = .true.
21  IF(ipair.LT.0)THEN
22  itev = .false.
23  ipair = -ipair
24  ENDIF
25  i1 = 2*ipair - 1
26  i = n + ipair
27  semi = -0.5*body(i)/h(ipair)
28 *
29 * Determine indices for primary & secondary star (donor & accretor).
30  j1 = 2*ipair - 1
31  j2 = i1 + 1
32  q = body(j1)/body(j2)
33  rl1 = rl(q)*semi
34 * Evaluate Roche radius for the second star.
35  q = 1.0/q
36  rl2 = rl(q)*semi
37 *
38 * Compare scaled Roche radii when choosing the primary.
39  IF (radius(j1)/rl1.LT.radius(j2)/rl2) THEN
40  j1 = i1 + 1
41  j2 = i1
42  rl1 = rl2
43  END IF
44 * IF(TEV(J1).LE.0.0)THEN
45 * WRITE(6,99)J1,NAME(J1),NAME(I),KSTAR(J1),TEV(J1)+TOFF,TTOT
46 *99 FORMAT(' TRFLOW WARNING! ',3I6,I4,2F10.3)
47 * ENDIF
48 *
49 * Exit with large interval if semi-major axis is negative.
50  IF(semi.LE.0.d0)THEN
51  dtr = 1.0d+10
52  goto 210
53 * Assign DTR for active Roche during coasting.
54  ELSEIF(kstar(i).GT.0.AND.mod(kstar(i),2).NE.0)THEN
55  IF(tev(i).GT.9.9d+09) tev(i) = tev0(j1)
56  dtr = tev(i) - time
57  goto 210
58  ENDIF
59 *
60 * Convert Roche radius, radius and initial & current mass to SU.
61  rls = rl1*su
62  m0 = body0(j1)*zmbar
63  m1 = body(j1)*zmbar
64  mc = 0.d0
65  kw = kstar(j1)
66  CALL star(kw,m0,m1,tm,tn,tscls,lums,gb,zpars)
67  age = tev0(j1)*tstar - epoch(j1)
68  CALL hrdiag(m0,age,m1,tm,tn,tscls,lums,gb,zpars,
69  & rsj,lum,kw,mc,rcc,menv,renv,k2)
70 *
71 * If the star already fills its Roche lobe exit with DTR = 0.
72 *
73  IF(rsj.GE.rls)THEN
74  dtr = 0.d0
75  radius(j1) = rsj/su
76  goto 210
77  ENDIF
78 *
79 * Exit with large interval if primary has terminated its evolution.
80  IF(kstar(j1).GE.10)THEN
81  dtr = 1.0d+10
82  goto 210
83  ENDIF
84 *
85 * If the star is on the MS then see if it will fill its Roche Lobe
86 * during this stage.
87 *
88  t = 2.0d+10
89  IF(kw.LE.1.OR.kw.EQ.7)THEN
90  aj = 0.99d0*tm
91  CALL hrdiag(m0,aj,m1,tm,tn,tscls,lums,gb,zpars,
92  & rhigh,lum,kw,mc,rcc,menv,renv,k2)
93  IF(rhigh.GT.rls.AND.aj.GT.age)THEN
94 * Find the overflow time using binary chopping.
95  tlow = age
96  rlow = rsj
97  thigh = aj
98  it = 0
99  10 t = 0.5d0*(thigh + tlow)
100  CALL hrdiag(m0,t,m1,tm,tn,tscls,lums,gb,zpars,
101  & rr,lum,kw,mc,rcc,menv,renv,k2)
102  IF(rr.LT.rls)THEN
103  tlow = t
104  rlow = rr
105  ELSE
106 * Stop when 0.1% accuracy is achieved.
107  it = it + 1
108  IF(it.EQ.25)THEN
109  IF(abs(rr-rls).GT.0.1)THEN
110  WRITE(38,*)' TRFLOW KW1: ',rr,rls,t,tm
111  ENDIF
112  ENDIF
113  IF(abs(rr - rls).LE.0.001*rr.OR.it.GT.25) goto 200
114  thigh = t
115  rhigh = rr
116  ENDIF
117  goto 10
118  ENDIF
119  ENDIF
120 *
121 * See if the star will fill its Roche Lobe on the HG if not already
122 * past this stage.
123 *
124  IF(kw.LE.2)THEN
125  aj = min(tn,tscls(1))*(1.d0-eps)
126  CALL hrdiag(m0,aj,m1,tm,tn,tscls,lums,gb,zpars,
127  & rhigh,lum,kw,mc,rcc,menv,renv,k2)
128  IF(rhigh.GT.rls.AND.aj.GT.age)THEN
129 * Solve for the overflow time.
130  rr = rtmsf(m0)
131  t = log(rls/rr)/log(rhigh/rr)
132  t = tm + t*(tscls(1) - tm)
133  goto 200
134  ENDIF
135  IF(tn.LE.tscls(1)) goto 200
136  ENDIF
137 *
138 * If the star is pre-helium ignition see if it will fill its Roche
139 * lobe before helium ignition.
140 *
141  IF(kw.LE.3)THEN
142  aj = min(tn,tscls(2))*(1.d0-eps)
143  CALL hrdiag(m0,aj,m1,tm,tn,tscls,lums,gb,zpars,
144  & rhigh,lum,kw,mc,rcc,menv,renv,k2)
145  IF(rhigh.GE.rls)THEN
146 * Solve for the luminosity corresponding to the Roche radius and
147 * then for the overflow time.
148  it = 0
149  lum = lums(3)
150  30 it = it + 1
151  del = rgbf(m1,lum) - rls
152  IF(it.EQ.25.AND.abs(del).GT.0.1)THEN
153  WRITE(38,*)' TRFLOW KW3: ',rgbf(m1,lum),rls,lum,lums(3)
154  ENDIF
155  IF(abs(del/rls).LE.tol.OR.it.EQ.25) goto 40
156  der = rgbdf(m1,lum)
157  lum = lum - del/der
158  goto 30
159  40 CONTINUE
160  IF(lum.LE.lums(6))THEN
161  t = tscls(4) - (1.d0/((gb(5)-1.d0)*gb(1)*gb(4)))*
162  & ((gb(4)/lum)**((gb(5)-1.d0)/gb(5)))
163  ELSE
164  t = tscls(5) - (1.d0/((gb(6)-1.d0)*gb(1)*gb(3)))*
165  & ((gb(3)/lum)**((gb(6)-1.d0)/gb(6)))
166  ENDIF
167  goto 200
168  ELSE
169 * If a low mass star has not yet filled its Roche Lobe then it will
170 * not do it until after the Helium Flash. Thus we don't let it go
171 * any further until it has actually become a type 4.
172  IF(m0.LT.zpars(2).OR.tn.LE.tscls(2)) goto 200
173  ENDIF
174  ENDIF
175 *
176 * Check for overflow during the CHeB stage.
177 *
178  IF(kw.EQ.4)THEN
179 *
180  IF(tn.LT.(tscls(2)+tscls(3)))THEN
181  t = max(tscls(2),age)
182  aj = t + 0.5d0*(tn - t)
183  CALL hrdiag(m0,aj,m1,tm,tn,tscls,lums,gb,zpars,
184  & rhigh,lum,kw,mc,rcc,menv,renv,k2)
185  IF(rhigh.LT.rls)THEN
186 * If the evolution is due to end during CHeB then the pertubation
187 * functions for small envelope mass will take effect at some point
188 * causing the radius to decrease. We assume that this point is after
189 * AJ and quit with T as an underestimate of the overflow time.
190  t = aj
191  goto 200
192  ENDIF
193  ELSE
194  aj = (tscls(2)+tscls(3))*(1.d0-eps)
195  CALL hrdiag(m0,aj,m1,tm,tn,tscls,lums,gb,zpars,
196  & rhigh,lum,kw,mc,rcc,menv,renv,k2)
197  IF(rhigh.LT.rls) goto 50
198  ENDIF
199 * Find the overflow time using binary chopping.
200  tlow = age
201  rlow = rsj
202  thigh = aj
203  it = 0
204  60 t = 0.5d0*(thigh + tlow)
205  CALL hrdiag(m0,t,m1,tm,tn,tscls,lums,gb,zpars,
206  & rr,lum,kw,mc,rcc,menv,renv,k2)
207  IF(rr.LT.rls)THEN
208  tlow = t
209  rlow = rr
210  ELSE
211 * Stop when 0.1% accuracy is achieved.
212  it = it + 1
213  IF(it.EQ.25)THEN
214  IF(abs(rr-rls).GT.0.1)THEN
215  WRITE(38,*)' TRFLOW KW4: ',rr,rls,tscls(2),t,tscls(3)
216  ENDIF
217  ENDIF
218  IF(abs(rr - rls).LE.0.001*rr.OR.it.GT.25) goto 200
219  thigh = t
220  rhigh = rr
221  ENDIF
222  goto 60
223  50 CONTINUE
224  ENDIF
225 *
226 * The star now has only until the end of the AGB phase to fill its lobe.
227 *
228  IF(kw.LE.6)THEN
229 * Solve for the luminosity corresponding to the Roche radius and
230 * then for the overflow time.
231  it = 0
232  IF(kw.EQ.6)THEN
233  lum = lums(8)
234  ELSE
235  lum = lums(7)
236  ENDIF
237  rr = ragbf(m1,lum,zpars(2))
238  IF(rr.GT.rls)THEN
239 * In this case the solution should already have been found. If it hasn't
240 * then most likely the envelope is small and the pertubation functions
241 * are taking effect.
242  t = tn
243  goto 200
244  ENDIF
245  70 it = it + 1
246  rr = ragbf(m1,lum,zpars(2))
247  del = rr - rls
248  IF(it.EQ.25.AND.abs(del).GT.0.1)THEN
249  WRITE(38,*)' TRFLOW KW6: ',rr,rls,lum,lums(7)
250  ENDIF
251  IF(abs(del/rls).LE.tol.OR.it.EQ.25) goto 80
252  der = ragbdf(m1,lum,zpars(2))
253  lum = lum - del/der
254  goto 70
255  80 CONTINUE
256  IF(lum.LE.lums(8))THEN
257  IF(lum.LE.lums(6))THEN
258  t = tscls(7) - (1.d0/((gb(5)-1.d0)*gb(8)*gb(4)))*
259  & ((gb(4)/lum)**((gb(5)-1.d0)/gb(5)))
260  ELSE
261  t = tscls(8) - (1.d0/((gb(6)-1.d0)*gb(8)*gb(3)))*
262  & ((gb(3)/lum)**((gb(6)-1.d0)/gb(6)))
263  ENDIF
264  ELSE
265  IF(lum.LE.lums(6))THEN
266  t = tscls(10) - (1.d0/((gb(5)-1.d0)*gb(2)*gb(4)))*
267  & ((gb(4)/lum)**((gb(5)-1.d0)/gb(5)))
268  ELSE
269  t = tscls(11) - (1.d0/((gb(6)-1.d0)*gb(2)*gb(3)))*
270  & ((gb(3)/lum)**((gb(6)-1.d0)/gb(6)))
271  ENDIF
272  ENDIF
273  goto 200
274  ENDIF
275 *
276  IF(kw.LE.9)THEN
277 *
278 * The star has until the end of the Helium GB to fill its lobe.
279 *
280  lum = (rls/0.08d0)**(4.0/3.0)
281  cm = 2.0d-03*m1**2.5/(2.d0 + m1**5)
282  IF(cm*lum.GE.100.d0)THEN
283  dtr = 1.0d+10
284  goto 210
285  ENDIF
286  rx = rzhef(m1)
287  rr = rx*(lum/lums(2))**0.2 +
288  & 0.02d0*(exp(cm*lum) - exp(cm*lums(2)))
289  IF(rr.LT.rls)THEN
290 * Find the overflow luminosity using binary chopping.
291  tlow = lum
292  rlow = rr
293  89 lum = 2.0*lum
294  rr = rx*(lum/lums(2))**0.2 +
295  & 0.02d0*(exp(cm*lum) - exp(cm*lums(2)))
296  IF(rr.LE.rls) goto 89
297  thigh = lum
298  rhigh = rr
299  it = 0
300  90 lum = 0.5d0*(thigh + tlow)
301  rr = rx*(lum/lums(2))**0.2 +
302  & 0.02d0*(exp(cm*lum) - exp(cm*lums(2)))
303  IF(rr.LT.rls)THEN
304  tlow = lum
305  rlow = rr
306  ELSE
307 * Stop when 0.1% accuracy is achieved.
308  it = it + 1
309  IF(it.EQ.25)THEN
310  IF(abs(rr-rls).GT.0.1)THEN
311  WRITE(38,*)' TRFLOW KW9: ',rr,rls,t,tm
312  ENDIF
313  ENDIF
314  IF(abs(rr - rls).LE.0.001*rr.OR.it.GT.25) goto 95
315  thigh = lum
316  rhigh = rr
317  ENDIF
318  goto 90
319  95 CONTINUE
320  ENDIF
321  IF(lum.LE.lums(6))THEN
322  t = tscls(4) - (1.d0/((gb(5)-1.d0)*gb(8)*gb(4)))*
323  & ((gb(4)/lum)**((gb(5)-1.d0)/gb(5)))
324  ELSE
325  t = tscls(5) - (1.d0/((gb(6)-1.d0)*gb(8)*gb(3)))*
326  & ((gb(3)/lum)**((gb(6)-1.d0)/gb(6)))
327  ENDIF
328  ENDIF
329 *
330  200 CONTINUE
331  IF(t.GE.tn)THEN
332  dtr = 1.0d+10
333  ELSE
334  tt = (t+epoch(j1))/tstar
335  dtr = tt - time
336  IF(dtr.LT.0.d0.AND.itev)THEN
337  tev(j1) = min(tev(j1),time)
338  tev(j2) = tev(j1)
339  ENDIF
340  dtr = max(dtr,0.d0)
341  IF(tev(j1).LT.tev0(j1))THEN
342  WRITE (6,101)j1,kstar(j1),tev(j1)+toff,ttot
343  101 FORMAT(' TRFLOW WARNING! TEV<TEV0 ',i6,i4,2f10.3)
344  ENDIF
345  ENDIF
346 *
347  210 RETURN
348  END
349 ***