Nbody6
 All Files Functions Variables
zfuncs.f
Go to the documentation of this file.
1 ***
2  real*8 FUNCTION lzamsf(m)
3  implicit none
4  real*8 m,mx,a(200)
5  common /mscff/ a
6 *
7 * A function to evaluate Lzams
8 * ( from Tout et al., 1996, MNRAS, 281, 257 ).
9 *
10  mx = sqrt(m)
11  lzamsf = (a(1)*m**5*mx + a(2)*m**11)/
12  & (a(3) + m**3 + a(4)*m**5 + a(5)*m**7 +
13  & a(6)*m**8 + a(7)*m**9*mx)
14 *
15  return
16  end
17 ***
18  real*8 FUNCTION rzamsf(m)
19  implicit none
20  real*8 m,mx,a(200)
21  common /mscff/ a
22 *
23 * A function to evaluate Rzams
24 * ( from Tout et al., 1996, MNRAS, 281, 257 ).
25 *
26  mx = sqrt(m)
27  rzamsf = ((a(8)*m**2 + a(9)*m**6)*mx + a(10)*m**11 +
28  & (a(11) + a(12)*mx)*m**19)/
29  & (a(13) + a(14)*m**2 +
30  & (a(15)*m**8 + m**18 + a(16)*m**19)*mx)
31 *
32  return
33  end
34 ***
35  real*8 FUNCTION tbgbf(m)
36  implicit none
37  real*8 m,a(200)
38  common /mscff/ a
39 *
40 * A function to evaluate the lifetime to the BGB or to
41 * Helium ignition if no FGB exists.
42 * (JH 24/11/97)
43 *
44  tbgbf = (a(17) + a(18)*m**4 + a(19)*m**(11.d0/2.d0) + m**7)/
45  & (a(20)*m**2 + a(21)*m**7)
46 *
47  return
48  end
49 ***
50  real*8 FUNCTION tbgbdf(m)
51  implicit none
52  real*8 m,mx,f,df,g,dg,a(200)
53  common /mscff/ a
54 *
55 * A function to evaluate the derivitive of the lifetime to the BGB
56 * (or to Helium ignition if no FGB exists) wrt mass.
57 * (JH 24/11/97)
58 *
59  mx = sqrt(m)
60  f = a(17) + a(18)*m**4 + a(19)*m**5*mx + m**7
61  df = 4.d0*a(18)*m**3 + 5.5d0*a(19)*m**4*mx + 7.d0*m**6
62  g = a(20)*m**2 + a(21)*m**7
63  dg = 2.d0*a(20)*m + 7.d0*a(21)*m**6
64  tbgbdf = (df*g - f*dg)/(g*g)
65 *
66  return
67  end
68 ***
69  real*8 FUNCTION tbgdzf(m)
70  implicit none
71  real*8 m,mx,f,df,g,dg,a(200)
72  common /mscff/ a
73 *
74 * A function to evaluate the derivitive of the lifetime to the BGB
75 * (or to Helium ignition if no FGB exists) wrt Z.
76 * (JH 14/12/98)
77 *
78  mx = m**5*sqrt(m)
79  f = a(17) + a(18)*m**4 + a(19)*mx + m**7
80  df = a(117) + a(118)*m**4 + a(119)*mx
81  g = a(20)*m**2 + a(21)*m**7
82  dg = a(120)*m**2
83  tbgdzf = (df*g - f*dg)/(g*g)
84 *
85  return
86  end
87 ***
88  real*8 FUNCTION thookf(m)
89  implicit none
90  real*8 m,a(200)
91  common /mscff/ a
92 *
93 * A function to evaluate the lifetime to the end of the MS
94 * hook ( for those models that have one ) as a fraction of
95 * the lifetime to the BGB
96 * Note that this function is only valid for M > Mhook.
97 * (JH 24/11/97)
98 *
99  thookf = 1.d0 - 0.01d0*max(a(22)/m**a(23),a(24)+a(25)/m**a(26))
100  thookf = max(thookf,0.5d0)
101 *
102  return
103  end
104 ***
105  real*8 FUNCTION ltmsf(m)
106  implicit none
107  real*8 m,a(200)
108  common /mscff/ a
109 *
110 * A function to evaluate the luminosity at the end of the MS
111 * (JH 24/11/97)
112 *
113  ltmsf = (a(27)*m**3 + a(28)*m**4 + a(29)*m**(a(32)+1.8d0))/
114  & (a(30) + a(31)*m**5 + m**a(32))
115 *
116  return
117  end
118 ***
119  real*8 FUNCTION lalphf(m)
120  implicit none
121  real*8 m,mcut,a(200)
122  common /mscff/ a
123 *
124 * A function to evaluate the Luminosity alpha coefficent.
125 * (JH 24/11/97)
126 *
127  mcut = 2.d0
128  if(m.ge.mcut)then
129  lalphf = (a(33) + a(34)*m**a(36))/(m**0.4d0 + a(35)*m**1.9d0)
130  else
131  if(m.le.0.5d0)then
132  lalphf = a(39)
133  elseif(m.le.0.7d0)then
134  lalphf = a(39) + ((0.3d0 - a(39))/0.2d0)*(m - 0.5d0)
135  elseif(m.le.a(37))then
136  lalphf = 0.3d0 + ((a(40)-0.3d0)/(a(37)-0.7d0))*(m - 0.7d0)
137  elseif(m.le.a(38))then
138  lalphf = a(40) + ((a(41)-a(40))/(a(38)-a(37)))*(m - a(37))
139  else
140  lalphf = a(41) + ((a(42)-a(41))/(mcut-a(38)))*(m - a(38))
141  endif
142  endif
143 *
144  return
145  end
146 ***
147  real*8 FUNCTION lbetaf(m)
148  implicit none
149  real*8 m,a1,a(200)
150  common /mscff/ a
151 *
152 * A function to evaluate the Luminosity beta coefficent.
153 * (JH 24/11/97)
154 *
155  lbetaf = a(43) - a(44)*m**a(45)
156  lbetaf = max(lbetaf,0.d0)
157  if(m.gt.a(46).and.lbetaf.gt.0.d0)then
158  a1 = a(43) - a(44)*a(46)**a(45)
159  lbetaf = a1 - 10.d0*a1*(m - a(46))
160  lbetaf = max(lbetaf,0.d0)
161  endif
162 *
163  return
164  end
165 ***
166  real*8 FUNCTION lnetaf(m)
167  implicit none
168  real*8 m,a(200)
169  common /mscff/ a
170 *
171 * A function to evaluate the Luminosity neta exponent.
172 * (JH 24/11/97)
173 *
174  if(m.le.1.d0)then
175  lnetaf = 10.d0
176  elseif(m.ge.1.1d0)then
177  lnetaf = 20.d0
178  else
179  lnetaf = 10.d0 + 100.d0*(m - 1.d0)
180  endif
181  lnetaf = min(lnetaf,a(97))
182 *
183  return
184  end
185 ***
186  real*8 FUNCTION lhookf(m,mhook)
187  implicit none
188  real*8 m,mhook,a2,a(200)
189  common /mscff/ a
190 *
191 * A function to evalute the luminosity at the start of
192 * the MS hook ( for those stars that have one ).
193 * Note that this function is only valid for M > Mhook.
194 * (JH 24/11/97)
195 *
196  if(m.le.mhook)then
197  lhookf = 0.d0
198  elseif(m.ge.a(51))then
199  lhookf = min(a(47)/m**a(48),a(49)/m**a(50))
200  else
201  a2 = min(a(47)/a(51)**a(48),a(49)/a(51)**a(50))
202  lhookf = a2*((m-mhook)/(a(51)-mhook))**0.4d0
203  endif
204 *
205  return
206  end
207 ***
208  real*8 FUNCTION rtmsf(m)
209  implicit none
210  real*8 m,m2,rchk,a(200)
211  common /mscff/ a
212  real*8 rzamsf
213  external rzamsf
214 *
215 * A function to evaluate the radius at the end of the MS
216 * Note that a safety check is added to ensure Rtms > Rzams
217 * when extrapolating the function to low masses.
218 * (JH 24/11/97)
219 *
220  m2 = a(62) + 0.1d0
221  if(m.le.a(62))then
222  rchk = 1.5d0*rzamsf(m)
223  rtmsf = max(rchk,(a(52) + a(53)*m**a(55))/(a(54) + m**a(56)))
224  elseif(m.ge.m2)then
225  rtmsf = (a(57)*m**3+a(58)*m**a(61)+a(59)*m**(a(61)+1.5d0))/
226  & (a(60) + m**5)
227  else
228  rtmsf = a(63) + ((a(64) - a(63))/0.1d0)*(m - a(62))
229  endif
230 *
231  return
232  end
233 ***
234  real*8 FUNCTION ralphf(m)
235  implicit none
236  real*8 m,a5,a(200)
237  common /mscff/ a
238 *
239 * A function to evaluate the radius alpha coefficent.
240 * (JH 24/11/97)
241 *
242  if(m.le.0.5d0)then
243  ralphf = a(73)
244  elseif(m.le.0.65d0)then
245  ralphf = a(73) + ((a(74) - a(73))/0.15d0)*(m - 0.5d0)
246  elseif(m.le.a(70))then
247  ralphf = a(74) + ((a(75)-a(74))/(a(70)-0.65d0))*(m - 0.65d0)
248  elseif(m.le.a(71))then
249  ralphf = a(75) + ((a(76) - a(75))/(a(71) - a(70)))*(m - a(70))
250  elseif(m.le.a(72))then
251  ralphf = (a(65)*m**a(67))/(a(66) + m**a(68))
252  else
253  a5 = (a(65)*a(72)**a(67))/(a(66) + a(72)**a(68))
254  ralphf = a5 + a(69)*(m - a(72))
255  endif
256 *
257  return
258  end
259 ***
260  real*8 FUNCTION rbetaf(m)
261  implicit none
262  real*8 m,m2,m3,b2,b3,a(200)
263  common /mscff/ a
264 *
265 * A function to evaluate the radius beta coefficent.
266 * (JH 24/11/97)
267 *
268  m2 = 2.d0
269  m3 = 16.d0
270  if(m.le.1.d0)then
271  rbetaf = 1.06d0
272  elseif(m.le.a(82))then
273  rbetaf = 1.06d0 + ((a(81)-1.06d0)/(a(82)-1.d0))*(m-1.d0)
274  elseif(m.le.m2)then
275  b2 = (a(77)*m2**(7.d0/2.d0))/(a(78) + m2**a(79))
276  rbetaf = a(81) + ((b2-a(81))/(m2-a(82)))*(m-a(82))
277  elseif(m.le.m3)then
278  rbetaf = (a(77)*m**(7.d0/2.d0))/(a(78) + m**a(79))
279  else
280  b3 = (a(77)*m3**(7.d0/2.d0))/(a(78) + m3**a(79))
281  rbetaf = b3 + a(80)*(m - m3)
282  endif
283  rbetaf = rbetaf - 1.d0
284 *
285  return
286  end
287 ***
288  real*8 FUNCTION rgammf(m)
289  implicit none
290  real*8 m,m1,b1,a(200)
291  common /mscff/ a
292 *
293 * A function to evaluate the radius gamma coefficent.
294 * (JH 24/11/97)
295 *
296  m1 = 1.d0
297  if(m.gt.(a(88)+0.1d0))then
298  rgammf = 0.d0
299  else
300  b1 = max(0.d0,a(83) + a(84)*(m1-a(85))**a(86))
301  if(m.le.m1)then
302  rgammf = a(83) + a(84)*abs(m-a(85))**a(86)
303  elseif(m.le.a(88))then
304  rgammf = b1 + (a(89) - b1)*((m - m1)/(a(88) - m1))**a(87)
305  else
306  if(a(88).gt.m1) b1 = a(89)
307  rgammf = b1 - 10.d0*b1*(m - a(88))
308  endif
309  rgammf = max(rgammf,0.d0)
310  endif
311 *
312  return
313  end
314 ***
315  real*8 FUNCTION rhookf(m,mhook)
316  implicit none
317  real*8 m,mhook,m2,b2,a(200)
318  common /mscff/ a
319 *
320 * A function to evalute the radius at the start of
321 * the MS hook ( for those stars that have one ).
322 * Note that this function is only valid for M > Mhook.
323 * (JH 24/11/97)
324 *
325  if(m.le.mhook)then
326  rhookf = 0.d0
327  elseif(m.le.a(94))then
328  rhookf = a(95)*sqrt((m-mhook)/(a(94)-mhook))
329  elseif(m.le.2.d0)then
330  m2 = 2.d0
331  b2 = (a(90) + a(91)*m2**(7.d0/2.d0))/
332  & (a(92)*m2**3 + m2**a(93)) - 1.d0
333  rhookf = a(95) + (b2-a(95))*((m-a(94))/(m2-a(94)))**a(96)
334  else
335  rhookf = (a(90) + a(91)*m**(7.d0/2.d0))/
336  & (a(92)*m**3 + m**a(93)) - 1.d0
337  endif
338 *
339  return
340  end
341 ***
342  real*8 FUNCTION lbgbf(m)
343  real*8 m,a(200)
344  common /gbcff/ a
345 *
346 * A function to evaluate the luminosity at the end of the
347 * FGB ( for those models that have one )
348 * Note that this function is only valid for LM & IM stars
349 * (JH 24/11/97)
350 *
351  lbgbf = (a(1)*m**a(5) + a(2)*m**a(8))/
352  & (a(3) + a(4)*m**a(7) + m**a(6))
353 *
354  return
355  end
356 ***
357  real*8 FUNCTION lbgbdf(m)
358  real*8 m,a(200)
359  real*8 f,df,g,dg
360  common /gbcff/ a
361 *
362 * A function to evaluate the derivitive of the Lbgb function.
363 * Note that this function is only valid for LM & IM stars
364 * (JH 24/11/97)
365 *
366  f = a(1)*m**a(5) + a(2)*m**a(8)
367  df = a(5)*a(1)*m**(a(5)-1.d0) + a(8)*a(2)*m**(a(8)-1.d0)
368  g = a(3) + a(4)*m**a(7) + m**a(6)
369  dg = a(7)*a(4)*m**(a(7)-1.d0) + a(6)*m**(a(6)-1.d0)
370 *
371  lbgbdf = (df*g - f*dg)/(g*g)
372 *
373  return
374  end
375 ***
376  real*8 FUNCTION lbagbf(m,mhefl)
377  implicit none
378  real*8 m,mhefl,a4,a(200)
379  common /gbcff/ a
380 *
381 * A function to evaluate the BAGB luminosity. (OP 21/04/98)
382 * Continuity between LM and IM functions is ensured by setting
383 * gbp(16) = lbagbf(mhefl,0.0) with gbp(16) = 1.0.
384 *
385  a4 = (a(9)*mhefl**a(10) - a(16))/(exp(mhefl*a(11))*a(16))
386 *
387  if(m.lt.mhefl)then
388  lbagbf = a(9)*m**a(10)/(1.d0 + a4*exp(m*a(11)))
389  else
390  lbagbf = (a(12) + a(13)*m**(a(15)+1.8d0))/(a(14) + m**a(15))
391  endif
392 *
393  return
394  end
395 ***
396  real*8 FUNCTION rgbf(m,lum)
397  implicit none
398  real*8 m,lum,a1,a(200)
399  common /gbcff/ a
400 *
401 * A function to evaluate radius on the GB.
402 * (JH 24/11/97)
403 *
404  a1 = min(a(20)/m**a(21),a(22)/m**a(23))
405  rgbf = a1*(lum**a(18) + a(17)*lum**a(19))
406 *
407  return
408  end
409 ***
410  real*8 FUNCTION rgbdf(m,lum)
411  implicit none
412  real*8 m,lum,a1,a(200)
413  common /gbcff/ a
414 *
415 * A function to evaluate radius derivitive on the GB (as f(L)).
416 * (JH 24/11/97)
417 *
418  a1 = min(a(20)/m**a(21),a(22)/m**a(23))
419  rgbdf = a1*(a(18)*lum**(a(18)-1.d0) +
420  & a(17)*a(19)*lum**(a(19)-1.d0))
421 *
422  return
423  end
424 ***
425  real*8 FUNCTION ragbf(m,lum,mhelf)
426  implicit none
427  real*8 m,lum,mhelf,m1,a1,a4,xx,a(200)
428  common /gbcff/ a
429 *
430 * A function to evaluate radius on the AGB.
431 * (JH 24/11/97)
432 *
433  m1 = mhelf - 0.2d0
434  if(m.ge.mhelf)then
435  xx = a(24)
436  elseif(m.ge.m1)then
437  xx = 1.d0 + 5.d0*(a(24)-1.d0)*(m-m1)
438  else
439  xx = 1.d0
440  endif
441  a4 = xx*a(19)
442  if(m.le.m1)then
443  a1 = a(29) + a(30)*m
444  elseif(m.ge.mhelf)then
445  a1 = min(a(25)/m**a(26),a(27)/m**a(28))
446  else
447  a1 = a(31) + 5.d0*(a(32)-a(31))*(m-m1)
448  endif
449 *
450  ragbf = a1*(lum**a(18) + a(17)*lum**a4)
451 *
452  return
453  end
454 ***
455  real*8 FUNCTION ragbdf(m,lum,mhelf)
456  implicit none
457  real*8 m,lum,mhelf,m1,a1,a4,xx,a(200)
458  common /gbcff/ a
459 *
460 * A function to evaluate radius derivitive on the AGB (as f(L)).
461 * (JH 24/11/97)
462 *
463  m1 = mhelf - 0.2d0
464  if(m.ge.mhelf)then
465  xx = a(24)
466  elseif(m.ge.m1)then
467  xx = 1.d0 + 5.d0*(a(24)-1.d0)*(m-m1)
468  else
469  xx = 1.d0
470  endif
471  a4 = xx*a(19)
472  if(m.le.m1)then
473  a1 = a(29) + a(30)*m
474  elseif(m.ge.mhelf)then
475  a1 = min(a(25)/m**a(26),a(27)/m**a(28))
476  else
477  a1 = a(31) + 5.d0*(a(32)-a(31))*(m-m1)
478  endif
479 *
480  ragbdf = a1*(a(18)*lum**(a(18)-1.d0) +
481  & a(17)*a4*lum**(a4-1.d0))
482 *
483  return
484  end
485 ***
486  real*8 FUNCTION mctmsf(m)
487  implicit none
488  real*8 m,m525
489 *
490 * A function to evaluate core mass at the end of the MS as a
491 * fraction of the BGB value, i.e. this must be multiplied by
492 * the BGB value (see below) to give the actual core mass (JH 5/9/99)
493 *
494  m525 = m**(21.d0/4.d0)
495  mctmsf = (1.586d0 + m525)/(2.434d0 + 1.02d0*m525)
496 *
497  return
498  end
499 ***
500  real*8 FUNCTION mcheif(m,mhefl,mchefl)
501  implicit none
502  real*8 m,mhefl,mchefl,mcbagb,a3,a(200)
503  common /gbcff/ a
504  real*8 mcagbf
505  external mcagbf
506 *
507 * A function to evaluate core mass at BGB or He ignition
508 * (depending on mchefl) for IM & HM stars (OP 25/11/97)
509 *
510  mcbagb = mcagbf(m)
511  a3 = mchefl**4 - a(33)*mhefl**a(34)
512  mcheif = min(0.95d0*mcbagb,(a3 + a(33)*m**a(34))**(1.d0/4.d0))
513 *
514  return
515  end
516 ***
517  real*8 FUNCTION mheif(mc,mhefl,mchefl)
518  implicit none
519  real*8 mc,mhefl,mchefl,m1,m2,a3,a(200)
520  common /gbcff/ a
521  real*8 mbagbf
522  external mbagbf
523 *
524 * A function to evaluate mass at BGB or He ignition
525 * (depending on mchefl) for IM & HM stars by inverting
526 * mcheif
527 *
528  m1 = mbagbf(mc/0.95d0)
529  a3 = mchefl**4 - a(33)*mhefl**a(34)
530  m2 = ((mc**4 - a3)/a(33))**(1.d0/a(34))
531  mheif = max(m1,m2)
532 *
533  return
534  end
535 ***
536  real*8 FUNCTION mcagbf(m)
537  implicit none
538  real*8 m,a(200)
539  common /gbcff/ a
540 *
541 * A function to evaluate core mass at the BAGB (OP 25/11/97)
542 *
543  mcagbf = (a(37) + a(35)*m**a(36))**(1.d0/4.d0)
544 *
545  return
546  end
547 ***
548  real*8 FUNCTION mbagbf(mc)
549  implicit none
550  real*8 mc,mc4,a(200)
551  common /gbcff/ a
552 *
553 * A function to evaluate mass at the BAGB by inverting mcagbf.
554 *
555  mc4 = mc**4
556  if(mc4.gt.a(37))then
557  mbagbf = ((mc4 - a(37))/a(35))**(1.d0/a(36))
558  else
559  mbagbf = 0.d0
560  endif
561 *
562  return
563  end
564 ***
565  real*8 FUNCTION mcgbtf(t,A,GB,tinf1,tinf2,tx)
566  implicit none
567  real*8 t,a,gb(10),tinf1,tinf2,tx
568 *
569 * A function to evaluate Mc given t for GB, AGB and NHe stars
570 *
571  if(t.le.tx)then
572  mcgbtf = ((gb(5)-1.d0)*a*gb(4)*(tinf1 - t))**
573  & (1.d0/(1.d0-gb(5)))
574  else
575  mcgbtf = ((gb(6)-1.d0)*a*gb(3)*(tinf2 - t))**
576  & (1.d0/(1.d0-gb(6)))
577  endif
578 *
579  return
580  end
581 ***
582  real*8 FUNCTION lgbtf(t,A,GB,tinf1,tinf2,tx)
583  implicit none
584  real*8 t,a,gb(10),tinf1,tinf2,tx
585 *
586 * A function to evaluate L given t for GB, AGB and NHe stars
587 *
588  if(t.le.tx)then
589  lgbtf = gb(4)*(((gb(5)-1.d0)*a*gb(4)*(tinf1 - t))**
590  & (gb(5)/(1.d0-gb(5))))
591  else
592  lgbtf = gb(3)*(((gb(6)-1.d0)*a*gb(3)*(tinf2 - t))**
593  & (gb(6)/(1.d0-gb(6))))
594  endif
595 *
596  return
597  end
598 ***
599  real*8 FUNCTION mcgbf(lum,GB,lx)
600  implicit none
601  real*8 lum,gb(10),lx
602 *
603 * A function to evaluate Mc given L for GB, AGB and NHe stars
604 *
605  if(lum.le.lx)then
606  mcgbf = (lum/gb(4))**(1.d0/gb(5))
607  else
608  mcgbf = (lum/gb(3))**(1.d0/gb(6))
609  endif
610 *
611  return
612  end
613 ***
614  real*8 FUNCTION lmcgbf(mc,GB)
615  implicit none
616  real*8 mc,gb(10)
617 *
618 * A function to evaluate L given Mc for GB, AGB and NHe stars
619 *
620  if(mc.le.gb(7))then
621  lmcgbf = gb(4)*(mc**gb(5))
622  else
623  lmcgbf = gb(3)*(mc**gb(6))
624  endif
625 *
626  return
627  end
628 ***
629  real*8 FUNCTION lheif(m,mhefl)
630  implicit none
631  real*8 m,mhefl,a(200)
632  common /gbcff/ a
633 *
634 * A function to evaluate He-ignition luminosity (OP 24/11/97)
635 * Continuity between the LM and IM functions is ensured with a first
636 * call setting lhefl = lHeIf(mhefl,0.0)
637 *
638  if(m.lt.mhefl)then
639  lheif = a(38)*m**a(39)/(1.d0 + a(41)*exp(m*a(40)))
640  else
641  lheif = (a(42) + a(43)*m**3.8d0)/(a(44) + m**2)
642  endif
643 *
644  return
645  end
646 ***
647  real*8 FUNCTION lhef(m)
648  implicit none
649  real*8 m,a(200)
650  common /gbcff/ a
651 *
652 * A function to evaluate the ratio LHe,min/LHeI (OP 20/11/97)
653 * Note that this function is everywhere <= 1, and is only valid
654 * for IM stars
655 *
656  lhef = (a(45) + a(46)*m**(a(48)+0.1d0))/(a(47) + m**a(48))
657 *
658  return
659  end
660 ***
661  real*8 FUNCTION rminf(m)
662  implicit none
663  real*8 m,mx,a(200)
664  common /gbcff/ a
665 *
666 * A function to evaluate the minimum radius during He-burning
667 * for IM & HM stars (OP 20/11/97)
668 *
669  mx = m**a(53)
670  rminf = (a(49)*m + (a(50)*m)**a(52)*mx)/(a(51) + mx)
671 *
672  return
673  end
674 ***
675  real*8 FUNCTION thef(m,mc,mhefl)
676  implicit none
677  real*8 m,mc,mhefl,mm,a(200)
678  common /gbcff/ a
679  real*8 themsf
680  external themsf
681 *
682 * A function to evaluate the He-burning lifetime. (OP 26/11/97)
683 * For IM & HM stars, tHef is relative to tBGB.
684 * Continuity between LM and IM stars is ensured by setting
685 * thefl = tHef(mhefl,0.0,,0.0), and the call to themsf ensures
686 * continuity between HB and NHe stars as Menv -> 0.
687 *
688  if(m.le.mhefl)then
689  mm = max((mhefl - m)/(mhefl - mc),1.0d-12)
690  thef = (a(54) + (themsf(mc) - a(54))*mm**a(55))*
691  & (1.d0 + a(57)*exp(m*a(56)))
692  else
693  mm = m**5
694  thef = (a(58)*m**a(61) + a(59)*mm)/(a(60) + mm)
695  endif
696 *
697  return
698  end
699 ***
700  real*8 FUNCTION tblf(m,mhefl,mfgb)
701  implicit none
702  real*8 m,mhefl,mfgb,mr,m1,m2,r1,a(200)
703  common /gbcff/ a
704  real*8 lheif,rminf,ragbf
705  external lheif,rminf,ragbf
706 *
707 * A function to evaluate the blue-loop fraction of the He-burning
708 * lifetime for IM & HM stars (OP 28/01/98)
709 *
710  mr = mhefl/mfgb
711  if(m.le.mfgb) then
712  m1 = m/mfgb
713  m2 = log10(m1)/log10(mr)
714  m2 = max(m2,1.0d-12)
715  tblf = a(64)*m1**a(63) + a(65)*m2**a(62)
716  else
717  r1 = 1.d0 - rminf(m)/ragbf(m,lheif(m,mhefl),mhefl)
718  r1 = max(r1,1.0d-12)
719  tblf = a(66)*m**a(67)*r1**a(68)
720  end if
721  tblf = min(1.d0,max(0.d0,tblf))
722  if(tblf.lt.1.0d-10) tblf = 0.d0
723 *
724  return
725  end
726 ***
727  real*8 FUNCTION lzahbf(m,mc,mhefl)
728  implicit none
729  real*8 m,mc,mhefl,mm,a4,a5,a(200)
730  common /gbcff/ a
731  real*8 lzhef
732  external lzhef
733 *
734 * A function to evaluate the ZAHB luminosity for LM stars. (OP 28/01/98)
735 * Continuity with LHe,min for IM stars is ensured by setting
736 * lx = lHeif(mhefl,z,0.0,1.0)*lHef(mhefl,z,mfgb), and the call to lzhef
737 * ensures continuity between the ZAHB and the NHe-ZAMS as Menv -> 0.
738 *
739  a5 = lzhef(mc)
740  a4 = (a(69) + a5 - a(74))/((a(74) - a5)*exp(a(71)*mhefl))
741  mm = max((m-mc)/(mhefl - mc),1.0d-12)
742  lzahbf = a5 + (1.d0 + a(72))*a(69)*mm**a(70)/
743  & ((1.d0 + a(72)*mm**a(73))*(1.d0 + a4*exp(m*a(71))))
744 *
745  return
746  end
747 ***
748  real*8 FUNCTION rzahbf(m,mc,mhefl)
749  implicit none
750  real*8 m,mc,mhefl,rx,ry,mm,f,a(200)
751  common /gbcff/ a
752  real*8 rzhef,rgbf,lzahbf
753 *
754 * A function to evaluate the ZAHB radius for LM stars. (OP 28/01/98)
755 * Continuity with R(LHe,min) for IM stars is ensured by setting
756 * lx = lHeif(mhefl,z,0.0,1.0)*lHef(mhefl,z,mfgb), and the call to rzhef
757 * ensures continuity between the ZAHB and the NHe-ZAMS as Menv -> 0.
758 *
759  rx = rzhef(mc)
760  ry = rgbf(m,lzahbf(m,mc,mhefl))
761  mm = max((m-mc)/(mhefl - mc),1.0d-12)
762  f = (1.d0 + a(76))*mm**a(75)/(1.d0 + a(76)*mm**a(77))
763  rzahbf = (1.d0 - f)*rx + f*ry
764 *
765  return
766  end
767 ***
768  real*8 FUNCTION lzhef(m)
769  implicit none
770  real*8 m,m15
771 *
772 * A function to evaluate Naked Helium star 'ZAMS' luminosity
773 *
774  m15 = m*sqrt(m)
775  lzhef = 1.5262d+04*m**(41.d0/4.d0)/
776  & (0.0469d0 + m**6*(31.18d0 + m15*(29.54d0 + m15)))
777 *
778  return
779  end
780 ***
781  real*8 FUNCTION rzhef(m)
782  implicit none
783  real*8 m
784 *
785 * A function to evaluate Helium star 'ZAMS' radius
786 *
787  rzhef = 0.2391d0*m**4.6d0/(0.0065d0 + (0.162d0 + m)*m**3)
788 *
789  return
790  end
791 ***
792  real*8 FUNCTION themsf(m)
793  implicit none
794  real*8 m
795 *
796 * A function to evaluate Helium star main sequence lifetime
797 *
798  themsf = (0.4129d0 + 18.81d0*m**4 + 1.853d0*m**6)/m**(13.d0/2.d0)
799 *
800  return
801  end
802 ***
803  real*8 FUNCTION rhehgf(m,lum,rx,lx)
804  implicit none
805  real*8 m,lum,rx,lx,cm
806 *
807 * A function to evaluate Helium star radius on the Hertzsprung gap
808 * from its mass and luminosity.
809 *
810  cm = 2.0d-03*m**(5.d0/2.d0)/(2.d0 + m**5)
811  rhehgf = rx*(lum/lx)**0.2d0 + 0.02d0*(exp(cm*lum) - exp(cm*lx))
812 *
813  return
814  end
815 ***
816  real*8 FUNCTION rhegbf(lum)
817  implicit none
818  real*8 lum
819 *
820 * A function to evaluate Helium star radius on the giant branch.
821 *
822  rhegbf = 0.08d0*lum**(3.d0/4.d0)
823 *
824  return
825  end
826 ***
827  real*8 FUNCTION lpertf(m,mew)
828  implicit none
829  real*8 m,mew
830  real*8 b,c
831 *
832 * A function to obtain the exponent that perturbs luminosity.
833 *
834  b = 0.002d0*max(1.d0,2.5d0/m)
835  c = 3.d0
836  lpertf = ((1.d0 + b**c)*((mew/b)**c))/(1.d0+(mew/b)**c)
837 *
838  return
839  end
840 ***
841  real*8 FUNCTION rpertf(m,mew,r,rc)
842  implicit none
843  real*8 m,mew,r,rc
844  real*8 a,b,c,q
845 *
846 * A function to obtain the exponent that perturbs radius.
847 *
848  a = 0.1d0
849  b = 0.006d0*max(1.d0,2.5d0/m)
850  c = 3.d0
851  q = log(r/rc)
852  rpertf = ((1.d0 + b**c)*((mew/b)**c)*(mew**(a/q)))/
853  & (1.d0+(mew/b)**c)
854 *
855  return
856  end
857 ***
858  real*8 FUNCTION vrotf(m)
859  implicit none
860  real*8 m
861 *
862  vrotf = 330.d0*m**3.3d0/(15.d0 + m**3.45d0)
863 *
864  return
865  end
866 ***
867  real*8 FUNCTION celamf(kw,m,lum,rad,rzams,menv,fac)
868  implicit none
869  integer kw
870  real*8 m,lum,rad,rzams,menv,fac
871  real*8 lam1,lam2,m1,logm,logl
872  real*8 aa,bb,cc,dd
873 *
874 * A function to estimate lambda for common-envelope.
875 *
876  if(fac.ge.0.d0)then
877 *
878 * No fits yet for naked He stars...
879 *
880  if(kw.gt.6)then
881  celamf = 0.5d0
882  goto 90
883  endif
884 *
885  if(menv.gt.0.d0)then
886 * Formulae for giant-like stars; also used for HG and CHeB stars close
887 * to the Hayashi track.
888  logl = log10(lum)
889  logm = log10(m)
890  if(kw.le.5)then
891  m1 = m
892  if(kw.gt.3) m1 = 100.d0
893  lam1 = 3.d0/(2.4d0 + 1.d0/m1**(3.d0/2.d0)) - 0.15d0*logl
894  lam1 = min(lam1,0.8d0)
895  else
896  lam1 = -3.5d0 - 0.75d0*logm + logl
897  endif
898  if(kw.gt.3)then
899  lam2 = min(0.9d0,0.58d0 + 0.75d0*logm) - 0.08d0*logl
900  if(kw.lt.6)then
901  lam1 = min(lam2,lam1)
902  else
903  lam1 = max(lam2,lam1)
904  lam1 = min(lam1,1.d0)
905  endif
906  endif
907  lam1 = 2.d0*lam1
908  if(fac.gt.0.d0)then
909 * Use a fraction FAC of the ionization energy in the energy balance.
910  if(kw.le.3)then
911  aa = min(1.2d0*(logm - 0.25d0)**2 - 0.7d0,-0.5d0)
912  else
913  aa = max(-0.2d0 - logm,-0.5d0)
914  endif
915  bb = max(3.d0 - 5.d0*logm,1.5d0)
916  cc = max(3.7d0 + 1.6d0*logm,3.3d0 + 2.1d0*logm)
917  lam2 = aa + atan(bb*(cc - logl))
918  if(kw.le.3)then
919  dd = max(0.d0,min(0.15d0,0.15d0 - 0.25d0*logm))
920  lam2 = lam2 + dd*(logl - 2.d0)
921  endif
922  lam2 = max(lam2,1.d-2)
923  lam2 = max(1.d0/lam2,lam1)
924  if(fac.ge.1.d0)then
925  lam1 = lam2
926  else
927  lam1 = lam1 + fac*(lam2 - lam1)
928  endif
929  endif
930  endif
931 *
932  if(menv.lt.1.d0)then
933 * Formula for HG stars; also reasonable for CHeB stars in blue loop.
934  lam2 = 0.42d0*(rzams/rad)**0.4d0
935 * Alternatively:
936 * lam2 = 0.3d0*(rtms/rad)**0.4d0
937  lam2 = 2.d0*lam2
938  endif
939 *
940  if(menv.le.0.d0)then
941  celamf = lam2
942  elseif(menv.ge.1.d0)then
943  celamf = lam1
944  else
945 * Interpolate between HG and GB values depending on conv. envelope mass.
946  celamf = lam2 + sqrt(menv)*(lam1 - lam2)
947  endif
948 *
949  90 continue
950 *
951  else
952  celamf = -1.d0*fac
953  endif
954 *
955  return
956  end
957 ***