Nbody6
 All Files Functions Variables
hrdiag.f
Go to the documentation of this file.
1 ***
2  SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars,
3  & r,lum,kw,mc,rc,menv,renv,k2)
4 *
5 *
6 * H-R diagram for population I stars.
7 * -----------------------------------
8 *
9 * Computes the new mass, luminosity, radius & stellar type.
10 * Input (MASS, AJ, TM, TN, LUMS & TSCLS) supplied by routine STAR.
11 * Ref: P.P. Eggleton, M.J. Fitchett & C.A. Tout (1989) Ap.J. 347, 998.
12 *
13 * Revised 27th March 1995 by C. A. Tout;
14 * 24th October 1995 to include metallicity;
15 * 14th November 1996 to include naked helium stars;
16 * 28th February 1997 to allow accretion induced supernovae.
17 *
18 * Revised 5th April 1997 by J. R. Hurley
19 * to include from Z=0.0001 to Z=0.03, convective overshooting,
20 * MS hook and more elaborate CHeB.
21 *
22 * Revised 21st January 2011 by A. D. Railton
23 * to include pre-mainsequence evolution for 0.1-8.0 Msun
24 * with solar metallicity. Note KW=-1 for preMS evolution.
25 * Use negative aj for preMS.
26 *
27  implicit none
28 *
29  integer kw,kwp
30  integer wdflag,nsflag
31  parameter(wdflag=1,nsflag=1)
32 *
33  real*8 mass,aj,mt,tm,tn,tscls(20),lums(10),gb(10),zpars(20)
34  real*8 r,lum,mc,rc,menv,renv,k2
35  real*8 mch,mlp,tiny
36  parameter(mch=1.44d0,mlp=12.d0,tiny=1.0d-14)
37  real*8 mxns,mxns0,mxns1
38  parameter(mxns0=1.8d0,mxns1=3.d0)
39  real*8 mass0,mt0,mtc
40 *
41  real*8 thook,thg,tbagb,tau,tloop,taul,tauh,tau1,tau2,dtau,texp
42  real*8 lx,ly,dell,alpha,beta,neta
43  real*8 rx,ry,delr,rzams,rtms,gamma,rmin,taumin,rg
44  parameter(taumin=5.0d-08)
45  real*8 mcmax,mcx,mcy,mcbagb,lambda
46  real*8 am,xx,fac,rdgen,mew,lum0,kap,zeta,ahe,aco
47  parameter(lum0=7.0d+04,kap=-0.5d0,ahe=4.d0,aco=16.d0)
48 *
49  real*8 tprems,pre1,pre2,pre3,pre4
50 *
51  real*8 thookf,tblf
52  real*8 lalphf,lbetaf,lnetaf,lhookf,lgbtf,lmcgbf,lzhef,lpertf
53  real*8 rzamsf,rtmsf,ralphf,rbetaf,rgammf,rhookf
54  real*8 rgbf,rminf,ragbf,rzahbf,rzhef,rhehgf,rhegbf,rpertf
55  real*8 mctmsf,mcgbtf,mcgbf,mcheif,mcagbf,lzahbf
56  external thookf,tblf
57  external lalphf,lbetaf,lnetaf,lhookf,lgbtf,lmcgbf,lzhef,lpertf
58  external rzamsf,rtmsf,ralphf,rbetaf,rgammf,rhookf
59  external rgbf,rminf,ragbf,rzahbf,rzhef,rhehgf,rhegbf,rpertf
60  external mctmsf,mcgbtf,mcgbf,mcheif,mcagbf,lzahbf
61 *
62 *
63 * ---------------------------------------------------------------------
64 * MASS Stellar mass in solar units (input: old; output: new value).
65 * AJ Current age in Myr.
66 * MT Current mass in solar units (used for R).
67 * TM Main sequence time.
68 * TN Nuclear burning time.
69 * TSCLS Time scale for different stages.
70 * LUMS Characteristic luminosity.
71 * GB Giant Branch parameters
72 * ZPARS Parameters for distinguishing various mass intervals.
73 * R Stellar radius in solar units.
74 * TE Effective temperature (suppressed).
75 * KW Classification type (-1 -> 15).
76 * MC Core mass.
77 * ---------------------------------------------------------------------
78 *
79 * Set maximum NS mass depending on which NS mass prescription is used.
80  mxns = mxns0
81  if(nsflag.eq.1) mxns = mxns1
82 *
83  mass0 = mass
84  if(mass0.gt.100.d0) mass = 100.d0
85  mt0 = mt
86  if(mt0.gt.100.d0) mt = 100.d0
87 *
88 * Make evolutionary changes to stars that have not reached KW > 5.
89 *
90  if(kw.gt.6) goto 90
91 *
92  tbagb = tscls(2) + tscls(3)
93  thg = tscls(1) - tm
94 *
95  rzams = rzamsf(mass)
96  rtms = rtmsf(mass)
97 *
98  if(aj.ge.0.d0.and.kw.eq.-1)then
99  if(mass.le.0.7d0)then
100  kw = 0
101  else
102  kw = 1
103  endif
104  endif
105 *
106  if(aj.lt.0.d0)then
107 *
108 * PreMS evolution (valid for 0.1<=M<=8.0).
109 *
110  kw = -1
111  tprems = -1.d0*aj/tscls(15)
112 * Note: tprems cannot exceed 1 - if it does, start at top of Hayashi track.
113  if(tprems.gt.1.d0)then
114  tprems = 1.d0
115  endif
116 *
117  if(mass.le.1.d0)then
118  pre1 = 0.d0
119  pre2 = 0.d0
120  pre3 = 7.432d-02 - 9.43d-02*mass + 7.439d-02*mass**2
121  endif
122  if(mass.gt.1.d0.and.mass.lt.2.d0)then
123  pre1 = -4.00772d0 + 4.00772d0*mass
124  pre2 = 8.5656d0 - 8.5656d0*mass
125  pre3 = -4.50678d0 + 4.56118d0*mass
126  endif
127  if(mass.ge.2.d0)then
128  pre1 = 1.60324d0 + 2.20401d0*mass - 0.60433d0*mass**2 +
129  & 5.172d-02*mass**3
130  pre2 = -4.56878d0 - 4.05305d0*mass + 1.24575*mass**2 -
131  & 0.10922d0*mass**3
132  pre3 = 3.01153 + 1.85745*mass -0.64290d0*mass**2 +
133  & 5.759d-02*mass**3
134  endif
135 *
136  rzams = rzamsf(mass)
137  r = rzams*10.d0**((pre1*tprems**3 + pre2*tprems**4 +
138  & pre3*tprems**5)/(1.05d0-tprems))
139 *
140  pre1 = -2.63181d0 + 3.16607d0*mass - 3.30223d0*mass**2 +
141  & 0.83556d0*mass**3 - 0.06356d0*mass**4
142  pre2 = -11.70230d0 + 16.60510d0*mass - 9.69755d0*mass**2 +
143  & 2.42426d0*mass**3 - 0.27213d0*mass**4 + 0.01134d0*mass**5
144  pre3 = 26.19360d0 - 35.09590d0*mass + 20.64280d0*mass**2 -
145  & 5.18601d0*mass**3 + 0.58360d0*mass**4 - 0.02434d0*mass**5
146  pre4 = -14.64590d0 + 18.55660d0*mass - 10.95070d0*mass**2 +
147  & 2.75979d0*mass**3 - 0.31103d0*mass**4 + 0.01298d0*mass**5
148 *
149  lum = lums(1)*10.d0**((exp(pre1*tprems**2) - 1.d0)*
150  & (pre2*tprems + pre3*tprems**2 + pre4*tprems**3)/
151  & (1.05d0-tprems))
152 *
153  elseif(aj.lt.tscls(1))then
154 *
155 * Either on MS or HG
156 *
157  rg = rgbf(mt,lums(3))
158 *
159  if(aj.lt.tm)then
160 *
161 * Main sequence star.
162 *
163  mc = 0.d0
164  tau = aj/tm
165  thook = thookf(mass)*tscls(1)
166  zeta = 0.01d0
167  tau1 = min(1.d0,aj/thook)
168  tau2 = max(0.d0,
169  & min(1.d0,(aj-(1.d0-zeta)*thook)/(zeta*thook)))
170 *
171  dell = lhookf(mass,zpars(1))
172  dtau = tau1**2 - tau2**2
173  alpha = lalphf(mass)
174  beta = lbetaf(mass)
175  neta = lnetaf(mass)
176  lx = log10(lums(2)/lums(1))
177  if(tau.gt.taumin)then
178  xx = alpha*tau + beta*tau**neta +
179  & (lx - alpha - beta)*tau**2 - dell*dtau
180  else
181  xx = alpha*tau + (lx - alpha)*tau**2 - dell*dtau
182  endif
183  lum = lums(1)*10.d0**xx
184 *
185  delr = rhookf(mass,zpars(1))
186  dtau = tau1**3 - tau2**3
187  alpha = ralphf(mass)
188  beta = rbetaf(mass)
189  gamma = rgammf(mass)
190  rx = log10(rtms/rzams)
191 * Note that the use of taumin is a slightly pedantic attempt to
192 * avoid floating point underflow. It IS overkill!
193  if(tau.gt.taumin)then
194  xx = alpha*tau + beta*tau**10 + gamma*tau**40 +
195  & (rx - alpha - beta - gamma)*tau**3 - delr*dtau
196  else
197  xx = alpha*tau + (rx - alpha)*tau**3 - delr*dtau
198  endif
199  r = rzams*10.d0**xx
200 *
201  if(mass.lt.(zpars(1)-0.3d0))then
202  kw = 0
203 * This following is given by Chris for low mass MS stars which will be
204 * substantially degenerate. We need the Hydrogen abundance, X, which we
205 * calculate from Z assuming that the helium abundance, Y, is calculated
206 * according to Y = 0.24 + 2*Z
207  rdgen = 0.0258d0*((1.d0+zpars(11))**(5.d0/3.d0))*
208  & (mass**(-1.d0/3.d0))
209  r = max(rdgen,r)
210  else
211  kw = 1
212  endif
213 *
214  else
215 *
216 * Star is on the HG
217 *
218  mcx = mc
219  if(mass.le.zpars(2))then
220  mc = mcgbf(lums(3),gb,lums(6))
221  elseif(mass.le.zpars(3))then
222  mc = mcheif(mass,zpars(2),zpars(9))
223  else
224  mc = mcheif(mass,zpars(2),zpars(10))
225  endif
226  neta = mctmsf(mass)
227  tau = (aj - tm)/thg
228  mc = ((1.d0 - tau)*neta + tau)*mc
229  mc = max(mc,mcx)
230 *
231 * Test whether core mass has reached total mass.
232 *
233  if(mc.ge.mt)then
234  aj = 0.d0
235  if(mass.gt.zpars(2))then
236 *
237 * Zero-age helium star
238 *
239  mc = 0.d0
240  mass = mt
241  kw = 7
242  CALL star(kw,mass,mt,tm,tn,tscls,lums,gb,zpars)
243  else
244 *
245 * Zero-age helium white dwarf.
246 *
247  mc = mt
248  mass = mt
249  kw = 10
250  endif
251  else
252  lum = lums(2)*(lums(3)/lums(2))**tau
253  if(mass.le.zpars(3))then
254  rx = rg
255  else
256 * He-ignition and end of HG occur at Rmin
257  rmin = rminf(mass)
258  ry = ragbf(mt,lums(4),zpars(2))
259  rx = min(rmin,ry)
260  if(mass.le.mlp)then
261  texp = log(mass/mlp)/log(zpars(3)/mlp)
262  rx = rg
263  rx = rmin*(rx/rmin)**texp
264  endif
265  tau2 = tblf(mass,zpars(2),zpars(3))
266  if(tau2.lt.tiny) rx = ry
267  endif
268  r = rtms*(rx/rtms)**tau
269  kw = 2
270  endif
271 *
272  endif
273 *
274 * Now the GB, CHeB and AGB evolution.
275 *
276  elseif(aj.lt.tscls(2))then
277 *
278 * Red Giant.
279 *
280  kw = 3
281  lum = lgbtf(aj,gb(1),gb,tscls(4),tscls(5),tscls(6))
282  if(mass.le.zpars(2))then
283 * Star has a degenerate He core which grows on the GB
284  mc = mcgbf(lum,gb,lums(6))
285  else
286 * Star has a non-degenerate He core which may grow, but
287 * only slightly, on the GB
288  tau = (aj - tscls(1))/(tscls(2) - tscls(1))
289  mcx = mcheif(mass,zpars(2),zpars(9))
290  mcy = mcheif(mass,zpars(2),zpars(10))
291  mc = mcx + (mcy - mcx)*tau
292  endif
293  r = rgbf(mt,lum)
294  rg = r
295  if(mc.ge.mt)then
296  aj = 0.d0
297  if(mass.gt.zpars(2))then
298 *
299 * Zero-age helium star
300 *
301  mc = 0.d0
302  mass = mt
303  kw = 7
304  CALL star(kw,mass,mt,tm,tn,tscls,lums,gb,zpars)
305  else
306 *
307 * Zero-age helium white dwarf.
308 *
309  mc = mt
310  mass = mt
311  kw = 10
312  endif
313  endif
314 *
315  elseif(aj.lt.tbagb)then
316 *
317 * Core helium burning star.
318 *
319  if(kw.eq.3.and.mass.le.zpars(2))then
320  mass = mt
321  CALL star(kw,mass,mt,tm,tn,tscls,lums,gb,zpars)
322  aj = tscls(2)
323  endif
324  if(mass.le.zpars(2))then
325  mcx = mcgbf(lums(4),gb,lums(6))
326  else
327  mcx = mcheif(mass,zpars(2),zpars(10))
328  end if
329  tau = (aj - tscls(2))/tscls(3)
330  mc = mcx + (mcagbf(mass) - mcx)*tau
331 *
332  if(mass.le.zpars(2))then
333  lx = lums(5)
334  ly = lums(7)
335  rx = rzahbf(mt,mc,zpars(2))
336  rg = rgbf(mt,lx)
337  rmin = rg*zpars(13)**(mass/zpars(2))
338  texp = min(max(0.4d0,rmin/rx),2.5d0)
339  ry = ragbf(mt,ly,zpars(2))
340  if(rmin.lt.rx)then
341  taul = (log(rx/rmin))**(1.d0/3.d0)
342  else
343  rmin = rx
344  taul = 0.d0
345  endif
346  tauh = (log(ry/rmin))**(1.d0/3.d0)
347  tau2 = taul*(tau - 1.d0) + tauh*tau
348  r = rmin*exp(abs(tau2)**3)
349  rg = rg + tau*(ry - rg)
350  lum = lx*(ly/lx)**(tau**texp)
351  elseif(mass.gt.zpars(3))then
352 *
353 * For HM stars He-ignition takes place at Rmin in the HG, and CHeB
354 * consists of a blue phase (before tloop) and a RG phase (after tloop).
355 *
356  tau2 = tblf(mass,zpars(2),zpars(3))
357  tloop = tscls(2) + tau2*tscls(3)
358  rmin = rminf(mass)
359  rg = rgbf(mt,lums(4))
360  rx = ragbf(mt,lums(4),zpars(2))
361  rmin = min(rmin, rx)
362  if(mass.le.mlp) then
363  texp = log(mass/mlp)/log(zpars(3)/mlp)
364  rx = rg
365  rx = rmin*(rx/rmin)**texp
366  else
367  rx = rmin
368  end if
369  texp = min(max(0.4d0,rmin/rx),2.5d0)
370  lum = lums(4)*(lums(7)/lums(4))**(tau**texp)
371  if(aj.lt.tloop)then
372  ly = lums(4)*(lums(7)/lums(4))**(tau2**texp)
373  ry = ragbf(mt,ly,zpars(2))
374  taul = 0.d0
375  if(abs(rmin-rx).gt.tiny)then
376  taul = (log(rx/rmin))**(1.d0/3.d0)
377  endif
378  tauh = 0.d0
379  if(ry.gt.rmin) tauh = (log(ry/rmin))**(1.d0/3.d0)
380  tau = (aj - tscls(2))/(tau2*tscls(3))
381  tau2 = taul*(tau - 1.d0) + tauh*tau
382  r = rmin*exp(abs(tau2)**3)
383  rg = rg + tau*(ry - rg)
384  else
385  r = ragbf(mt,lum,zpars(2))
386  rg = r
387  end if
388  else
389 *
390 * For IM stars CHeB consists of a RG phase (before tloop) and a blue
391 * loop (after tloop).
392 *
393  tau2 = 1.d0 - tblf(mass,zpars(2),zpars(3))
394  tloop = tscls(2) + tau2*tscls(3)
395  if(aj.lt.tloop)then
396  tau = (tloop - aj)/(tau2*tscls(3))
397  lum = lums(5)*(lums(4)/lums(5))**(tau**3)
398  r = rgbf(mt,lum)
399  rg = r
400  else
401  lx = lums(5)
402  ly = lums(7)
403  rx = rgbf(mt,lx)
404  rmin = rminf(mt)
405  texp = min(max(0.4d0,rmin/rx),2.5d0)
406  ry = ragbf(mt,ly,zpars(2))
407  if(rmin.lt.rx)then
408  taul = (log(rx/rmin))**(1.d0/3.d0)
409  else
410  rmin = rx
411  taul = 0.d0
412  endif
413  tauh = (log(ry/rmin))**(1.d0/3.d0)
414  tau = (aj - tloop)/(tscls(3) - (tloop - tscls(2)))
415  tau2 = taul*(tau - 1.d0) + tauh*tau
416  r = rmin*exp(abs(tau2)**3)
417  rg = rx + tau*(ry - rx)
418  lum = lx*(ly/lx)**(tau**texp)
419  endif
420  endif
421 *
422 * Test whether core mass exceeds total mass.
423 *
424  if(mc.ge.mt)then
425 *
426 * Evolved MS naked helium star.
427 *
428  kw = 7
429  xx = (aj - tscls(2))/tscls(3)
430  mass = mt
431  CALL star(kw,mass,mt,tm,tn,tscls,lums,gb,zpars)
432  aj = xx*tm
433  else
434  kw = 4
435  endif
436 *
437  else
438 *
439 * Asymptotic Red Giant.
440 *
441 * On the AGB the He core mass remains constant until at Ltp it
442 * is caught by the C core mass and they grow together.
443 *
444  mcbagb = mcagbf(mass)
445  mcx = mcgbtf(tbagb,gb(8),gb,tscls(7),tscls(8),tscls(9))
446  mcmax = max(max(mch,0.773d0*mcbagb-0.35d0),1.05d0*mcx)
447 *
448  if(aj.lt.tscls(13))then
449  mcx = mcgbtf(aj,gb(8),gb,tscls(7),tscls(8),tscls(9))
450  mc = mcbagb
451  lum = lmcgbf(mcx,gb)
452  if(mt.le.mc)then
453 *
454 * Evolved naked helium star as the envelope is lost but the
455 * star has not completed its interior burning. The star becomes
456 * a post-HeMS star.
457 *
458  kw = 9
459  mt = mc
460  mass = mt
461  mc = mcx
462  CALL star(kw,mass,mt,tm,tn,tscls,lums,gb,zpars)
463  if(mc.le.gb(7))then
464  aj = tscls(4) - (1.d0/((gb(5)-1.d0)*gb(8)*gb(4)))*
465  & (mc**(1.d0-gb(5)))
466  else
467  aj = tscls(5) - (1.d0/((gb(6)-1.d0)*gb(8)*gb(3)))*
468  & (mc**(1.d0-gb(6)))
469  endif
470  aj = max(aj,tm)
471  goto 90
472  else
473  kw = 5
474  endif
475  else
476  kw = 6
477  mc = mcgbtf(aj,gb(2),gb,tscls(10),tscls(11),tscls(12))
478  lum = lmcgbf(mc,gb)
479 *
480 * Approximate 3rd Dredge-up on AGB by limiting Mc.
481 *
482  lambda = min(0.9d0,0.3d0+0.001d0*mass**5)
483  tau = tscls(13)
484  mcx = mcgbtf(tau,gb(2),gb,tscls(10),tscls(11),tscls(12))
485  mcy = mc
486  mc = mc - lambda*(mcy-mcx)
487  mcx = mc
488  mcmax = min(mt,mcmax)
489  endif
490  r = ragbf(mt,lum,zpars(2))
491  rg = r
492 *
493 * Mc,x represents the C core mass and we now test whether it
494 * exceeds either the total mass or the maximum allowed core mass.
495 *
496  if(mcmax-mcx.lt.tiny)then
497  aj = 0.d0
498  mc = mcmax
499  if(mc.lt.mch)then
500  mt = mc
501  if(mcbagb.lt.1.6d0)then
502 *
503 * Zero-age Carbon/Oxygen White Dwarf
504 *
505  kw = 11
506  else
507 *
508 * Zero-age Oxygen/Neon White Dwarf
509 *
510  kw = 12
511  endif
512  mass = mt
513  else
514  if(mcbagb.lt.1.6d0)then
515 *
516 * Star is not massive enough to ignite C burning.
517 * so no remnant is left after the SN
518 *
519  kw = 15
520  aj = 0.d0
521  mt = 0.d0
522  lum = 1.0d-10
523  r = 1.0d-10
524  else
525  if(nsflag.eq.0)then
526  mt = 1.17d0 + 0.09d0*mc
527  elseif(nsflag.ge.1)then
528 c*
529 c* Use NS/BH mass given by Belczynski et al. 2002, ApJ, 572, 407.
530 c*
531 c if(mc.lt.2.5d0)then
532 c mcx = 0.161767d0*mc + 1.067055d0
533 c else
534 c mcx = 0.314154d0*mc + 0.686088d0
535 c endif
536 c if(mc.le.5.d0)then
537 c mt = mcx
538 c elseif(mc.lt.7.6d0)then
539 c mt = mcx + (mc - 5.d0)*(mt - mcx)/2.6d0
540 c endif
541 * New formula for remnant masses from Eldridge & Tout (MNRAS, 2004).
542  if(mc.lt.6d0)then
543  mcx=1.44d0
544  else
545  mcx=1.4512017*mc -6.5913737d-3*mc*mc -6.1073371
546  endif
547  mt=mcx
548 
549  endif
550  mc = mt
551  if(mt.le.mxns)then
552 *
553 * Zero-age Neutron star
554 *
555  kw = 13
556  else
557 *
558 * Zero-age Black hole
559 *
560  kw = 14
561  endif
562  endif
563  endif
564  endif
565 *
566  endif
567 *
568  90 continue
569 *
570  if(kw.ge.7.and.kw.le.9)then
571 *
572 * Naked Helium Star
573 *
574  rzams = rzhef(mt)
575  rx = rzams
576  if(aj.lt.tm)then
577 *
578 * Main Sequence
579 *
580  kw = 7
581  tau = aj/tm
582  am = max(0.d0,0.85d0-0.08d0*mass)
583  lum = lums(1)*(1.d0+0.45d0*tau+am*tau**2)
584  am = max(0.d0,0.4d0-0.22d0*log10(mt))
585  r = rx*(1.d0+am*(tau-tau**6))
586  rg = rx
587 * Star has no core mass and hence no memory of its past
588 * which is why we subject mass and mt to mass loss for
589 * this phase.
590  mc = 0.d0
591  if(mt.lt.zpars(10)) kw = 10
592  else
593 *
594 * Helium Shell Burning
595 *
596  kw = 8
597  lum = lgbtf(aj,gb(8),gb,tscls(4),tscls(5),tscls(6))
598  r = rhehgf(mt,lum,rx,lums(2))
599  rg = rhegbf(lum)
600  if(r.ge.rg)then
601  kw = 9
602  r = rg
603  endif
604  mc = mcgbf(lum,gb,lums(6))
605  mtc = min(mt,1.45d0*mt-0.31d0)
606  mcmax = min(mtc,max(mch,0.773d0*mass-0.35d0))
607  if(mcmax-mc.lt.tiny)then
608  aj = 0.d0
609  mc = mcmax
610  if(mc.lt.mch)then
611  if(mass.lt.1.6d0)then
612 *
613 * Zero-age Carbon/Oxygen White Dwarf
614 *
615  mt = max(mc,(mc+0.31d0)/1.45d0)
616  kw = 11
617  else
618 *
619 * Zero-age Oxygen/Neon White Dwarf
620 *
621  mt = mc
622  kw = 12
623  endif
624  mass = mt
625  else
626  if(mass.lt.1.6d0)then
627 *
628 * Star is not massive enough to ignite C burning.
629 * so no remnant is left after the SN
630 *
631  kw = 15
632  aj = 0.d0
633  mt = 0.d0
634  lum = 1.0d-10
635  r = 1.0d-10
636  else
637  if(nsflag.eq.0)then
638  mt = 1.17d0 + 0.09d0*mc
639  elseif(nsflag.ge.1)then
640  if(mc.lt.2.5d0)then
641  mcx = 0.161767d0*mc + 1.067055d0
642  else
643  mcx = 0.314154d0*mc + 0.686088d0
644  endif
645  if(mc.le.5.d0)then
646  mt = mcx
647  elseif(mc.lt.7.6d0)then
648  mt = mcx + (mc - 5.d0)*(mt - mcx)/2.6d0
649  endif
650  endif
651  mc = mt
652  if(mt.le.mxns)then
653 *
654 * Zero-age Neutron star
655 *
656  kw = 13
657  else
658 *
659 * Zero-age Black hole
660 *
661  kw = 14
662  endif
663  endif
664  endif
665  endif
666  endif
667  endif
668 *
669  if(kw.ge.10.and.kw.le.12)then
670 *
671 * White dwarf.
672 *
673  mc = mt
674  if(mc.ge.mch.and.kw.eq.12)then
675 *
676 * Accretion induced supernova with no remnant
677 * unless WD is ONe.
678 *
679  kw = 13
680  mt = 1.3d0
681 * mt = 1.17d0 + 0.09d0*mc
682  elseif(mc.ge.(mch-0.06d0).and.kw.le.11)then
683  kw = 15
684  aj = 0.d0
685  mt = 0.d0
686  lum = 1.0d-10
687  r = 1.0d-10
688  else
689 *
690  if(kw.eq.10)then
691  xx = ahe
692  else
693  xx = aco
694  endif
695 *
696  if(wdflag.eq.0)then
697 *
698 * Mestel cooling
699 *
700  lum = 635.d0*mt*zpars(14)/(xx*(aj+0.1d0))**1.4d0
701 *
702  elseif(wdflag.ge.1)then
703 *
704 * modified-Mestel cooling
705 *
706  if(aj.lt.9000.0)then
707  lum = 300.d0*mt*zpars(14)/(xx*(aj+0.1d0))**1.18d0
708  else
709  fac = (9000.1d0*xx)**5.3d0
710  lum = 300.d0*fac*mt*zpars(14)/(xx*(aj+0.1d0))**6.48d0
711  endif
712 *
713  endif
714 *
715  if(mt.lt.0.000005d0)then
716  r = 0.009d0
717  elseif(mt.lt.0.0005d0)then
718  r = 0.09d0
719  else
720  r = 0.0115d0*sqrt(max(1.48204d-06,(mch/mt)**(2.d0/3.d0)
721  & - (mt/mch)**(2.d0/3.d0)))
722  r = min(0.1d0,r)
723  endif
724 *
725  endif
726  endif
727 *
728  if(kw.eq.13)then
729 *
730 * Neutron Star.
731 *
732  mc = mt
733  if(mc.gt.mxns)then
734 *
735 * Accretion induced Black Hole?
736 *
737  kw = 14
738  aj = 0.d0
739  else
740  lum = 0.02d0*(mt**0.67d0)/(max(aj,0.1d0))**2
741  r = 1.4d-05
742  endif
743  endif
744 *
745  if(kw.eq.14)then
746 *
747 * Black hole
748 *
749  mc = mt
750  lum = 1.0d-10
751  r = 4.24d-06*mt
752  endif
753 *
754 * Calculate the core radius and the luminosity and radius of the
755 * remnant that the star will become.
756 *
757  tau = 0.d0
758  if(kw.le.1.or.kw.eq.7)then
759  rc = 0.d0
760  elseif(kw.le.3)then
761  if(mass.gt.zpars(2))then
762  lx = lzhef(mc)
763  rx = rzhef(mc)
764  rc = rx
765  else
766  if(wdflag.eq.0)then
767  lx = 635.d0*mc*zpars(14)/((ahe*0.1d0)**1.4d0)
768  elseif(wdflag.ge.1)then
769  lx = 300.d0*mc*zpars(14)/((ahe*0.1d0)**1.18d0)
770  endif
771  rx = 0.0115d0*sqrt(max(1.48204d-06,
772  & (mch/mc)**(2.d0/3.d0)-(mc/mch)**(2.d0/3.d0)))
773  rc = 5.d0*rx
774  endif
775  elseif(kw.eq.4)then
776  tau = (aj - tscls(2))/tscls(3)
777  kwp = 7
778  CALL star(kwp,mc,mc,tm,tn,tscls,lums,gb,zpars)
779  am = max(0.d0,0.85d0-0.08d0*mc)
780  lx = lums(1)*(1.d0+0.45d0*tau+am*tau**2)
781  rx = rzhef(mc)
782  am = max(0.d0,0.4d0-0.22d0*log10(mc))
783  rx = rx*(1.d0+am*(tau-tau**6))
784  CALL star(kw,mass,mt,tm,tn,tscls,lums,gb,zpars)
785  rc = rx
786  elseif(kw.eq.5)then
787  kwp = 9
788  if(tn.gt.tbagb) tau = 3.d0*(aj-tbagb)/(tn-tbagb)
789  CALL star(kwp,mc,mc,tm,tn,tscls,lums,gb,zpars)
790  lx = lmcgbf(mcx,gb)
791  if(tau.lt.1.d0) lx = lums(2)*(lx/lums(2))**tau
792  rx = rzhef(mc)
793  rx = min(rhehgf(mc,lx,rx,lums(2)),rhegbf(lx))
794  CALL star(kw,mass,mt,tm,tn,tscls,lums,gb,zpars)
795  rc = rx
796  elseif(kw.le.9)then
797  if(wdflag.eq.0)then
798  lx = 635.d0*mc*zpars(14)/((aco*0.1d0)**1.4d0)
799  elseif(wdflag.ge.1)then
800  lx = 300.d0*mc*zpars(14)/((aco*0.1d0)**1.18d0)
801  endif
802  rx = 0.0115d0*sqrt(max(1.48204d-06,
803  & (mch/mc)**(2.d0/3.d0) - (mc/mch)**(2.d0/3.d0)))
804  rc = 5.d0*rx
805  else
806  rc = r
807  menv = 1.0d-10
808  renv = 1.0d-10
809  k2 = 0.21d0
810  endif
811 *
812 * Perturb the luminosity and radius due to small envelope mass.
813 *
814  if(kw.ge.2.and.kw.le.9.and.kw.ne.7)then
815  mew = ((mt-mc)/mt)*min(5.d0,max(1.2d0,(lum/lum0)**kap))
816  if(kw.ge.8) mew = ((mtc-mc)/mtc)*5.d0
817  if(mew.lt.1.d0)then
818  xx = lpertf(mt,mew)
819  lum = lx*(lum/lx)**xx
820  if(r.le.rx)then
821  xx = 0.d0
822  else
823  xx = rpertf(mt,mew,r,rx)
824  endif
825  r = rx*(r/rx)**xx
826  endif
827  rc = min(rc,r)
828  endif
829 *
830 * Calculate mass and radius of convective envelope, and envelope
831 * gyration radius.
832 *
833  if(kw.ge.0.and.kw.lt.10)then
834  CALL mrenv(kw,mass,mt,mc,lum,r,rc,aj,tm,lums(2),lums(3),
835  & lums(4),rzams,rtms,rg,menv,renv,k2)
836  endif
837 *
838  if(mass.gt.99.99d0)then
839  mass = mass0
840  endif
841  if(mt.gt.99.99d0)then
842  mt = mt0
843  endif
844 *
845  return
846  end
847 ***