cc// deltat.f *** SUBROUTINE deltat(kw,age,tm,tn,tscls,dt,dtr) implicit none * INTEGER kw REAL*8 age,tm,tn,tscls(20) REAL*8 dt,dtr REAL*8 pts1,pts2,pts3 COMMON /POINTS/ pts1,pts2,pts3 * * Base new time scale for changes in radius & mass on stellar type. * if(kw.le.1)then dt = pts1*tm dtr = tm - age elseif(kw.eq.2)then dt = pts1*(tscls(1) - tm) dtr = tscls(1) - age elseif(kw.eq.3)then if(age.lt.tscls(6))then dt = pts2*(tscls(4) - age) else dt = pts2*(tscls(5) - age) endif dtr = MIN(tscls(2),tn) - age elseif(kw.eq.4)then dt = pts2*tscls(3) dtr = MIN(tn,tscls(2) + tscls(3)) - age elseif(kw.eq.5)then if(age.lt.tscls(9))then dt = pts3*(tscls(7) - age) else dt = pts3*(tscls(8) - age) endif dtr = MIN(tn,tscls(13)) - age elseif(kw.eq.6)then if(age.lt.tscls(12))then dt = pts3*(tscls(10) - age) else dt = pts3*(tscls(11) - age) endif dt = MIN(dt,0.005d0) dtr = tn - age elseif(kw.eq.7)then dt = pts1*tm dtr = tm - age elseif(kw.eq.8.or.kw.eq.9)then if(age.lt.tscls(6))then dt = pts2*(tscls(4) - age) else dt = pts2*(tscls(5) - age) endif dtr = tn - age else * dt = MAX(0.1d0,age*10.d0) dt = MAX(0.1d0,dt*10.d0) dt = MIN(dt,5.0d+02) dtr = dt endif * RETURN END *** cc//evolv1.f *** SUBROUTINE evolv1(kw,mass,mt,r,lum,mc,rc,menv,renv,ospin, & epoch,tm,tphys,tphysf,dtp,z,zpars,vkick,vs) c-------------------------------------------------------------c c c Evolves a single star. c Mass loss is an option. c The timestep is not constant but determined by certain criteria. c Plots the HRD and variables as a function of time. c c Written by Jarrod Hurley 26/08/97 at the Institute of c Astronomy, Cambridge. c c-------------------------------------------------------------c c c STELLAR TYPES - KW c c 0 - deeply or fully convective low mass MS star c 1 - Main Sequence star c 2 - Hertzsprung Gap c 3 - First Giant Branch c 4 - Core Helium Burning c 5 - First Asymptotic Giant Branch c 6 - Second Asymptotic Giant Branch c 7 - Main Sequence Naked Helium star c 8 - Hertzsprung Gap Naked Helium star c 9 - Giant Branch Naked Helium star c 10 - Helium White Dwarf c 11 - Carbon/Oxygen White Dwarf c 12 - Oxygen/Neon White Dwarf c 13 - Neutron Star c 14 - Black Hole c 15 - Massless Supernova c c-------------------------------------------------------------c implicit none * integer kw,it,ip,jp,j,kwold,rflag integer nv parameter(nv=50000) * real*8 mass,z,aj real*8 epoch,tphys,tphys2,tmold,tbgold real*8 mt,tm,tn,tphysf,dtp,tsave real*8 tscls(20),lums(10),GB(10),zpars(20) real*8 r,lum,mc,teff,rc,menv,renv,vs(3) real*8 ospin,jspin,djt,djmb,k2,k3,vkick parameter(k3=0.21d0) real*8 m0,r1,lum1,mc1,rc1,menv1,renv1,k21 real*8 dt,dtm,dtr,dr,dtdr,dms,dml,mt2,rl real*8 tol,tiny parameter(tol=1.0d-10,tiny=1.0d-14) real*8 ajhold,rm0,eps,alpha2 parameter(eps=1.0d-06,alpha2=0.09d0) real*8 mlwind,vrotf external mlwind,vrotf logical iplot,isave REAL*8 neta,bwind,hewind,mxns COMMON /VALUE1/ neta,bwind,hewind,mxns REAL*8 pts1,pts2,pts3 COMMON /POINTS/ pts1,pts2,pts3 REAL scm(50000,14),spp(20,3) COMMON /SINGLE/ scm,spp * dtm = 0.d0 r = 0.d0 lum = 0.d0 mc = 0.d0 mc1 = 0.d0 rc = 0.d0 rl = 0.d0 if(ospin.le.0.d0)then ospin = 1.0d-10 jspin = 1.0d-10 endif k2 = 0.15d0 rflag = 0 * * Setup variables which control the output (if it is required). * ip = 0 jp = 0 tsave = tphys isave = .true. iplot = .false. if(dtp.le.0.d0)then iplot = .true. isave = .false. tsave = tphysf elseif(dtp.gt.tphysf)then isave = .false. tsave = tphysf endif * do 10 , j = 1,nv * if(neta.gt.tiny.and.j.gt.1)then * * Calculate mass loss from the previous timestep. * dt = 1.0d+06*dtm dms = mlwind(kw,lum,r,mt,mc,rl,z)*dt if(kw.lt.10)then dml = mt - mc if(dml.lt.dms)then dtm = (dml/dms)*dtm dms = dml endif endif else dms = 0.d0 endif * * Limit to 1% mass loss. * if(dms.gt.0.01d0*mt)then dtm = 0.01d0*mt*dtm/dms dms = 0.01d0*mt endif * * Calculate the rate of angular momentum loss due to magnetic braking * and/or mass loss. * if(j.gt.1)then djt = (2.d0/3.d0)*(dms/(1.0d+06*dtm))*r*r*ospin if(mt.gt.0.35d0.and.kw.lt.10)then djmb = 5.83d-16*menv*(r*ospin)**3/mt djt = djt + djmb endif endif * * Update mass and time and reset epoch for a MS (and possibly a HG) star. * if(dms.gt.0.d0)then mt = mt - dms if(kw.le.2.or.kw.eq.7)then m0 = mass mc1 = mc mass = mt tmold = tm tbgold = tscls(1) CALL star(kw,mass,mt,tm,tn,tscls,lums,GB,zpars) if(kw.eq.2)then if(GB(9).lt.mc1.or.m0.gt.zpars(3))then mass = m0 else epoch = tm + (tscls(1) - tm)*(ajhold-tmold)/ & (tbgold - tmold) epoch = tphys - epoch endif else epoch = tphys - ajhold*tm/tmold endif endif endif tphys2 = tphys tphys = tphys + dtm * * Find the landmark luminosities and timescales as well as setting * the GB parameters. * aj = tphys - epoch CALL star(kw,mass,mt,tm,tn,tscls,lums,GB,zpars) * * Find the current radius, luminosity, core mass and stellar type * given the initial mass, current mass, metallicity and age * kwold = kw CALL hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars, & r,lum,kw,mc,rc,menv,renv,k2) * * If mass loss has occurred and no type change then check that we * have indeed limited the radius change to 10%. * if(kw.eq.kwold.and.dms.gt.0.d0.and.rflag.ne.0)then mt2 = mt + dms dml = dms/dtm it = 0 20 dr = r - rm0 if(ABS(dr).gt.0.1d0*rm0)then it = it + 1 if(it.eq.20.and.kw.eq.4) goto 30 if(it.gt.30)then WRITE(99,*)' DANGER1! ',it,kw,mass,dr,rm0 WRITE(*,*)' STOP: EVOLV1 FATAL ERROR ' CALL exit(0) STOP endif dtdr = dtm/ABS(dr) dtm = alpha2*MAX(rm0,r)*dtdr if(it.ge.20) dtm = 0.5d0*dtm if(dtm.lt.1.0d-07*aj) goto 30 dms = dtm*dml mt = mt2 - dms if(kw.le.2.or.kw.eq.7)then mass = mt CALL star(kw,mass,mt,tm,tn,tscls,lums,GB,zpars) if(kw.eq.2)then if(GB(9).lt.mc1.or.m0.gt.zpars(3))then mass = m0 else epoch = tm + (tscls(1) - tm)*(ajhold-tmold)/ & (tbgold - tmold) epoch = tphys2 - epoch endif else epoch = tphys2 - ajhold*tm/tmold endif endif tphys = tphys2 + dtm aj = tphys - epoch mc = mc1 CALL star(kw,mass,mt,tm,tn,tscls,lums,GB,zpars) CALL hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars, & r,lum,kw,mc,rc,menv,renv,k2) goto 20 endif 30 continue endif * * Initialize or adjust the spin of the star. * if(j.eq.1)then if(tphys.lt.tiny.and.ospin.lt.0.001d0)then ospin = 45.35d0*vrotf(mt)/r endif jspin = ospin*(k2*r*r*(mt-mc)+k3*rc*rc*mc) else jspin = MAX(1.0d-10,jspin - djt*1.0d+06*dtm) ospin = jspin/(k2*r*r*(mt-mc)+k3*rc*rc*mc) endif * * Test for changes in evolution type. * if(j.eq.1.or.kw.ne.kwold)then * * Force new NS or BH to have a one second period. * if(kw.eq.13.or.kw.eq.14)then ospin = 2.0d+08 jspin = k3*rc*rc*mc*ospin CALL kick(kw,mass,mt,0.d0,0.d0,-1.d0,0.d0,vs) vkick = dsqrt(vs(1)*vs(1)+vs(2)*vs(2)+vs(3)*vs(3)) endif jp = jp + 1 spp(jp,1) = tphys spp(jp,2) = float(kw) if(kw.eq.15)then spp(jp,3) = mass goto 90 else spp(jp,3) = mt endif endif * * Record values for plotting and reset epoch. * epoch = tphys - aj if((isave.and.tphys.ge.tsave).or.iplot)then ip = ip + 1 scm(ip,1) = tphys scm(ip,2) = float(kw) scm(ip,3) = mass scm(ip,4) = mt scm(ip,5) = log10(lum) scm(ip,6) = log10(r) teff = 1000.d0*((1130.d0*lum/(r**2.d0))**(1.d0/4.d0)) scm(ip,7) = log10(teff) scm(ip,8) = mc scm(ip,9) = rc scm(ip,10) = menv scm(ip,11) = renv scm(ip,12) = epoch scm(ip,13) = ospin if(isave) tsave = tsave + dtp if(tphysf.lt.tiny)then ip = ip + 1 do 35 , it = 1,13 scm(ip,it) = scm(ip-1,it) 35 continue endif endif * if(tphys.ge.tphysf)then jp = jp + 1 spp(jp,1) = tphys spp(jp,2) = float(kw) spp(jp,3) = mt goto 90 endif * * Record radius and current age. * rm0 = r ajhold = aj if(kw.ne.kwold) kwold = kw CALL deltat(kw,aj,tm,tn,tscls,dtm,dtr) * * Check for type change. * it = 0 m0 = mass if((dtr-dtm).le.tol.and.kw.le.9)then * * Check final radius for too large a jump. * aj = MAX(aj,aj*(1.d0-eps)+dtr) mc1 = mc CALL hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars, & r1,lum1,kw,mc1,rc1,menv1,renv1,k21) dr = r1 - rm0 if(ABS(dr).gt.0.1d0*rm0)then dtm = dtr - ajhold*eps dtdr = dtm/ABS(dr) dtm = alpha2*MAX(r1,rm0)*dtdr goto 40 else dtm = dtr goto 50 endif endif * * Limit to a 10% increase in radius assuming no further mass loss * and thus that the pertubation functions due to small envelope mass * will not change the radius. * 40 aj = ajhold + dtm mc1 = mc CALL hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars, & r1,lum1,kw,mc1,rc1,menv1,renv1,k21) dr = r1 - rm0 it = it + 1 if(it.eq.20.and.kw.eq.4) goto 50 if(it.gt.30)then WRITE(99,*)' DANGER2! ',it,kw,mass,dr,rm0 WRITE(*,*)' STOP: EVOLV1 FATAL ERROR ' CALL exit(0) STOP endif if(ABS(dr).gt.0.1d0*rm0)then dtdr = dtm/ABS(dr) dtm = alpha2*MAX(rm0,r1)*dtdr if(it.ge.20) dtm = 0.5d0*dtm goto 40 endif * 50 continue * * Ensure that change of type has not occurred during radius check. * This is rare but may occur for HG stars of ZAMS mass > 50 Msun. * if(kw.ne.kwold)then kw = kwold mass = m0 CALL star(kw,mass,mt,tm,tn,tscls,lums,GB,zpars) endif * * Choose minimum of time-scale and remaining interval (> 100 yrs). * dtm = MAX(dtm,1.0d-07*aj) dtm = MIN(dtm,tsave-tphys) * 10 continue * 90 continue * tphysf = tphys scm(ip+1,1) = -1.0 spp(jp+1,1) = -1.0 if(ip.ge.nv)then WRITE(99,*)' EVOLV1 ARRAY ERROR ',mass WRITE(*,*)' STOP: EVOLV1 ARRAY ERROR ' CALL exit(0) STOP endif * RETURN END *** cc//hrdiag.f *** SUBROUTINE hrdiag(mass,aj,mt,tm,tn,tscls,lums,GB,zpars, & r,lum,kw,mc,rc,menv,renv,k2) * * * H-R diagram for population I stars. * ----------------------------------- * * Computes the new mass, luminosity, radius & stellar type. * Input (MASS, AJ, TM, TN, LUMS & TSCLS) supplied by routine STAR. * Ref: P.P. Eggleton, M.J. Fitchett & C.A. Tout (1989) Ap.J. 347, 998. * * Revised 27th March 1995 by C. A. Tout; * 24th October 1995 to include metallicity; * 14th November 1996 to include naked helium stars; * 28th February 1997 to allow accretion induced supernovae. * * Revised 5th April 1997 by J. R. Hurley * to include Z=0.001 as well as Z=0.02, convective overshooting, * MS hook and more elaborate CHeB * implicit none * integer kw,kwp INTEGER ceflag,tflag,ifflag,nsflag,wdflag COMMON /FLAGS/ ceflag,tflag,ifflag,nsflag,wdflag * CCC added by Long Wang & Peter Berczik 2014.10. real*8 fallback common /fall/fallback real*8 mass,aj,mt,tm,tn,tscls(20),lums(10),GB(10),zpars(20) real*8 r,lum,mc,rc,menv,renv,k2 real*8 mch,mlp,tiny parameter(mch=1.44d0,mlp=12.d0,tiny=1.0d-14) real*8 mass0,mt0,mtc REAL*8 neta,bwind,hewind,mxns COMMON /VALUE1/ neta,bwind,hewind,mxns * real*8 thook,thg,tbagb,tau,tloop,taul,tauh,tau1,tau2,dtau,texp real*8 lx,ly,dell,alpha,beta,eta real*8 rx,ry,delr,rzams,rtms,gamma,rmin,taumin,rg parameter(taumin=5.0d-08) real*8 mcmax,mcx,mcy,mcbagb,lambda real*8 am,xx,fac,rdgen,mew,lum0,kap,zeta,ahe,aco parameter(lum0=7.0d+04,kap=-0.5d0,ahe=4.d0,aco=16.d0) * real*8 thookf,tblf real*8 lalphf,lbetaf,lnetaf,lhookf,lgbtf,lmcgbf,lzhef,lpertf real*8 rzamsf,rtmsf,ralphf,rbetaf,rgammf,rhookf real*8 rgbf,rminf,ragbf,rzahbf,rzhef,rhehgf,rhegbf,rpertf real*8 mctmsf,mcgbtf,mcgbf,mcheif,mcagbf,lzahbf external thookf,tblf external lalphf,lbetaf,lnetaf,lhookf,lgbtf,lmcgbf,lzhef,lpertf external rzamsf,rtmsf,ralphf,rbetaf,rgammf,rhookf external rgbf,rminf,ragbf,rzahbf,rzhef,rhehgf,rhegbf,rpertf external mctmsf,mcgbtf,mcgbf,mcheif,mcagbf,lzahbf * * * --------------------------------------------------------------------- * MASS Stellar mass in solar units (input: old; output: new value). * AJ Current age in Myr. * MT Current mass in solar units (used for R). * TM Main sequence time. * TN Nuclear burning time. * TSCLS Time scale for different stages. * LUMS Characteristic luminosity. * GB Giant Branch parameters * ZPARS Parameters for distinguishing various mass intervals. * R Stellar radius in solar units. * TE Effective temperature (suppressed). * KW Classification type (0 - 15). * MC Core mass. * --------------------------------------------------------------------- * * * Make evolutionary changes to stars that have not reached KW > 5. * CCC added by Long Wang & Peter Berczik 2014.10. fallback = 0.0d0 mass0 = mass if(mass0.gt.100.d0) mass = 100.d0 mt0 = mt if(mt0.gt.100.d0) mt = 100.d0 * if(kw.gt.6) goto 90 * tbagb = tscls(2) + tscls(3) thg = tscls(1) - tm * rzams = rzamsf(mass) rtms = rtmsf(mass) * if(aj.lt.tscls(1))then * * Either on MS or HG * rg = rgbf(mt,lums(3)) * if(aj.lt.tm)then * * Main sequence star. * mc = 0.d0 tau = aj/tm thook = thookf(mass)*tscls(1) zeta = 0.01d0 tau1 = MIN(1.d0,aj/thook) tau2 = MAX(0.d0, & MIN(1.d0,(aj-(1.d0-zeta)*thook)/(zeta*thook))) * dell = lhookf(mass,zpars(1)) dtau = tau1**2 - tau2**2 alpha = lalphf(mass) beta = lbetaf(mass) eta = lnetaf(mass) lx = LOG10(lums(2)/lums(1)) if(tau.gt.taumin)then xx = alpha*tau + beta*tau**eta + & (lx - alpha - beta)*tau**2 - dell*dtau else xx = alpha*tau + (lx - alpha)*tau**2 - dell*dtau endif lum = lums(1)*10.d0**xx * delr = rhookf(mass,zpars(1)) dtau = tau1**3 - tau2**3 alpha = ralphf(mass) beta = rbetaf(mass) gamma = rgammf(mass) rx = LOG10(rtms/rzams) * Note that the use of taumin is a slightly pedantic attempt to * avoid floating point underflow. It IS overkill! if(tau.gt.taumin)then xx = alpha*tau + beta*tau**10 + gamma*tau**40 + & (rx - alpha - beta - gamma)*tau**3 - delr*dtau else xx = alpha*tau + (rx - alpha)*tau**3 - delr*dtau endif r = rzams*10.d0**xx * if(mass.lt.(zpars(1)-0.3d0))then kw = 0 * This following is given by Chris for low mass MS stars which will be * substantially degenerate. We need the Hydrogen abundance, X, which we * calculate from Z assuming that the helium abundance, Y, is calculated * according to Y = 0.24 + 2*Z rdgen = 0.0258d0*((1.d0+zpars(11))**(5.d0/3.d0))* & (mass**(-1.d0/3.d0)) r = MAX(rdgen,r) else kw = 1 endif * else * * Star is on the HG * mcx = mc if(mass.le.zpars(2))then mc = mcgbf(lums(3),GB,lums(6)) elseif(mass.le.zpars(3))then mc = mcheif(mass,zpars(2),zpars(9)) else mc = mcheif(mass,zpars(2),zpars(10)) endif eta = mctmsf(mass) tau = (aj - tm)/thg mc = ((1.d0 - tau)*eta + tau)*mc mc = MAX(mc,mcx) * * Test whether core mass has reached total mass. * if(mc.ge.mt)then aj = 0.d0 if(mass.gt.zpars(2))then * * Zero-age helium star * mc = 0.d0 mass = mt kw = 7 CALL star(kw,mass,mt,tm,tn,tscls,lums,GB,zpars) else * * Zero-age helium white dwarf. * mc = mt mass = mt kw = 10 endif else lum = lums(2)*(lums(3)/lums(2))**tau if(mass.le.zpars(3))then rx = rg else * He-ignition and end of HG occur at Rmin rmin = rminf(mass) ry = ragbf(mt,lums(4),zpars(2)) rx = MIN(rmin,ry) if(mass.le.mlp)then texp = log(mass/mlp)/log(zpars(3)/mlp) rx = rg rx = rmin*(rx/rmin)**texp endif tau2 = tblf(mass,zpars(2),zpars(3)) if(tau2.lt.tiny) rx = ry endif r = rtms*(rx/rtms)**tau kw = 2 endif * endif * * Now the GB, CHeB and AGB evolution. * elseif(aj.lt.tscls(2))then * * Red Giant. * kw = 3 lum = lgbtf(aj,GB(1),GB,tscls(4),tscls(5),tscls(6)) if(mass.le.zpars(2))then * Star has a degenerate He core which grows on the GB mc = mcgbf(lum,GB,lums(6)) else * Star has a non-degenerate He core which may grow, but * only slightly, on the GB tau = (aj - tscls(1))/(tscls(2) - tscls(1)) mcx = mcheif(mass,zpars(2),zpars(9)) mcy = mcheif(mass,zpars(2),zpars(10)) mc = mcx + (mcy - mcx)*tau endif r = rgbf(mt,lum) rg = r if(mc.ge.mt)then aj = 0.d0 if(mass.gt.zpars(2))then * * Zero-age helium star * mc = 0.d0 mass = mt kw = 7 CALL star(kw,mass,mt,tm,tn,tscls,lums,GB,zpars) else * * Zero-age helium white dwarf. * mc = mt mass = mt kw = 10 endif endif * elseif(aj.lt.tbagb)then * * Core helium burning star. * if(kw.eq.3.and.mass.le.zpars(2))then mass = mt CALL star(kw,mass,mt,tm,tn,tscls,lums,GB,zpars) aj = tscls(2) endif if(mass.le.zpars(2))then mcx = mcgbf(lums(4),GB,lums(6)) else mcx = mcheif(mass,zpars(2),zpars(10)) endif tau = (aj - tscls(2))/tscls(3) mc = mcx + (mcagbf(mass) - mcx)*tau * if(mass.le.zpars(2))then lx = lums(5) ly = lums(7) rx = rzahbf(mt,mc,zpars(2)) rg = rgbf(mt,lx) rmin = rg*zpars(13)**(mass/zpars(2)) texp = MIN(MAX(0.4d0,rmin/rx),2.5d0) ry = ragbf(mt,ly,zpars(2)) if(rmin.lt.rx)then taul = (log(rx/rmin))**(1.d0/3.d0) else rmin = rx taul = 0.d0 endif tauh = (log(ry/rmin))**(1.d0/3.d0) tau2 = taul*(tau - 1.d0) + tauh*tau r = rmin*exp(abs(tau2)**3) rg = rg + tau*(ry - rg) lum = lx*(ly/lx)**(tau**texp) elseif(mass.gt.zpars(3))then * * For HM stars He-ignition takes place at Rmin in the HG, and CHeB * consists of a blue phase (before tloop) and a RG phase (after tloop). * tau2 = tblf(mass,zpars(2),zpars(3)) tloop = tscls(2) + tau2*tscls(3) rmin = rminf(mass) rg = rgbf(mt,lums(4)) rx = ragbf(mt,lums(4),zpars(2)) rmin = MIN(rmin, rx) if(mass.le.mlp) then texp = log(mass/mlp)/log(zpars(3)/mlp) rx = rg rx = rmin*(rx/rmin)**texp else rx = rmin end if texp = MIN(MAX(0.4d0,rmin/rx),2.5d0) lum = lums(4)*(lums(7)/lums(4))**(tau**texp) if(aj.lt.tloop)then ly = lums(4)*(lums(7)/lums(4))**(tau2**texp) ry = ragbf(mt,ly,zpars(2)) taul = 0.d0 if(ABS(rmin-rx).gt.tiny)then taul = (log(rx/rmin))**(1.d0/3.d0) endif tauh = 0.d0 if(ry.gt.rmin) tauh = (log(ry/rmin))**(1.d0/3.d0) tau = (aj - tscls(2))/(tau2*tscls(3)) tau2 = taul*(tau - 1.d0) + tauh*tau r = rmin*exp(abs(tau2)**3) rg = rg + tau*(ry - rg) else r = ragbf(mt,lum,zpars(2)) rg = r end if else * * For IM stars CHeB consists of a RG phase (before tloop) and a blue * loop (after tloop). * tau2 = 1.d0 - tblf(mass,zpars(2),zpars(3)) tloop = tscls(2) + tau2*tscls(3) if(aj.lt.tloop)then tau = (tloop - aj)/(tau2*tscls(3)) lum = lums(5)*(lums(4)/lums(5))**(tau**3) r = rgbf(mt,lum) rg = r else lx = lums(5) ly = lums(7) rx = rgbf(mt,lx) rmin = rminf(mt) texp = MIN(MAX(0.4d0,rmin/rx),2.5d0) ry = ragbf(mt,ly,zpars(2)) if(rmin.lt.rx)then taul = (log(rx/rmin))**(1.d0/3.d0) else rmin = rx taul = 0.d0 endif tauh = (log(ry/rmin))**(1.d0/3.d0) tau = (aj - tloop)/(tscls(3) - (tloop - tscls(2))) tau2 = taul*(tau - 1.d0) + tauh*tau r = rmin*exp(abs(tau2)**3) rg = rx + tau*(ry - rx) lum = lx*(ly/lx)**(tau**texp) endif endif * * Test whether core mass exceeds total mass. * if(mc.ge.mt)then * * Evolved MS naked helium star. * kw = 7 xx = (aj - tscls(2))/tscls(3) mass = mt CALL star(kw,mass,mt,tm,tn,tscls,lums,GB,zpars) aj = xx*tm else kw = 4 endif * else * * Asymptotic Red Giant. * * On the AGB the He core mass remains constant until at Ltp it * is caught by the C core mass and they grow together. * mcbagb = mcagbf(mass) mcx = mcgbtf(tbagb,GB(8),GB,tscls(7),tscls(8),tscls(9)) mcmax = MAX(MAX(mch,0.773d0*mcbagb-0.35d0),1.05d0*mcx) * if(aj.lt.tscls(13))then mcx = mcgbtf(aj,GB(8),GB,tscls(7),tscls(8),tscls(9)) mc = mcbagb lum = lmcgbf(mcx,GB) if(mt.le.mc)then * * Evolved naked helium star as the envelope is lost but the * star has not completed its interior burning. The star becomes * a post-HeMS star. * kw = 9 mt = mc mass = mt mc = mcx CALL star(kw,mass,mt,tm,tn,tscls,lums,GB,zpars) if(mc.le.GB(7))then aj = tscls(4) - (1.d0/((GB(5)-1.d0)*GB(8)*GB(4)))* & (mc**(1.d0-GB(5))) else aj = tscls(5) - (1.d0/((GB(6)-1.d0)*GB(8)*GB(3)))* & (mc**(1.d0-GB(6))) endif aj = MAX(aj,tm) goto 90 else kw = 5 endif else kw = 6 mc = mcgbtf(aj,GB(2),GB,tscls(10),tscls(11),tscls(12)) lum = lmcgbf(mc,GB) * * Approximate 3rd Dredge-up on AGB by limiting Mc. * lambda = MIN(0.9d0,0.3d0+0.001d0*mass**5) tau = tscls(13) mcx = mcgbtf(tau,GB(2),GB,tscls(10),tscls(11),tscls(12)) mcy = mc mc = mc - lambda*(mcy-mcx) mcx = mc mcmax = MIN(mt,mcmax) endif r = ragbf(mt,lum,zpars(2)) rg = r * * Mc,x represents the C core mass and we now test whether it * exceeds either the total mass or the maximum allowed core mass. * if(mcmax-mcx.lt.tiny)then aj = 0.d0 mc = mcmax if(mc.lt.mch)then if(ifflag.ge.1)then * * Invoke WD IFMR from HPE, 1995, MNRAS, 272, 800. * if(zpars(14).ge.1.0d-08)then mc = MIN(0.36d0+0.104d0*mass,0.58d0+0.061d0*mass) mc = MAX(0.54d0+0.042d0*mass,mc) if(mass.lt.1.d0) mc = 0.46d0 else mc = MIN(0.29d0+0.178d0*mass,0.65d0+0.062d0*mass) mc = MAX(0.54d0+0.073d0*mass,mc) endif mc = MIN(mch,mc) endif * mt = mc if(mcbagb.lt.1.6d0)then * * Zero-age Carbon/Oxygen White Dwarf * kw = 11 else * * Zero-age Oxygen/Neon White Dwarf * kw = 12 endif mass = mt * else if(mcbagb.lt.1.6d0)then * * Star is not massive enough to ignite C burning. * so no remnant is left after the SN * kw = 15 aj = 0.d0 mt = 0.d0 lum = 1.0d-10 r = 1.0d-10 else if(nsflag.eq.0)then mt = 1.17d0 + 0.09d0*mc elseif(nsflag.ge.1)then * * Use NS/BH mass given by Belczynski et al. 2002, ApJ, 572, 407. * if(mc.lt.2.5d0) then mcx = 0.161767d0*mc + 1.067055d0 else mcx = 0.314154d0*mc + 0.686088d0 endif if(mc.le.5.d0) then mt = mcx elseif(mc.lt.7.6d0) then mt = mcx + (mc - 5.d0)*(mt - mcx)/2.6d0 CCC added by Long Wang & Peter Berczik 2014.10. fallback = (mc - 5.0d0)/2.6d0 endif CCC added by Long Wang & Peter Berczik 2014.10. if(mc.gt.7.6d0) fallback = 1.0d0 endif mc = mt if(mt.le.mxns)then * * Zero-age Neutron star * kw = 13 else * * Zero-age Black hole * kw = 14 endif endif endif endif * endif * 90 continue * if(kw.ge.7.and.kw.le.9)then * * Naked Helium Star * rzams = rzhef(mt) rx = rzams if(aj.lt.tm)then * * Main Sequence * kw = 7 tau = aj/tm am = MAX(0.d0,0.85d0-0.08d0*mass) lum = lums(1)*(1.d0+0.45d0*tau+am*tau**2) am = MAX(0.d0,0.4d0-0.22d0*LOG10(mt)) r = rx*(1.d0+am*(tau-tau**6)) rg = rx * Star has no core mass and hence no memory of its past * which is why we subject mass and mt to mass loss for * this phase. mc = 0.d0 if(mt.lt.zpars(10)) kw = 10 else * * Helium Shell Burning * kw = 8 lum = lgbtf(aj,GB(8),GB,tscls(4),tscls(5),tscls(6)) r = rhehgf(mt,lum,rx,lums(2)) rg = rhegbf(lum) if(r.ge.rg)then kw = 9 r = rg endif mc = mcgbf(lum,GB,lums(6)) mtc = MIN(mt,1.45d0*mt-0.31d0) mcmax = MIN(mtc,MAX(mch,0.773d0*mass-0.35d0)) if(mcmax-mc.lt.tiny)then aj = 0.d0 mc = mcmax if(mc.lt.mch)then if(mass.lt.1.6d0)then * * Zero-age Carbon/Oxygen White Dwarf * mt = MAX(mc,(mc+0.31d0)/1.45d0) kw = 11 else * * Zero-age Oxygen/Neon White Dwarf * mt = mc kw = 12 endif mass = mt else if(mass.lt.1.6d0)then * * Star is not massive enough to ignite C burning. * so no remnant is left after the SN * kw = 15 aj = 0.d0 mt = 0.d0 lum = 1.0d-10 r = 1.0d-10 else if(nsflag.eq.0)then mt = 1.17d0 + 0.09d0*mc elseif(nsflag.ge.1)then if(mc.lt.2.5d0)then mcx = 0.161767d0*mc + 1.067055d0 else mcx = 0.314154d0*mc + 0.686088d0 endif if(mc.le.5.d0)then mt = mcx elseif(mc.lt.7.6d0)then mt = mcx + (mc - 5.d0)*(mt - mcx)/2.6d0 endif endif mc = mt if(mt.le.mxns)then * * Zero-age Neutron star * kw = 13 else * * Zero-age Black hole * kw = 14 endif endif endif endif endif endif * if(kw.ge.10.and.kw.le.12)then * * White dwarf. * mc = mt if(mc.ge.mch)then * * Accretion induced supernova with no remnant * unless WD is ONe in which case we assume a NS * of minimum mass is the remnant. * if(kw.eq.12)then kw = 13 aj = 0.d0 mt = 1.3d0 else kw = 15 aj = 0.d0 mt = 0.d0 lum = 1.0d-10 r = 1.0d-10 endif else * if(kw.eq.10)then xx = ahe else xx = aco endif * if(wdflag.eq.0)then * * Mestel cooling * lum = 635.d0*mt*zpars(14)/(xx*(aj+0.1d0))**1.4d0 * elseif(wdflag.ge.1)then * * modified-Mestel cooling * if(aj.lt.9000.0)then lum = 300.d0*mt*zpars(14)/(xx*(aj+0.1d0))**1.18d0 else fac = (9000.1d0*xx)**5.3d0 lum = 300.d0*fac*mt*zpars(14)/(xx*(aj+0.1d0))**6.48d0 endif * endif * r = 0.0115d0*SQRT(MAX(1.48204d-06,(mch/mt)**(2.d0/3.d0) & - (mt/mch)**(2.d0/3.d0))) r = MIN(0.1d0,r) if(mt.lt.0.0005d0) r = 0.09d0 if(mt.lt.0.000005d0) r = 0.009d0 * endif endif * if(kw.eq.13)then * * Neutron Star. * mc = mt if(mc.gt.mxns)then * * Accretion induced Black Hole? * kw = 14 aj = 0.d0 else lum = 0.02d0*(mt**0.67d0)/(MAX(aj,0.1d0))**2 r = 1.4d-05 endif endif * if(kw.eq.14)then * * Black hole * mc = mt lum = 1.0d-10 r = 4.24d-06*mt endif * * Calculate the core radius and the luminosity and radius of the * remnant that the star will become. * tau = 0.d0 if(kw.le.1.or.kw.eq.7)then rc = 0.d0 elseif(kw.le.3)then if(mass.gt.zpars(2))then lx = lzhef(mc) rx = rzhef(mc) rc = rx else if(wdflag.eq.0)then lx = 635.d0*mc*zpars(14)/((ahe*0.1d0)**1.4d0) elseif(wdflag.ge.1)then lx = 300.d0*mc*zpars(14)/((ahe*0.1d0)**1.18d0) endif rx = 0.0115d0*SQRT(MAX(1.48204d-06, & (mch/mc)**(2.d0/3.d0)-(mc/mch)**(2.d0/3.d0))) rc = 5.d0*rx endif elseif(kw.eq.4)then tau = (aj - tscls(2))/tscls(3) kwp = 7 CALL star(kwp,mc,mc,tm,tn,tscls,lums,GB,zpars) am = MAX(0.d0,0.85d0-0.08d0*mc) lx = lums(1)*(1.d0+0.45d0*tau+am*tau**2) rx = rzhef(mc) am = MAX(0.d0,0.4d0-0.22d0*LOG10(mc)) rx = rx*(1.d0+am*(tau-tau**6)) CALL star(kw,mass,mt,tm,tn,tscls,lums,GB,zpars) rc = rx elseif(kw.eq.5)then kwp = 9 if(tn.gt.tbagb) tau = 3.d0*(aj-tbagb)/(tn-tbagb) CALL star(kwp,mc,mc,tm,tn,tscls,lums,GB,zpars) lx = lmcgbf(mcx,GB) if(tau.lt.1.d0) lx = lums(2)*(lx/lums(2))**tau rx = rzhef(mc) rx = MIN(rhehgf(mc,lx,rx,lums(2)),rhegbf(lx)) CALL star(kw,mass,mt,tm,tn,tscls,lums,GB,zpars) rc = rx elseif(kw.le.9)then if(wdflag.eq.0)then lx = 635.d0*mc*zpars(14)/((aco*0.1d0)**1.4d0) elseif(wdflag.ge.1)then lx = 300.d0*mc*zpars(14)/((aco*0.1d0)**1.18d0) endif rx = 0.0115d0*SQRT(MAX(1.48204d-06, & (mch/mc)**(2.d0/3.d0) - (mc/mch)**(2.d0/3.d0))) rc = 5.d0*rx else rc = r menv = 1.0d-10 renv = 1.0d-10 k2 = 0.21d0 endif * * Perturb the luminosity and radius due to small envelope mass. * if(kw.ge.2.and.kw.le.9.and.kw.ne.7)then mew = ((mt-mc)/mt)*MIN(5.d0,MAX(1.2d0,(lum/lum0)**kap)) if(kw.ge.8) mew = ((mtc-mc)/mtc)*5.d0 if(mew.lt.1.d0)then xx = lpertf(mt,mew) lum = lx*(lum/lx)**xx if(r.le.rx)then xx = 0.d0 else xx = rpertf(mt,mew,r,rx) endif r = rx*(r/rx)**xx endif rc = MIN(rc,r) endif * * Calculate mass and radius of convective envelope, and envelope * gyration radius. * if(kw.lt.10)then CALL mrenv(kw,mass,mt,mc,lum,r,rc,aj,tm,lums(2),lums(3), & lums(4),rzams,rtms,rg,menv,renv,k2) endif * if(mass.gt.99.99d0)then mass = mass0 endif if(mt.gt.99.99d0)then mt = mt0 endif * return end *** cc//kick.f *** SUBROUTINE kick(kw,m1,m1n,m2,ecc,sep,jorb,vs) implicit none * integer kw,k INTEGER idum COMMON /VALUE3/ idum INTEGER idum2,iy,ir(32) COMMON /RAND3/ idum2,iy,ir integer bhflag real*8 m1,m2,m1n,ecc,sep,jorb,ecc2 real*8 pi,twopi,gmrkm,yearsc,rsunkm parameter(yearsc=3.1557d+07,rsunkm=6.96d+05) real*8 mm,em,dif,der,del,r real*8 u1,u2,vk,v(4),s,theta,phi real*8 sphi,cphi,stheta,ctheta,salpha,calpha real*8 vr,vr2,vk2,vn2,hn2 real*8 mu,cmu,vs(3),v1,v2,mx1,mx2 real*8 sigma COMMON /VALUE4/ sigma,bhflag real ran3,xx external ran3 CCC added by Long Wang & Peter Berczik 2014.10. real*8 fallback common /fall/fallback * do k = 1,3 vs(k) = 0.d0 enddo * if(kw.eq.14.and.bhflag.eq.0) goto 95 * pi = ACOS(-1.d0) twopi = 2.d0*pi * Conversion factor to ensure velocities are in km/s using mass and * radius in solar units. gmrkm = 1.906125d+05 * * Find the initial separation by randomly choosing a mean anomaly. if(sep.gt.0.d0.and.ecc.ge.0.d0)then xx = RAN3(idum) mm = xx*twopi em = mm 2 dif = em - ecc*SIN(em) - mm if(ABS(dif/mm).le.1.0d-04) goto 3 der = 1.d0 - ecc*COS(em) del = dif/der em = em - del goto 2 3 continue r = sep*(1.d0 - ecc*COS(em)) * * Find the initial relative velocity vector. salpha = SQRT((sep*sep*(1.d0-ecc*ecc))/(r*(2.d0*sep-r))) calpha = (-1.d0*ecc*SIN(em))/SQRT(1.d0-ecc*ecc*COS(em)*COS(em)) vr2 = gmrkm*(m1+m2)*(2.d0/r - 1.d0/sep) vr = SQRT(vr2) else vr = 0.d0 vr2 = 0.d0 salpha = 0.d0 calpha = 0.d0 endif * * Generate Kick Velocity using Maxwellian Distribution (Phinney 1992). * Use Henon's method for pairwise components (Douglas Heggie 22/5/97). do 20 k = 1,2 u1 = RAN3(idum) u2 = RAN3(idum) * Generate two velocities from polar coordinates S & THETA. s = sigma*SQRT(-2.d0*LOG(1.d0 - u1)) theta = twopi*u2 v(2*k-1) = s*COS(theta) v(2*k) = s*SIN(theta) 20 continue vk2 = v(1)**2 + v(2)**2 + v(3)**2 vk = SQRT(vk2) CCC added by Long Wang & Peter Berczik 2014.10. if(kw.eq.14) then vk = vk*(1.0d0-fallback) vk2 = vk*vk endif if((kw.eq.14.and.bhflag.eq.0).or.kw.lt.0)then vk2 = 0.d0 vk = 0.d0 if(kw.lt.0) kw = 13 endif sphi = -1.d0 + 2.d0*u1 phi = ASIN(sphi) cphi = COS(phi) stheta = SIN(theta) ctheta = COS(theta) * WRITE(66,*)' KICK VK PHI THETA ',vk,phi,theta if(sep.le.0.d0.or.ecc.lt.0.d0) goto 90 * * Determine the magnitude of the new relative velocity. vn2 = vk2+vr2-2.d0*vk*vr*(ctheta*cphi*salpha-stheta*cphi*calpha) * Calculate the new semi-major axis. sep = 2.d0/r - vn2/(gmrkm*(m1n+m2)) sep = 1.d0/sep * if(sep.le.0.d0)then * ecc = 1.1d0 * goto 90 * endif * Determine the magnitude of the cross product of the separation vector * and the new relative velocity. v1 = vk2*sphi*sphi v2 = (vk*ctheta*cphi-vr*salpha)**2 hn2 = r*r*(v1 + v2) * Calculate the new eccentricity. ecc2 = 1.d0 - hn2/(gmrkm*sep*(m1n+m2)) ecc2 = MAX(ecc2,0.d0) ecc = SQRT(ecc2) * Calculate the new orbital angular momentum taking care to convert * hn to units of Rsun^2/yr. jorb = (m1n*m2/(m1n+m2))*SQRT(hn2)*(yearsc/rsunkm) * Determine the angle between the new and old orbital angular * momentum vectors. cmu = (vr*salpha-vk*ctheta*cphi)/SQRT(v1 + v2) mu = ACOS(cmu) * Calculate the components of the velocity of the new centre-of-mass. 90 continue if(ecc.le.1.0)then * Calculate the components of the velocity of the new centre-of-mass. mx1 = vk*m1n/(m1n+m2) mx2 = vr*(m1-m1n)*m2/((m1n+m2)*(m1+m2)) vs(1) = mx1*ctheta*cphi + mx2*salpha vs(2) = mx1*stheta*cphi + mx2*calpha vs(3) = mx1*sphi else * Calculate the relative hyperbolic velocity at infinity (simple method). sep = r/(ecc-1.d0) * cmu = SQRT(ecc-1.d0) * mu = ATAN(cmu) mu = ACOS(1.d0/ecc) vr2 = gmrkm*(m1n+m2)/sep vr = SQRT(vr2) vs(1) = vr*SIN(mu) vs(2) = vr*COS(mu) vs(3) = 0.d0 ecc = MIN(ecc,99.99d0) endif * 95 continue * RETURN END *** cc//mlwind.f *** real*8 FUNCTION mlwind(kw,lum,r,mt,mc,rl,z) implicit none integer kw real*8 lum,r,mt,mc,rl,z real*8 dml,dms,dmt,p0,x,mew,lum0,kap,neta,bwind,hewind,mxns parameter(lum0=7.0d+04,kap=-0.5d0) common /value1/ neta,bwind,hewind,mxns * * Calculate stellar wind mass loss. * * Apply mass loss of Nieuwenhuijzen & de Jager, A&A, 1990, 231, 134, * for massive stars over the entire HRD. dms = 0.d0 if(lum.gt.4000.d0)then x = MIN(1.d0,(lum-4000.d0)/500.d0) dms = 9.6d-15*x*(r**0.81d0)*(lum**1.24d0)*(mt**0.16d0) dms = dms*(z/0.02d0)**(1.d0/2.d0) endif if(kw.ge.2.and.kw.le.9)then * 'Reimers' mass loss dml = neta*4.0d-13*r*lum/mt if(rl.gt.0.d0) dml = dml*(1.d0 + bwind*(MIN(0.5d0,(r/rl)))**6) * Apply mass loss of Vassiliadis & Wood, ApJ, 1993, 413, 641, * for high pulsation periods on AGB. if(kw.eq.5.or.kw.eq.6)then p0 = -2.07d0 - 0.9d0*log10(mt) + 1.94d0*log10(r) p0 = 10.d0**p0 p0 = MIN(p0,2000.d0) dmt = -11.4d0+0.0125d0*(p0-100.d0*MAX(mt-2.5d0,0.d0)) dmt = 10.d0**dmt dmt = 1.d0*MIN(dmt,1.36d-09*lum) dml = MAX(dml,dmt) endif if(kw.gt.6)then dms = MAX(dml,1.0d-13*hewind*lum**(3.d0/2.d0)) else dms = MAX(dml,dms) mew = ((mt-mc)/mt)*MIN(5.d0,MAX(1.2d0,(lum/lum0)**kap)) * reduced WR-like mass loss for small H-envelope mass if(mew.lt.1.d0)then dml = 1.0d-13*lum**(3.d0/2.d0)*(1.d0 - mew) dms = MAX(dml,dms) end if * LBV-like mass loss beyond the Humphreys-Davidson limit. x = 1.0d-5*r*sqrt(lum) if(lum.gt.6.0d+05.and.x.gt.1.d0)then dml = 0.1d0*(x-1.d0)**3*(lum/6.0d+05-1.d0) dms = dms + dml endif endif endif * mlwind = dms * return end *** cc//mrenv.f *** SUBROUTINE mrenv(kw,mass,mt,mc,lum,rad,rc,aj,tm,ltms,lbgb,lhei, & rzams,rtms,rg,menv,renv,k2e) implicit none integer kw real*8 mass,mt,mc,lum,rad,rc,aj,tm real*8 k2e,menv,menvg,menvt,menvz,renv,renvg,renvt,renvz real*8 A,B,C,D,E,F,x,y real*8 k2bgb,k2g,k2z,logm,logmt,lbgb,ltms,lhei,rg,rtms,rzams real*8 teff,tebgb,tetms,tau,tauenv,tautms * * A function to estimate the mass and radius of the convective envelope, * as well as the gyration radius of the envelope. * N.B. Valid only for Z=0.02! * * The following input is needed from HRDIAG: * kw = stellar type * mass = zero-age stellar mass * mt = actual mass * mc = core mass (not really needed, can also be done outside subroutine) * lum = luminosity * rad = radius * rc = core radius (not really needed...) * aj = age * tm = main-sequence lifetime * ltms = luminosity at TMS, lums(2) * lbgb = luminosity at BGB, lums(3) * lhei = luminosity at He ignition, lums(4) * rzams = radius at ZAMS * rtms = radius at TMS * rg = giant branch or Hayashi track radius, approporaite for the type. * For kw=1 or 2 this is radius at BGB, and for kw=4 either GB or * AGB radius at present luminosity. * logm = log10(mass) A = MIN(0.81d0,MAX(0.68d0,0.68d0+0.4d0*logm)) C = MAX(-2.5d0,MIN(-1.5d0,-2.5d0+5.d0*logm)) D = -0.1d0 E = 0.025d0 * * Zero-age and BGB values of k^2. * k2z = MIN(0.21d0,MAX(0.09d0-0.27d0*logm,0.037d0+0.033d0*logm)) if(logm.gt.1.3d0) k2z = k2z - 0.055d0*(logm-1.3d0)**2 k2bgb = MIN(0.15d0,MIN(0.147d0+0.03d0*logm,0.162d0-0.04d0*logm)) * if(kw.ge.3.and.kw.le.6)then * * Envelope k^2 for giant-like stars; this will be modified for non-giant * CHeB stars or small envelope mass below. * Formula is fairly accurate for both FGB and AGB stars if M <= 10, and * gives reasonable values for higher masses. Mass dependence is on actual * rather than ZA mass, expected to work for mass-losing stars (but not * tested!). The slightly complex appearance is to insure continuity at * the BGB, which depends on the ZA mass. * logmt = log10(mt) F = 0.208d0 + 0.125d0*logmt - 0.035d0*logmt**2 B = 1.0d+04*mt**(3.d0/2.d0)/(1.d0+0.1d0*mt**(3.d0/2.d0)) x = ((lum-lbgb)/B)**2 y = (F - 0.033d0*log10(lbgb))/k2bgb - 1.d0 k2g = (F - 0.033d0*log10(lum) + 0.4d0*x)/(1.d0+y*(lbgb/lum)+x) elseif(kw.eq.9)then * * Rough fit for for HeGB stars... * B = 3.0d+04*mt**(3.d0/2.d0) x = (MAX(0.d0,lum/B-0.5d0))**2 k2g = (k2bgb + 0.4d0*x)/(1.d0 + 0.4d0*x) else k2g = k2bgb endif * if(kw.le.2)then menvg = 0.5d0 renvg = 0.65d0 elseif(kw.eq.3.and.lum.lt.3.d0*lbgb)then * * FGB stars still close to the BGB do not yet have a fully developed CE. * x = MIN(3.d0,lhei/lbgb) tau = MAX(0.d0,MIN(1.d0,(x-lum/lbgb)/(x-1.d0))) menvg = 1.d0 - 0.5d0*tau**2 renvg = 1.d0 - 0.35d0*tau**2 else menvg = 1.d0 renvg = 1.d0 endif * if(rad.lt.rg)then * * Stars not on the Hayashi track: MS and HG stars, non-giant CHeB stars, * HeMS and HeHG stars, as well as giants with very small envelope mass. * if(kw.le.6)then * * Envelope k^2 fitted for MS and HG stars. * Again, pretty accurate for M <= 10 but less so for larger masses. * [Note that this represents the whole star on the MS, so there is a * discontinuity in stellar k^2 between MS and HG - okay for stars with a * MS hook but low-mass stars should preferably be continous...] * * For other types of star not on the Hayashi track we use the same fit as * for HG stars, this is not very accurate but has the correct qualitative * behaviour. For CheB stars this is an overestimate because they appear * to have a more centrally concentrated envelope than HG stars. * k2e = (k2z-E)*(rad/rzams)**C + E*(rad/rzams)**D elseif(kw.eq.7)then * Rough fit for naked He MS stars. tau = aj/tm k2e = 0.08d0 - 0.03d0*tau elseif(kw.le.9)then * Rough fit for HeHG stars. k2e = 0.08d0*rzams/rad endif * * tauenv measures proximity to the Hayashi track in terms of Teff. * If tauenv>0 then an appreciable convective envelope is present, and * k^2 needs to be modified. * if(kw.le.2)then teff = sqrt(sqrt(lum)/rad) tebgb = sqrt(sqrt(lbgb)/rg) tauenv = MAX(0.d0,MIN(1.d0,(tebgb/teff-A)/(1.d0-A))) else tauenv = MAX(0.d0,MIN(1.d0,(sqrt(rad/rg)-A)/(1.d0-A))) endif * if(tauenv.gt.0.d0)then menv = menvg*tauenv**5 renv = renvg*tauenv**(5.d0/4.d0) if(kw.le.1)then * Zero-age values for CE mass and radius. x = MAX(0.d0,MIN(1.d0,(0.1d0-logm)/0.55d0)) menvz = 0.18d0*x + 0.82d0*x**5 renvz = 0.4d0*x**(1.d0/4.d0) + 0.6d0*x**10 y = 2.d0 + 8.d0*x * Values for CE mass and radius at start of the HG. tetms = sqrt(sqrt(ltms)/rtms) tautms = MAX(0.d0,MIN(1.d0,(tebgb/tetms-A)/(1.d0-A))) menvt = menvg*tautms**5 renvt = renvg*tautms**(5.d0/4.d0) * Modified expressions during MS evolution. tau = aj/tm if(tautms.gt.0.d0)then menv = menvz + tau**y*menv*(menvt - menvz)/menvt renv = renvz + tau**y*renv*(renvt - renvz)/renvt else menv = 0.d0 renv = 0.d0 endif k2e = k2e + tau**y*tauenv**3*(k2g - k2e) else k2e = k2e + tauenv**3*(k2g - k2e) endif else menv = 0.d0 renv = 0.d0 endif else * * All other stars should be true giants. * menv = menvg renv = renvg k2e = k2g endif * menv = menv*(mt - mc) renv = renv*(rad - rc) menv = MAX(menv,1.0d-10) renv = MAX(renv,1.0d-10) * return end *** cc//ran3.f *** REAL FUNCTION ran3(IDUM) * * Random number generator from Numerical Recipes, Press et al. pg 272. * IMPLICIT NONE INTEGER j,k,im1,im2,imm1,ia1,ia2,iq1,iq2,ir1,ir2,ntab,ndiv PARAMETER(im1=2147483563,im2=2147483399,ia1=40014,ia2=40692) PARAMETER(iq1=53668,iq2=52774,ir1=12211,ir2=3791,ntab=32) INTEGER idum INTEGER idum2,iy,ir(ntab) COMMON /RAND3/ idum2,iy,ir DATA idum2/123456789/, iy/0/, ir/ntab*0/ REAL am * am = 1.0/float(im1) imm1 = im1 - 1 ndiv = 1 + imm1/ntab * if(idum.le.0)then idum = MAX(-idum,1) idum2 = idum do 11 , j = ntab+8,1,-1 k = idum/iq1 idum = ia1*(idum-k*iq1)-k*ir1 if(idum.lt.0) idum = idum + im1 if(j.le.ntab) ir(j) = idum 11 continue iy = ir(1) endif k = idum/iq1 idum = ia1*(idum-k*iq1)-k*ir1 if(idum.lt.0) idum = idum + im1 k = idum2/iq2 idum2 = ia2*(idum2-k*iq2)-k*ir2 if(idum2.lt.0) idum2 = idum2 + im2 j = 1 + iy/ndiv iy = ir(j) - idum2 ir(j) = idum if(iy.lt.1) iy = iy + imm1 ran3 = am*iy * RETURN END *** cc//star.f *** SUBROUTINE star(kw,mass,mt,tm,tn,tscls,lums,GB,zpars) * * * Stellar luminosity & evolution time. * ------------------------------------ * implicit none * integer kw * real*8 mass,mt,tm,tn,tscls(20),lums(10),GB(10),zpars(20) real*8 tgb,tbagb,mch,mcmax,mc1,mc2,mcbagb,dx,am real*8 lambda,tau,mtc,mass0 parameter(mch=1.44d0) * real*8 lzamsf,lzahbf,lzhef real*8 tbgbf,thookf,tHef,themsf,mcgbf,mcagbf,mcheif,mcgbtf real*8 ltmsf,lbgbf,lHeIf,lHef,lbagbf,lmcgbf external lzamsf,lzahbf,lzhef external tbgbf,thookf,tHef,themsf,mcgbf,mcagbf,mcheif,mcgbtf external ltmsf,lbgbf,lHeIf,lHef,lbagbf,lmcgbf * * Computes the characteristic luminosities at different stages (LUMS), * and various timescales (TSCLS). * Ref: P.P. Eggleton, M.J. Fitchett & C.A. Tout (1989) Ap.J. 347, 998. * * Revised 27th March 1995 by C. A. Tout * and 24th October 1995 to include metallicity * and 13th December 1996 to include naked helium stars * * Revised 5th April 1997 by J. R. Hurley * to include Z=0.001 as well as Z=0.02, convective overshooting, * MS hook and more elaborate CHeB. It now also sets the Giant * Branch parameters relevant to the mass of the star. * * ------------------------------------------------------------ * Times: 1; BGB 2; He ignition 3; He burning * 4; Giant t(inf1) 5; Giant t(inf2) 6; Giant t(Mx) * 7; FAGB t(inf1) 8; FAGB t(inf2) 9; FAGB t(Mx) * 10; SAGB t(inf1) 11; SAGB t(inf2) 12; SAGB t(Mx) * 13; TP 14; t(Mcmax) * * LUMS: 1; ZAMS 2; End MS 3; BGB * 4; He ignition 5; He burning 6; L(Mx) * 7; BAGB 8; TP * * GB: 1; effective A(H) 2; A(H,He) 3; B * 4; D 5; p 6; q * 7; Mx 8; A(He) 9; Mc,BGB * * ------------------------------------------------------------ * * mass0 = mass if(mass0.gt.100.d0) mass = 100.d0 * if(kw.ge.7.and.kw.le.9) goto 90 if(kw.ge.10) goto 95 * * MS and BGB times * tscls(1) = tbgbf(mass) tm = MAX(zpars(8),thookf(mass))*tscls(1) * * Zero- and terminal age main sequence luminosity * lums(1) = lzamsf(mass) lums(2) = ltmsf(mass) * * Set the GB parameters * GB(1) = MAX(-4.8d0,MIN(-5.7d0+0.8d0*mass,-4.1d0+0.14d0*mass)) GB(1) = 10.d0**GB(1) GB(2) = 1.27d-05 GB(8) = 8.0d-05 GB(3) = MAX(3.0d+04,500.d0 + 1.75d+04*mass**0.6d0) if(mass.le.2.0)then GB(4) = zpars(6) GB(5) = 6.d0 GB(6) = 3.d0 elseif(mass.lt.2.5)then dx = zpars(6) - (0.975d0*zpars(6) - 0.18d0*2.5d0) GB(4) = zpars(6) - dx*(mass - 2.d0)/(0.5d0) GB(5) = 6.d0 - (mass - 2.d0)/(0.5d0) GB(6) = 3.d0 - (mass - 2.d0)/(0.5d0) else GB(4) = MAX(-1.d0,0.5d0*zpars(6) - 0.06d0*mass) GB(4) = MAX(GB(4),0.975d0*zpars(6) - 0.18d0*mass) GB(5) = 5.d0 GB(6) = 2.d0 endif GB(4) = 10.d0**GB(4) GB(7) = (GB(3)/GB(4))**(1.d0/(GB(5)-GB(6))) * * Change in slope of giant L-Mc relation. lums(6) = GB(4)*GB(7)**GB(5) * * HeI ignition luminosity lums(4) = lHeIf(mass,zpars(2)) lums(7) = lbagbf(mass,zpars(2)) * if(mass.lt.0.1d0.and.kw.le.1)then tscls(2) = 1.1d0*tscls(1) tscls(3) = 0.1d0*tscls(1) lums(3) = lbgbf(mass) goto 96 endif * if(mass.le.zpars(3))then * Base of the giant branch luminosity lums(3) = lbgbf(mass) * Set GB timescales tscls(4) = tscls(1) + (1.d0/((GB(5)-1.d0)*GB(1)*GB(4)))* & ((GB(4)/lums(3))**((GB(5)-1.d0)/GB(5))) tscls(6) = tscls(4) - (tscls(4) - tscls(1))*((lums(3)/lums(6)) & **((GB(5)-1.d0)/GB(5))) tscls(5) = tscls(6) + (1.d0/((GB(6)-1.d0)*GB(1)*GB(3)))* & ((GB(3)/lums(6))**((GB(6)-1.d0)/GB(6))) * Set Helium ignition time if(lums(4).le.lums(6))then tscls(2) = tscls(4) - (1.d0/((GB(5)-1.d0)*GB(1)*GB(4)))* & ((GB(4)/lums(4))**((GB(5)-1.d0)/GB(5))) else tscls(2) = tscls(5) - (1.d0/((GB(6)-1.d0)*GB(1)*GB(3)))* & ((GB(3)/lums(4))**((GB(6)-1.d0)/GB(6))) endif tgb = tscls(2) - tscls(1) if(mass.le.zpars(2))then mc1 = mcgbf(lums(4),GB,lums(6)) mc2 = mcagbf(mass) lums(5) = lzahbf(mass,mc1,zpars(2)) tscls(3) = tHef(mass,mc1,zpars(2)) else lums(5) = lHef(mass)*lums(4) tscls(3) = tHef(mass,1.d0,zpars(2))*tscls(1) endif else * Note that for M>zpars(3) there is no GB as the star goes from * HG -> CHeB -> AGB. So in effect tscls(1) refers to the time of * Helium ignition and not the BGB. tscls(2) = tscls(1) tscls(3) = tHef(mass,1.d0,zpars(2))*tscls(1) * This now represents the luminosity at the end of CHeB, ie. BAGB lums(5) = lums(7) * We set lums(3) to be the luminosity at the end of the HG lums(3) = lums(4) endif * * Set the core mass at the BGB. * if(mass.le.zpars(2))then GB(9) = mcgbf(lums(3),GB,lums(6)) elseif(mass.le.zpars(3))then GB(9) = mcheif(mass,zpars(2),zpars(9)) else GB(9) = mcheif(mass,zpars(2),zpars(10)) endif * * FAGB time parameters * tbagb = tscls(2) + tscls(3) tscls(7) = tbagb + (1.d0/((GB(5)-1.d0)*GB(8)*GB(4)))* & ((GB(4)/lums(7))**((GB(5)-1.d0)/GB(5))) tscls(9) = tscls(7) - (tscls(7) - tbagb)*((lums(7)/lums(6)) & **((GB(5)-1.d0)/GB(5))) tscls(8) = tscls(9) + (1.d0/((GB(6)-1.d0)*GB(8)*GB(3)))* & ((GB(3)/lums(6))**((GB(6)-1.d0)/GB(6))) * * Now to find Ltp and ttp using Mc,He,tp * mcbagb = mcagbf(mass) mc1 = mcbagb if(mc1.ge.0.8d0.and.mc1.lt.2.25d0)then * The star undergoes dredge-up at Ltp causing a decrease in Mc,He mc1 = 0.44d0*mc1 + 0.448d0 endif lums(8) = lmcgbf(mc1,GB) if(mc1.le.GB(7))then tscls(13) = tscls(7) - (1.d0/((GB(5)-1.d0)*GB(8)*GB(4)))* & (mc1**(1.d0-GB(5))) else tscls(13) = tscls(8) - (1.d0/((GB(6)-1.d0)*GB(8)*GB(3)))* & (mc1**(1.d0-GB(6))) endif * * SAGB time parameters * if(mc1.le.GB(7))then tscls(10) = tscls(13) + (1.d0/((GB(5)-1.d0)*GB(2)*GB(4)))* & ((GB(4)/lums(8))**((GB(5)-1.d0)/GB(5))) tscls(12) = tscls(10) - (tscls(10) - tscls(13))* & ((lums(8)/lums(6))**((GB(5)-1.d0)/GB(5))) tscls(11) = tscls(12) + (1.d0/((GB(6)-1.d0)*GB(2)*GB(3)))* & ((GB(3)/lums(6))**((GB(6)-1.d0)/GB(6))) else tscls(10) = tscls(7) tscls(12) = tscls(9) tscls(11) = tscls(13) + (1.d0/((GB(6)-1.d0)*GB(2)*GB(3)))* & ((GB(3)/lums(8))**((GB(6)-1.d0)/GB(6))) endif * * Get an idea of when Mc,C = Mc,C,max on the AGB tau = tscls(2) + tscls(3) mc2 = mcgbtf(tau,GB(8),GB,tscls(7),tscls(8),tscls(9)) mcmax = MAX(MAX(mch,0.773d0*mcbagb - 0.35d0),1.05d0*mc2) * if(mcmax.le.mc1)then if(mcmax.le.GB(7))then tscls(14) = tscls(7) - (1.d0/((GB(5)-1.d0)*GB(8)*GB(4)))* & (mcmax**(1.d0-GB(5))) else tscls(14) = tscls(8) - (1.d0/((GB(6)-1.d0)*GB(8)*GB(3)))* & (mcmax**(1.d0-GB(6))) endif else * Star is on SAGB and we need to increase mcmax if any 3rd * dredge-up has occurred. lambda = MIN(0.9d0,0.3d0+0.001d0*mass**5) mcmax = (mcmax - lambda*mc1)/(1.d0 - lambda) if(mcmax.le.GB(7))then tscls(14) = tscls(10) - (1.d0/((GB(5)-1.d0)*GB(2)*GB(4)))* & (mcmax**(1.d0-GB(5))) else tscls(14) = tscls(11) - (1.d0/((GB(6)-1.d0)*GB(2)*GB(3)))* & (mcmax**(1.d0-GB(6))) endif endif tscls(14) = MAX(tbagb,tscls(14)) if(mass.ge.100.d0)then tn = tscls(2) goto 100 endif * * Calculate the nuclear timescale - the time of exhausting * nuclear fuel without further mass loss. * This means we want to find when Mc = Mt which defines Tn and will * be used in determining the timestep required. Note that after some * stars reach Mc = Mt there will be a Naked Helium Star lifetime * which is also a nuclear burning period but is not included in Tn. * if(ABS(mt-mcbagb).lt.1.0d-14.and.kw.lt.5)then tn = tbagb else * Note that the only occurence of Mc being double-valued is for stars * that have a dredge-up. If Mt = Mc where Mc could be the value taken * from CHeB or from the AGB we need to check the current stellar type. if(mt.gt.mcbagb.or.(mt.ge.mc1.and.kw.gt.4))then if(kw.eq.6)then lambda = MIN(0.9d0,0.3d0+0.001d0*mass**5) mc1 = (mt - lambda*mc1)/(1.d0 - lambda) else mc1 = mt endif if(mc1.le.GB(7))then tn = tscls(10) - (1.d0/((GB(5)-1.d0)*GB(2)*GB(4)))* & (mc1**(1.d0-GB(5))) else tn = tscls(11) - (1.d0/((GB(6)-1.d0)*GB(2)*GB(3)))* & (mc1**(1.d0-GB(6))) endif else if(mass.gt.zpars(3))then mc1 = mcheif(mass,zpars(2),zpars(10)) if(mt.le.mc1)then tn = tscls(2) else tn = tscls(2) + tscls(3)*((mt - mc1)/(mcbagb - mc1)) endif elseif(mass.le.zpars(2))then mc1 = mcgbf(lums(3),GB,lums(6)) mc2 = mcgbf(lums(4),GB,lums(6)) if(mt.le.mc1)then tn = tscls(1) elseif(mt.le.mc2)then if(mt.le.GB(7))then tn = tscls(4) - (1.d0/((GB(5)-1.d0)*GB(1)*GB(4)))* & (mt**(1.d0-GB(5))) else tn = tscls(5) - (1.d0/((GB(6)-1.d0)*GB(1)*GB(3)))* & (mt**(1.d0-GB(6))) endif else tn = tscls(2) + tscls(3)*((mt - mc2)/(mcbagb - mc2)) endif else mc1 = mcheif(mass,zpars(2),zpars(9)) mc2 = mcheif(mass,zpars(2),zpars(10)) if(mt.le.mc1)then tn = tscls(1) elseif(mt.le.mc2)then tn = tscls(1) + tgb*((mt - mc1)/(mc2 - mc1)) else tn = tscls(2) + tscls(3)*((mt - mc2)/(mcbagb - mc2)) endif endif endif endif tn = MIN(tn,tscls(14)) * goto 100 * 90 continue * * Calculate Helium star Main Sequence lifetime. * tm = themsf(mass) tscls(1) = tm * * Zero- and terminal age Helium star main sequence luminosity * lums(1) = lzhef(mass) am = MAX(0.d0,0.85d0-0.08d0*mass) lums(2) = lums(1)*(1.d0+0.45d0+am) * * Set the Helium star GB parameters * GB(8) = 8.0d-05 GB(3) = 4.1d+04 GB(4) = 5.5d+04/(1.d0+0.4d0*mass**4) GB(5) = 5.d0 GB(6) = 3.d0 GB(7) = (GB(3)/GB(4))**(1.d0/(GB(5)-GB(6))) * Change in slope of giant L-Mc relation. lums(6) = GB(4)*GB(7)**GB(5) * *** Set Helium star GB timescales * mc1 = mcgbf(lums(2),GB,lums(6)) tscls(4) = tm + (1.d0/((GB(5)-1.d0)*GB(8)*GB(4)))* & mc1**(1.d0-GB(5)) tscls(6) = tscls(4) - (tscls(4) - tm)*((GB(7)/mc1) & **(1.d0-GB(5))) tscls(5) = tscls(6) + (1.d0/((GB(6)-1.d0)*GB(8)*GB(3)))* & GB(7)**(1.d0-GB(6)) * * Get an idea of when Mc = MIN(Mt,Mc,C,max) on the GB mtc = MIN(mt,1.45d0*mt-0.31d0) if(mtc.le.0.d0) mtc = mt mcmax = MIN(mtc,MAX(mch,0.773d0*mass-0.35d0)) if(mcmax.le.GB(7))then tscls(14) = tscls(4) - (1.d0/((GB(5)-1.d0)*GB(8)*GB(4)))* & (mcmax**(1.d0-GB(5))) else tscls(14) = tscls(5) - (1.d0/((GB(6)-1.d0)*GB(8)*GB(3)))* & (mcmax**(1.d0-GB(6))) endif tscls(14) = MAX(tscls(14),tm) tn = tscls(14) * goto 100 * 95 continue tm = 1.0d+10 tscls(1) = tm 96 continue tn = 1.0d+10 * 100 continue mass = mass0 * return end *** cc//zcnsts.f *** SUBROUTINE zcnsts(z,zpars) * implicit none integer kw * real*8 z,zpars(20) real*8 tm,tn,tscls(20),lums(10),GB(10) real*8 lzs,dlzs,lz,lzd,dum1,m1,m2,rr,rb,mhefl,lhefl,thefl,lx real*8 tbgbf,thef,lbagbf,lheif,lhef,lzahbf real*8 rgbf,ragbf,rminf,mcgbf external tbgbf,thef,lbagbf,lheif,lhef,lzahbf external rgbf,ragbf,rminf,mcgbf * cc// include 'zdata.h' * * Initialization data set for fitting formulae. * ----------------------------------------------------------- * real*8 xz(76),xt(31),xl(72),xr(119),xg(112),xh(99) * * data for Lzams(1->35) and Rzams(36->76) * data xz / & 3.970417d-01, -3.2913574d-01, 3.4776688d-01, 3.7470851d-01, & 9.011915d-02, & 8.527626d+00,-2.441225973d+01, 5.643597107d+01, 3.706152575d+01, & 5.4562406d+00, & 2.5546d-04, -1.23461d-03, -2.3246d-04, 4.5519d-04, & 1.6176d-04, & 5.432889d+00, -8.62157806d+00, 1.344202049d+01, & 1.451584135d+01, 3.39793084d+00, & 5.563579d+00,-1.032345224d+01, 1.944322980d+01, & 1.897361347d+01, 4.16903097d+00, & 7.8866060d-01, -2.90870942d+00, 6.54713531d+00, & 4.05606657d+00, 5.3287322d-01, & 5.86685d-03, -1.704237d-02, 3.872348d-02, 2.570041d-02, & 3.83376d-03, & 1.715359d+00, 6.2246212d-01, -9.2557761d-01, -1.16996966d+00, & -3.0631491d-01, & 6.597788d+00, -4.2450044d-01,-1.213339427d+01,-1.073509484d+01, & -2.51487077d+00, & 1.008855000d+01, -7.11727086d+00,-3.167119479d+01, & -2.424848322d+01,-5.33608972d+00, & 1.012495d+00, 3.2699690d-01, -9.23418d-03, -3.876858d-02, & -4.12750d-03, & 7.490166d-02, 2.410413d-02, 7.233664d-02, 3.040467d-02, & 1.97741d-03, 1.077422d-02, & 3.082234d+00, 9.447205d-01, -2.15200882d+00, -2.49219496d+00, & -6.3848738d-01, & 1.784778d+01, -7.4534569d+00,-4.896066856d+01,-4.005386135d+01, & -9.09331816d+00, & 2.2582d-04, -1.86899d-03, 3.88783d-03, 1.42402d-03,-7.671d-05/ * * data for Tbgb(1->17) and Thook(18->31) * data xt /1.593890d+03, 2.053038d+03, 1.231226d+03, & 2.327785d+02, 2.706708d+03, 1.483131d+03, & 5.772723d+02, 7.411230d+01, 1.466143d+02, & -1.048442d+02,-6.795374d+01,-1.391127d+01, & 4.141960d-02, 4.564888d-02, 2.958542d-02, & 5.571483d-03, 3.426349d-01, & 1.949814d+01, 1.758178d+00,-6.008212d+00, & -4.470533d+00, 4.903830d+00, 5.212154d-02, & 3.166411d-02,-2.750074d-03,-2.271549d-03, & 1.312179d+00,-3.294936d-01, 9.231860d-02, & 2.610989d-02, 8.073972d-01/ * * data for Ltms(1->27), Lalpha(28->43), Lbeta(44->56) and Lhook(57->72) * data xl /1.031538d+00,-2.434480d-01, 7.732821d+00, & 6.460705d+00, 1.374484d+00, 1.043715d+00, & -1.577474d+00,-5.168234d+00,-5.596506d+00, & -1.299394d+00, 7.859573d+02,-8.542048d+00, & -2.642511d+01,-9.585707d+00, 3.858911d+03, & 2.459681d+03,-7.630093d+01,-3.486057d+02, & -4.861703d+01, 2.888720d+02, 2.952979d+02, & 1.850341d+02, 3.797254d+01, 7.196580d+00, & 5.613746d-01, 3.805871d-01, 8.398728d-02, & 2.321400d-01, 1.828075d-03,-2.232007d-02, & -3.378734d-03, 1.163659d-02, 3.427682d-03, & 1.421393d-03,-3.710666d-03, 1.048020d-02, & -1.231921d-02,-1.686860d-02,-4.234354d-03, & 1.555590d+00,-3.223927d-01,-5.197429d-01, & -1.066441d-01, & 3.855707d-01,-6.104166d-01, 5.676742d+00, & 1.060894d+01, 5.284014d+00, 3.579064d-01, & -6.442936d-01, 5.494644d+00, 1.054952d+01, & 5.280991d+00, 9.587587d-01, 8.777464d-01, & 2.017321d-01, & 1.910302d-01, 1.158624d-01, 3.348990d-02, & 2.599706d-03, 3.931056d-01, 7.277637d-02, & -1.366593d-01,-4.508946d-02, 3.267776d-01, & 1.204424d-01, 9.988332d-02, 2.455361d-02, & 5.990212d-01, 5.570264d-02, 6.207626d-02, & 1.777283d-02/ * * data for Rtms(1->40), Ralpha(41->64), Rbeta(65->83), Rgamma(84->103) * and Rhook(104->119) * data xr /2.187715d-01,-2.154437d+00,-3.768678d+00, & -1.975518d+00,-3.021475d-01, 1.466440d+00, & 1.839725d+00, 6.442199d+00, 4.023635d+00, & 6.957529d-01, 2.652091d+01, 8.178458d+01, & 1.156058d+02, 7.633811d+01, 1.950698d+01, & 1.472103d+00,-2.947609d+00,-3.312828d+00, & -9.945065d-01, 3.071048d+00,-5.679941d+00, & -9.745523d+00,-3.594543d+00,-8.672073d-02, & 2.617890d+00, 1.019135d+00,-3.292551d-02, & -7.445123d-02, 1.075567d-02, 1.773287d-02, & 9.610479d-03, 1.732469d-03, 1.476246d+00, & 1.899331d+00, 1.195010d+00, 3.035051d-01, & 5.502535d+00,-6.601663d-02, 9.968707d-02, & 3.599801d-02, & 4.907546d-01,-1.683928d-01,-3.108742d-01, & -7.202918d-02, 4.537070d+00,-4.465455d+00, & -1.612690d+00,-1.623246d+00, 1.796220d+00, & 2.814020d-01, 1.423325d+00, 3.421036d-01, & 2.256216d+00, 3.773400d-01, 1.537867d+00, & 4.396373d-01, 1.564231d-03, 1.653042d-03, & -4.439786d-03,-4.951011d-03,-1.216530d-03, & 5.210157d+00,-4.143695d+00,-2.120870d+00, & 1.071489d+00,-1.164852d-01,-8.623831d-02, & -1.582349d-02, 7.108492d-01, 7.935927d-01, & 3.926983d-01, 3.622146d-02, 3.478514d+00, & -2.585474d-02,-1.512955d-02,-2.833691d-03, & 3.969331d-03, 4.539076d-03, 1.720906d-03, & 1.897857d-04, 9.132108d-01,-1.653695d-01, & 3.636784d-02, & 1.192334d-02, 1.083057d-02, 1.230969d+00, & 1.551656d+00,-1.668868d-01, 5.818123d-01, & -1.105027d+01,-1.668070d+01, 7.615495d-01, & 1.068243d-01,-2.011333d-01,-9.371415d-02, & -1.015564d-01,-2.161264d-01,-5.182516d-02, & -3.868776d-01,-5.457078d-01,-1.463472d-01, & 9.409838d+00, 1.522928d+00, & 7.330122d-01, 5.192827d-01, 2.316416d-01, & 8.346941d-03, 1.172768d+00,-1.209262d-01, & -1.193023d-01,-2.859837d-02, 3.982622d-01, & -2.296279d-01,-2.262539d-01,-5.219837d-02, & 3.571038d+00,-2.223625d-02,-2.611794d-02, & -6.359648d-03/ * * data for Lbgb(1->24), Lbagb(25->44), Rgb(45->66), Ragb(67->100), * Mchei(101->102) and Mcbagb(103->112) * data xg /9.511033d+01, 6.819618d+01,-1.045625d+01, & -1.474939d+01, 3.113458d+01, 1.012033d+01, & -4.650511d+00,-2.463185d+00, 1.413057d+00, & 4.578814d-01,-6.850581d-02,-5.588658d-02, & 3.910862d+01, 5.196646d+01, 2.264970d+01, & 2.873680d+00, 4.597479d+00,-2.855179d-01, & 2.709724d-01, 6.682518d+00, 2.827718d-01, & -7.294429d-02, 4.637345d+00, 9.301992d+00, & 1.626062d+02,-1.168838d+01,-5.498343d+00, & 3.336833d-01,-1.458043d-01,-2.011751d-02, & 7.425137d+01, 1.790236d+01, 3.033910d+01, & 1.018259d+01, 9.268325d+02,-9.739859d+01, & -7.702152d+01,-3.158268d+01, 1.127018d+01, & 1.622158d+00,-1.443664d+00,-9.474699d-01, & 2.474401d+00, 3.892972d-01, & 9.960283d-01, 8.164393d-01, 2.383830d+00, & 2.223436d+00, 8.638115d-01, 1.231572d-01, & 2.561062d-01, 7.072646d-02,-5.444596d-02, & -5.798167d-02,-1.349129d-02, 1.157338d+00, & 1.467883d+00, 4.299661d+00, 3.130500d+00, & 6.992080d-01, 1.640687d-02, 4.022765d-01, & 3.050010d-01, 9.962137d-01, 7.914079d-01, & 1.728098d-01, & 1.125124d+00, 1.306486d+00, 3.622359d+00, & 2.601976d+00, 3.031270d-01,-1.343798d-01, & 3.349489d-01, 4.531269d-03, 1.131793d-01, & 2.300156d-01, 7.632745d-02, 1.467794d+00, & 2.798142d+00, 9.455580d+00, 8.963904d+00, & 3.339719d+00, 4.426929d-01, 4.658512d-01, & 2.597451d-01, 9.048179d-01, 7.394505d-01, & 1.607092d-01, 1.110866d+00, 9.623856d-01, & 2.735487d+00, 2.445602d+00, 8.826352d-01, & 1.140142d-01,-1.584333d-01,-1.728865d-01, & -4.461431d-01,-3.925259d-01,-1.276203d-01, & -1.308728d-02, & 9.796164d-02, 1.350554d+00, & 1.445216d-01,-6.180219d-02, 3.093878d-02, & 1.567090d-02, 1.304129d+00, 1.395919d-01, & 4.142455d-03,-9.732503d-03, 5.114149d-01, & -1.160850d-02/ * * data for Lhei(1->14), Lhe(15->25) Rmin(26->43), The(44->65), * Tbl(66->79), Lzahb(80->87) and Rzahb(88->99) * data xh /2.751631d+03, 3.557098d+02, & -3.820831d-02, 5.872664d-02, 1.5d+01, & 1.071738d+02,-8.970339d+01,-3.949739d+01, & 7.348793d+02,-1.531020d+02,-3.793700d+01, & 9.219293d+00,-2.005865d+00,-5.561309d-01, & 2.917412d+00, 1.575290d+00, 5.751814d-01, & 6.371760d-01, 3.880523d-01, 4.916389d+00, & 2.862149d+00, 7.844850d-01, 3.629118d+00, & -9.112722d-01, 1.042291d+00, & 1.609901d+01, 7.391573d+00, 2.277010d+01, & 8.334227d+00, 1.747500d-01, 6.271202d-02, & -2.324229d-02,-1.844559d-02, 2.752869d+00, & 2.729201d-02, 4.996927d-01, 2.496551d-01, & -9.138012d-02,-3.671407d-01, 3.518506d+00, & 1.112440d+00,-4.556216d-01,-2.179426d-01, & 1.314955d+02, 2.009258d+01,-5.143082d-01, & -1.379140d+00, 1.823973d+01,-3.074559d+00, & -4.307878d+00, 1.5d1, & 2.327037d+00, 2.403445d+00, 1.208407d+00, & 2.087263d-01, 1.079113d-01, 1.762409d-02, & 1.096601d-02, 3.058818d-03, 2.327409d+00, & 6.901582d-01,-2.158431d-01,-1.084117d-01, & 1.997378d+00,-8.126205d-01, & 2.214315d+00,-1.975747d+00, 4.805428d-01, & 2.471620d+00,-5.401682d+00, 3.247361d+00, & 5.072525d+00, 1.146189d+01, 6.961724d+00, & 1.316965d+00, 5.139740d+00, 1.127733d+00, & 2.344416d-01,-3.793726d-01, & 5.496045d+01,-1.289968d+01, 6.385758d+00, & 1.832694d+00,-5.766608d-02, 5.696128d-02, & 1.211104d+02, 1.647903d+00, & 2.063983d+00, 7.363827d-01, 2.654323d-01, & -6.140719d-02, 2.214088d+02, 2.187113d+02, & 1.170177d+01,-2.635340d+01, 2.003160d+00, & 9.388871d-01, 9.656450d-01, 2.362266d-01/ * *** real*8 msp(200),gbp(200),c(5) common /MSCFF/ msp common /GBCFF/ gbp data c /3.040581d-01, 8.049509d-02, 8.967485d-02, & 8.780198d-02, 2.219170d-02/ * * ------------------------------------------------------------ * * zpars: 1; M below which hook doesn't appear on MS, Mhook. * 2; M above which He ignition occurs non-degenerately, Mhef. * 3; M above which He ignition occurs on the HG, Mfgb. * 4; M below which C/O ignition doesn't occur, Mup. * 5; M above which C ignites in the centre, Mec. * 6; value of log D for M<= zpars(3) * 7; value of x for Rgb propto M^(-x) * 8; value of x for tMS = MAX(tHOOK,x*tBGB) * 9; constant for McHeIf when computing Mc,BGB, mchefl. * 10; constant for McHeIf when computing Mc,HeI, mchefl. * 11; hydrogen abundance. * 12; helium abundance. * 13; constant x in rmin = rgb*x**y used by LM CHeB. * 14; z**0.4 to be used for WD L formula. * * ------------------------------------------------------------ * lzs = log10(z/0.02d0) dlzs = 1.d0/(z*log(10.d0)) lz = log10(z) lzd = lzs + 1.d0 * zpars(1) = 1.0185d0 + lzs*(0.16015d0 + lzs*0.0892d0) zpars(2) = 1.995d0 + lzs*(0.25d0 + lzs*0.087d0) zpars(3) = 16.5d0*z**0.06d0/(1.d0 + (1.0d-04/z)**1.27d0) zpars(4) = MAX(6.11044d0 + 1.02167d0*lzs, 5.d0) zpars(5) = zpars(4) + 1.8d0 zpars(6) = 5.37d0 + lzs*0.135d0 zpars(7) = c(1) + lzs*(c(2) + lzs*(c(3) + lzs*(c(4) + lzs*c(5)))) zpars(8) = MAX(0.95d0,MAX(0.95d0-(10.d0/3.d0)*(z-0.01d0), & MIN(0.99d0,0.98d0-(100.d0/7.d0)*(z-0.001d0)))) *** * Lzams msp(1) = xz(1)+lzs*(xz(2)+lzs*(xz(3)+lzs*(xz(4)+lzs*xz(5)))) msp(2) = xz(6)+lzs*(xz(7)+lzs*(xz(8)+lzs*(xz(9)+lzs*xz(10)))) msp(3) = xz(11)+lzs*(xz(12)+lzs*(xz(13)+lzs*(xz(14)+lzs*xz(15)))) msp(4) = xz(16)+lzs*(xz(17)+lzs*(xz(18)+lzs*(xz(19)+lzs*xz(20)))) msp(5) = xz(21)+lzs*(xz(22)+lzs*(xz(23)+lzs*(xz(24)+lzs*xz(25)))) msp(6) = xz(26)+lzs*(xz(27)+lzs*(xz(28)+lzs*(xz(29)+lzs*xz(30)))) msp(7) = xz(31)+lzs*(xz(32)+lzs*(xz(33)+lzs*(xz(34)+lzs*xz(35)))) * Rzams msp(8) = xz(36)+lzs*(xz(37)+lzs*(xz(38)+lzs*(xz(39)+lzs*xz(40)))) msp(9) = xz(41)+lzs*(xz(42)+lzs*(xz(43)+lzs*(xz(44)+lzs*xz(45)))) msp(10) = xz(46)+lzs*(xz(47)+lzs*(xz(48)+lzs*(xz(49)+lzs*xz(50)))) msp(11) = xz(51)+lzs*(xz(52)+lzs*(xz(53)+lzs*(xz(54)+lzs*xz(55)))) msp(12) = xz(56)+lzs*(xz(57)+lzs*(xz(58)+lzs*(xz(59)+lzs*xz(60)))) msp(13) = xz(61) msp(14) = xz(62)+lzs*(xz(63)+lzs*(xz(64)+lzs*(xz(65)+lzs*xz(66)))) msp(15) = xz(67)+lzs*(xz(68)+lzs*(xz(69)+lzs*(xz(70)+lzs*xz(71)))) msp(16) = xz(72)+lzs*(xz(73)+lzs*(xz(74)+lzs*(xz(75)+lzs*xz(76)))) * Tbgb msp(17) = xt(1)+lzs*(xt(2)+lzs*(xt(3)+lzs*xt(4))) msp(18) = xt(5)+lzs*(xt(6)+lzs*(xt(7)+lzs*xt(8))) msp(19) = xt(9)+lzs*(xt(10)+lzs*(xt(11)+lzs*xt(12))) msp(20) = xt(13)+lzs*(xt(14)+lzs*(xt(15)+lzs*xt(16))) msp(21) = xt(17) * dTbgb/dz msp(117) = dlzs*(xt(2)+lzs*(2.d0*xt(3)+3.d0*lzs*xt(4))) msp(118) = dlzs*(xt(6)+lzs*(2.d0*xt(7)+3.d0*lzs*xt(8))) msp(119) = dlzs*(xt(10)+lzs*(2.d0*xt(11)+3.d0*lzs*xt(12))) msp(120) = dlzs*(xt(14)+lzs*(2.d0*xt(15)+3.d0*lzs*xt(16))) * Thook msp(22) = xt(18)+lzs*(xt(19)+lzs*(xt(20)+lzs*xt(21))) msp(23) = xt(22) msp(24) = xt(23)+lzs*(xt(24)+lzs*(xt(25)+lzs*xt(26))) msp(25) = xt(27)+lzs*(xt(28)+lzs*(xt(29)+lzs*xt(30))) msp(26) = xt(31) * Ltms msp(27) = xl(1)+lzs*(xl(2)+lzs*(xl(3)+lzs*(xl(4)+lzs*xl(5)))) msp(28) = xl(6)+lzs*(xl(7)+lzs*(xl(8)+lzs*(xl(9)+lzs*xl(10)))) msp(29) = xl(11)+lzs*(xl(12)+lzs*(xl(13)+lzs*xl(14))) msp(30) = xl(15)+lzs*(xl(16)+lzs*(xl(17)+lzs*(xl(18)+lzs*xl(19)))) msp(27) = msp(27)*msp(30) msp(28) = msp(28)*msp(30) msp(31) = xl(20)+lzs*(xl(21)+lzs*(xl(22)+lzs*xl(23))) msp(32) = xl(24)+lzs*(xl(25)+lzs*(xl(26)+lzs*xl(27))) * Lalpha m2 = 2.d0 msp(33) = xl(28)+lzs*(xl(29)+lzs*(xl(30)+lzs*xl(31))) msp(34) = xl(32)+lzs*(xl(33)+lzs*(xl(34)+lzs*xl(35))) msp(35) = xl(36)+lzs*(xl(37)+lzs*(xl(38)+lzs*xl(39))) msp(36) = xl(40)+lzs*(xl(41)+lzs*(xl(42)+lzs*xl(43))) msp(37) = MAX(0.9d0,1.1064d0+lzs*(0.415d0+0.18d0*lzs)) msp(38) = MAX(1.d0,1.19d0+lzs*(0.377d0+0.176d0*lzs)) if(z.gt.0.01d0)then msp(37) = MIN(msp(37),1.d0) msp(38) = MIN(msp(38),1.1d0) endif msp(39) = MAX(0.145d0,0.0977d0-lzs*(0.231d0+0.0753d0*lzs)) msp(40) = MIN(0.24d0+lzs*(0.18d0+0.595d0*lzs),0.306d0+0.053d0*lzs) msp(41) = MIN(0.33d0+lzs*(0.132d0+0.218d0*lzs), & 0.3625d0+0.062d0*lzs) msp(42) = (msp(33)+msp(34)*m2**msp(36))/ & (m2**0.4d0+msp(35)*m2**1.9d0) * Lbeta msp(43) = xl(44)+lzs*(xl(45)+lzs*(xl(46)+lzs*(xl(47)+lzs*xl(48)))) msp(44) = xl(49)+lzs*(xl(50)+lzs*(xl(51)+lzs*(xl(52)+lzs*xl(53)))) msp(45) = xl(54)+lzs*(xl(55)+lzs*xl(56)) msp(46) = MIN(1.4d0,1.5135d0+0.3769d0*lzs) msp(46) = MAX(0.6355d0-0.4192d0*lzs,MAX(1.25d0,msp(46))) * Lhook msp(47) = xl(57)+lzs*(xl(58)+lzs*(xl(59)+lzs*xl(60))) msp(48) = xl(61)+lzs*(xl(62)+lzs*(xl(63)+lzs*xl(64))) msp(49) = xl(65)+lzs*(xl(66)+lzs*(xl(67)+lzs*xl(68))) msp(50) = xl(69)+lzs*(xl(70)+lzs*(xl(71)+lzs*xl(72))) msp(51) = MIN(1.4d0,1.5135d0+0.3769d0*lzs) msp(51) = MAX(0.6355d0-0.4192d0*lzs,MAX(1.25d0,msp(51))) * Rtms msp(52) = xr(1)+lzs*(xr(2)+lzs*(xr(3)+lzs*(xr(4)+lzs*xr(5)))) msp(53) = xr(6)+lzs*(xr(7)+lzs*(xr(8)+lzs*(xr(9)+lzs*xr(10)))) msp(54) = xr(11)+lzs*(xr(12)+lzs*(xr(13)+lzs*(xr(14)+lzs*xr(15)))) msp(55) = xr(16)+lzs*(xr(17)+lzs*(xr(18)+lzs*xr(19))) msp(56) = xr(20)+lzs*(xr(21)+lzs*(xr(22)+lzs*xr(23))) msp(52) = msp(52)*msp(54) msp(53) = msp(53)*msp(54) msp(57) = xr(24) msp(58) = xr(25)+lzs*(xr(26)+lzs*(xr(27)+lzs*xr(28))) msp(59) = xr(29)+lzs*(xr(30)+lzs*(xr(31)+lzs*xr(32))) msp(60) = xr(33)+lzs*(xr(34)+lzs*(xr(35)+lzs*xr(36))) msp(61) = xr(37)+lzs*(xr(38)+lzs*(xr(39)+lzs*xr(40))) * msp(62) = MAX(0.097d0-0.1072d0*(lz+3.d0),MAX(0.097d0,MIN(0.1461d0, & 0.1461d0+0.1237d0*(lz+2.d0)))) msp(62) = 10.d0**msp(62) m2 = msp(62) + 0.1d0 msp(63) = (msp(52)+msp(53)*msp(62)**msp(55))/ & (msp(54)+msp(62)**msp(56)) msp(64) = (msp(57)*m2**3+msp(58)*m2**msp(61)+ & msp(59)*m2**(msp(61)+1.5d0))/(msp(60)+m2**5) * Ralpha msp(65) = xr(41)+lzs*(xr(42)+lzs*(xr(43)+lzs*xr(44))) msp(66) = xr(45)+lzs*(xr(46)+lzs*(xr(47)+lzs*xr(48))) msp(67) = xr(49)+lzs*(xr(50)+lzs*(xr(51)+lzs*xr(52))) msp(68) = xr(53)+lzs*(xr(54)+lzs*(xr(55)+lzs*xr(56))) msp(69) = xr(57)+lzs*(xr(58)+lzs*(xr(59)+lzs*(xr(60)+lzs*xr(61)))) msp(70) = MAX(0.9d0,MIN(1.d0,1.116d0+0.166d0*lzs)) msp(71) = MAX(1.477d0+0.296d0*lzs,MIN(1.6d0,-0.308d0-1.046d0*lzs)) msp(71) = MAX(0.8d0,MIN(0.8d0-2.d0*lzs,msp(71))) msp(72) = xr(62)+lzs*(xr(63)+lzs*xr(64)) msp(73) = MAX(0.065d0,0.0843d0-lzs*(0.0475d0+0.0352d0*lzs)) msp(74) = 0.0736d0+lzs*(0.0749d0+0.04426d0*lzs) if(z.lt.0.004d0) msp(74) = MIN(0.055d0,msp(74)) msp(75) = MAX(0.091d0,MIN(0.121d0,0.136d0+0.0352d0*lzs)) msp(76) = (msp(65)*msp(71)**msp(67))/(msp(66) + msp(71)**msp(68)) if(msp(70).gt.msp(71))then msp(70) = msp(71) msp(75) = msp(76) endif * Rbeta msp(77) = xr(65)+lzs*(xr(66)+lzs*(xr(67)+lzs*xr(68))) msp(78) = xr(69)+lzs*(xr(70)+lzs*(xr(71)+lzs*xr(72))) msp(79) = xr(73)+lzs*(xr(74)+lzs*(xr(75)+lzs*xr(76))) msp(80) = xr(77)+lzs*(xr(78)+lzs*(xr(79)+lzs*xr(80))) msp(81) = xr(81)+lzs*(xr(82)+lzs*lzs*xr(83)) if(z.gt.0.01d0) msp(81) = MAX(msp(81),0.95d0) msp(82) = MAX(1.4d0,MIN(1.6d0,1.6d0+lzs*(0.764d0+0.3322d0*lzs))) * Rgamma msp(83) = MAX(xr(84)+lzs*(xr(85)+lzs*(xr(86)+lzs*xr(87))), & xr(96)+lzs*(xr(97)+lzs*xr(98))) msp(84) = MIN(0.d0,xr(88)+lzs*(xr(89)+lzs*(xr(90)+lzs*xr(91)))) msp(84) = MAX(msp(84),xr(99)+lzs*(xr(100)+lzs*xr(101))) msp(85) = xr(92)+lzs*(xr(93)+lzs*(xr(94)+lzs*xr(95))) msp(85) = MAX(0.d0,MIN(msp(85),7.454d0+9.046d0*lzs)) msp(86) = MIN(xr(102)+lzs*xr(103),MAX(2.d0,-13.3d0-18.6d0*lzs)) msp(87) = MIN(1.5d0,MAX(0.4d0,2.493d0+1.1475d0*lzs)) msp(88) = MAX(1.d0,MIN(1.27d0,0.8109d0-0.6282d0*lzs)) msp(88) = MAX(msp(88),0.6355d0-0.4192d0*lzs) msp(89) = MAX(5.855420d-02,-0.2711d0-lzs*(0.5756d0+0.0838d0*lzs)) * Rhook msp(90) = xr(104)+lzs*(xr(105)+lzs*(xr(106)+lzs*xr(107))) msp(91) = xr(108)+lzs*(xr(109)+lzs*(xr(110)+lzs*xr(111))) msp(92) = xr(112)+lzs*(xr(113)+lzs*(xr(114)+lzs*xr(115))) msp(93) = xr(116)+lzs*(xr(117)+lzs*(xr(118)+lzs*xr(119))) msp(94) = MIN(1.25d0, & MAX(1.1d0,1.9848d0+lzs*(1.1386d0+0.3564d0*lzs))) msp(95) = 0.063d0 + lzs*(0.0481d0 + 0.00984d0*lzs) msp(96) = MIN(1.3d0,MAX(0.45d0,1.2d0+2.45d0*lzs)) * Lneta if(z.gt.0.0009d0)then msp(97) = 10.d0 else msp(97) = 20.d0 endif * Lbgb gbp(1) = xg(1)+lzs*(xg(2)+lzs*(xg(3)+lzs*xg(4))) gbp(2) = xg(5)+lzs*(xg(6)+lzs*(xg(7)+lzs*xg(8))) gbp(3) = xg(9)+lzs*(xg(10)+lzs*(xg(11)+lzs*xg(12))) gbp(4) = xg(13)+lzs*(xg(14)+lzs*(xg(15)+lzs*xg(16))) gbp(5) = xg(17)+lzs*(xg(18)+lzs*xg(19)) gbp(6) = xg(20)+lzs*(xg(21)+lzs*xg(22)) gbp(3) = gbp(3)**gbp(6) gbp(7) = xg(23) gbp(8) = xg(24) * Lbagb * set gbp(16) = 1.d0 until it is reset later with an initial * call to Lbagbf using mass = zpars(2) and mhefl = 0.0 gbp(9) = xg(25) + lzs*(xg(26) + lzs*xg(27)) gbp(10) = xg(28) + lzs*(xg(29) + lzs*xg(30)) gbp(11) = 15.d0 gbp(12) = xg(31)+lzs*(xg(32)+lzs*(xg(33)+lzs*xg(34))) gbp(13) = xg(35)+lzs*(xg(36)+lzs*(xg(37)+lzs*xg(38))) gbp(14) = xg(39)+lzs*(xg(40)+lzs*(xg(41)+lzs*xg(42))) gbp(15) = xg(43)+lzs*xg(44) gbp(12) = gbp(12)**gbp(15) gbp(14) = gbp(14)**gbp(15) gbp(16) = 1.d0 * Rgb gbp(17) = -4.6739d0-0.9394d0*lz gbp(17) = 10.d0**gbp(17) gbp(17) = MAX(gbp(17),-0.04167d0+55.67d0*z) gbp(17) = MIN(gbp(17),0.4771d0-9329.21d0*z**2.94d0) gbp(18) = MIN(0.54d0,0.397d0+lzs*(0.28826d0+0.5293d0*lzs)) gbp(19) = MAX(-0.1451d0,-2.2794d0-lz*(1.5175d0+0.254d0*lz)) gbp(19) = 10.d0**gbp(19) if(z.gt.0.004d0)then gbp(19) = MAX(gbp(19),0.7307d0+14265.1d0*z**3.395d0) endif gbp(20) = xg(45)+lzs*(xg(46)+lzs*(xg(47)+lzs*(xg(48)+ & lzs*(xg(49)+lzs*xg(50))))) gbp(21) = xg(51)+lzs*(xg(52)+lzs*(xg(53)+lzs*(xg(54)+lzs*xg(55)))) gbp(22) = xg(56)+lzs*(xg(57)+lzs*(xg(58)+lzs*(xg(59)+ & lzs*(xg(60)+lzs*xg(61))))) gbp(23) = xg(62)+lzs*(xg(63)+lzs*(xg(64)+lzs*(xg(65)+lzs*xg(66)))) * Ragb gbp(24) = MIN(0.99164d0-743.123d0*z**2.83d0, & 1.0422d0+lzs*(0.13156d0+0.045d0*lzs)) gbp(25) = xg(67)+lzs*(xg(68)+lzs*(xg(69)+lzs*(xg(70)+ & lzs*(xg(71)+lzs*xg(72))))) gbp(26) = xg(73)+lzs*(xg(74)+lzs*(xg(75)+lzs*(xg(76)+lzs*xg(77)))) gbp(27) = xg(78)+lzs*(xg(79)+lzs*(xg(80)+lzs*(xg(81)+ & lzs*(xg(82)+lzs*xg(83))))) gbp(28) = xg(84)+lzs*(xg(85)+lzs*(xg(86)+lzs*(xg(87)+lzs*xg(88)))) gbp(29) = xg(89)+lzs*(xg(90)+lzs*(xg(91)+lzs*(xg(92)+ & lzs*(xg(93)+lzs*xg(94))))) gbp(30) = xg(95)+lzs*(xg(96)+lzs*(xg(97)+lzs*(xg(98)+ & lzs*(xg(99)+lzs*xg(100))))) m1 = zpars(2) - 0.2d0 gbp(31) = gbp(29) + gbp(30)*m1 gbp(32) = MIN(gbp(25)/zpars(2)**gbp(26),gbp(27)/zpars(2)**gbp(28)) * Mchei gbp(33) = xg(101)**4 gbp(34) = xg(102)*4.d0 * Mcagb gbp(35) = xg(103)+lzs*(xg(104)+lzs*(xg(105)+lzs*xg(106))) gbp(36) = xg(107)+lzs*(xg(108)+lzs*(xg(109)+lzs*xg(110))) gbp(37) = xg(111)+lzs*xg(112) gbp(35) = gbp(35)**4 gbp(36) = gbp(36)*4.d0 gbp(37) = gbp(37)**4 * Lhei * set gbp(41) = -1.d0 until it is reset later with an initial * call to Lheif using mass = zpars(2) and mhefl = 0.0 gbp(38) = xh(1)+lzs*xh(2) gbp(39) = xh(3)+lzs*xh(4) gbp(40) = xh(5) gbp(41) = -1.d0 gbp(42) = xh(6)+lzs*(xh(7)+lzs*xh(8)) gbp(43) = xh(9)+lzs*(xh(10)+lzs*xh(11)) gbp(44) = xh(12)+lzs*(xh(13)+lzs*xh(14)) gbp(42) = gbp(42)**2 gbp(44) = gbp(44)**2 * Lhe gbp(45) = xh(15)+lzs*(xh(16)+lzs*xh(17)) if(lzs.gt.-1.d0)then gbp(46) = 1.d0 - xh(19)*(lzs+1.d0)**xh(18) else gbp(46) = 1.d0 endif gbp(47) = xh(20)+lzs*(xh(21)+lzs*xh(22)) gbp(48) = xh(23)+lzs*(xh(24)+lzs*xh(25)) gbp(45) = gbp(45)**gbp(48) gbp(47) = gbp(47)**gbp(48) gbp(46) = gbp(46)/zpars(3)**0.1d0+(gbp(46)*gbp(47)-gbp(45))/ & zpars(3)**(gbp(48)+0.1d0) * Rmin gbp(49) = xh(26)+lzs*(xh(27)+lzs*(xh(28)+lzs*xh(29))) gbp(50) = xh(30)+lzs*(xh(31)+lzs*(xh(32)+lzs*xh(33))) gbp(51) = xh(34)+lzs*(xh(35)+lzs*(xh(36)+lzs*xh(37))) gbp(52) = 5.d0+xh(38)*z**xh(39) gbp(53) = xh(40)+lzs*(xh(41)+lzs*(xh(42)+lzs*xh(43))) gbp(49) = gbp(49)**gbp(53) gbp(51) = gbp(51)**(2.d0*gbp(53)) * The * set gbp(57) = -1.d0 until it is reset later with an initial * call to Thef using mass = zpars(2), mc = 0.0 and mhefl = 0.0 gbp(54) = xh(44)+lzs*(xh(45)+lzs*(xh(46)+lzs*xh(47))) gbp(55) = xh(48)+lzs*(xh(49)+lzs*xh(50)) gbp(55) = MAX(gbp(55),1.d0) gbp(56) = xh(51) gbp(57) = -1.d0 gbp(58) = xh(52)+lzs*(xh(53)+lzs*(xh(54)+lzs*xh(55))) gbp(59) = xh(56)+lzs*(xh(57)+lzs*(xh(58)+lzs*xh(59))) gbp(60) = xh(60)+lzs*(xh(61)+lzs*(xh(62)+lzs*xh(63))) gbp(61) = xh(64)+lzs*xh(65) gbp(58) = gbp(58)**gbp(61) gbp(60) = gbp(60)**5 * Tbl dum1 = zpars(2)/zpars(3) gbp(62) = xh(66)+lzs*xh(67) gbp(62) = -gbp(62)*log10(dum1) gbp(63) = xh(68) if(lzd.gt.0.d0) then gbp(64) = 1.d0-lzd*(xh(69)+lzd*(xh(70)+lzd*xh(71))) else gbp(64) = 1.d0 end if gbp(65) = 1.d0-gbp(64)*dum1**gbp(63) gbp(66) = 1.d0 - lzd*(xh(77) + lzd*(xh(78) + lzd*xh(79))) gbp(67) = xh(72) + lzs*(xh(73) + lzs*(xh(74) + lzs*xh(75))) gbp(68) = xh(76) * Lzahb gbp(69) = xh(80) + lzs*(xh(81) + lzs*xh(82)) gbp(70) = xh(83) + lzs*(xh(84) + lzs*xh(85)) gbp(71) = 15.d0 gbp(72) = xh(86) gbp(73) = xh(87) * Rzahb gbp(75) = xh(88) + lzs*(xh(89) + lzs*(xh(90) + lzs*xh(91))) gbp(76) = xh(92) + lzs*(xh(93) + lzs*(xh(94) + lzs*xh(95))) gbp(77) = xh(96) + lzs*(xh(97) + lzs*(xh(98) + lzs*xh(99))) *** * finish Lbagb mhefl = 0.d0 lx = lbagbf(zpars(2),mhefl) gbp(16) = lx * finish LHeI dum1 = 0.d0 lhefl = lheif(zpars(2),mhefl) gbp(41) = (gbp(38)*zpars(2)**gbp(39)-lhefl)/ & (EXP(zpars(2)*gbp(40))*lhefl) * finish THe thefl = thef(zpars(2),dum1,mhefl)*tbgbf(zpars(2)) gbp(57) = (thefl-gbp(54))/(gbp(54)*EXP(gbp(56)*zpars(2))) * finish Tblf rb = ragbf(zpars(3),lheif(zpars(3),zpars(2)),mhefl) rr = 1.d0 - rminf(zpars(3))/rb rr = MAX(rr,1.0d-12) gbp(66) = gbp(66)/(zpars(3)**gbp(67)*rr**gbp(68)) * finish Lzahb gbp(74) = lhefl*lHef(zpars(2)) *** kw = 0 tm = 0.d0 tn = 0.d0 CALL star(kw,zpars(2),zpars(2),tm,tn,tscls,lums,GB,zpars) zpars(9) = mcgbf(lums(3),GB,lums(6)) zpars(10) = mcgbf(lums(4),GB,lums(6)) * set the hydrogen and helium abundances zpars(11) = 0.76d0 - 3.d0*z zpars(12) = 0.24d0 + 2.d0*z * set constant for low-mass CHeB stars zpars(13) = rminf(zpars(2))/ & rgbf(zpars(2),lzahbf(zpars(2),zpars(9),zpars(2))) * zpars(14) = z**0.4d0 * return end *** cc//zfuncs.f *** real*8 FUNCTION lzamsf(m) implicit none real*8 m,mx,a(200) common /MSCFF/ a * * A function to evaluate Lzams * ( from Tout et al., 1996, MNRAS, 281, 257 ). * mx = SQRT(m) lzamsf = (a(1)*m**5*mx + a(2)*m**11)/ & (a(3) + m**3 + a(4)*m**5 + a(5)*m**7 + & a(6)*m**8 + a(7)*m**9*mx) * return end *** real*8 FUNCTION rzamsf(m) implicit none real*8 m,mx,a(200) common /MSCFF/ a * * A function to evaluate Rzams * ( from Tout et al., 1996, MNRAS, 281, 257 ). * mx = SQRT(m) rzamsf = ((a(8)*m**2 + a(9)*m**6)*mx + a(10)*m**11 + & (a(11) + a(12)*mx)*m**19)/ & (a(13) + a(14)*m**2 + & (a(15)*m**8 + m**18 + a(16)*m**19)*mx) * return end *** real*8 FUNCTION tbgbf(m) implicit none real*8 m,a(200) common /MSCFF/ a * * A function to evaluate the lifetime to the BGB or to * Helium ignition if no FGB exists. * (JH 24/11/97) * tbgbf = (a(17) + a(18)*m**4 + a(19)*m**(11.d0/2.d0) + m**7)/ & (a(20)*m**2 + a(21)*m**7) * return end *** real*8 FUNCTION tbgbdf(m) implicit none real*8 m,mx,f,df,g,dg,a(200) common /MSCFF/ a * * A function to evaluate the derivitive of the lifetime to the BGB * (or to Helium ignition if no FGB exists) wrt mass. * (JH 24/11/97) * mx = SQRT(m) f = a(17) + a(18)*m**4 + a(19)*m**5*mx + m**7 df = 4.d0*a(18)*m**3 + 5.5d0*a(19)*m**4*mx + 7.d0*m**6 g = a(20)*m**2 + a(21)*m**7 dg = 2.d0*a(20)*m + 7.d0*a(21)*m**6 tbgbdf = (df*g - f*dg)/(g*g) * return end *** real*8 FUNCTION tbgdzf(m) implicit none real*8 m,mx,f,df,g,dg,a(200) common /MSCFF/ a * * A function to evaluate the derivitive of the lifetime to the BGB * (or to Helium ignition if no FGB exists) wrt Z. * (JH 14/12/98) * mx = m**5*SQRT(m) f = a(17) + a(18)*m**4 + a(19)*mx + m**7 df = a(117) + a(118)*m**4 + a(119)*mx g = a(20)*m**2 + a(21)*m**7 dg = a(120)*m**2 tbgdzf = (df*g - f*dg)/(g*g) * return end *** real*8 FUNCTION thookf(m) implicit none real*8 m,a(200) common /MSCFF/ a * * A function to evaluate the lifetime to the end of the MS * hook ( for those models that have one ) as a fraction of * the lifetime to the BGB * Note that this function is only valid for M > Mhook. * (JH 24/11/97) * thookf = 1.d0 - 0.01d0*MAX(a(22)/m**a(23),a(24)+a(25)/m**a(26)) thookf = MAX(thookf,0.5d0) * return end *** real*8 FUNCTION ltmsf(m) implicit none real*8 m,a(200) common /MSCFF/ a * * A function to evaluate the luminosity at the end of the MS * (JH 24/11/97) * ltmsf = (a(27)*m**3 + a(28)*m**4 + a(29)*m**(a(32)+1.8d0))/ & (a(30) + a(31)*m**5 + m**a(32)) * return end *** real*8 FUNCTION lalphf(m) implicit none real*8 m,mcut,a(200) common /MSCFF/ a * * A function to evaluate the Luminosity alpha coefficent. * (JH 24/11/97) * mcut = 2.d0 if(m.ge.mcut)then lalphf = (a(33) + a(34)*m**a(36))/(m**0.4d0 + a(35)*m**1.9d0) else if(m.le.0.5d0)then lalphf = a(39) elseif(m.le.0.7d0)then lalphf = a(39) + ((0.3d0 - a(39))/0.2d0)*(m - 0.5d0) elseif(m.le.a(37))then lalphf = 0.3d0 + ((a(40)-0.3d0)/(a(37)-0.7d0))*(m - 0.7d0) elseif(m.le.a(38))then lalphf = a(40) + ((a(41)-a(40))/(a(38)-a(37)))*(m - a(37)) else lalphf = a(41) + ((a(42)-a(41))/(mcut-a(38)))*(m - a(38)) endif endif * return end *** real*8 FUNCTION lbetaf(m) implicit none real*8 m,a1,a(200) common /MSCFF/ a * * A function to evaluate the Luminosity beta coefficent. * (JH 24/11/97) * lbetaf = a(43) - a(44)*m**a(45) lbetaf = MAX(lbetaf,0.d0) if(m.gt.a(46).and.lbetaf.gt.0.d0)then a1 = a(43) - a(44)*a(46)**a(45) lbetaf = a1 - 10.d0*a1*(m - a(46)) lbetaf = MAX(lbetaf,0.d0) endif * return end *** real*8 FUNCTION lnetaf(m) implicit none real*8 m,a(200) common /MSCFF/ a * * A function to evaluate the Luminosity neta exponent. * (JH 24/11/97) * if(m.le.1.d0)then lnetaf = 10.d0 elseif(m.ge.1.1d0)then lnetaf = 20.d0 else lnetaf = 10.d0 + 100.d0*(m - 1.d0) endif lnetaf = MIN(lnetaf,a(97)) * return end *** real*8 FUNCTION lhookf(m,mhook) implicit none real*8 m,mhook,a2,a(200) common /MSCFF/ a * * A function to evalute the luminosity at the start of * the MS hook ( for those stars that have one ). * Note that this function is only valid for M > Mhook. * (JH 24/11/97) * if(m.le.mhook)then lhookf = 0.d0 elseif(m.ge.a(51))then lhookf = MIN(a(47)/m**a(48),a(49)/m**a(50)) else a2 = MIN(a(47)/a(51)**a(48),a(49)/a(51)**a(50)) lhookf = a2*((m-mhook)/(a(51)-mhook))**0.4d0 endif * return end *** real*8 FUNCTION rtmsf(m) implicit none real*8 m,m2,rchk,a(200) common /MSCFF/ a real*8 rzamsf external rzamsf * * A function to evaluate the radius at the end of the MS * Note that a safety check is added to ensure Rtms > Rzams * when extrapolating the function to low masses. * (JH 24/11/97) * m2 = a(62) + 0.1d0 if(m.le.a(62))then rchk = 1.5d0*rzamsf(m) rtmsf = MAX(rchk,(a(52) + a(53)*m**a(55))/(a(54) + m**a(56))) elseif(m.ge.m2)then rtmsf = (a(57)*m**3+a(58)*m**a(61)+a(59)*m**(a(61)+1.5d0))/ & (a(60) + m**5) else rtmsf = a(63) + ((a(64) - a(63))/0.1d0)*(m - a(62)) endif * return end *** real*8 FUNCTION ralphf(m) implicit none real*8 m,a5,a(200) common /MSCFF/ a * * A function to evaluate the radius alpha coefficent. * (JH 24/11/97) * if(m.le.0.5d0)then ralphf = a(73) elseif(m.le.0.65d0)then ralphf = a(73) + ((a(74) - a(73))/0.15d0)*(m - 0.5d0) elseif(m.le.a(70))then ralphf = a(74) + ((a(75)-a(74))/(a(70)-0.65d0))*(m - 0.65d0) elseif(m.le.a(71))then ralphf = a(75) + ((a(76) - a(75))/(a(71) - a(70)))*(m - a(70)) elseif(m.le.a(72))then ralphf = (a(65)*m**a(67))/(a(66) + m**a(68)) else a5 = (a(65)*a(72)**a(67))/(a(66) + a(72)**a(68)) ralphf = a5 + a(69)*(m - a(72)) endif * return end *** real*8 FUNCTION rbetaf(m) implicit none real*8 m,m2,m3,b2,b3,a(200) common /MSCFF/ a * * A function to evaluate the radius beta coefficent. * (JH 24/11/97) * m2 = 2.d0 m3 = 16.d0 if(m.le.1.d0)then rbetaf = 1.06d0 elseif(m.le.a(82))then rbetaf = 1.06d0 + ((a(81)-1.06d0)/(a(82)-1.d0))*(m-1.d0) elseif(m.le.m2)then b2 = (a(77)*m2**(7.d0/2.d0))/(a(78) + m2**a(79)) rbetaf = a(81) + ((b2-a(81))/(m2-a(82)))*(m-a(82)) elseif(m.le.m3)then rbetaf = (a(77)*m**(7.d0/2.d0))/(a(78) + m**a(79)) else b3 = (a(77)*m3**(7.d0/2.d0))/(a(78) + m3**a(79)) rbetaf = b3 + a(80)*(m - m3) endif rbetaf = rbetaf - 1.d0 * return end *** real*8 FUNCTION rgammf(m) implicit none real*8 m,m1,b1,a(200) common /MSCFF/ a * * A function to evaluate the radius gamma coefficent. * (JH 24/11/97) * m1 = 1.d0 if(m.gt.(a(88)+0.1d0))then rgammf = 0.d0 else b1 = MAX(0.d0,a(83) + a(84)*(m1-a(85))**a(86)) if(m.le.m1)then rgammf = a(83) + a(84)*ABS(m-a(85))**a(86) elseif(m.le.a(88))then rgammf = b1 + (a(89) - b1)*((m - m1)/(a(88) - m1))**a(87) else if(a(88).gt.m1) b1 = a(89) rgammf = b1 - 10.d0*b1*(m - a(88)) endif rgammf = MAX(rgammf,0.d0) endif * return end *** real*8 FUNCTION rhookf(m,mhook) implicit none real*8 m,mhook,m2,b2,a(200) common /MSCFF/ a * * A function to evalute the radius at the start of * the MS hook ( for those stars that have one ). * Note that this function is only valid for M > Mhook. * (JH 24/11/97) * if(m.le.mhook)then rhookf = 0.d0 elseif(m.le.a(94))then rhookf = a(95)*SQRT((m-mhook)/(a(94)-mhook)) elseif(m.le.2.d0)then m2 = 2.d0 b2 = (a(90) + a(91)*m2**(7.d0/2.d0))/ & (a(92)*m2**3 + m2**a(93)) - 1.d0 rhookf = a(95) + (b2-a(95))*((m-a(94))/(m2-a(94)))**a(96) else rhookf = (a(90) + a(91)*m**(7.d0/2.d0))/ & (a(92)*m**3 + m**a(93)) - 1.d0 endif * return end *** real*8 FUNCTION lbgbf(m) real*8 m,a(200) common /GBCFF/ a * * A function to evaluate the luminosity at the end of the * FGB ( for those models that have one ) * Note that this function is only valid for LM & IM stars * (JH 24/11/97) * lbgbf = (a(1)*m**a(5) + a(2)*m**a(8))/ & (a(3) + a(4)*m**a(7) + m**a(6)) * return end *** real*8 FUNCTION lbgbdf(m) real*8 m,a(200) real*8 f,df,g,dg common /GBCFF/ a * * A function to evaluate the derivitive of the Lbgb function. * Note that this function is only valid for LM & IM stars * (JH 24/11/97) * f = a(1)*m**a(5) + a(2)*m**a(8) df = a(5)*a(1)*m**(a(5)-1.d0) + a(8)*a(2)*m**(a(8)-1.d0) g = a(3) + a(4)*m**a(7) + m**a(6) dg = a(7)*a(4)*m**(a(7)-1.d0) + a(6)*m**(a(6)-1.d0) * lbgbdf = (df*g - f*dg)/(g*g) * return end *** real*8 FUNCTION lbagbf(m,mhefl) implicit none real*8 m,mhefl,a4,a(200) common /GBCFF/ a * * A function to evaluate the BAGB luminosity. (OP 21/04/98) * Continuity between LM and IM functions is ensured by setting * gbp(16) = lbagbf(mhefl,0.0) with gbp(16) = 1.0. * a4 = (a(9)*mhefl**a(10) - a(16))/(exp(mhefl*a(11))*a(16)) * if(m.lt.mhefl)then lbagbf = a(9)*m**a(10)/(1.d0 + a4*exp(m*a(11))) else lbagbf = (a(12) + a(13)*m**(a(15)+1.8d0))/(a(14) + m**a(15)) endif * return end *** real*8 FUNCTION rgbf(m,lum) implicit none real*8 m,lum,a1,a(200) common /GBCFF/ a * * A function to evaluate radius on the GB. * (JH 24/11/97) * a1 = MIN(a(20)/m**a(21),a(22)/m**a(23)) rgbf = a1*(lum**a(18) + a(17)*lum**a(19)) * return end *** real*8 FUNCTION rgbdf(m,lum) implicit none real*8 m,lum,a1,a(200) common /GBCFF/ a * * A function to evaluate radius derivitive on the GB (as f(L)). * (JH 24/11/97) * a1 = MIN(a(20)/m**a(21),a(22)/m**a(23)) rgbdf = a1*(a(18)*lum**(a(18)-1.d0) + & a(17)*a(19)*lum**(a(19)-1.d0)) * return end *** real*8 FUNCTION ragbf(m,lum,mhelf) implicit none real*8 m,lum,mhelf,m1,a1,a4,xx,a(200) common /GBCFF/ a * * A function to evaluate radius on the AGB. * (JH 24/11/97) * m1 = mhelf - 0.2d0 if(m.ge.mhelf)then xx = a(24) elseif(m.ge.m1)then xx = 1.d0 + 5.d0*(a(24)-1.d0)*(m-m1) else xx = 1.d0 endif a4 = xx*a(19) if(m.le.m1)then a1 = a(29) + a(30)*m elseif(m.ge.mhelf)then a1 = MIN(a(25)/m**a(26),a(27)/m**a(28)) else a1 = a(31) + 5.d0*(a(32)-a(31))*(m-m1) endif * ragbf = a1*(lum**a(18) + a(17)*lum**a4) * return end *** real*8 FUNCTION ragbdf(m,lum,mhelf) implicit none real*8 m,lum,mhelf,m1,a1,a4,xx,a(200) common /GBCFF/ a * * A function to evaluate radius derivitive on the AGB (as f(L)). * (JH 24/11/97) * m1 = mhelf - 0.2d0 if(m.ge.mhelf)then xx = a(24) elseif(m.ge.m1)then xx = 1.d0 + 5.d0*(a(24)-1.d0)*(m-m1) else xx = 1.d0 endif a4 = xx*a(19) if(m.le.m1)then a1 = a(29) + a(30)*m elseif(m.ge.mhelf)then a1 = MIN(a(25)/m**a(26),a(27)/m**a(28)) else a1 = a(31) + 5.d0*(a(32)-a(31))*(m-m1) endif * ragbdf = a1*(a(18)*lum**(a(18)-1.d0) + & a(17)*a4*lum**(a4-1.d0)) * return end *** real*8 FUNCTION mctmsf(m) implicit none real*8 m,m525 * * A function to evaluate core mass at the end of the MS as a * fraction of the BGB value, i.e. this must be multiplied by * the BGB value (see below) to give the actual core mass (JH 5/9/99) * m525 = m**(21.d0/4.d0) mctmsf = (1.586d0 + m525)/(2.434d0 + 1.02d0*m525) * return end *** real*8 FUNCTION mcheif(m,mhefl,mchefl) implicit none real*8 m,mhefl,mchefl,mcbagb,a3,a(200) common /GBCFF/ a real*8 mcagbf external mcagbf * * A function to evaluate core mass at BGB or He ignition * (depending on mchefl) for IM & HM stars (OP 25/11/97) * mcbagb = mcagbf(m) a3 = mchefl**4 - a(33)*mhefl**a(34) mcheif = MIN(0.95d0*mcbagb,(a3 + a(33)*m**a(34))**(1.d0/4.d0)) * return end *** real*8 FUNCTION mheif(mc,mhefl,mchefl) implicit none real*8 mc,mhefl,mchefl,m1,m2,a3,a(200) common /GBCFF/ a real*8 mbagbf external mbagbf * * A function to evaluate mass at BGB or He ignition * (depending on mchefl) for IM & HM stars by inverting * mcheif * m1 = mbagbf(mc/0.95d0) a3 = mchefl**4 - a(33)*mhefl**a(34) m2 = ((mc**4 - a3)/a(33))**(1.d0/a(34)) mheif = MAX(m1,m2) * return end *** real*8 FUNCTION mcagbf(m) implicit none real*8 m,a(200) common /GBCFF/ a * * A function to evaluate core mass at the BAGB (OP 25/11/97) * mcagbf = (a(37) + a(35)*m**a(36))**(1.d0/4.d0) * return end *** real*8 FUNCTION mbagbf(mc) implicit none real*8 mc,mc4,a(200) common /GBCFF/ a * * A function to evaluate mass at the BAGB by inverting mcagbf. * mc4 = mc**4 if(mc4.gt.a(37))then mbagbf = ((mc4 - a(37))/a(35))**(1.d0/a(36)) else mbagbf = 0.d0 endif * return end *** real*8 FUNCTION mcgbtf(t,A,GB,tinf1,tinf2,tx) implicit none real*8 t,A,GB(10),tinf1,tinf2,tx * * A function to evaluate Mc given t for GB, AGB and NHe stars * if(t.le.tx)then mcgbtf = ((GB(5)-1.d0)*A*GB(4)*(tinf1 - t))** & (1.d0/(1.d0-GB(5))) else mcgbtf = ((GB(6)-1.d0)*A*GB(3)*(tinf2 - t))** & (1.d0/(1.d0-GB(6))) endif * return end *** real*8 FUNCTION lgbtf(t,A,GB,tinf1,tinf2,tx) implicit none real*8 t,A,GB(10),tinf1,tinf2,tx * * A function to evaluate L given t for GB, AGB and NHe stars * if(t.le.tx)then lgbtf = GB(4)*(((GB(5)-1.d0)*A*GB(4)*(tinf1 - t))** & (GB(5)/(1.d0-GB(5)))) else lgbtf = GB(3)*(((GB(6)-1.d0)*A*GB(3)*(tinf2 - t))** & (GB(6)/(1.d0-GB(6)))) endif * return end *** real*8 FUNCTION mcgbf(lum,GB,lx) implicit none real*8 lum,GB(10),lx * * A function to evaluate Mc given L for GB, AGB and NHe stars * if(lum.le.lx)then mcgbf = (lum/GB(4))**(1.d0/GB(5)) else mcgbf = (lum/GB(3))**(1.d0/GB(6)) endif * return end *** real*8 FUNCTION lmcgbf(mc,GB) implicit none real*8 mc,GB(10) * * A function to evaluate L given Mc for GB, AGB and NHe stars * if(mc.le.GB(7))then lmcgbf = GB(4)*(mc**GB(5)) else lmcgbf = GB(3)*(mc**GB(6)) endif * return end *** real*8 FUNCTION lHeIf(m,mhefl) implicit none real*8 m,mhefl,a(200) common /GBCFF/ a * * A function to evaluate He-ignition luminosity (OP 24/11/97) * Continuity between the LM and IM functions is ensured with a first * call setting lhefl = lHeIf(mhefl,0.0) * if(m.lt.mhefl)then lHeIf = a(38)*m**a(39)/(1.d0 + a(41)*EXP(m*a(40))) else lHeIf = (a(42) + a(43)*m**3.8d0)/(a(44) + m**2) endif * return end *** real*8 FUNCTION lHef(m) implicit none real*8 m,a(200) common /GBCFF/ a * * A function to evaluate the ratio LHe,min/LHeI (OP 20/11/97) * Note that this function is everywhere <= 1, and is only valid * for IM stars * lHef = (a(45) + a(46)*m**(a(48)+0.1d0))/(a(47) + m**a(48)) * return end *** real*8 FUNCTION rminf(m) implicit none real*8 m,mx,a(200) common /GBCFF/ a * * A function to evaluate the minimum radius during He-burning * for IM & HM stars (OP 20/11/97) * mx = m**a(53) rminf = (a(49)*m + (a(50)*m)**a(52)*mx)/(a(51) + mx) * return end *** real*8 FUNCTION tHef(m,mc,mhefl) implicit none real*8 m,mc,mhefl,mm,a(200) common /GBCFF/ a real*8 themsf external themsf * * A function to evaluate the He-burning lifetime. (OP 26/11/97) * For IM & HM stars, tHef is relative to tBGB. * Continuity between LM and IM stars is ensured by setting * thefl = tHef(mhefl,0.0,,0.0), and the call to themsf ensures * continuity between HB and NHe stars as Menv -> 0. * if(m.le.mhefl)then mm = MAX((mhefl - m)/(mhefl - mc),1.0d-12) tHef = (a(54) + (themsf(mc) - a(54))*mm**a(55))* & (1.d0 + a(57)*EXP(m*a(56))) else mm = m**5 tHef = (a(58)*m**a(61) + a(59)*mm)/(a(60) + mm) endif * return end *** real*8 FUNCTION tblf(m,mhefl,mfgb) implicit none real*8 m,mhefl,mfgb,mr,m1,m2,r1,a(200) common /GBCFF/ a real*8 lheif,rminf,ragbf external lheif,rminf,ragbf * * A function to evaluate the blue-loop fraction of the He-burning * lifetime for IM & HM stars (OP 28/01/98) * mr = mhefl/mfgb if(m.le.mfgb) then m1 = m/mfgb m2 = log10(m1)/log10(mr) m2 = max(m2,1.0d-12) tblf = a(64)*m1**a(63) + a(65)*m2**a(62) else r1 = 1.d0 - rminf(m)/ragbf(m,lheif(m,mhefl),mhefl) r1 = max(r1,1.0d-12) tblf = a(66)*m**a(67)*r1**a(68) end if tblf = MIN(1.d0,MAX(0.d0,tblf)) if(tblf.lt.1.0d-10) tblf = 0.d0 * return end *** real*8 FUNCTION lzahbf(m,mc,mhefl) implicit none real*8 m,mc,mhefl,mm,a4,a5,a(200) common /GBCFF/ a real*8 lzhef external lzhef * * A function to evaluate the ZAHB luminosity for LM stars. (OP 28/01/98) * Continuity with LHe,min for IM stars is ensured by setting * lx = lHeif(mhefl,z,0.0,1.0)*lHef(mhefl,z,mfgb), and the call to lzhef * ensures continuity between the ZAHB and the NHe-ZAMS as Menv -> 0. * a5 = lzhef(mc) a4 = (a(69) + a5 - a(74))/((a(74) - a5)*exp(a(71)*mhefl)) mm = MAX((m-mc)/(mhefl - mc),1.0d-12) lzahbf = a5 + (1.d0 + a(72))*a(69)*mm**a(70)/ & ((1.d0 + a(72)*mm**a(73))*(1.d0 + a4*EXP(m*a(71)))) * return end *** real*8 FUNCTION rzahbf(m,mc,mhefl) implicit none real*8 m,mc,mhefl,rx,ry,mm,f,a(200) common /GBCFF/ a real*8 rzhef,rgbf,lzahbf * * A function to evaluate the ZAHB radius for LM stars. (OP 28/01/98) * Continuity with R(LHe,min) for IM stars is ensured by setting * lx = lHeif(mhefl,z,0.0,1.0)*lHef(mhefl,z,mfgb), and the call to rzhef * ensures continuity between the ZAHB and the NHe-ZAMS as Menv -> 0. * rx = rzhef(mc) ry = rgbf(m,lzahbf(m,mc,mhefl)) mm = MAX((m-mc)/(mhefl - mc),1.0d-12) f = (1.d0 + a(76))*mm**a(75)/(1.d0 + a(76)*mm**a(77)) rzahbf = (1.d0 - f)*rx + f*ry * return end *** real*8 FUNCTION lzhef(m) implicit none real*8 m,m15 * * A function to evaluate Naked Helium star 'ZAMS' luminosity * m15 = m*SQRT(m) lzhef = 1.5262d+04*m**(41.d0/4.d0)/ & (0.0469d0 + m**6*(31.18d0 + m15*(29.54d0 + m15))) * return end *** real*8 FUNCTION rzhef(m) implicit none real*8 m * * A function to evaluate Helium star 'ZAMS' radius * rzhef = 0.2391d0*m**4.6d0/(0.0065d0 + (0.162d0 + m)*m**3) * return end *** real*8 FUNCTION themsf(m) implicit none real*8 m * * A function to evaluate Helium star main sequence lifetime * themsf = (0.4129d0 + 18.81d0*m**4 + 1.853d0*m**6)/m**(13.d0/2.d0) * return end *** real*8 FUNCTION rhehgf(m,lum,rx,lx) implicit none real*8 m,lum,rx,lx,cm * * A function to evaluate Helium star radius on the Hertzsprung gap * from its mass and luminosity. * cm = 2.0d-03*m**(5.d0/2.d0)/(2.d0 + m**5) rhehgf = rx*(lum/lx)**0.2d0 + 0.02d0*(EXP(cm*lum) - EXP(cm*lx)) * return end *** real*8 FUNCTION rhegbf(lum) implicit none real*8 lum * * A function to evaluate Helium star radius on the giant branch. * rhegbf = 0.08d0*lum**(3.d0/4.d0) * return end *** real*8 FUNCTION lpertf(m,mew) implicit none real*8 m,mew real*8 b,c * * A function to obtain the exponent that perturbs luminosity. * b = 0.002d0*MAX(1.d0,2.5d0/m) c = 3.d0 lpertf = ((1.d0 + b**c)*((mew/b)**c))/(1.d0+(mew/b)**c) * return end *** real*8 FUNCTION rpertf(m,mew,r,rc) implicit none real*8 m,mew,r,rc real*8 a,b,c,q,fac,facmax * * A function to obtain the exponent that perturbs radius. * if(mew.le.0.d0)then rpertf = 0.d0 else a = 0.1d0 b = 0.006d0*MAX(1.d0,2.5d0/m) c = 3.d0 q = log(r/rc) fac = a/q facmax = -14.d0/log10(mew) fac = MIN(fac,facmax) rpertf = ((1.d0 + b**c)*((mew/b)**c)*(mew**fac))/ & (1.d0+(mew/b)**c) endif * return end *** real*8 FUNCTION vrotf(m) implicit none real*8 m * vrotf = 330.d0*m**3.3d0/(15.d0 + m**3.45d0) * return end *** real*8 FUNCTION celamf(kw,m,lum,rad,rzams,menv,fac) implicit none integer kw real*8 m,lum,rad,rzams,menv,fac real*8 lam1,lam2,m1,logm,logl real*8 aa,bb,cc,dd * * A function to estimate lambda for common-envelope. * if(fac.ge.0.d0)then * * No fits yet for naked He stars... * if(kw.gt.6)then celamf = 0.5d0 goto 90 endif * if(menv.gt.0.d0)then * Formulae for giant-like stars; also used for HG and CHeB stars close * to the Hayashi track. logl = log10(lum) logm = log10(m) if(kw.le.5)then m1 = m if(kw.gt.3) m1 = 100.d0 lam1 = 3.d0/(2.4d0 + 1.d0/m1**(3.d0/2.d0)) - 0.15d0*logl lam1 = MIN(lam1,0.8d0) else lam1 = -3.5d0 - 0.75d0*logm + logl endif if(kw.gt.3)then lam2 = MIN(0.9d0,0.58d0 + 0.75d0*logm) - 0.08d0*logl if(kw.lt.6)then lam1 = MIN(lam2,lam1) else lam1 = MAX(lam2,lam1) lam1 = MIN(lam1,1.d0) endif endif lam1 = 2.d0*lam1 if(fac.gt.0.d0)then * Use a fraction FAC of the ionization energy in the energy balance. if(kw.le.3)then aa = MIN(1.2d0*(logm - 0.25d0)**2 - 0.7d0,-0.5d0) else aa = MAX(-0.2d0 - logm,-0.5d0) endif bb = MAX(3.d0 - 5.d0*logm,1.5d0) cc = MAX(3.7d0 + 1.6d0*logm,3.3d0 + 2.1d0*logm) lam2 = aa + ATAN(bb*(cc - logl)) if(kw.le.3)then dd = MAX(0.d0,MIN(0.15d0,0.15d0 - 0.25d0*logm)) lam2 = lam2 + dd*(logl - 2.d0) endif lam2 = MAX(lam2,1.d-2) lam2 = MAX(1.d0/lam2,lam1) if(fac.ge.1.d0)then lam1 = lam2 else lam1 = lam1 + fac*(lam2 - lam1) endif endif endif * if(menv.lt.1.d0)then * Formula for HG stars; also reasonable for CHeB stars in blue loop. lam2 = 0.42d0*(rzams/rad)**0.4d0 * Alternatively: * lam2 = 0.3d0*(rtms/rad)**0.4d0 lam2 = 2.d0*lam2 endif * if(menv.le.0.d0)then celamf = lam2 elseif(menv.ge.1.d0)then celamf = lam1 else * Interpolate between HG and GB values depending on conv. envelope mass. celamf = lam2 + sqrt(menv)*(lam1 - lam2) endif * 90 continue * else celamf = -1.d0*fac endif * return end *** cc//comenv.f *** SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, & M02,M2,MC2,AJ2,JSPIN2,KW2, & ZPARS,ECC,SEP,JORB,COEL) * * Common Envelope Evolution. * * Author : C. A. Tout * Date : 18th September 1996 * * Redone : J. R. Hurley * Date : 7th July 1998 * IMPLICIT NONE * INTEGER KW1,KW2,KW INTEGER KTYPE(0:14,0:14) COMMON /TYPES/ KTYPE INTEGER ceflag,tflag,ifflag,nsflag,wdflag COMMON /FLAGS/ ceflag,tflag,ifflag,nsflag,wdflag * REAL*8 M01,M1,MC1,AJ1,JSPIN1,R1,L1,K21 REAL*8 M02,M2,MC2,AJ2,JSPIN2,R2,L2,K22,MC22 REAL*8 TSCLS1(20),TSCLS2(20),LUMS(10),GB(10),TM1,TM2,TN,ZPARS(20) REAL*8 EBINDI,EBINDF,EORBI,EORBF,ECIRC,SEPF,SEPL,MF,XX REAL*8 CONST,DELY,DERI,DELMF,MC3,FAGE1,FAGE2 REAL*8 ECC,SEP,JORB,TB,OORB,OSPIN1,OSPIN2,TWOPI REAL*8 RC1,RC2,Q1,Q2,RL1,RL2,LAMB1,LAMB2 REAL*8 MENV,RENV,MENVD,RZAMS,VS(3) REAL*8 AURSUN,K3,ALPHA1,LAMBDA PARAMETER (AURSUN = 214.95D0,K3 = 0.21D0) COMMON /VALUE2/ ALPHA1,LAMBDA LOGICAL COEL REAL*8 CELAMF,RL,RZAMSF EXTERNAL CELAMF,RL,RZAMSF * * Common envelope evolution - entered only when KW1 = 2, 3, 4, 5, 6, 8 or 9. * * For simplicity energies are divided by -G. * TWOPI = 2.D0*ACOS(-1.D0) COEL = .FALSE. * * Obtain the core masses and radii. * KW = KW1 CALL star(KW1,M01,M1,TM1,TN,TSCLS1,LUMS,GB,ZPARS) CALL hrdiag(M01,AJ1,M1,TM1,TN,TSCLS1,LUMS,GB,ZPARS, & R1,L1,KW1,MC1,RC1,MENV,RENV,K21) OSPIN1 = JSPIN1/(K21*R1*R1*(M1-MC1)+K3*RC1*RC1*MC1) MENVD = MENV/(M1-MC1) RZAMS = RZAMSF(M01) LAMB1 = CELAMF(KW,M01,L1,R1,RZAMS,MENVD,LAMBDA) KW = KW2 CALL star(KW2,M02,M2,TM2,TN,TSCLS2,LUMS,GB,ZPARS) CALL hrdiag(M02,AJ2,M2,TM2,TN,TSCLS2,LUMS,GB,ZPARS, & R2,L2,KW2,MC2,RC2,MENV,RENV,K22) OSPIN2 = JSPIN2/(K22*R2*R2*(M2-MC2)+K3*RC2*RC2*MC2) * * Calculate the binding energy of the giant envelope (multiplied by lambda). * EBINDI = M1*(M1-MC1)/(LAMB1*R1) * * If the secondary star is also giant-like add its envelopes's energy. * EORBI = M1*M2/(2.D0*SEP) IF(KW2.GE.2.AND.KW2.LE.9.AND.KW2.NE.7)THEN MENVD = MENV/(M2-MC2) RZAMS = RZAMSF(M02) LAMB2 = CELAMF(KW,M02,L2,R2,RZAMS,MENVD,LAMBDA) EBINDI = EBINDI + M2*(M2-MC2)/(LAMB2*R2) * * Calculate the initial orbital energy * IF(CEFLAG.NE.3) EORBI = MC1*MC2/(2.D0*SEP) ELSE IF(CEFLAG.NE.3) EORBI = MC1*M2/(2.D0*SEP) ENDIF * * Allow for an eccentric orbit. * ECIRC = EORBI/(1.D0 - ECC*ECC) * * Calculate the final orbital energy without coalescence. * EORBF = EORBI + EBINDI/ALPHA1 * * If the secondary is on the main sequence see if it fills its Roche lobe. * IF(KW2.LE.1.OR.KW2.EQ.7)THEN SEPF = MC1*M2/(2.D0*EORBF) Q1 = MC1/M2 Q2 = 1.D0/Q1 RL1 = RL(Q1) RL2 = RL(Q2) IF(RC1/RL1.GE.R2/RL2)THEN * * The helium core of a very massive star of type 4 may actually fill * its Roche lobe in a wider orbit with a very low-mass secondary. * IF(RC1.GT.RL1*SEPF)THEN COEL = .TRUE. SEPL = RC1/RL1 ENDIF ELSE IF(R2.GT.RL2*SEPF)THEN COEL = .TRUE. SEPL = R2/RL2 ENDIF ENDIF IF(COEL)THEN * KW = KTYPE(KW1,KW2) - 100 MC3 = MC1 IF(KW2.EQ.7.AND.KW.EQ.4) MC3 = MC3 + M2 * * Coalescence - calculate final binding energy. * EORBF = MAX(MC1*M2/(2.D0*SEPL),EORBI) EBINDF = EBINDI - ALPHA1*(EORBF - EORBI) ELSE * * Primary becomes a black hole, neutron star, white dwarf or helium star. * MF = M1 M1 = MC1 CALL star(KW1,M01,M1,TM1,TN,TSCLS1,LUMS,GB,ZPARS) CALL hrdiag(M01,AJ1,M1,TM1,TN,TSCLS1,LUMS,GB,ZPARS, & R1,L1,KW1,MC1,RC1,MENV,RENV,K21) IF(KW1.GE.13)THEN CALL kick(KW1,MF,M1,M2,ECC,SEPF,JORB,VS) IF(ECC.GT.1.D0) GOTO 30 ENDIF ENDIF ELSE * * Degenerate or giant secondary. Check if the least massive core fills its * Roche lobe. * SEPF = MC1*MC2/(2.D0*EORBF) Q1 = MC1/MC2 Q2 = 1.D0/Q1 RL1 = RL(Q1) RL2 = RL(Q2) IF(RC1/RL1.GE.RC2/RL2)THEN IF(RC1.GT.RL1*SEPF)THEN COEL = .TRUE. SEPL = RC1/RL1 ENDIF ELSE IF(RC2.GT.RL2*SEPF)THEN COEL = .TRUE. SEPL = RC2/RL2 ENDIF ENDIF * IF(COEL)THEN * * If the secondary was a neutron star or black hole the outcome * is an unstable Thorne-Zytkow object that leaves only the core. * SEPF = 0.D0 IF(KW2.GE.13)THEN MC1 = MC2 M1 = MC1 MC2 = 0.D0 M2 = 0.D0 KW1 = KW2 KW2 = 15 AJ1 = 0.D0 * * The envelope mass is not required in this case. * GOTO 30 ENDIF * KW = KTYPE(KW1,KW2) - 100 MC3 = MC1 + MC2 * * Calculate the final envelope binding energy. * EORBF = MAX(MC1*MC2/(2.D0*SEPL),EORBI) EBINDF = EBINDI - ALPHA1*(EORBF - EORBI) * * Check if we have the merging of two degenerate cores and if so * then see if the resulting core will survive or change form. * IF(KW1.EQ.6.AND.(KW2.EQ.6.OR.KW2.GE.11))THEN CALL dgcore(KW1,KW2,KW,MC1,MC2,MC3,EBINDF) ENDIF IF(KW1.LE.3.AND.M01.LE.ZPARS(2))THEN IF((KW2.GE.2.AND.KW2.LE.3.AND.M02.LE.ZPARS(2)).OR. & KW2.EQ.10)THEN CALL dgcore(KW1,KW2,KW,MC1,MC2,MC3,EBINDF) IF(KW.GE.10)THEN KW1 = KW M1 = MC3 MC1 = MC3 IF(KW.LT.15) M01 = MC3 AJ1 = 0.D0 MC2 = 0.D0 M2 = 0.D0 KW2 = 15 GOTO 30 ENDIF ENDIF ENDIF * ELSE * * The cores do not coalesce - assign the correct masses and ages. * MF = M1 M1 = MC1 CALL star(KW1,M01,M1,TM1,TN,TSCLS1,LUMS,GB,ZPARS) CALL hrdiag(M01,AJ1,M1,TM1,TN,TSCLS1,LUMS,GB,ZPARS, & R1,L1,KW1,MC1,RC1,MENV,RENV,K21) IF(KW1.GE.13)THEN CALL kick(KW1,MF,M1,M2,ECC,SEPF,JORB,VS) IF(ECC.GT.1.D0) GOTO 30 ENDIF MF = M2 KW = KW2 M2 = MC2 CALL star(KW2,M02,M2,TM2,TN,TSCLS2,LUMS,GB,ZPARS) CALL hrdiag(M02,AJ2,M2,TM2,TN,TSCLS2,LUMS,GB,ZPARS, & R2,L2,KW2,MC2,RC2,MENV,RENV,K22) IF(KW2.GE.13.AND.KW.LT.13)THEN CALL kick(KW2,MF,M2,M1,ECC,SEPF,JORB,VS) IF(ECC.GT.1.D0) GOTO 30 ENDIF ENDIF ENDIF * IF(COEL)THEN MC22 = MC2 IF(KW.EQ.4.OR.KW.EQ.7)THEN * If making a helium burning star calculate the fractional age * depending on the amount of helium that has burnt. IF(KW1.LE.3)THEN FAGE1 = 0.D0 ELSEIF(KW1.GE.6)THEN FAGE1 = 1.D0 ELSE FAGE1 = (AJ1 - TSCLS1(2))/(TSCLS1(13) - TSCLS1(2)) ENDIF IF(KW2.LE.3.OR.KW2.EQ.10)THEN FAGE2 = 0.D0 ELSEIF(KW2.EQ.7)THEN FAGE2 = AJ2/TM2 MC22 = M2 ELSEIF(KW2.GE.6)THEN FAGE2 = 1.D0 ELSE FAGE2 = (AJ2 - TSCLS2(2))/(TSCLS2(13) - TSCLS2(2)) ENDIF ENDIF ENDIF * * Now calculate the final mass following coelescence. This requires a * Newton-Raphson iteration. * IF(COEL)THEN * * Calculate the orbital spin just before coalescence. * TB = (SEPL/AURSUN)*SQRT(SEPL/(AURSUN*(MC1+MC2))) OORB = TWOPI/TB * XX = 1.D0 + ZPARS(7) IF(EBINDF.LE.0.D0)THEN MF = MC3 GOTO 20 ELSE CONST = ((M1+M2)**XX)*(M1-MC1+M2-MC22)*EBINDF/EBINDI ENDIF * * Initial Guess. * MF = MAX(MC1 + MC22,(M1 + M2)*(EBINDF/EBINDI)**(1.D0/XX)) 10 DELY = (MF**XX)*(MF - MC1 - MC22) - CONST * IF(ABS(DELY/MF**(1.D0+XX)).LE.1.0D-02) GOTO 20 IF(ABS(DELY/MF).LE.1.0D-03) GOTO 20 DERI = MF**ZPARS(7)*((1.D0+XX)*MF - XX*(MC1 + MC22)) DELMF = DELY/DERI MF = MF - DELMF GOTO 10 * * Set the masses and separation. * 20 IF(MC22.EQ.0.D0) MF = MAX(MF,MC1+M2) M2 = 0.D0 M1 = MF KW2 = 15 * * Combine the core masses. * IF(KW.EQ.2)THEN CALL star(KW,M1,M1,TM2,TN,TSCLS2,LUMS,GB,ZPARS) IF(GB(9).GE.MC1)THEN M01 = M1 AJ1 = TM2 + (TSCLS2(1) - TM2)*(AJ1-TM1)/(TSCLS1(1) - TM1) CALL star(KW,M01,M1,TM1,TN,TSCLS1,LUMS,GB,ZPARS) ENDIF ELSEIF(KW.EQ.7)THEN M01 = M1 CALL star(KW,M01,M1,TM1,TN,TSCLS1,LUMS,GB,ZPARS) AJ1 = TM1*(FAGE1*MC1 + FAGE2*MC22)/(MC1 + MC22) ELSEIF(KW.EQ.4.OR.MC2.GT.0.D0.OR.KW.NE.KW1)THEN IF(KW.EQ.4) AJ1 = (FAGE1*MC1 + FAGE2*MC22)/(MC1 + MC22) MC1 = MC1 + MC2 MC2 = 0.D0 * * Obtain a new age for the giant. * CALL gntage(MC1,M1,KW,ZPARS,M01,AJ1) CALL star(KW,M01,M1,TM1,TN,TSCLS1,LUMS,GB,ZPARS) ENDIF CALL hrdiag(M01,AJ1,M1,TM1,TN,TSCLS1,LUMS,GB,ZPARS, & R1,L1,KW,MC1,RC1,MENV,RENV,K21) JSPIN1 = OORB*(K21*R1*R1*(M1-MC1)+K3*RC1*RC1*MC1) KW1 = KW ECC = 0.D0 ELSE * * Check if any eccentricity remains in the orbit by first using * energy to circularise the orbit before removing angular momentum. * (note this should not be done in case of CE SN ... fix). * IF(EORBF.LT.ECIRC)THEN ECC = SQRT(1.D0 - EORBF/ECIRC) ELSE ECC = 0.D0 ENDIF * * Set both cores in co-rotation with the orbit on exit of CE, * TB = (SEPF/AURSUN)*SQRT(SEPF/(AURSUN*(M1+M2))) OORB = TWOPI/TB JORB = M1*M2/(M1+M2)*SQRT(1.D0-ECC*ECC)*SEPF*SEPF*OORB * JSPIN1 = OORB*(K21*R1*R1*(M1-MC1)+K3*RC1*RC1*MC1) * JSPIN2 = OORB*(K22*R2*R2*(M2-MC2)+K3*RC2*RC2*MC2) * * or, leave the spins of the cores as they were on entry. * Tides will deal with any synchronization later. * JSPIN1 = OSPIN1*(K21*R1*R1*(M1-MC1)+K3*RC1*RC1*MC1) JSPIN2 = OSPIN2*(K22*R2*R2*(M2-MC2)+K3*RC2*RC2*MC2) ENDIF 30 SEP = SEPF RETURN END *** cc//corerd.f *** REAL*8 FUNCTION CORERD(KW,MC,M0,MFLASH) * * A function to determine the radius of the core of a giant-like star. * NOTE: this is out of date so rc should be obtained using HRDIAG! * It is still OK to use but bear in mind that the core radius calculated * for non-degenerate giant cores is only a rough estimate. * * Author : C. A. Tout * Date : 26th February 1997 * Updated 6/1/98 by J. Hurley * IMPLICIT NONE INTEGER KW REAL*8 MC,MCH,M0,MFLASH PARAMETER (MCH = 1.44d0) * * First do the black holes and neutron stars. * IF(KW.EQ.14)THEN CORERD = 4.24d-06*MC ELSEIF(KW.EQ.13)THEN CORERD = 1.4d-05 * * Main sequence stars. * ELSEIF(KW.LE.1.OR.KW.EQ.7)THEN CORERD = 0.d0 * * Core-helium-burning stars, FAGB stars and non-degenerate giant cores. * ELSEIF(KW.EQ.4.OR.KW.EQ.5.OR.(KW.LE.3.AND.M0.GT.MFLASH))THEN CORERD = 0.2239d0*MC**0.62d0 * * The degenerate giants and white dwarfs. * ELSE CORERD = 0.0115d0*SQRT(MAX(1.48204d-06,(MCH/MC)**(2.d0/3.d0) & - (MC/MCH)**(2.d0/3.d0))) * * Degenerate giants have hot subdwarf cores. * IF(KW.LE.9) CORERD = 5.d0*CORERD ENDIF * RETURN END *** cc//dgcore.f *** SUBROUTINE dgcore(kw1,kw2,kw3,m1,m2,m3,ebinde) * * A routine to determine the outcome of a collision or coalescence * of two degenerate cores. * Entered with kw1,kw2 = 2 or 3 with M <= Mflash, 6, 10, 11 or 12 * implicit none * integer kw1,kw2,kw3 * real*8 m1,m2,m3,ebinde real*8 r1,r2,r3,mhe,mc,mne,ebindi,ebindf,deleb,de,enuc real*8 temp,x,y,m0,mflash real*8 cvhe,cvc,cvne parameter(cvhe=3.1d+07,cvc=8.27d+06,cvne=7.44d+06) real*8 ehe,ec,ene parameter(ehe=5.812d+17,ec=2.21d+17,ene=2.06d+17) real*8 the,tc,gmr,mch parameter(the=1.0d+08,tc=1.0d+09,gmr=1.906d+15,mch=1.44d0) * real*8 corerd external corerd * * Calculate the core radii setting m0 < mflash using dummy values as we * know it to be true if kw = 2 or 3. m0 = 1.d0 mflash = 2.d0 r1 = corerd(kw1,m1,m0,mflash) r2 = corerd(kw2,m2,m0,mflash) r3 = corerd(kw3,m3,m0,mflash) * Calculate the initial binding energy of the two seperate cores and the * difference between this and the final binding energy of the new core. ebindi = m1*m1/r1 + m2*m2/r2 ebindf = m3*m3/r3 deleb = ABS(ebindi - ebindf) * If an envelope is present reduce its binding energy by the amount * of energy liberated by the coalescence. ebinde = MAX(0.d0,ebinde - deleb) if(kw1.gt.3) goto 90 * Distribute the mass into core mass groups where mhe represents Helium * core mass, mc represents Carbon core mass and mne represents a Neon * core mass or any mass that is all converted Carbon. mhe = 0.d0 mc = 0.d0 mne = 0.d0 if(kw1.le.3.or.kw1.eq.10)then mhe = mhe + m1 elseif(kw1.eq.12)then mne = mne + m1 else mc = mc + m1 endif if(kw2.le.3.or.kw2.eq.10)then mhe = mhe + m2 elseif(kw2.eq.12)then mne = mne + m2 else mc = mc + m2 endif * Calculate the temperature generated by the merging of the cores. temp = (deleb/(cvhe*mhe+cvc*mc+cvne*mne))*gmr * * To decide if He is converted to C we use: * 3He4 -> C12 , T > 10^8 K , 7.274 Mev released, * to decide if C is converted further we use: * 2C12 -> Ne20 + alpha , T > 10^9 K , 4.616 Mev released. * and to decide if O is converted further we use: * 2O16 -> P31 + p , T > 10^9 K , 7.677 Mev released. * To obtain the heat capacity of an O/Ne WD and to gain an idea of the * energy released from the further processing of an O/Ne WD we use: * 2Ne20 + gamma -> O16 + Mg24 +gamma , T > 10^9 K , 4.583 Mev released. * For a CO core the composition is assumed to be 20% C, 80% O and for * an ONe core 80% O, 20% Ne. * * Decide if conversion of elements can take place. * if(temp.gt.the)then x = 1.d0 * else * x = 0.d0 * endif * if(temp.gt.tc)then * y = 1.d0 * else y = 0.d0 * endif * Calculate the nuclear energy generated from any element conversion. enuc = (x*ehe*mhe + y*(ec*mc + ene*mne))/gmr * Calculate the difference between the binding energy of the star * (core + envelope) and the nuclear energy. The new star will be * destroyed in a SN if dE < 0. de = (ebindf + ebinde) - enuc * If the star survives and an envelope is present then reduce the * envelope binding energy by the amount of liberated nuclear energy. * The envelope will not survive if its binding energy is reduced to <= 0. if(de.ge.0.d0) ebinde = MAX(0.d0,ebinde - enuc) * Now determine the final evolution state of the merged core. if(de.lt.0.d0) kw3 = 15 if(kw3.eq.3)then if(x.gt.0.5d0)then kw3 = 6 elseif(ebinde.le.0.d0)then kw3 = 10 endif elseif(kw3.eq.4)then if(x.gt.0.5d0)then kw3 = 6 elseif(ebinde.le.0.d0)then kw3 = 7 endif endif if(kw3.eq.6.and.y.lt.0.5d0)then if(ebinde.le.0.d0) kw3 = 11 elseif(kw3.eq.6.and.y.gt.0.5d0)then if(ebinde.le.0.d0) kw3 = 12 endif if(kw3.eq.10.and.x.gt.0.5d0) kw3 = 11 if(kw3.eq.11.and.y.gt.0.5d0) kw3 = 12 if(kw3.ge.10.and.kw3.le.12.and.m3.ge.mch) kw3 = 15 * if(kw3.eq.15) m3 = 0.d0 90 continue * return end *** cc//evolv2.f *** SUBROUTINE evolv2(kstar,mass0,mass,rad,lumin,massc,radc, & menv,renv,ospin,epoch,tms, & tphys,tphysf,dtp,z,zpars,tb,ecc,vkick) implicit none *** * * B I N A R Y * *********** * * Roche lobe overflow. * -------------------- * * Developed by Jarrod Hurley, IOA, Cambridge. * ......................................................... * * Advice by Christopher Tout, Onno Pols & Sverre Aarseth. * ++++++++++++++++++++++++++++++++++++++++++++++++++ * * Adapted from Aarseth's code 21st September 1996. * Fully revised on 27th November 1996 to remove vestiges of N-body code and * incorporate corrections. * Fully revised on 1st April 1998 to include new stellar evolution formulae * and associated binary evolution changes. * Fully revised on 4th July 1998 to include eccentricity, tidal * circularization, wind accretion, velocity kicks for supernovae and all * associated orbital momentum changes. * Revised on 31st October 2000 to upgrade restrictions imposed on the * timestep owing to magnetic braking and orbital angular momentum changes. * *** * * See Tout et al., 1997, MNRAS, 291, 732 for a description of many of the * processes in this code as well as the relevant references mentioned * within the code. * * Reference for the stellar evolution formulae is Hurley, Pols & Tout, * 2000, MNRAS, 315, 543 (SSE paper). * Reference for the binary evolution algorithm is Hurley, Tout & Pols, * 2002, MNRAS, 329, 897 (BSE paper). * *** * * March 2001 * * Changes since version 3, i.e. since production of Paper3: * * 1) The Eddington limit flag (on/off) has been replaced by an * Eddington limit multiplicative factor (eddfac). So if you * want to neglect the Eddington limit you would set eddfac * to a large value. * * 2) To determine whether material transferred during RLOF forms * an accretion disk around the secondary or hits the secondary * in a direct stream we calculate a minimum radial distance, rmin, * of the mass stream from the secondary. This is taken from eq.(1) * of Ulrich & Burger (1976, ApJ, 206, 509) which they fitted to * the calculations of Lubow & Shu (1974, ApJ, 198, 383). * If rmin is less than the radius of the secondary then an * accretion disk is not formed. * Note that the formula for rmin given by Ulrich & Burger is valid * for all q whereas that given by Nelemans et al. (2001, A&A, * submitted) in their eq.(6) is only valid for q < 1 where * they define q = Mdonor/Maccretor, i.e. DD systems. * * 3) The changes to orbital and spin angular momentum owing to * RLOF mass transfer have been improved, and an new input option * now exists. * When mass is lost from the system during RLOF there are now * three choices as to how the orbital angular momentum is * affected: a) the lost material carries with it a fraction * gamma of the orbital angular momentum, i.e. * dJorb = gamma*dm*a^2*omega_orb; b) the material carries with it * the specific angular momentum of the primary, i.e. * dJorb = dm*a_1^2*omega_orb; or c) assume the material is lost * from the system as if a wind from the secondary, i.e. * dJorb = dm*a_2^2*omega_orb. * The parameter gamma is an input option. * Choice c) is used if the mass transfer is super-Eddington * or the system is experiencing novae eruptions. * In all other cases choice a) is used if gamma > 0.0, b) if * gamma = -1.0 and c) is used if gamma = -2.0. * The primary spin angular momentum is reduced by an amount * dm1*r_1^2*omega_1 when an amount of mass dm1 is transferred * from the primary. * If the secondary accretes through a disk then its spin * angular momentum is altered by assuming that the material * falls onto the star from the inner edge of a Keplerian * disk and that the system is in a steady state, i.e. * an amount dm2*SQRT(G*m_2*r_2). * If there is no accretion disk then we calculate the angular * momentum of the transferred material by using the radius at * at which the disk would have formed (rdisk = 1.7*rmin, see * Ulrich & Burger 1976) if allowed, i.e. the angular momentum * of the inner Lagrangian point, and add this directly to * the secondary, i.e. an amount dm2*SQRT(G*m_2*rdisk). * Total angular momentum is conserved in this model. * * 4) Now using q_crit = 3.0 for MS-MS Roche systems (previously we * had nothing). This corresponds roughly to R proportional to M^5 * which should be true for the majority of the MS (varies from * (M^17 -> M^2). If q > q_crit then contact occurs. * For CHeB primaries we also take q_crit = 3.0 and allow * common-envelope to occur if this is exceeded. * * 5) The value of lambda used in calculations of the envelope binding * energy for giants in common-envelope is now variable (see function * in zfuncs). The lambda function has been fitted by Onno to detailed * models ... he will write about this soon! * * 6) Note that eq.42 in the paper is missing a SQRT around the * MR^2/a^5 part. This needs to be corrected in any code update * paper with a thanks to Jeremy Sepinsky (student at NorthWestern). * It is ok in the code. * * March 2003 * * New input options added: * * ifflag - for the mass of a WD you can choose to use the mass that * results from the evolution algorithm (basically a competition * between core-mass growth and envelope mass-loss) or use the IFMR * proposed by Han, Podsiadlowski & Eggleton, 1995, MNRAS, 272, 800 * [>0 activates HPE IFMR]. * * wdflag - for the cooling of WDs you can choose to use either the standard * Mestel cooling law (see SSE paper) or a modified-Mestel law that * is better matched to detailed models (provided by Brad Hansen * ... see Hurley & Shara, 2003, ApJ, May 20, in press) * [>0 activates modified-Mestel]. * * bhflag - choose whether or not black holes should get velocity kicks * at formation * [0= no kick; >0 kick]. * * nsflag - for the mass of neutron stars and black holes you can use either * the SSE prescription or the prescription presented by * Belczynski et al. 2002, ApJ, 572, 407 who found that SSE was * underestimating the masses of these stars. In either case you also * need to set the maximum NS mass (mxns) for the prescription * [0= SSE, mxns=1.8; >0 Belczynski, mxns=3.0]. * * Sept 2004 * * Input options added/changed: * * ceflag - set to 3 this uses de Kool (or Podsiadlowski) CE prescription, * other options, such as Yungelson, could be added as well. * * hewind - factor to control the amount of He star mass-loss, i.e. * 1.0e-13*hewind*L^(2/3) gives He star mass-loss. * * * ++++++++++++++++++++++++++++++++++++++++++++++++++ *** * INTEGER loop,iter,intpol,k,ip,jp,j1,j2 PARAMETER(loop=20000) INTEGER kstar(2),kw,kst,kw1,kw2,kmin,kmax INTEGER ktype(0:14,0:14) COMMON /TYPES/ ktype INTEGER ceflag,tflag,ifflag,nsflag,wdflag COMMON /FLAGS/ ceflag,tflag,ifflag,nsflag,wdflag * REAL*8 km,km0,tphys,tphys0,dtm0,tphys00 REAL*8 tphysf,dtp,tsave REAL*8 aj(2),aj0(2),epoch(2),tms(2),tbgb(2),tkh(2),dtmi(2) REAL*8 mass0(2),mass(2),massc(2),menv(2),mass00(2),mcxx(2) REAL*8 rad(2),rol(2),rol0(2),rdot(2),radc(2),renv(2),radx(2) REAL*8 lumin(2),k2str(2),q(2),dms(2),dmr(2),dmt(2),vkick(2) REAL*8 dml,vorb2,vwind2,omv2,ivsqm,lacc,vs(3) REAL*8 sep,dr,tb,dme,tdyn,taum,dm1,dm2,dmchk,qc,dt,pd,rlperi REAL*8 m1ce,m2ce,mch,tmsnew,dm22,mew PARAMETER(mch=1.44d0) REAL*8 yeardy,aursun PARAMETER(yeardy=365.24d0,aursun=214.95d0) REAL*8 acc1,tiny PARAMETER(acc1=3.920659d+08,tiny=1.0d-14) REAL*8 ecc,ecc1,tc,tcirc,ttid,ecc2,omecc2,sqome2,sqome3,sqome5 REAL*8 f1,f2,f3,f4,f5,f,raa2,raa6,eqspin,rg2,tcqr REAL*8 k3,mr23yr,twopi PARAMETER(k3=0.21d0,mr23yr=0.4311d0) REAL*8 jspin(2),ospin(2),jorb,oorb,jspbru,ospbru REAL*8 delet,delet1,dspint(2),djspint(2),djtx(2) REAL*8 dtj,djorb,djgr,djmb,djt,djtt,rmin,rdisk REAL*8 neta,bwind,hewind,mxns COMMON /VALUE1/ neta,bwind,hewind,mxns REAL*8 beta,xi,acc2,epsnov,eddfac,gamma COMMON /VALUE5/ beta,xi,acc2,epsnov,eddfac,gamma * REAL*8 z,tm,tn,m0,mt,rm,lum,mc,rc,me,re,k2,age,dtm,dtr REAL*8 tscls(20),lums(10),GB(10),zpars(20) REAL*8 zero,ngtv,ngtv2,mt2,rrl1,rrl2,mcx,teff1,teff2 REAL*8 mass1i,mass2i,tbi,ecci LOGICAL coel,com,prec,inttry,change,snova,sgl,bsymb,esymb,bss LOGICAL supedd,novae,disk LOGICAL isave,iplot REAL*8 rl,mlwind,vrotf,corerd EXTERNAL rl,mlwind,vrotf,corerd REAL bcm(50000,34),bpp(80,10) COMMON /BINARY/ bcm,bpp * * Save the initial state. * mass1i = mass0(1) mass2i = mass0(2) tbi = tb ecci = ecc * zero = 0.d0 ngtv = -1.d0 ngtv2 = -2.d0 twopi = 2.d0*ACOS(-1.d0) * * Initialize the parameters. * kmin = 1 kmax = 2 sgl = .false. mt2 = MIN(mass(1),mass(2)) kst = 0 * if(mt2.lt.tiny.or.tb.le.0.d0)then sgl = .true. if(mt2.lt.tiny)then mt2 = 0.d0 if(mass(1).lt.tiny)then if(tphys.lt.tiny)then mass0(1) = 0.01d0 mass(1) = mass0(1) kst = 1 else kmin = 2 lumin(1) = 1.0d-10 rad(1) = 1.0d-10 massc(1) = 0.d0 dmt(1) = 0.d0 dmr(1) = 0.d0 endif ospin(1) = 1.0d-10 jspin(1) = 1.0d-10 else if(tphys.lt.tiny)then mass0(2) = 0.01d0 mass(2) = mass0(2) kst = 2 else kmax = 1 lumin(2) = 1.0d-10 rad(2) = 1.0d-10 massc(2) = 0.d0 dmt(2) = 0.d0 dmr(2) = 0.d0 endif ospin(2) = 1.0d-10 jspin(2) = 1.0d-10 endif endif ecc = -1.d0 tb = 0.d0 sep = 1.0d+10 oorb = 0.d0 jorb = 0.d0 if(ospin(1).lt.0.0) ospin(1) = 1.0d-10 if(ospin(2).lt.0.0) ospin(2) = 1.0d-10 q(1) = 1.0d+10 q(2) = 1.0d+10 rol(1) = 1.0d+10 rol(2) = 1.0d+10 else tb = tb/yeardy sep = aursun*(tb*tb*(mass(1) + mass(2)))**(1.d0/3.d0) oorb = twopi/tb jorb = mass(1)*mass(2)/(mass(1)+mass(2)) & *SQRT(1.d0-ecc*ecc)*sep*sep*oorb if(ospin(1).lt.0.d0) ospin(1) = oorb if(ospin(2).lt.0.d0) ospin(2) = oorb endif * do 500 , k = kmin,kmax age = tphys - epoch(k) CALL star(kstar(k),mass0(k),mass(k),tm,tn,tscls,lums,GB,zpars) CALL hrdiag(mass0(k),age,mass(k),tm,tn,tscls,lums,GB,zpars, & rm,lum,kstar(k),mc,rc,me,re,k2) aj(k) = age epoch(k) = tphys - age rad(k) = rm lumin(k) = lum massc(k) = mc radc(k) = rc menv(k) = me renv(k) = re k2str(k) = k2 tms(k) = tm tbgb(k) = tscls(1) * if(tphys.lt.tiny.and.ospin(k).le.0.001d0)then ospin(k) = 45.35d0*vrotf(mass(k))/rm endif jspin(k) = ospin(k)*(k2*rm*rm*(mass(k)-mc)+k3*rc*rc*mc) if(.not.sgl)then q(k) = mass(k)/mass(3-k) rol(k) = rl(q(k))*sep endif rol0(k) = rol(k) dmr(k) = 0.d0 dmt(k) = 0.d0 djspint(k) = 0.d0 dtmi(k) = 1.0d+06 * 500 continue * if(mt2.lt.tiny)then sep = 0.d0 if(kst.gt.0)then mass0(kst) = 0.d0 mass(kst) = 0.d0 kmin = 3 - kst kmax = kmin endif endif * * On the first entry the previous timestep is zero to prevent mass loss. * dtm = 0.d0 delet = 0.d0 djorb = 0.d0 bss = .false. * * Setup variables which control the output (if it is required). * ip = 0 jp = 0 tsave = tphys isave = .true. iplot = .false. if(dtp.le.0.d0)then iplot = .true. isave = .false. tsave = tphysf elseif(dtp.gt.tphysf)then isave = .false. tsave = tphysf endif if(tphys.ge.tphysf) goto 140 * 4 iter = 0 intpol = 0 inttry = .false. change = .false. prec = .false. snova = .false. coel = .false. com = .false. bsymb = .false. esymb = .false. tphys0 = tphys ecc1 = ecc j1 = 1 j2 = 2 if(kstar(1).ge.10.and.kstar(1).le.14) dtmi(1) = 0.01d0 if(kstar(2).ge.10.and.kstar(2).le.14) dtmi(2) = 0.01d0 dm1 = 0.d0 dm2 = 0.d0 * 5 kw1 = kstar(1) kw2 = kstar(2) * dt = 1.0d+06*dtm eqspin = 0.d0 djtt = 0.d0 * if(intpol.eq.0.and.ABS(dtm).gt.tiny.and..not.sgl)then vorb2 = acc1*(mass(1)+mass(2))/sep ivsqm = 1.d0/SQRT(1.d0-ecc*ecc) do 501 , k = 1,2 * * Calculate wind mass loss from the previous timestep. * if(neta.gt.tiny)then rlperi = rol(k)*(1.d0-ecc) dmr(k) = mlwind(kstar(k),lumin(k),rad(k),mass(k), & massc(k),rlperi,z) * * Calculate how much of wind mass loss from companion will be * accreted (Boffin & Jorissen, A&A 1988, 205, 155). * vwind2 = 2.d0*beta*acc1*mass(k)/rad(k) omv2 = (1.d0 + vorb2/vwind2)**(3.d0/2.d0) dmt(3-k) = ivsqm*acc2*dmr(k)*((acc1*mass(3-k)/vwind2)**2) & /(2.d0*sep*sep*omv2) dmt(3-k) = MIN(dmt(3-k),0.8d0*dmr(k)) else dmr(k) = 0.d0 dmt(3-k) = 0.d0 endif 501 continue * * Diagnostic for Symbiotic-type stars. * if(neta.gt.tiny.and..not.esymb)then lacc = 3.14d+07*mass(j2)*dmt(j2)/rad(j2) lacc = lacc/lumin(j1) if((lacc.gt.0.01d0.and..not.bsymb).or. & (lacc.lt.0.01d0.and.bsymb))then jp = MIN(80,jp + 1) bpp(jp,1) = tphys bpp(jp,2) = mass(1) bpp(jp,3) = mass(2) bpp(jp,4) = float(kstar(1)) bpp(jp,5) = float(kstar(2)) bpp(jp,6) = sep bpp(jp,7) = ecc bpp(jp,8) = rad(1)/rol(1) bpp(jp,9) = rad(2)/rol(2) if(bsymb)then bpp(jp,10) = 13.0 esymb = .true. else bpp(jp,10) = 12.0 bsymb = .true. endif endif endif * * Calculate orbital angular momentum change due to wind mass loss. * ecc2 = ecc*ecc omecc2 = 1.d0 - ecc2 sqome2 = SQRT(omecc2) * djorb = ((dmr(1)+q(1)*dmt(1))*mass(2)*mass(2) + & (dmr(2)+q(2)*dmt(2))*mass(1)*mass(1))* & sep*sep*sqome2*oorb/(mass(1)+mass(2))**2 delet = ecc*(dmt(1)*(0.5d0/mass(1) + 1.d0/(mass(1)+mass(2))) + & dmt(2)*(0.5d0/mass(2) + 1.d0/(mass(1)+mass(2)))) * * For very close systems include angular momentum loss owing to * gravitational radiation. * if(sep.le.10.d0)then djgr = 8.315d-10*mass(1)*mass(2)*(mass(1)+mass(2))/ & (sep*sep*sep*sep) f1 = (19.d0/6.d0) + (121.d0/96.d0)*ecc2 sqome5 = sqome2**5 delet1 = djgr*ecc*f1/sqome5 djgr = djgr*jorb*(1.d0+0.875d0*ecc2)/sqome5 djorb = djorb + djgr delet = delet + delet1 endif * do 502 , k = 1,2 * * Calculate change in the intrinsic spin of the star. * djtx(k) = (2.d0/3.d0)*xi*dmt(k)*rad(3-k)*rad(3-k)*ospin(3-k) djspint(k) = (2.d0/3.d0)*(dmr(k)*rad(k)*rad(k)*ospin(k)) - & djtx(k) * * Include magnetic braking for stars that have appreciable convective * envelopes. This includes MS stars with M < 1.25, HG stars near the GB * and giants. MB is not allowed for fully convective MS stars. * if(mass(k).gt.0.35d0.and.kstar(k).lt.10)then djmb = 5.83d-16*menv(k)*(rad(k)*ospin(k))**3/mass(k) djspint(k) = djspint(k) + djmb * * Limit to a 3% angular momentum change for the star owing to MB. * This is found to work best with the maximum iteration of 20000, * i.e. does not create an excessive number of iterations, while not * affecting the evolution outcome when compared with a 2% restriction. * if(djmb.gt.tiny)then dtj = 0.03d0*jspin(k)/ABS(djmb) dt = MIN(dt,dtj) endif endif * * Calculate circularization, orbital shrinkage and spin up. * dspint(k) = 0.d0 if(((kstar(k).le.9.and.rad(k).ge.0.01d0*rol(k)).or. & (kstar(k).ge.10.and.k.eq.j1)).and.tflag.gt.0)then * raa2 = (rad(k)/sep)**2 raa6 = raa2**3 * * Hut's polynomials. * f5 = 1.d0+ecc2*(3.d0+ecc2*0.375d0) f4 = 1.d0+ecc2*(1.5d0+ecc2*0.125d0) f3 = 1.d0+ecc2*(3.75d0+ecc2*(1.875d0+ecc2*7.8125d-02)) f2 = 1.d0+ecc2*(7.5d0+ecc2*(5.625d0+ecc2*0.3125d0)) f1 = 1.d0+ecc2*(15.5d0+ecc2*(31.875d0+ecc2*(11.5625d0 & +ecc2*0.390625d0))) * if((kstar(k).eq.1.and.mass(k).ge.1.25d0).or. & kstar(k).eq.4.or.kstar(k).eq.7)then * * Radiative damping (Zahn, 1977, A&A, 57, 383 and 1975, A&A, 41, 329). * tc = 1.592d-09*(mass(k)**2.84d0) f = 1.9782d+04*SQRT((mass(k)*rad(k)*rad(k))/sep**5)* & tc*(1.d0+q(3-k))**(5.d0/6.d0) tcqr = f*q(3-k)*raa6 rg2 = k2str(k) elseif(kstar(k).le.9)then * * Convective damping (Hut, 1981, A&A, 99, 126). * tc = mr23yr*(menv(k)*renv(k)*(rad(k)-0.5d0*renv(k))/ & (3.d0*lumin(k)))**(1.d0/3.d0) ttid = twopi/(1.0d-10 + ABS(oorb - ospin(k))) f = MIN(1.d0,(ttid/(2.d0*tc)**2)) tcqr = 2.d0*f*q(3-k)*raa6*menv(k)/(21.d0*tc*mass(k)) rg2 = (k2str(k)*(mass(k)-massc(k)))/mass(k) else * * Degenerate damping (Campbell, 1984, MNRAS, 207, 433) * f = 7.33d-09*(lumin(k)/mass(k))**(5.d0/7.d0) tcqr = f*q(3-k)*q(3-k)*raa2*raa2/(1.d0+q(3-k)) rg2 = k3 endif * * Circularization. * sqome3 = sqome2**3 delet1 = 27.d0*tcqr*(1.d0+q(3-k))*raa2*(ecc/sqome2**13)* & (f3 - (11.d0/18.d0)*sqome3*f4*ospin(k)/oorb) tcirc = ecc/(ABS(delet1) + 1.0d-20) delet = delet + delet1 * * Spin up of star. * dspint(k) = (3.d0*q(3-k)*tcqr/(rg2*omecc2**6))* & (f2*oorb - sqome3*f5*ospin(k)) * * Calculate the equilibrium spin at which no angular momentum * can be transferred. * eqspin = oorb*f2/(sqome3*f5) * * Calculate angular momentum change for the star owing to tides. * djt = (k2str(k)*(mass(k)-massc(k))*rad(k)*rad(k) + & k3*massc(k)*radc(k)*radc(k))*dspint(k) if(kstar(k).le.6.or.ABS(djt)/jspin(k).gt.0.1d0)then djtt = djtt + djt endif endif 502 continue * * Limit to 2% orbital angular momentum change. * djtt = djtt + djorb if(ABS(djtt).gt.tiny)then dtj = 0.02d0*jorb/ABS(djtt) dt = MIN(dt,dtj) endif dtm = dt/1.0d+06 * elseif(ABS(dtm).gt.tiny.and.sgl)then do 503 , k = kmin,kmax if(neta.gt.tiny)then rlperi = 0.d0 dmr(k) = mlwind(kstar(k),lumin(k),rad(k),mass(k), & massc(k),rlperi,z) else dmr(k) = 0.d0 endif dmt(k) = 0.d0 djspint(k) = (2.d0/3.d0)*dmr(k)*rad(k)*rad(k)*ospin(k) if(mass(k).gt.0.35d0.and.kstar(k).lt.10)then djmb = 5.83d-16*menv(k)*(rad(k)*ospin(k))**3/mass(k) djspint(k) = djspint(k) + djmb if(djmb.gt.tiny)then dtj = 0.03d0*jspin(k)/ABS(djmb) dt = MIN(dt,dtj) endif endif 503 continue dtm = dt/1.0d+06 endif * do 504 , k = kmin,kmax * dms(k) = (dmr(k) - dmt(k))*dt if(kstar(k).lt.10)then dml = mass(k) - massc(k) if(dml.lt.dms(k))then dml = MAX(dml,2.d0*tiny) dtm = (dml/dms(k))*dtm if(k.eq.2) dms(1) = dms(1)*dml/dms(2) dms(k) = dml dt = 1.0d+06*dtm endif * * Limit to 1% mass loss. * if(dms(k).gt.0.01d0*mass(k))then dtm = 0.01d0*mass(k)*dtm/dms(k) if(k.eq.2) dms(1) = dms(1)*0.01d0*mass(2)/dms(2) dms(k) = 0.01d0*mass(k) dt = 1.0d+06*dtm endif endif * 504 continue * * Update mass and intrinsic spin (checking that the star is not spun * past the equilibrium) and reset epoch for a MS (and possibly a HG) star. * do 505 , k = kmin,kmax * if(eqspin.gt.0.d0.and.ABS(dspint(k)).gt.tiny)then if(intpol.eq.0)then if(dspint(k).ge.0.d0)then dspint(k) = MIN(dspint(k),(eqspin-ospin(k))/dt) else dspint(k) = MAX(dspint(k),(eqspin-ospin(k))/dt) endif djt = (k2str(k)*(mass(k)-massc(k))*rad(k)*rad(k) + & k3*massc(k)*radc(k)*radc(k))*dspint(k) djorb = djorb + djt djspint(k) = djspint(k) - djt endif endif * jspin(k) = MAX(1.0d-10,jspin(k) - djspint(k)*dt) * * Ensure that the star does not spin up beyond break-up. * ospbru = twopi*SQRT(mass(k)*aursun**3/rad(k)**3) jspbru = (k2str(k)*(mass(k)-massc(k))*rad(k)*rad(k) + & k3*massc(k)*radc(k)*radc(k))*ospbru if(jspin(k).gt.jspbru.and.ABS(dtm).gt.tiny)then mew = 1.d0 if(djtx(k).gt.0.d0)then mew = MIN(mew,(jspin(k) - jspbru)/djtx(k)) endif jspin(k) = jspbru * If excess material should not be accreted, activate next line. * dms(k) = dms(k) + (1.d0 - mew)*dmt(k)*dt endif * if(ABS(dms(k)).gt.tiny)then mass(k) = mass(k) - dms(k) if(kstar(k).le.2.or.kstar(k).eq.7)then m0 = mass0(k) mass0(k) = mass(k) CALL star(kstar(k),mass0(k),mass(k),tm,tn,tscls, & lums,GB,zpars) if(kstar(k).eq.2)then if(GB(9).lt.massc(k).or.m0.gt.zpars(3))then mass0(k) = m0 else epoch(k) = tm + (tscls(1) - tm)*(aj(k)-tms(k))/ & (tbgb(k) - tms(k)) epoch(k) = tphys - epoch(k) endif else epoch(k) = tphys - aj(k)*tm/tms(k) endif endif endif * 505 continue * if(.not.sgl)then * ecc1 = ecc1 - delet*dt ecc = MAX(ecc1,0.d0) if(ecc.lt.1.0d-10) ecc = 0.d0 * if(ecc.ge.1.d0) goto 135 * jorb = jorb - djorb*dt sep = (mass(1) + mass(2))*jorb*jorb/ & ((mass(1)*mass(2)*twopi)**2*aursun**3*(1.d0-ecc*ecc)) tb = (sep/aursun)*SQRT(sep/(aursun*(mass(1)+mass(2)))) oorb = twopi/tb endif * * Advance the time. * if(intpol.eq.0)then tphys0 = tphys dtm0 = dtm endif tphys = tphys + dtm * do 6 , k = kmin,kmax * * Acquire stellar parameters (M, R, L, Mc & K*) at apparent evolution age. * age = tphys - epoch(k) aj0(k) = age kw = kstar(k) m0 = mass0(k) mt = mass(k) mc = massc(k) if(intpol.eq.0) mcxx(k) = mc if(intpol.gt.0) mc = mcxx(k) mass00(k) = m0 * * Masses over 100Msun should probably not be trusted in the * evolution formulae. * if(mt.gt.100.d0)then WRITE(99,*)' MASS EXCEEDED ',mass1i,mass2i,tbi,ecci,mt goto 140 endif * CALL star(kw,m0,mt,tm,tn,tscls,lums,GB,zpars) CALL hrdiag(m0,age,mt,tm,tn,tscls,lums,GB,zpars, & rm,lum,kw,mc,rc,me,re,k2) * if(kw.ne.15)then ospin(k) = jspin(k)/(k2*(mt-mc)*rm*rm+k3*mc*rc*rc) endif * * At this point there may have been a supernova. * if(kw.ne.kstar(k).and.kstar(k).le.12.and. & (kw.eq.13.or.kw.eq.14))then if(sgl)then CALL kick(kw,mass(k),mt,0.d0,0.d0,-1.d0,0.d0,vs) vkick(k) = dsqrt(vs(1)*vs(1)+vs(2)*vs(2)+vs(3)*vs(3)) else CALL kick(kw,mass(k),mt,mass(3-k),ecc,sep,jorb,vs) vkick(k) = dsqrt(vs(1)*vs(1)+vs(2)*vs(2)+vs(3)*vs(3)) if(ecc.gt.1.d0)then kstar(k) = kw mass(k) = mt epoch(k) = tphys - age goto 135 endif tb = (sep/aursun)*SQRT(sep/(aursun*(mt+mass(3-k)))) oorb = twopi/tb endif snova = .true. endif * if(kw.ne.kstar(k))then change = .true. mass(k) = mt dtmi(k) = 0.01d0 if(kw.eq.15)then kstar(k) = kw goto 135 endif mass0(k) = m0 epoch(k) = tphys - age if(kw.gt.6.and.kstar(k).le.6)then bsymb = .false. esymb = .false. endif endif * * Force new NS or BH to have a second period. * if(kstar(k).eq.13.or.kstar(k).eq.14)then if(tphys-epoch(k).lt.tiny)then ospin(k) = 2.0d+08 jspin(k) = k3*rc*rc*mc*ospin(k) endif endif * * Set radius derivative for later interpolation. * if(ABS(dtm).gt.tiny)then rdot(k) = ABS(rm - rad(k))/dtm else rdot(k) = 0.d0 endif * * Base new time scale for changes in radius & mass on stellar type. * dt = dtmi(k) CALL deltat(kw,age,tm,tn,tscls,dt,dtr) * * Choose minimum of time-scale and remaining interval. * dtmi(k) = MIN(dt,dtr) * * Save relevent solar quantities. * aj(k) = age kstar(k) = kw rad(k) = rm lumin(k) = lum massc(k) = mc radc(k) = rc menv(k) = me renv(k) = re k2str(k) = k2 tms(k) = tm tbgb(k) = tscls(1) * * Check for blue straggler formation. * if(kw.le.1.and.tm.lt.tphys.and..not.bss)then bss = .true. jp = MIN(80,jp + 1) bpp(jp,1) = tphys bpp(jp,2) = mass(1) bpp(jp,3) = mass(2) bpp(jp,4) = float(kstar(1)) bpp(jp,5) = float(kstar(2)) bpp(jp,6) = sep bpp(jp,7) = ecc bpp(jp,8) = rad(1)/rol(1) bpp(jp,9) = rad(2)/rol(2) bpp(jp,10) = 14.0 endif * 6 continue * if(.not.sgl)then * * Determine the mass ratios. * do 506 , k = 1,2 q(k) = mass(k)/mass(3-k) 506 continue * * Determine the Roche lobe radii and adjust the radius derivative. * do 507 , k = 1,2 rol(k) = rl(q(k))*sep if(ABS(dtm).gt.tiny)then rdot(k) = rdot(k) + (rol(k) - rol0(k))/dtm rol0(k) = rol(k) endif 507 continue else do 508 , k = kmin,kmax rol(k) = 10000.d0*rad(k) 508 continue endif * if((tphys.lt.tiny.and.ABS(dtm).lt.tiny.and. & (mass2i.lt.0.1d0.or..not.sgl)).or.snova)then jp = MIN(80,jp + 1) bpp(jp,1) = tphys bpp(jp,2) = mass(1) bpp(jp,3) = mass(2) bpp(jp,4) = float(kstar(1)) bpp(jp,5) = float(kstar(2)) bpp(jp,6) = sep bpp(jp,7) = ecc bpp(jp,8) = rad(1)/rol(1) bpp(jp,9) = rad(2)/rol(2) bpp(jp,10) = 1.0 if(snova)then bpp(jp,10) = 2.0 dtm = 0.d0 goto 4 endif endif * if((isave.and.tphys.ge.tsave).or.iplot)then if(sgl.or.(rad(1).lt.rol(1).and.rad(2).lt.rol(2)). & or.tphys.lt.tiny)then ip = ip + 1 bcm(ip,1) = tphys bcm(ip,2) = float(kstar(1)) bcm(ip,3) = mass0(1) bcm(ip,4) = mass(1) bcm(ip,5) = log10(lumin(1)) bcm(ip,6) = log10(rad(1)) teff1 = 1000.d0*((1130.d0*lumin(1)/ & (rad(1)**2.d0))**(1.d0/4.d0)) bcm(ip,7) = log10(teff1) bcm(ip,8) = massc(1) bcm(ip,9) = radc(1) bcm(ip,10) = menv(1) bcm(ip,11) = renv(1) bcm(ip,12) = epoch(1) bcm(ip,13) = ospin(1) bcm(ip,14) = dmt(1) - dmr(1) bcm(ip,15) = rad(1)/rol(1) bcm(ip,16) = float(kstar(2)) bcm(ip,17) = mass0(2) bcm(ip,18) = mass(2) bcm(ip,19) = log10(lumin(2)) bcm(ip,20) = log10(rad(2)) teff2 = 1000.d0*((1130.d0*lumin(2)/ & (rad(2)**2.d0))**(1.d0/4.d0)) bcm(ip,21) = log10(teff2) bcm(ip,22) = massc(2) bcm(ip,23) = radc(2) bcm(ip,24) = menv(2) bcm(ip,25) = renv(2) bcm(ip,26) = epoch(2) bcm(ip,27) = ospin(2) bcm(ip,28) = dmt(2) - dmr(2) bcm(ip,29) = rad(2)/rol(2) bcm(ip,30) = tb bcm(ip,31) = sep bcm(ip,32) = ecc if(isave) tsave = tsave + dtp endif endif * * If not interpolating set the next timestep. * if(intpol.eq.0)then dtm = MAX(1.0d-07*tphys,MIN(dtmi(1),dtmi(2))) dtm = MIN(dtm,tsave-tphys) if(iter.eq.0) dtm0 = dtm endif if(sgl) goto 98 * * Set j1 to the donor - the primary * and j2 to the accretor - the secondary. * if(intpol.eq.0)then if(rad(1)/rol(1).ge.rad(2)/rol(2))then j1 = 1 j2 = 2 else j1 = 2 j2 = 1 endif endif * * Test whether Roche lobe overflow has begun. * if(rad(j1).gt.rol(j1))then * * Interpolate back until the primary is just filling its Roche lobe. * if(rad(j1).ge.1.002d0*rol(j1))then if(intpol.eq.0) tphys00 = tphys intpol = intpol + 1 if(iter.eq.0) goto 7 if(inttry) goto 7 if(intpol.ge.100)then WRITE(99,*)' INTPOL EXCEEDED ',mass1i,mass2i,tbi,ecci goto 140 endif dr = rad(j1) - 1.001d0*rol(j1) if(ABS(rdot(j1)).lt.tiny.or.prec)then goto 7 endif dtm = -dr/ABS(rdot(j1)) if(ABS(tphys0-tphys).gt.tiny) dtm = MAX(dtm,tphys0-tphys) if(kstar(1).ne.kw1)then kstar(1) = kw1 mass0(1) = mass00(1) epoch(1) = tphys - aj0(1) endif if(kstar(2).ne.kw2)then kstar(2) = kw2 mass0(2) = mass00(2) epoch(2) = tphys - aj0(2) endif change = .false. else * * Enter Roche lobe overflow * if(tphys.ge.tphysf) goto 140 goto 7 endif else * * Check if already interpolating. * if(intpol.gt.0)then intpol = intpol + 1 if(intpol.ge.80)then inttry = .true. endif if(ABS(rdot(j1)).lt.tiny)then prec = .true. dtm = 1.0d-07*tphys else dr = rad(j1) - 1.001d0*rol(j1) dtm = -dr/ABS(rdot(j1)) endif if((tphys+dtm).ge.tphys00)then * * If this occurs then most likely the star is a high mass type 4 * where the radius can change very sharply or possibly there is a * discontinuity in the radius as a function of time and HRDIAG * needs to be checked! * dtm = 0.5d0*(tphys00 - tphys0) dtm = MAX(dtm,1.0d-10) prec = .true. endif tphys0 = tphys endif endif * * Check for collision at periastron. * pd = sep*(1.d0 - ecc) if(pd.lt.(rad(1)+rad(2)).and.intpol.eq.0) goto 130 * * Go back for the next step or interpolation. * 98 continue if(tphys.ge.tphysf.and.intpol.eq.0) goto 140 if(change)then change = .false. jp = MIN(80,jp + 1) bpp(jp,1) = tphys bpp(jp,2) = mass(1) bpp(jp,3) = mass(2) bpp(jp,4) = float(kstar(1)) bpp(jp,5) = float(kstar(2)) bpp(jp,6) = sep bpp(jp,7) = ecc bpp(jp,8) = rad(1)/rol(1) bpp(jp,9) = rad(2)/rol(2) bpp(jp,10) = 2.0 endif * iter = iter + 1 * if(iter.ge.loop)then WRITE(99,*)' MAXIMUM ITER EXCEEDED ',mass1i,mass2i,tbi,ecci goto 140 endif goto 5 * * Set the nuclear timescale in years and slow-down factor. * 7 km0 = dtm0*1.0d+03/tb if(km0.lt.tiny) km0 = 0.5d0 * * Force co-rotation of primary and orbit to ensure that the tides do not * lead to unstable Roche (not currently used). * * if(ospin(j1).gt.1.05d0*oorb)then * ospin(j1) = oorb * jspin(j1) = (k2str(j1)*rad(j1)*rad(j1)*(mass(j1)-massc(j1))+ * & k3*radc(j1)*radc(j1)*massc(j1))*ospin(j1) * endif * iter = 0 coel = .false. change = .false. radx(j1) = MAX(radc(j1),rol(j1)) radx(j2) = rad(j2) * jp = MIN(80,jp + 1) bpp(jp,1) = tphys bpp(jp,2) = mass(1) bpp(jp,3) = mass(2) bpp(jp,4) = float(kstar(1)) bpp(jp,5) = float(kstar(2)) bpp(jp,6) = sep bpp(jp,7) = ecc bpp(jp,8) = rad(1)/rol(1) bpp(jp,9) = rad(2)/rol(2) bpp(jp,10) = 3.0 * if(iplot.and.tphys.gt.tiny)then ip = ip + 1 bcm(ip,1) = tphys bcm(ip,2) = float(kstar(1)) bcm(ip,3) = mass0(1) bcm(ip,4) = mass(1) bcm(ip,5) = log10(lumin(1)) bcm(ip,6) = log10(rad(1)) teff1 = 1000.d0*((1130.d0*lumin(1)/ & (rad(1)**2.d0))**(1.d0/4.d0)) bcm(ip,7) = log10(teff1) bcm(ip,8) = massc(1) bcm(ip,9) = radc(1) bcm(ip,10) = menv(1) bcm(ip,11) = renv(1) bcm(ip,12) = epoch(1) bcm(ip,13) = ospin(1) bcm(ip,14) = 0.0 bcm(ip,15) = rad(1)/rol(1) bcm(ip,16) = float(kstar(2)) bcm(ip,17) = mass0(2) bcm(ip,18) = mass(2) bcm(ip,19) = log10(lumin(2)) bcm(ip,20) = log10(rad(2)) teff2 = 1000.d0*((1130.d0*lumin(2)/ & (rad(2)**2.d0))**(1.d0/4.d0)) bcm(ip,21) = log10(teff2) bcm(ip,22) = massc(2) bcm(ip,23) = radc(2) bcm(ip,24) = menv(2) bcm(ip,25) = renv(2) bcm(ip,26) = epoch(2) bcm(ip,27) = ospin(2) bcm(ip,28) = 0.0 bcm(ip,29) = rad(2)/rol(2) bcm(ip,30) = tb bcm(ip,31) = sep bcm(ip,32) = ecc endif * * Eddington limit for accretion on to the secondary in one orbit. * 8 dme = 2.08d-03*eddfac*(1.d0/(1.d0 + zpars(11)))*rad(j2)*tb supedd = .false. novae = .false. disk = .false. * * Determine whether the transferred material forms an accretion * disk around the secondary or hits the secondary in a direct * stream, by using eq.(1) of Ulrich & Burger (1976, ApJ, 206, 509) * fitted to the calculations of Lubow & Shu (1974, ApJ, 198, 383). * * if(kstar(j2).ge.10) disk = .true. rmin = 0.0425d0*sep*(q(j2)*(1.d0+q(j2)))**(1.d0/4.d0) if(rmin.gt.rad(j2)) disk = .true. * * Kelvin-Helmholtz time from the modified classical expression. * do 13 , k = 1,2 tkh(k) = 1.0d+07*mass(k)/(rad(k)*lumin(k)) if(kstar(k).le.1.or.kstar(k).eq.7.or.kstar(k).ge.10)then tkh(k) = tkh(k)*mass(k) else tkh(k) = tkh(k)*(mass(k) - massc(k)) endif 13 continue * * Dynamical timescale for the primary. * tdyn = 5.05d-05*SQRT(rad(j1)**3/mass(j1)) * * Identify special cases. * if(kstar(j1).eq.2)then qc = 4.d0 elseif(kstar(j1).eq.3.or.kstar(j1).eq.5.or.kstar(j1).eq.6)then * qc = (1.67d0-zpars(7)+2.d0*(massc(j1)/mass(j1))**5)/2.13d0 * Alternatively use condition of Hjellming & Webbink, 1987, ApJ, 318, 794. qc = 0.362 + 1.0/(3.0*(1.0 - massc(j1)/mass(j1))) * Or allow all cases to avoid common-envelope. * qc = 100.d0 elseif(kstar(j1).eq.8.or.kstar(j1).eq.9)then qc = 0.784d0 else qc = 3.d0 endif * if(kstar(j1).eq.0.and.q(j1).gt.0.695d0)then * * This will be dynamical mass transfer of a similar nature to * common-envelope evolution. The result is always a single * star placed in *2. * taum = SQRT(tkh(j1)*tdyn) dm1 = mass(j1) if(kstar(j2).le.1)then * * Restrict accretion to thermal timescale of secondary. * dm2 = taum/tkh(j2)*dm1 mass(j2) = mass(j2) + dm2 * * Rejuvenate if the star is still on the main sequence. * mass0(j2) = mass(j2) CALL star(kstar(j2),mass0(j2),mass(j2),tmsnew,tn, & tscls,lums,GB,zpars) * If the star has no convective core then the effective age decreases, * otherwise it will become younger still. if(mass(j2).lt.0.35d0.or.mass(j2).gt.1.25d0)then aj(j2) = tmsnew/tms(j2)*aj(j2)*(mass(j2) - dm2)/mass(j2) else aj(j2) = tmsnew/tms(j2)*aj(j2) endif epoch(j2) = tphys - aj(j2) elseif(kstar(j2).le.6)then * * Add all the material to the giant's envelope. * dm2 = dm1 mass(j2) = mass(j2) + dm2 if(kstar(j2).eq.2)then mass0(j2) = mass(j2) CALL star(kstar(j2),mass0(j2),mass(j2),tmsnew,tn,tscls, & lums,GB,zpars) aj(j2) = tmsnew + tscls(1)*(aj(j2)-tms(j2))/tbgb(j2) epoch(j2) = tphys - aj(j2) endif elseif(kstar(j2).le.12)then * * Form a new giant envelope. * dm2 = dm1 kst = ktype(kstar(j1),kstar(j2)) if(kst.gt.100) kst = kst - 100 if(kst.eq.4)then aj(j2) = aj(j2)/tms(j2) massc(j2) = mass(j2) endif * * Check for planets or low-mass WDs. * if((kstar(j2).eq.10.and.mass(j2).lt.0.05d0).or. & (kstar(j2).ge.11.and.mass(j2).lt.0.5d0))then kst = kstar(j1) mass(j1) = mass(j2) + dm2 mass(j2) = 0.d0 else mass(j2) = mass(j2) + dm2 CALL gntage(massc(j2),mass(j2),kst,zpars, & mass0(j2),aj(j2)) epoch(j2) = tphys - aj(j2) endif kstar(j2) = kst else * * The neutron star or black hole simply accretes at the Eddington rate. * dm2 = MIN(dme*taum/tb,dm1) if(dm2.lt.dm1) supedd = .true. mass(j2) = mass(j2) + dm2 endif coel = .true. if(mass(j2).gt.0.d0)then mass(j1) = 0.d0 kstar(j1) = 15 else kstar(j1) = kstar(j2) kstar(j2) = 15 endif goto 135 elseif(((ABS(ABS(2*kstar(j1)-11)-3).eq.2.or.kstar(j1).eq.9). & and.(q(j1).gt.qc.or.radx(j1).le.radc(j1))).or. & (kstar(j1).eq.2.and.q(j1).gt.qc).or. & (kstar(j1).eq.4.and.q(j1).gt.qc))then * * Common-envelope evolution. * m1ce = mass(j1) m2ce = mass(j2) CALL comenv(mass0(j1),mass(j1),massc(j1),aj(j1),jspin(j1), & kstar(j1),mass0(j2),mass(j2),massc(j2),aj(j2), & jspin(j2),kstar(j2),zpars,ecc,sep,jorb,coel) * jp = MIN(80,jp + 1) bpp(jp,1) = tphys bpp(jp,2) = mass(1) if(kstar(1).eq.15) bpp(jp,2) = mass0(1) bpp(jp,3) = mass(2) if(kstar(2).eq.15) bpp(jp,3) = mass0(2) bpp(jp,4) = float(kstar(1)) bpp(jp,5) = float(kstar(2)) bpp(jp,6) = sep bpp(jp,7) = ecc bpp(jp,8) = rad(1)/rol(1) bpp(jp,9) = rad(2)/rol(2) bpp(jp,10) = 7.0 * epoch(j1) = tphys - aj(j1) if(coel)then com = .true. goto 135 endif epoch(j2) = tphys - aj(j2) if(ecc.gt.1.d0)then if(kstar(1).ge.13)then rc = corerd(kstar(1),mass(1),mass(1),zpars(2)) ospin(1) = jspin(1)/(k3*rc*rc*mass(1)) endif if(kstar(2).ge.13)then rc = corerd(kstar(2),mass(2),mass(2),zpars(2)) ospin(2) = jspin(2)/(k3*rc*rc*mass(2)) endif goto 135 endif * * Next step should be made without changing the time. * dm1 = m1ce - mass(j1) dm2 = mass(j2) - m2ce dm22 = dm2 dtm = 0.d0 * * Reset orbital parameters as separation may have changed. * tb = (sep/aursun)*SQRT(sep/(aursun*(mass(1)+mass(2)))) oorb = twopi/tb elseif(kstar(j1).ge.10.and.kstar(j1).le.12.and. & q(j1).gt.0.628d0)then * * Dynamic transfer from a white dwarf. Secondary will have KW > 9. * taum = SQRT(tkh(j1)*tdyn) dm1 = mass(j1) if(eddfac.lt.10.d0)then dm2 = MIN(dme*taum/tb,dm1) if(dm2.lt.dm1) supedd = .true. else dm2 = dm1 endif mass(j2) = mass(j2) + dm2 * if(kstar(j1).eq.10.and.kstar(j2).eq.10)then * * Assume the energy released by ignition of the triple-alpha reaction * is enough to destroy the star. * kstar(j2) = 15 mass(j2) = 0.d0 elseif(kstar(j1).eq.10.or.kstar(j2).eq.10)then * * Should be helium overflowing onto a CO or ONe core in which case the * helium swells up to form a giant envelope so a HeGB star is formed. * Allowance for the rare case of CO or ONe flowing onto He is made. * kst = 9 if(kstar(j2).eq.10) massc(j2) = dm2 CALL gntage(massc(j2),mass(j2),kst,zpars,mass0(j2),aj(j2)) kstar(j2) = kst epoch(j2) = tphys - aj(j2) elseif(kstar(j2).le.12)then mass0(j2) = mass(j2) if(kstar(j1).eq.12.and.kstar(j2).eq.11)then * * Mixture of ONe and CO will result in an ONe product. * kstar(j2) = 12 endif endif kstar(j1) = 15 mass(j1) = 0.d0 * * Might be a supernova that destroys the system. * if(kstar(j2).le.11.and.mass(j2).gt.mch)then kstar(j2) = 15 mass(j2) = 0.d0 endif coel = .true. goto 135 elseif(kstar(j1).eq.13)then * * Gamma ray burster? * dm1 = mass(j1) mass(j1) = 0.d0 kstar(j1) = 15 dm2 = dm1 mass(j2) = mass(j2) + dm2 kstar(j2) = 14 coel = .true. goto 135 elseif(kstar(j1).eq.14)then * * Both stars are black holes. Let them merge quietly. * dm1 = mass(j1) mass(j1) = 0.d0 kstar(j1) = 15 dm2 = dm1 mass(j2) = mass(j2) + dm2 coel = .true. goto 135 else * * Mass transfer in one Kepler orbit. * dm1 = 3.0d-06*tb*(LOG(rad(j1)/rol(j1))**3)* & MIN(mass(j1),5.d0)**2 if(kstar(j1).eq.2)then mew = (mass(j1) - massc(j1))/mass(j1) dm1 = MAX(mew,0.01d0)*dm1 elseif(kstar(j1).ge.10)then * dm1 = dm1*1.0d+03/MAX(rad(j1),1.0d-04) dm1 = dm1*1.0d+03*mass(j1)/MAX(rad(j1),1.0d-04) endif kst = kstar(j2) * * Possibly mass transfer needs to be reduced if primary is rotating * faster than the orbit (not currently implemented). * * spnfac = MIN(3.d0,MAX(ospin(j1)/oorb,1.d0)) * dm1 = dm1/spnfac**2 * * Limit mass transfer to the thermal rate for remaining giant-like stars * and to the dynamical rate for all others. * if(kstar(j1).ge.2.and.kstar(j1).le.9.and.kstar(j1).ne.7)then *** * JH_temp ... this may be good for HG RLOF?? * if(kstar(j1).eq.2)then * mew = rad(j1)/rol(j1) - 1.d0 * mew = 2.d0*mew * dm1 = dm1*10.d0**mew * endif *** dm1 = MIN(dm1,mass(j1)*tb/tkh(j1)) elseif(rad(j1).gt.10.d0*rol(j1).or.(kstar(j1).le.1.and. & kstar(j2).le.1.and.q(j1).gt.qc))then * * Allow the stars to merge with the product in *1. * m1ce = mass(j1) m2ce = mass(j2) CALL mix(mass0,mass,aj,kstar,zpars) dm1 = m1ce - mass(j1) dm2 = mass(j2) - m2ce * * Next step should be made without changing the time. * dtm = 0.d0 epoch(1) = tphys - aj(1) coel = .true. goto 135 else dm1 = MIN(dm1,mass(j1)*tb/tdyn) endif * * Calculate wind mass loss from the stars during one orbit. * vorb2 = acc1*(mass(1)+mass(2))/sep ivsqm = 1.d0/SQRT(1.d0-ecc*ecc) do 14 , k = 1,2 if(neta.gt.tiny)then rlperi = rol(k)*(1.d0-ecc) dmr(k) = mlwind(kstar(k),lumin(k),radx(k), & mass(k),massc(k),rlperi,z) vwind2 = 2.d0*beta*acc1*mass(k)/radx(k) omv2 = (1.d0 + vorb2/vwind2)**(3.d0/2.d0) dmt(3-k) = ivsqm*acc2*dmr(k)*((acc1*mass(3-k)/vwind2)**2) & /(2.d0*sep*sep*omv2) dmt(3-k) = MIN(dmt(3-k),dmr(k)) else dmr(k) = 0.d0 dmt(3-k) = 0.d0 endif 14 continue * do 15 , k = 1,2 dms(k) = (dmr(k)-dmt(k))*tb 15 continue * * Increase time-scale to relative mass loss of 0.5% but not more than twice. * KM is the number of orbits for the timestep. * km = MIN(2.d0*km0,5.0d-03/ & MAX(ABS(dm1+dms(j1))/mass(j1),dms(j2)/mass(j2))) km0 = km * * Modify time-step & mass loss terms by speed-up factor. * dt = km*tb dtm = dt/1.0d+06 * * Take the stellar evolution timestep into account but don't let it * be overly restrictive for long lived phases. * if(iter.le.1000) dtm = MIN(dtm,dtmi(1),dtmi(2)) dtm = MIN(dtm,tsave-tphys) dt = dtm*1.0d+06 km = dt/tb * * Decide between accreted mass by secondary and/or system mass loss. * taum = mass(j2)/dm1*tb if(kstar(j2).le.2.or.kstar(j2).eq.4)then * * Limit according to the thermal timescale of the secondary. * dm2 = MIN(1.d0,10.d0*taum/tkh(j2))*dm1 elseif(kstar(j2).ge.7.and.kstar(j2).le.9)then * * Naked helium star secondary swells up to a core helium burning star * or SAGB star unless the primary is also a helium star. * if(kstar(j1).ge.7)then dm2 = MIN(1.d0,10.d0*taum/tkh(j2))*dm1 else dm2 = dm1 dmchk = dm2 - 1.05d0*dms(j2) if(dmchk.gt.0.d0.and.dm2/mass(j2).gt.1.0d-04)then kst = MIN(6,2*kstar(j2)-10) if(kst.eq.4)then aj(j2) = aj(j2)/tms(j2) mcx = mass(j2) else mcx = massc(j2) endif mt2 = mass(j2) + km*(dm2 - dms(j2)) CALL gntage(mcx,mt2,kst,zpars,mass0(j2),aj(j2)) epoch(j2) = tphys + dtm - aj(j2) * jp = MIN(80,jp + 1) bpp(jp,1) = tphys bpp(jp,2) = mass(j1) bpp(jp,3) = mt2 bpp(jp,4) = float(kstar(j1)) bpp(jp,5) = float(kst) bpp(jp,6) = sep bpp(jp,7) = ecc bpp(jp,8) = rad(1)/rol(1) bpp(jp,9) = rad(2)/rol(2) bpp(jp,10) = 8.0 if(j1.eq.2)then bpp(jp,2) = mt2 bpp(jp,3) = mass(j1) bpp(jp,4) = float(kst) bpp(jp,5) = float(kstar(j1)) endif * endif endif elseif(kstar(j1).le.6.and. & (kstar(j2).ge.10.and.kstar(j2).le.12))then * * White dwarf secondary. * if(dm1/tb.lt.2.71d-07)then if(dm1/tb.lt.1.03d-07)then * * Accrete until a nova explosion blows away most of the accreted material. * novae = .true. dm2 = MIN(dm1,dme) if(dm2.lt.dm1) supedd = .true. dm22 = epsnov*dm2 else * * Steady burning at the surface * dm2 = dm1 endif else * * Make a new giant envelope. * dm2 = dm1 * * Check for planets or low-mass WDs. * if((kstar(j2).eq.10.and.mass(j2).lt.0.05d0).or. & (kstar(j2).ge.11.and.mass(j2).lt.0.5d0))then kst = kstar(j2) else kst = MIN(6,3*kstar(j2)-27) mt2 = mass(j2) + km*(dm2 - dms(j2)) CALL gntage(massc(j2),mt2,kst,zpars,mass0(j2),aj(j2)) epoch(j2) = tphys + dtm - aj(j2) * jp = MIN(80,jp + 1) bpp(jp,1) = tphys bpp(jp,2) = mass(j1) bpp(jp,3) = mt2 bpp(jp,4) = float(kstar(j1)) bpp(jp,5) = float(kst) bpp(jp,6) = sep bpp(jp,7) = ecc bpp(jp,8) = rad(1)/rol(1) bpp(jp,9) = rad(2)/rol(2) bpp(jp,10) = 8.0 if(j1.eq.2)then bpp(jp,2) = mt2 bpp(jp,3) = mass(j1) bpp(jp,4) = float(kst) bpp(jp,5) = float(kstar(j1)) endif * endif * endif elseif(kstar(j2).ge.10)then * * Impose the Eddington limit. * dm2 = MIN(dm1,dme) if(dm2.lt.dm1) supedd = .true. * else * * We have a giant whose envelope can absorb any transferred material. * dm2 = dm1 endif if(.not.novae) dm22 = dm2 * if(kst.ge.10.and.kst.le.12)then mt2 = mass(j2) + km*(dm22 - dms(j2)) if(kstar(j1).le.10.and.kst.eq.10.and.mt2.ge.0.7d0)then * * HeWD can only accrete helium-rich material up to a mass of 0.7 when * it is destroyed in a possible Type 1a SN. * mass(j1) = mass(j1) - km*(dm1 + dms(j1)) mass(j2) = 0.d0 kstar(j2) = 15 goto 135 elseif(kstar(j1).le.10.and.kst.ge.11)then * * CO and ONeWDs accrete helium-rich material until the accumulated * material exceeds a mass of 0.15 when it ignites. For a COWD with * mass less than 0.95 the system will be destroyed as an ELD in a * possible Type 1a SN. COWDs with mass greater than 0.95 and ONeWDs * will survive with all the material converted to ONe (JH 30/09/99). * ** Now changed to an ELD for all COWDs when 0.15 accreted (JH 11/01/00). * if((mt2-mass0(j2)).ge.0.15d0)then if(kst.eq.11)then mass(j1) = mass(j1) - km*(dm1 + dms(j1)) mass(j2) = 0.d0 kstar(j2) = 15 goto 135 endif mass0(j2) = mt2 endif else mass0(j2) = mt2 endif * * If the Chandrasekhar limit is exceeded for a white dwarf then destroy * the white dwarf in a supernova. If the WD is ONe then a neutron star * will survive the supernova and we let HRDIAG take care of this when * the stars are next updated. * if(kst.eq.10.or.kst.eq.11)then if(mt2.ge.mch)then dm1 = mch - mass(j2) + km*dms(j2) mass(j1) = mass(j1) - dm1 - km*dms(j1) mass(j2) = 0.d0 kstar(j2) = 15 goto 135 endif endif endif * * Modify mass loss terms by speed-up factor. * dm1 = km*dm1 dm2 = km*dm2 dm22 = km*dm22 dme = km*dme * * Calculate orbital angular momentum change due to system mass loss. * djorb = ((dmr(1)+q(1)*dmt(1))*mass(2)*mass(2) + & (dmr(2)+q(2)*dmt(2))*mass(1)*mass(1))/ & (mass(1)+mass(2))**2 djorb = djorb*dt * * For super-Eddington mass transfer rates, for gamma = -2.0, * and for novae systems, assume that material is lost from * the system as if a wind from the secondary. * If gamma = -1.0 then assume the lost material carries with it * the specific angular momentum of the primary and for all * gamma > 0.0 assume that it takes away a fraction gamma of * the orbital angular momentum. * if(supedd.or.novae.or.gamma.lt.-1.5d0)then djorb = djorb + (dm1 - dm22)*mass(j1)*mass(j1)/ & (mass(1)+mass(2))**2 elseif(gamma.ge.0.d0)then djorb = djorb + gamma*(dm1 - dm2) else djorb = djorb + (dm1 - dm2)*mass(j2)*mass(j2)/ & (mass(1)+mass(2))**2 endif * ecc2 = ecc*ecc omecc2 = 1.d0 - ecc2 sqome2 = SQRT(omecc2) * djorb = djorb*sep*sep*sqome2*oorb delet = 0.d0 * * For very close systems include angular momentum loss mechanisms. * if(sep.le.10.d0)then djgr = 8.315d-10*mass(1)*mass(2)*(mass(1)+mass(2))/ & (sep*sep*sep*sep) f1 = (19.d0/6.d0) + (121.d0/96.d0)*ecc2 sqome5 = sqome2**5 delet1 = djgr*ecc*f1/sqome5 djgr = djgr*jorb*(1.d0+0.875d0*ecc2)/sqome5 djorb = djorb + djgr*dt delet = delet + delet1*dt endif * do 602 , k = 1,2 * dms(k) = km*dms(k) if(kstar(k).lt.10) dms(k) = MIN(dms(k),mass(k) - massc(k)) * * Calculate change in the intrinsic spin of the star. * djspint(k) = (2.d0/3.d0)*(dmr(k)*radx(k)*radx(k)*ospin(k) - & xi*dmt(k)*radx(3-k)*radx(3-k)*ospin(3-k)) djspint(k) = djspint(k)*dt * if(mass(k).gt.0.35d0.and.kstar(k).lt.10)then djmb = 5.83d-16*menv(k)*(rad(k)*ospin(k))**3/mass(k) djspint(k) = djspint(k) + djmb*dt endif * 602 continue * * Adjust the spin angular momentum of each star owing to mass transfer * and conserve total angular momentum. * djt = dm1*radx(j1)*radx(j1)*ospin(j1) djspint(j1) = djspint(j1) + djt djorb = djorb - djt if(disk)then * * Alter spin of the degenerate secondary by assuming that material * falls onto the star from the inner edge of a Keplerian accretion * disk and that the system is in a steady state. * djt = dm2*twopi*aursun*SQRT(aursun*mass(j2)*radx(j2)) djspint(j2) = djspint(j2) - djt djorb = djorb + djt * else * * No accretion disk. * Calculate the angular momentum of the transferred material by * using the radius of the disk (see Ulrich & Burger) that would * have formed if allowed. * rdisk = 1.7d0*rmin djt = dm2*twopi*aursun*SQRT(aursun*mass(j2)*rdisk) djspint(j2) = djspint(j2) - djt djorb = djorb + djt * endif djtx(2) = djt * * Adjust the secondary spin if a nova eruption has occurred. * if(novae)then djt = (dm2 - dm22)*radx(j2)*radx(j2)*ospin(j2) djspint(j2) = djspint(j2) + djt djtx(2) = djtx(2) - djt endif * * Calculate circularization, orbital shrinkage and spin up. * do 603 , k = 1,2 * dspint(k) = 0.d0 if(((kstar(k).le.9.and.rad(k).ge.0.01d0*rol(k)).or. & (kstar(k).ge.10.and.k.eq.j1)).and.tflag.gt.0)then * raa2 = (radx(k)/sep)**2 raa6 = raa2**3 * f5 = 1.d0+ecc2*(3.d0+ecc2*0.375d0) f4 = 1.d0+ecc2*(1.5d0+ecc2*0.125d0) f3 = 1.d0+ecc2*(3.75d0+ecc2*(1.875d0+ecc2*7.8125d-02)) f2 = 1.d0+ecc2*(7.5d0+ecc2*(5.625d0+ecc2*0.3125d0)) f1 = 1.d0+ecc2*(15.5d0+ecc2*(31.875d0+ecc2*(11.5625d0 & +ecc2*0.390625d0))) * if((kstar(k).eq.1.and.mass(k).ge.1.25d0).or. & kstar(k).eq.4.or.kstar(k).eq.7)then tc = 1.592d-09*(mass(k)**2.84d0) f = 1.9782d+04*SQRT((mass(k)*radx(k)*radx(k))/sep**5)* & tc*(1.d0+q(3-k))**(5.d0/6.d0) tcqr = f*q(3-k)*raa6 rg2 = k2str(k) elseif(kstar(k).le.9)then renv(k) = MIN(renv(k),radx(k)-radc(k)) renv(k) = MAX(renv(k),1.0d-10) tc = mr23yr*(menv(k)*renv(k)*(radx(k)-0.5d0*renv(k))/ & (3.d0*lumin(k)))**(1.d0/3.d0) ttid = twopi/(1.0d-10 + ABS(oorb - ospin(k))) f = MIN(1.d0,(ttid/(2.d0*tc)**2)) tcqr = 2.d0*f*q(3-k)*raa6*menv(k)/(21.d0*tc*mass(k)) rg2 = (k2str(k)*(mass(k)-massc(k)))/mass(k) else f = 7.33d-09*(lumin(k)/mass(k))**(5.d0/7.d0) tcqr = f*q(3-k)*q(3-k)*raa2*raa2/(1.d0+q(3-k)) rg2 = k3 endif sqome3 = sqome2**3 delet1 = 27.d0*tcqr*(1.d0+q(3-k))*raa2*(ecc/sqome2**13)* & (f3 - (11.d0/18.d0)*sqome3*f4*ospin(k)/oorb) tcirc = ecc/(ABS(delet1) + 1.0d-20) delet = delet + delet1*dt dspint(k) = (3.d0*q(3-k)*tcqr/(rg2*omecc2**6))* & (f2*oorb - sqome3*f5*ospin(k)) eqspin = oorb*f2/(sqome3*f5) if(dt.gt.0.d0)then if(dspint(k).ge.0.d0)then dspint(k) = MIN(dt*dspint(k),eqspin-ospin(k))/dt else dspint(k) = MAX(dt*dspint(k),eqspin-ospin(k))/dt endif endif djt = (k2str(k)*(mass(k)-massc(k))*radx(k)*radx(k) + & k3*massc(k)*radc(k)*radc(k))*dspint(k) djorb = djorb + djt*dt djspint(k) = djspint(k) - djt*dt * endif * jspin(k) = MAX(1.0d-10,jspin(k) - djspint(k)) * * Ensure that the star does not spin up beyond break-up, and transfer * the excess angular momentum back to the orbit. * ospbru = twopi*SQRT(mass(k)*aursun**3/radx(k)**3) jspbru = (k2str(k)*(mass(k)-massc(k))*radx(k)*radx(k) + & k3*massc(k)*radc(k)*radc(k))*ospbru if(jspin(k).gt.jspbru)then mew = 1.d0 if(djtx(2).gt.0.d0)then mew = MIN(mew,(jspin(k) - jspbru)/djtx(2)) endif djorb = djorb - (jspin(k) - jspbru) jspin(k) = jspbru * If excess material should not be accreted, activate next line. * dm22 = (1.d0 - mew)*dm22 endif * 603 continue * * Update the masses. * kstar(j2) = kst mass(j1) = mass(j1) - dm1 - dms(j1) if(kstar(j1).le.1.or.kstar(j1).eq.7) mass0(j1) = mass(j1) mass(j2) = mass(j2) + dm22 - dms(j2) if(kstar(j2).le.1.or.kstar(j2).eq.7) mass0(j2) = mass(j2) * * For a HG star check if the initial mass can be reduced. * if(kstar(j1).eq.2.and.mass0(j1).le.zpars(3))then m0 = mass0(j1) mass0(j1) = mass(j1) CALL star(kstar(j1),mass0(j1),mass(j1),tmsnew,tn,tscls, & lums,GB,zpars) if(GB(9).lt.massc(j1))then mass0(j1) = m0 endif endif if(kstar(j2).eq.2.and.mass0(j2).le.zpars(3))then m0 = mass0(j2) mass0(j2) = mass(j2) CALL star(kstar(j2),mass0(j2),mass(j2),tmsnew,tn,tscls, & lums,GB,zpars) if(GB(9).lt.massc(j2))then mass0(j2) = m0 endif endif * ecc = ecc - delet ecc = MAX(ecc,0.d0) if(ecc.lt.1.0d-10) ecc = 0.d0 * if(ecc.ge.1.d0) goto 135 * * Ensure that Jorb does not become negative which could happen if the * primary overfills its Roche lobe initially. In this case we simply * allow contact to occur. * jorb = MAX(1.d0,jorb - djorb) sep = (mass(1) + mass(2))*jorb*jorb/ & ((mass(1)*mass(2)*twopi)**2*aursun**3*(1.d0-ecc*ecc)) tb = (sep/aursun)*SQRT(sep/(aursun*(mass(1)+mass(2)))) oorb = twopi/tb * endif * * Always rejuvenate the secondary and age the primary if they are on * the main sequence. * if(kstar(j1).le.2.or.kstar(j1).eq.7)then CALL star(kstar(j1),mass0(j1),mass(j1),tmsnew,tn,tscls, & lums,GB,zpars) if(kstar(j1).eq.2)then aj(j1) = tmsnew + (tscls(1) - tmsnew)*(aj(j1)-tms(j1))/ & (tbgb(j1) - tms(j1)) else aj(j1) = tmsnew/tms(j1)*aj(j1) endif epoch(j1) = tphys - aj(j1) endif * if(kstar(j2).le.2.or.kstar(j2).eq.7)then CALL star(kstar(j2),mass0(j2),mass(j2),tmsnew,tn,tscls, & lums,GB,zpars) if(kstar(j2).eq.2)then aj(j2) = tmsnew + (tscls(1) - tmsnew)*(aj(j2)-tms(j2))/ & (tbgb(j2) - tms(j2)) elseif((mass(j2).lt.0.35d0.or.mass(j2).gt.1.25d0). & and.kstar(j2).ne.7)then aj(j2) = tmsnew/tms(j2)*aj(j2)*(mass(j2) - dm22)/mass(j2) else aj(j2) = tmsnew/tms(j2)*aj(j2) endif epoch(j2) = tphys - aj(j2) endif * * Obtain the stellar parameters for the next step. * tphys = tphys + dtm do 90 , k = 1,2 age = tphys - epoch(k) m0 = mass0(k) mt = mass(k) mc = massc(k) * * Masses over 100Msun should probably not be trusted in the * evolution formulae. * if(mt.gt.100.d0)then WRITE(99,*)' MASS EXCEEDED ',mass1i,mass2i,tbi,ecci,mt goto 140 endif kw = kstar(k) CALL star(kw,m0,mt,tm,tn,tscls,lums,GB,zpars) CALL hrdiag(m0,age,mt,tm,tn,tscls,lums,GB,zpars, & rm,lum,kw,mc,rc,me,re,k2) * * Check for a supernova and correct the semi-major axis if so. * if(kw.ne.kstar(k).and.kstar(k).le.12.and. & (kw.eq.13.or.kw.eq.14))then dms(k) = mass(k) - mt CALL kick(kw,mass(k),mt,mass(3-k),ecc,sep,jorb,vs) vkick(k) = dsqrt(vs(1)*vs(1)+vs(2)*vs(2)+vs(3)*vs(3)) if(ecc.gt.1.d0)then kstar(k) = kw mass(k) = mt epoch(k) = tphys - age goto 135 endif tb = (sep/aursun)*SQRT(sep/(aursun*(mt+mass(3-k)))) oorb = twopi/tb endif if(kw.ne.kstar(k))then change = .true. if((kw.eq.13.or.kw.eq.14).and.kstar(k).le.12)then snova = .true. endif mass(k) = mt if(kw.eq.15)then kstar(k) = kw goto 135 endif mass0(k) = m0 epoch(k) = tphys - age endif * * Determine stellar evolution timescale for nuclear burning types. * if(kw.le.9)then CALL deltat(kw,age,tm,tn,tscls,dt,dtr) dtmi(k) = MIN(dt,dtr) * dtmi(k) = dtr dtmi(k) = MAX(1.0d-07,dtmi(k)) else dtmi(k) = 1.0d+10 endif * dtmi(k) = MAX((tn-age),1.0d-07) * * Save relevent solar quantities. * aj(k) = age kstar(k) = kw rad(k) = rm radx(k) = rm lumin(k) = lum massc(k) = mc radc(k) = rc menv(k) = me renv(k) = re k2str(k) = k2 tms(k) = tm tbgb(k) = tscls(1) * * Check for blue straggler formation. * if(kw.le.1.and.tm.lt.tphys.and..not.bss)then bss = .true. jp = MIN(80,jp + 1) bpp(jp,1) = tphys bpp(jp,2) = mass(1) bpp(jp,3) = mass(2) bpp(jp,4) = float(kstar(1)) bpp(jp,5) = float(kstar(2)) bpp(jp,6) = sep bpp(jp,7) = ecc bpp(jp,8) = rad(1)/rol(1) bpp(jp,9) = rad(2)/rol(2) bpp(jp,10) = 14.0 endif * 90 continue * do 100 , k = 1,2 q(k) = mass(k)/mass(3-k) rol(k) = rl(q(k))*sep 100 continue if(rad(j1).gt.rol(j1)) radx(j1) = MAX(radc(j1),rol(j1)) do 110 , k = 1,2 ospin(k) = jspin(k)/(k2str(k)*(mass(k)-massc(k))*radx(k)* & radx(k) + k3*massc(k)*radc(k)*radc(k)) 110 continue * if((isave.and.tphys.ge.tsave).or.iplot)then ip = ip + 1 bcm(ip,1) = tphys bcm(ip,2) = float(kstar(1)) bcm(ip,3) = mass0(1) bcm(ip,4) = mass(1) bcm(ip,5) = log10(lumin(1)) bcm(ip,6) = log10(rad(1)) teff1 = 1000.d0*((1130.d0*lumin(1)/ & (rad(1)**2.d0))**(1.d0/4.d0)) bcm(ip,7) = log10(teff1) bcm(ip,8) = massc(1) bcm(ip,9) = radc(1) bcm(ip,10) = menv(1) bcm(ip,11) = renv(1) bcm(ip,12) = epoch(1) bcm(ip,13) = ospin(1) bcm(ip,15) = rad(1)/rol(1) bcm(ip,16) = float(kstar(2)) bcm(ip,17) = mass0(2) bcm(ip,18) = mass(2) bcm(ip,19) = log10(lumin(2)) bcm(ip,20) = log10(rad(2)) teff2 = 1000.d0*((1130.d0*lumin(2)/ & (rad(2)**2.d0))**(1.d0/4.d0)) bcm(ip,21) = log10(teff2) bcm(ip,22) = massc(2) bcm(ip,23) = radc(2) bcm(ip,24) = menv(2) bcm(ip,25) = renv(2) bcm(ip,26) = epoch(2) bcm(ip,27) = ospin(2) bcm(ip,29) = rad(2)/rol(2) bcm(ip,30) = tb bcm(ip,31) = sep bcm(ip,32) = ecc dt = MAX(dtm,1.0d-12)*1.0d+06 if(j1.eq.1)then bcm(ip,14) = (-1.0*dm1 - dms(1))/dt bcm(ip,28) = (dm2 - dms(2))/dt else bcm(ip,14) = (dm2 - dms(1))/dt bcm(ip,28) = (-1.0*dm1 - dms(2))/dt endif if(isave) tsave = tsave + dtp endif * if(tphys.ge.tphysf) goto 140 * if(change)then change = .false. jp = MIN(80,jp + 1) bpp(jp,1) = tphys bpp(jp,2) = mass(1) bpp(jp,3) = mass(2) bpp(jp,4) = float(kstar(1)) bpp(jp,5) = float(kstar(2)) bpp(jp,6) = sep bpp(jp,7) = ecc bpp(jp,8) = rad(1)/rol(1) bpp(jp,9) = rad(2)/rol(2) bpp(jp,10) = 2.0 endif * * Test whether the primary still fills its Roche lobe. * if(rad(j1).gt.rol(j1).and..not.snova)then * * Test for a contact system * if(rad(j2).gt.rol(j2)) goto 130 iter = iter + 1 goto 8 else jp = MIN(80,jp + 1) bpp(jp,1) = tphys bpp(jp,2) = mass(1) bpp(jp,3) = mass(2) bpp(jp,4) = float(kstar(1)) bpp(jp,5) = float(kstar(2)) bpp(jp,6) = sep bpp(jp,7) = ecc bpp(jp,8) = rad(1)/rol(1) bpp(jp,9) = rad(2)/rol(2) bpp(jp,10) = 4.0 dtm = 0.d0 goto 4 endif * 130 continue * * Contact system. * coel = .true. * * If *1 or *2 is giant-like this will be common-envelope evolution. * m1ce = mass(j1) m2ce = mass(j2) rrl1 = MIN(999.999d0,rad(1)/rol(1)) rrl2 = MIN(999.999d0,rad(2)/rol(2)) * jp = MIN(80,jp + 1) bpp(jp,1) = tphys bpp(jp,2) = mass(1) bpp(jp,3) = mass(2) bpp(jp,4) = float(kstar(1)) bpp(jp,5) = float(kstar(2)) bpp(jp,6) = sep bpp(jp,7) = ecc bpp(jp,8) = rrl1 bpp(jp,9) = rrl2 bpp(jp,10) = 5.0 * if(kstar(j1).ge.2.and.kstar(j1).le.9.and.kstar(j1).ne.7)then CALL comenv(mass0(j1),mass(j1),massc(j1),aj(j1),jspin(j1), & kstar(j1),mass0(j2),mass(j2),massc(j2),aj(j2), & jspin(j2),kstar(j2),zpars,ecc,sep,jorb,coel) com = .true. elseif(kstar(j2).ge.2.and.kstar(j2).le.9.and.kstar(j2).ne.7)then CALL comenv(mass0(j2),mass(j2),massc(j2),aj(j2),jspin(j2), & kstar(j2),mass0(j1),mass(j1),massc(j1),aj(j1), & jspin(j1),kstar(j1),zpars,ecc,sep,jorb,coel) com = .true. else CALL mix(mass0,mass,aj,kstar,zpars) endif if(com)then jp = MIN(80,jp + 1) bpp(jp,1) = tphys bpp(jp,2) = mass(1) if(kstar(1).eq.15) bpp(jp,2) = mass0(1) bpp(jp,3) = mass(2) if(kstar(2).eq.15) bpp(jp,3) = mass0(2) bpp(jp,4) = float(kstar(1)) bpp(jp,5) = float(kstar(2)) bpp(jp,6) = sep bpp(jp,7) = ecc rrl1 = MIN(rrl1,0.99d0) rrl2 = MIN(rrl2,0.99d0) bpp(jp,8) = rrl1 bpp(jp,9) = rrl2 bpp(jp,10) = 7.0 endif epoch(1) = tphys - aj(1) epoch(2) = tphys - aj(2) if(.not.coel)then * * Next step should be made without changing the time. * if(ecc.gt.1.d0)then if(kstar(1).ge.13)then rc = corerd(kstar(1),mass(1),mass(1),zpars(2)) ospin(1) = jspin(1)/(k3*rc*rc*mass(1)) endif if(kstar(2).ge.13)then rc = corerd(kstar(2),mass(2),mass(2),zpars(2)) ospin(2) = jspin(2)/(k3*rc*rc*mass(2)) endif goto 135 endif dtm = 0.d0 * * Reset orbital parameters as separation may have changed. * tb = (sep/aursun)*SQRT(sep/(aursun*(mass(1)+mass(2)))) oorb = twopi/tb goto 4 endif * 135 continue * sgl = .true. if(kstar(1).ne.15.or.kstar(2).ne.15)then if(com)then com = .false. else jp = MIN(80,jp + 1) bpp(jp,1) = tphys bpp(jp,2) = mass(1) if(kstar(1).eq.15) bpp(jp,2) = mass0(1) bpp(jp,3) = mass(2) if(kstar(2).eq.15) bpp(jp,3) = mass0(2) bpp(jp,4) = float(kstar(1)) bpp(jp,5) = float(kstar(2)) bpp(jp,6) = zero bpp(jp,7) = zero bpp(jp,8) = zero bpp(jp,9) = ngtv if(coel)then bpp(jp,10) = 6.0 elseif(ecc.gt.1.d0)then * * Binary dissolved by a supernova or tides. * bpp(jp,6) = sep bpp(jp,7) = ecc bpp(jp,9) = ngtv2 bpp(jp,10) = 11.0 else bpp(jp,10) = 9.0 endif endif if(kstar(2).eq.15)then kmax = 1 rol(2) = -1.d0*rad(2) dtmi(2) = tphysf elseif(kstar(1).eq.15)then kmin = 2 rol(1) = -1.d0*rad(1) dtmi(1) = tphysf endif ecc = -1.d0 sep = 0.d0 dtm = 0.d0 coel = .false. goto 4 endif * 140 continue * if(com)then com = .false. else jp = MIN(80,jp + 1) bpp(jp,1) = tphys bpp(jp,2) = mass(1) if(kstar(1).eq.15.and.bpp(jp-1,4).lt.15.0)then bpp(jp,2) = mass0(1) endif bpp(jp,3) = mass(2) if(kstar(2).eq.15.and.bpp(jp-1,5).lt.15.0)then bpp(jp,3) = mass0(2) endif bpp(jp,4) = float(kstar(1)) bpp(jp,5) = float(kstar(2)) bpp(jp,6) = zero bpp(jp,7) = zero bpp(jp,8) = zero if(coel)then bpp(jp,9) = ngtv bpp(jp,10) = 6.0 elseif(kstar(1).eq.15.and.kstar(2).eq.15)then * * Cases of accretion induced supernova or single star supernova. * No remnant is left in either case. * bpp(jp,9) = ngtv2 bpp(jp,10) = 9.0 else bpp(jp,6) = sep bpp(jp,7) = ecc bpp(jp,8) = rad(1)/rol(1) bpp(jp,9) = rad(2)/rol(2) bpp(jp,10) = 10.0 endif endif * if((isave.and.tphys.ge.tsave).or.iplot)then ip = ip + 1 bcm(ip,1) = tphys bcm(ip,2) = float(kstar(1)) bcm(ip,3) = mass0(1) bcm(ip,4) = mass(1) bcm(ip,5) = log10(lumin(1)) bcm(ip,6) = log10(rad(1)) teff1 = 1000.d0*((1130.d0*lumin(1)/ & (rad(1)**2.d0))**(1.d0/4.d0)) bcm(ip,7) = log10(teff1) bcm(ip,8) = massc(1) bcm(ip,9) = radc(1) bcm(ip,10) = menv(1) bcm(ip,11) = renv(1) bcm(ip,12) = epoch(1) bcm(ip,13) = ospin(1) bcm(ip,15) = rad(1)/rol(1) bcm(ip,16) = float(kstar(2)) bcm(ip,17) = mass0(2) bcm(ip,18) = mass(2) bcm(ip,19) = log10(lumin(2)) bcm(ip,20) = log10(rad(2)) teff2 = 1000.d0*((1130.d0*lumin(2)/ & (rad(2)**2.d0))**(1.d0/4.d0)) bcm(ip,21) = log10(teff2) bcm(ip,22) = massc(2) bcm(ip,23) = radc(2) bcm(ip,24) = menv(2) bcm(ip,25) = renv(2) bcm(ip,26) = epoch(2) bcm(ip,27) = ospin(2) bcm(ip,29) = rad(2)/rol(2) bcm(ip,30) = tb bcm(ip,31) = sep bcm(ip,32) = ecc dt = MAX(dtm,1.0d-12)*1.0d+06 if(j1.eq.1)then bcm(ip,14) = (-1.0*dm1 - dms(1))/dt bcm(ip,28) = (dm2 - dms(2))/dt else bcm(ip,14) = (dm2 - dms(1))/dt bcm(ip,28) = (-1.0*dm1 - dms(2))/dt endif if(isave) tsave = tsave + dtp if(tphysf.le.0.d0)then ip = ip + 1 do 145 , k = 1,32 bcm(ip,k) = bcm(ip-1,k) 145 continue endif endif * tphysf = tphys if(sgl)then if(ecc.ge.0.d0.and.ecc.le.1.d0) ecc = -1.d0 tb = -1.d0 endif tb = tb*yeardy if(jp.ge.80)then WRITE(99,*)' EVOLV2 ARRAY ERROR ',mass1i,mass2i,tbi,ecci WRITE(*,*)' STOP: EVOLV2 ARRAY ERROR ' CALL exit(0) STOP elseif(jp.ge.40)then WRITE(99,*)' EVOLV2 ARRAY WARNING ',mass1i,mass2i,tbi,ecci,jp endif bcm(ip+1,1) = -1.0 bpp(jp+1,1) = -1.0 * RETURN END *** cc//gntage.f *** SUBROUTINE gntage(mc,mt,kw,zpars,m0,aj) * * A routine to determine the age of a giant from its core mass and type. * * Author : C. A. Tout * Date : 24th September 1996 * Revised: 21st February 1997 to include core-helium-burning stars * * Rewritten: 2nd January 1998 by J. R. Hurley to be compatible with * the new evolution routines and to include new stellar * types. * implicit none * integer kw integer j,jmax parameter(jmax=30) * real*8 mc,mt,m0,aj,tm,tn real*8 tscls(20),lums(10),GB(10),zpars(20) real*8 mmin,mmax,mmid,dm,f,fmid,dell,derl,lum real*8 macc,lacc,tiny parameter(macc=0.00001d0,lacc=0.0001d0,tiny=1.0d-14) real*8 mcx,mcy * real*8 mcheif,mcagbf,mheif,mbagbf,mcgbf,lmcgbf,lbgbf,lbgbdf external mcheif,mcagbf,mheif,mbagbf,mcgbf,lmcgbf,lbgbf,lbgbdf * * This should only be entered with KW = 3, 4, 5, 6 or 9 * * First we check that we don't have a CheB star * with too small a core mass. if(kw.eq.4)then * Set the minimum CHeB core mass using M = Mflash mcy = mcheif(zpars(2),zpars(2),zpars(10)) if(mc.le.mcy) kw = 3 * if(mc.le.mcy) WRITE(66,*)' GNTAGE4: changed to 3' endif * * Next we check that we don't have a GB star for M => Mfgb if(kw.eq.3)then * Set the maximum GB core mass using M = Mfgb mcy = mcheif(zpars(3),zpars(2),zpars(9)) if(mc.ge.mcy)then kw = 4 aj = 0.d0 * WRITE(66,*)' GNTAGE3: changed to 4' endif endif * if(kw.eq.6)then * * We try to start the star from the start of the SAGB by * setting Mc = Mc,TP. * mcy = 0.44d0*2.25d0 + 0.448d0 if(mc.gt.mcy)then * A type 6 with this sized core mass cannot exist as it should * already have become a NS or BH as a type 5. * We set it up so that it will. mcx = (mc + 0.35d0)/0.773d0 elseif(mc.ge.0.8d0)then mcx = (mc - 0.448d0)/0.44d0 else mcx = mc endif m0 = mbagbf(mcx) if(m0.lt.tiny)then * Carbon core mass is less then the minimum for the start of SAGB. * This must be the case of a low-mass C/O or O/Ne WD with only a * very small envelope added or possibly the merger of a helium star * with a main sequence star. We will set m0 = mt and then reset the * core mass to allow for some helium to be added to the C/O core. kw = 14 * WRITE(66,*)' GNTAGE6: changed to 4' else CALL star(kw,m0,mt,tm,tn,tscls,lums,GB,zpars) aj = tscls(13) endif endif * if(kw.eq.5)then * * We fit a Helium core mass at the base of the AGB. * m0 = mbagbf(mc) if(m0.lt.tiny)then * Helium core mass is less then the BAGB minimum. kw = 14 * WRITE(66,*)' GNTAGE5: changed to 4' else CALL star(kw,m0,mt,tm,tn,tscls,lums,GB,zpars) aj = tscls(2) + tscls(3) endif endif * * if(kw.eq.4)then * * The supplied age is actually the fractional age, fage, of CHeB lifetime * that has been completed, ie. 0 <= aj <= 1. * if(aj.lt.0.d0.or.aj.gt.1.d0)then * WRITE(99,*)' FATAL ERROR! GNTAGE4: fage out of bounds ' * WRITE(99,*)' FAGE ',aj * WRITE(*,*)' STOP: FATAL ERROR ' * CALL exit(0) * STOP aj = 0.d0 endif * Get the minimum, fage=1, and maximum, fage=0, allowable masses mcy = mcagbf(zpars(2)) if(mc.ge.mcy)then mmin = mbagbf(mc) else mmin = zpars(2) endif mmax = mheif(mc,zpars(2),zpars(10)) if(aj.lt.tiny)then m0 = mmax goto 20 elseif(aj.ge.1.d0)then m0 = mmin goto 20 endif * Use the bisection method to find m0 fmid = (1.d0-aj)*mcheif(mmax,zpars(2),zpars(10)) + & aj*mcagbf(mmax) - mc f = (1.d0-aj)*mcheif(mmin,zpars(2),zpars(10)) + & aj*mcagbf(mmin) - mc if(f*fmid.ge.0.d0)then * This will probably occur if mc is just greater than the minimum * allowed mass for a CHeB star and fage > 0. kw = 3 * WRITE(66,*)' GNTAGE4: changed to 3' goto 90 endif m0 = mmin dm = mmax - mmin do 10 , j = 1,jmax dm = 0.5d0*dm mmid = m0 + dm fmid = (1.d0-aj)*mcheif(mmid,zpars(2),zpars(10)) + & aj*mcagbf(mmid) - mc if(fmid.lt.0.d0) m0 = mmid if(ABS(dm).lt.macc.or.ABS(fmid).lt.tiny) goto 20 if(j.eq.jmax)then * WRITE(99,*)' FATAL ERROR! GNTAGE4: root not found ' * WRITE(*,*)' STOP: FATAL ERROR ' * CALL exit(0) * STOP m0 = mt aj = 0.d0 endif 10 continue 20 continue * CALL star(kw,m0,mt,tm,tn,tscls,lums,GB,zpars) aj = tscls(2) + aj*tscls(3) * endif * 90 continue * if(kw.eq.3)then * * First we double check that we don't have a GB star for M => Mfgb mcy = mcheif(zpars(3),zpars(2),zpars(9)) if(mc.ge.mcy)then * WRITE(99,*)' GNTAGE3: star too big for GB ' * WRITE(*,*)' STOP: FATAL ERROR ' * CALL exit(0) * STOP mc = 0.99d0*mcy endif * Next we find an m0 so as to place the star at the BGB mcx = mcheif(zpars(2),zpars(2),zpars(9)) if(mc.gt.mcx)then m0 = mheif(mc,zpars(2),zpars(9)) else * Use Newton-Raphson to find m0 from Lbgb m0 = zpars(2) CALL star(kw,m0,mt,tm,tn,tscls,lums,GB,zpars) lum = lmcgbf(mc,GB) j = 0 30 continue dell = lbgbf(m0) - lum if(ABS(dell/lum).le.lacc) goto 40 derl = lbgbdf(m0) m0 = m0 - dell/derl j = j + 1 if(j.eq.jmax)then * WRITE(99,*)' FATAL ERROR! GNTAGE3: root not found ' * WRITE(*,*)' STOP: FATAL ERROR ' * CALL exit(0) * STOP m0 = zpars(2) m0 = MAX(m0,mt) goto 40 endif goto 30 40 continue endif CALL star(kw,m0,mt,tm,tn,tscls,lums,GB,zpars) aj = tscls(1) + 1.0d-06*(tscls(2) - tscls(1)) * endif * if(kw.eq.8.or.kw.eq.9)then * * We make a post-MS naked helium star. * To make things easier we put the star at the TMS point * so it actually begins as type 8. * kw = 8 mmin = mc CALL star(kw,mmin,mc,tm,tn,tscls,lums,GB,zpars) mcx = mcgbf(lums(2),GB,lums(6)) if(mcx.ge.mc)then * WRITE(99,*)' FATAL ERROR! GNTAGE9: mmin too big ' * WRITE(*,*)' STOP: FATAL ERROR ' * CALL exit(0) * STOP m0 = mt goto 80 endif f = mcx - mc mmax = mt do 50 , j = 1,jmax CALL star(kw,mmax,mc,tm,tn,tscls,lums,GB,zpars) mcy = mcgbf(lums(2),GB,lums(6)) if(mcy.gt.mc) goto 60 mmax = 2.d0*mmax if(j.eq.jmax)then * WRITE(99,*)' FATAL ERROR! GNTAGE9: mmax not found ' * WRITE(*,*)' STOP: FATAL ERROR ' * CALL exit(0) * STOP m0 = mt goto 80 endif 50 continue 60 continue fmid = mcy - mc * Use the bisection method to find m0 if(f*fmid.ge.0.d0)then * WRITE(99,*)' FATAL ERROR! GNTAGE9: root not bracketed ' * WRITE(*,*)' STOP: FATAL ERROR ' * CALL exit(0) * STOP m0 = mt goto 80 endif m0 = mmin dm = mmax - mmin do 70 , j = 1,jmax dm = 0.5d0*dm mmid = m0 + dm CALL star(kw,mmid,mc,tm,tn,tscls,lums,GB,zpars) mcy = mcgbf(lums(2),GB,lums(6)) fmid = mcy - mc if(fmid.lt.0.d0) m0 = mmid if(ABS(dm).lt.macc.or.ABS(fmid).lt.tiny) goto 80 if(j.eq.jmax)then * WRITE(99,*)' FATAL ERROR! GNTAGE9: root not found ' * WRITE(*,*)' STOP: FATAL ERROR ' * CALL exit(0) * STOP m0 = mt goto 80 endif 70 continue 80 continue * CALL star(kw,m0,mt,tm,tn,tscls,lums,GB,zpars) aj = tm + 1.0d-10*tm * endif * if(kw.eq.14)then * kw = 4 m0 = mt mcy = mcagbf(m0) aj = mc/mcy CALL star(kw,m0,mt,tm,tn,tscls,lums,GB,zpars) if(m0.le.zpars(2))then mcx = mcgbf(lums(4),GB,lums(6)) else mcx = mcheif(m0,zpars(2),zpars(10)) end if mc = mcx + (mcy - mcx)*aj aj = tscls(2) + aj*tscls(3) endif * RETURN END *** cc//instar.f *** SUBROUTINE instar * * * Initialization of collision matrix. * ------------------------ * implicit none integer i,j,ktype(0:14,0:14) common /TYPES/ ktype * * Initialize stellar collision matrix. * ktype(0,0) = 1 do 10 , j = 1,6 ktype(0,j) = j ktype(1,j) = j 10 continue ktype(0,7) = 4 ktype(1,7) = 4 do 15 , j = 8,12 if(j.ne.10)then ktype(0,j) = 6 else ktype(0,j) = 3 endif ktype(1,j) = ktype(0,j) 15 continue ktype(2,2) = 3 do 20 , i = 3,14 ktype(i,i) = i 20 continue ktype(5,5) = 4 ktype(7,7) = 1 ktype(10,10) = 15 ktype(13,13) = 14 do 25 , i = 2,5 do 30 j = i+1,12 ktype(i,j) = 4 30 continue 25 continue ktype(2,3) = 3 ktype(2,6) = 5 ktype(2,10) = 3 ktype(2,11) = 5 ktype(2,12) = 5 ktype(3,6) = 5 ktype(3,10) = 3 ktype(3,11) = 5 ktype(3,12) = 5 ktype(6,7) = 4 ktype(6,8) = 6 ktype(6,9) = 6 ktype(6,10) = 5 ktype(6,11) = 6 ktype(6,12) = 6 ktype(7,8) = 8 ktype(7,9) = 9 ktype(7,10) = 7 ktype(7,11) = 9 ktype(7,12) = 9 ktype(8,9) = 9 ktype(8,10) = 7 ktype(8,11) = 9 ktype(8,12) = 9 ktype(9,10) = 7 ktype(9,11) = 9 ktype(9,12) = 9 ktype(10,11) = 9 ktype(10,12) = 9 ktype(11,12) = 12 do 35 , i = 0,12 ktype(i,13) = 13 ktype(i,14) = 14 35 continue ktype(13,14) = 14 * * Increase common-envelope cases by 100. do 40 , i = 0,9 do 45 , j = i,14 if(i.le.1.or.i.eq.7)then if(j.ge.2.and.j.le.9.and.j.ne.7)then ktype(i,j) = ktype(i,j) + 100 endif else ktype(i,j) = ktype(i,j) + 100 endif 45 continue 40 continue * * Assign the remaining values by symmetry. do 50 , i = 1,14 do 55 , j = 0,i-1 ktype(i,j) = ktype(j,i) 55 continue 50 continue * return end *** cc//mix.f *** SUBROUTINE MIX(M0,M,AJ,KS,ZPARS) * * Author : J. R. Hurley * Date : 7th July 1998 * * Evolution parameters for mixed star. * ------------------------------------ * implicit none * INTEGER KS(2),I1,I2,K1,K2,KW,ICASE INTEGER KTYPE(0:14,0:14) COMMON /TYPES/ KTYPE REAL*8 M0(2),M(2),AJ(2),ZPARS(20) REAL*8 TSCLS(20),LUMS(10),GB(10),TMS1,TMS2,TMS3,TN REAL*8 M01,M02,M03,M1,M2,M3,AGE1,AGE2,AGE3,MC3,MCH PARAMETER(MCH=1.44D0) REAL*8 NETA,BWIND,HEWIND,MXNS COMMON /VALUE1/ NETA,BWIND,HEWIND,MXNS * * * Define global indices with body #I1 being most evolved. IF(KS(1).GE.KS(2))THEN I1 = 1 I2 = 2 ELSE I1 = 2 I2 = 1 END IF * * Specify case index for collision treatment. K1 = KS(I1) K2 = KS(I2) ICASE = KTYPE(K1,K2) * if(icase.gt.100) WRITE(66,*)' MIX ERROR ICASE>100 ',icase,k1,k2 * * Determine evolution time scales for first star. M01 = M0(I1) M1 = M(I1) AGE1 = AJ(I1) CALL star(K1,M01,M1,TMS1,TN,TSCLS,LUMS,GB,ZPARS) * * Obtain time scales for second star. M02 = M0(I2) M2 = M(I2) AGE2 = AJ(I2) CALL star(K2,M02,M2,TMS2,TN,TSCLS,LUMS,GB,ZPARS) * * Check for planetary systems - defined as HeWDs and low-mass WDs! IF(K1.EQ.10.AND.M1.LT.0.05)THEN ICASE = K2 IF(K2.LE.1)THEN ICASE = 1 AGE1 = 0.D0 ENDIF ELSEIF(K1.GE.11.AND.M1.LT.0.5.AND.ICASE.EQ.6)THEN ICASE = 9 ENDIF IF(K2.EQ.10.AND.M2.LT.0.05)THEN ICASE = K1 IF(K1.LE.1)THEN ICASE = 1 AGE2 = 0.D0 ENDIF ENDIF * * Specify total mass. M3 = M1 + M2 M03 = M01 + M02 KW = ICASE AGE3 = 0.d0 * * Restrict merged stars to masses less than 100 Msun. IF(M3.GE.100.D0)THEN M3 = 99.D0 M03 = MIN(M03,M3) ENDIF * * Evaluate apparent age and other parameters. * IF(ICASE.EQ.1)THEN * Specify new age based on complete mixing. IF(K1.EQ.7) KW = 7 CALL star(KW,M03,M3,TMS3,TN,TSCLS,LUMS,GB,ZPARS) AGE3 = 0.1d0*TMS3*(AGE1*M1/TMS1 + AGE2*M2/TMS2)/M3 ELSEIF(ICASE.EQ.3.OR.ICASE.EQ.6.OR.ICASE.EQ.9)THEN MC3 = M1 CALL gntage(MC3,M3,KW,ZPARS,M03,AGE3) ELSEIF(ICASE.EQ.4)THEN MC3 = M1 AGE3 = AGE1/TMS1 CALL gntage(MC3,M3,KW,ZPARS,M03,AGE3) ELSEIF(ICASE.EQ.7)THEN CALL star(KW,M03,M3,TMS3,TN,TSCLS,LUMS,GB,ZPARS) AGE3 = TMS3*(AGE2*M2/TMS2)/M3 ELSEIF(ICASE.LE.12)THEN * Ensure that a new WD has the initial mass set correctly. M03 = M3 IF(ICASE.LT.12.AND.M3.GE.MCH)THEN M3 = 0.D0 KW = 15 ENDIF ELSEIF(ICASE.EQ.13.OR.ICASE.EQ.14)THEN * Set unstable Thorne-Zytkow object with fast mass loss of envelope * unless the less evolved star is a WD, NS or BH. IF(K2.LT.10)THEN M03 = M1 M3 = M1 ENDIF IF(ICASE.EQ.13.AND.M3.GT.MXNS) KW = 14 ELSEIF(ICASE.EQ.15)THEN M3 = 0.D0 ELSEIF(ICASE.GT.100)THEN * Common envelope case which should only be used after COMENV. KW = K1 AGE3 = AGE1 M3 = M1 M03 = M01 ELSE * This should not be reached. KW = 1 M03 = M3 ENDIF * * Put the result in *1. * KS(1) = KW KS(2) = 15 M(1) = M3 M(2) = 0.D0 M0(1) = M03 AJ(1) = AGE3 * RETURN END *** cc//rl.f *** REAL*8 FUNCTION RL(Q) IMPLICIT NONE REAL*8 Q,P * * A function to evaluate R_L/a(q), Eggleton 1983. * P = Q**(1.d0/3.d0) RL = 0.49d0*P*P/(0.6d0*P*P + LOG(1.d0+P)) * RETURN END ***