Nbody6
 All Files Functions Variables
star.f
Go to the documentation of this file.
1 ***
2  SUBROUTINE star(kw,mass,mt,tm,tn,tscls,lums,GB,zpars)
3 *
4 *
5 * Stellar luminosity & evolution time.
6 * ------------------------------------
7 *
8  implicit none
9 *
10  integer kw
11 *
12  real*8 mass,mt,tm,tn,tscls(20),lums(10),gb(10),zpars(20)
13  real*8 tgb,tbagb,mch,mcmax,mc1,mc2,mcbagb,dx,am
14  real*8 lambda,tau,mtc,mass0
15  parameter(mch=1.44d0)
16 *
17  real*8 lzamsf,lzahbf,lzhef
18  real*8 tbgbf,thookf,thef,themsf,mcgbf,mcagbf,mcheif,mcgbtf
19  real*8 ltmsf,lbgbf,lheif,lhef,lbagbf,lmcgbf
20  external lzamsf,lzahbf,lzhef
21  external tbgbf,thookf,thef,themsf,mcgbf,mcagbf,mcheif,mcgbtf
22  external ltmsf,lbgbf,lheif,lhef,lbagbf,lmcgbf
23 *
24 * Computes the characteristic luminosities at different stages (LUMS),
25 * and various timescales (TSCLS).
26 * Ref: P.P. Eggleton, M.J. Fitchett & C.A. Tout (1989) Ap.J. 347, 998.
27 *
28 * Revised 27th March 1995 by C. A. Tout
29 * and 24th October 1995 to include metallicity
30 * and 13th December 1996 to include naked helium stars
31 *
32 * Revised 5th April 1997 by J. R. Hurley
33 * to include Z=0.001 as well as Z=0.02, convective overshooting,
34 * MS hook and more elaborate CHeB. It now also sets the Giant
35 * Branch parameters relevant to the mass of the star.
36 *
37 * Revised 21st January 2011 by A. D. Railton
38 * to include the pre-mainsequence evolution for 0.1-8 solar masses
39 * with solar metallicity. New timescale (15) added.
40 * KW=-1 for preMS evolution.
41 *
42 * ------------------------------------------------------------
43 * Times: 1; BGB 2; He ignition 3; He burning
44 * 4; Giant t(inf1) 5; Giant t(inf2) 6; Giant t(Mx)
45 * 7; FAGB t(inf1) 8; FAGB t(inf2) 9; FAGB t(Mx)
46 * 10; SAGB t(inf1) 11; SAGB t(inf2) 12; SAGB t(Mx)
47 * 13; TP 14; t(Mcmax) 15; PreMS
48 *
49 * LUMS: 1; ZAMS 2; End MS 3; BGB
50 * 4; He ignition 5; He burning 6; L(Mx)
51 * 7; BAGB 8; TP
52 *
53 * GB: 1; effective A(H) 2; A(H,He) 3; B
54 * 4; D 5; p 6; q
55 * 7; Mx 8; A(He) 9; Mc,BGB
56 *
57 * ------------------------------------------------------------
58 *
59 *
60  mass0 = mass
61  if(mass0.gt.100.d0) mass = 100.d0
62 *
63  if(kw.ge.7.and.kw.le.9) goto 90
64  if(kw.ge.10) goto 95
65 *
66 * Mass should be in the range [0.1,8] (Extrapolate at your own risk!)
67 * PreMS timescale
68 *
69  if(mass.gt.0.01.and.mass.lt.8.0)then
70  tscls(15)=43.628d0-(35.835d0*mass**1.5029d-2)*
71  & exp(mass*3.9608d-3)
72  tscls(15)=10.d0**tscls(15)/1.0d+06
73  endif
74 *
75 * MS and BGB times
76 *
77  tscls(1) = tbgbf(mass)
78  tm = max(zpars(8),thookf(mass))*tscls(1)
79 *
80 * Zero- and terminal age main sequence luminosity
81 *
82  lums(1) = lzamsf(mass)
83  lums(2) = ltmsf(mass)
84 *
85 * Set the GB parameters
86 *
87  gb(1) = max(-4.8d0,min(-5.7d0+0.8d0*mass,-4.1d0+0.14d0*mass))
88  gb(1) = 10.d0**gb(1)
89  gb(2) = 1.27d-05
90  gb(8) = 8.0d-05
91  gb(3) = max(3.0d+04,500.d0 + 1.75d+04*mass**0.6d0)
92  if(mass.le.2.0)then
93  gb(4) = zpars(6)
94  gb(5) = 6.d0
95  gb(6) = 3.d0
96  elseif(mass.lt.2.5)then
97  dx = zpars(6) - (0.975d0*zpars(6) - 0.18d0*2.5d0)
98  gb(4) = zpars(6) - dx*(mass - 2.d0)/(0.5d0)
99  gb(5) = 6.d0 - (mass - 2.d0)/(0.5d0)
100  gb(6) = 3.d0 - (mass - 2.d0)/(0.5d0)
101  else
102  gb(4) = max(-1.d0,0.5d0*zpars(6) - 0.06d0*mass)
103  gb(4) = max(gb(4),0.975d0*zpars(6) - 0.18d0*mass)
104  gb(5) = 5.d0
105  gb(6) = 2.d0
106  endif
107  gb(4) = 10.d0**gb(4)
108  gb(7) = (gb(3)/gb(4))**(1.d0/(gb(5)-gb(6)))
109 *
110 * Change in slope of giant L-Mc relation.
111  lums(6) = gb(4)*gb(7)**gb(5)
112 *
113 * HeI ignition luminosity
114  lums(4) = lheif(mass,zpars(2))
115  lums(7) = lbagbf(mass,zpars(2))
116 *
117  if(mass.lt.0.1d0.and.kw.le.1)then
118  tscls(2) = 1.1d0*tscls(1)
119  tscls(3) = 0.1d0*tscls(1)
120  lums(3) = lbgbf(mass)
121  goto 96
122  endif
123 *
124  if(mass.le.zpars(3))then
125 * Base of the giant branch luminosity
126  lums(3) = lbgbf(mass)
127 * Set GB timescales
128  tscls(4) = tscls(1) + (1.d0/((gb(5)-1.d0)*gb(1)*gb(4)))*
129  & ((gb(4)/lums(3))**((gb(5)-1.d0)/gb(5)))
130  tscls(6) = tscls(4) - (tscls(4) - tscls(1))*((lums(3)/lums(6))
131  & **((gb(5)-1.d0)/gb(5)))
132  tscls(5) = tscls(6) + (1.d0/((gb(6)-1.d0)*gb(1)*gb(3)))*
133  & ((gb(3)/lums(6))**((gb(6)-1.d0)/gb(6)))
134 * Set Helium ignition time
135  if(lums(4).le.lums(6))then
136  tscls(2) = tscls(4) - (1.d0/((gb(5)-1.d0)*gb(1)*gb(4)))*
137  & ((gb(4)/lums(4))**((gb(5)-1.d0)/gb(5)))
138  else
139  tscls(2) = tscls(5) - (1.d0/((gb(6)-1.d0)*gb(1)*gb(3)))*
140  & ((gb(3)/lums(4))**((gb(6)-1.d0)/gb(6)))
141  endif
142  tgb = tscls(2) - tscls(1)
143  if(mass.le.zpars(2))then
144  mc1 = mcgbf(lums(4),gb,lums(6))
145  mc2 = mcagbf(mass)
146  lums(5) = lzahbf(mass,mc1,zpars(2))
147  tscls(3) = thef(mass,mc1,zpars(2))
148  else
149  lums(5) = lhef(mass)*lums(4)
150  tscls(3) = thef(mass,1.d0,zpars(2))*tscls(1)
151  endif
152  else
153 * Note that for M>zpars(3) there is no GB as the star goes from
154 * HG -> CHeB -> AGB. So in effect tscls(1) refers to the time of
155 * Helium ignition and not the BGB.
156  tscls(2) = tscls(1)
157  tscls(3) = thef(mass,1.d0,zpars(2))*tscls(1)
158 * This now represents the luminosity at the end of CHeB, ie. BAGB
159  lums(5) = lums(7)
160 * We set lums(3) to be the luminosity at the end of the HG
161  lums(3) = lums(4)
162  endif
163 *
164 * Set the core mass at the BGB.
165 *
166  if(mass.le.zpars(2))then
167  gb(9) = mcgbf(lums(3),gb,lums(6))
168  elseif(mass.le.zpars(3))then
169  gb(9) = mcheif(mass,zpars(2),zpars(9))
170  else
171  gb(9) = mcheif(mass,zpars(2),zpars(10))
172  endif
173 *
174 * FAGB time parameters
175 *
176  tbagb = tscls(2) + tscls(3)
177  tscls(7) = tbagb + (1.d0/((gb(5)-1.d0)*gb(8)*gb(4)))*
178  & ((gb(4)/lums(7))**((gb(5)-1.d0)/gb(5)))
179  tscls(9) = tscls(7) - (tscls(7) - tbagb)*((lums(7)/lums(6))
180  & **((gb(5)-1.d0)/gb(5)))
181  tscls(8) = tscls(9) + (1.d0/((gb(6)-1.d0)*gb(8)*gb(3)))*
182  & ((gb(3)/lums(6))**((gb(6)-1.d0)/gb(6)))
183 *
184 * Now to find Ltp and ttp using Mc,He,tp
185 *
186  mcbagb = mcagbf(mass)
187  mc1 = mcbagb
188  if(mc1.ge.0.8d0.and.mc1.lt.2.25d0)then
189 * The star undergoes dredge-up at Ltp causing a decrease in Mc,He
190  mc1 = 0.44d0*mc1 + 0.448d0
191  endif
192  lums(8) = lmcgbf(mc1,gb)
193  if(mc1.le.gb(7))then
194  tscls(13) = tscls(7) - (1.d0/((gb(5)-1.d0)*gb(8)*gb(4)))*
195  & (mc1**(1.d0-gb(5)))
196  else
197  tscls(13) = tscls(8) - (1.d0/((gb(6)-1.d0)*gb(8)*gb(3)))*
198  & (mc1**(1.d0-gb(6)))
199  endif
200 *
201 * SAGB time parameters
202 *
203  if(mc1.le.gb(7))then
204  tscls(10) = tscls(13) + (1.d0/((gb(5)-1.d0)*gb(2)*gb(4)))*
205  & ((gb(4)/lums(8))**((gb(5)-1.d0)/gb(5)))
206  tscls(12) = tscls(10) - (tscls(10) - tscls(13))*
207  & ((lums(8)/lums(6))**((gb(5)-1.d0)/gb(5)))
208  tscls(11) = tscls(12) + (1.d0/((gb(6)-1.d0)*gb(2)*gb(3)))*
209  & ((gb(3)/lums(6))**((gb(6)-1.d0)/gb(6)))
210  else
211  tscls(10) = tscls(7)
212  tscls(12) = tscls(9)
213  tscls(11) = tscls(13) + (1.d0/((gb(6)-1.d0)*gb(2)*gb(3)))*
214  & ((gb(3)/lums(8))**((gb(6)-1.d0)/gb(6)))
215  endif
216 *
217 * Get an idea of when Mc,C = Mc,C,max on the AGB
218  tau = tscls(2) + tscls(3)
219  mc2 = mcgbtf(tau,gb(8),gb,tscls(7),tscls(8),tscls(9))
220  mcmax = max(max(mch,0.773d0*mcbagb - 0.35d0),1.05d0*mc2)
221 *
222  if(mcmax.le.mc1)then
223  if(mcmax.le.gb(7))then
224  tscls(14) = tscls(7) - (1.d0/((gb(5)-1.d0)*gb(8)*gb(4)))*
225  & (mcmax**(1.d0-gb(5)))
226  else
227  tscls(14) = tscls(8) - (1.d0/((gb(6)-1.d0)*gb(8)*gb(3)))*
228  & (mcmax**(1.d0-gb(6)))
229  endif
230  else
231 * Star is on SAGB and we need to increase mcmax if any 3rd
232 * dredge-up has occurred.
233  lambda = min(0.9d0,0.3d0+0.001d0*mass**5)
234  mcmax = (mcmax - lambda*mc1)/(1.d0 - lambda)
235  if(mcmax.le.gb(7))then
236  tscls(14) = tscls(10) - (1.d0/((gb(5)-1.d0)*gb(2)*gb(4)))*
237  & (mcmax**(1.d0-gb(5)))
238  else
239  tscls(14) = tscls(11) - (1.d0/((gb(6)-1.d0)*gb(2)*gb(3)))*
240  & (mcmax**(1.d0-gb(6)))
241  endif
242  endif
243  tscls(14) = max(tbagb,tscls(14))
244  if(mass.ge.100.d0)then
245  tn = tscls(2)
246 * goto 100
247  endif
248 *
249 * Calculate the nuclear timescale - the time of exhausting
250 * nuclear fuel without further mass loss.
251 * This means we want to find when Mc = Mt which defines Tn and will
252 * be used in determining the timestep required. Note that after some
253 * stars reach Mc = Mt there will be a Naked Helium Star lifetime
254 * which is also a nuclear burning period but is not included in Tn.
255 *
256  if(abs(mt-mcbagb).lt.1.0d-14.and.kw.lt.5)then
257  tn = tbagb
258  else
259 * Note that the only occurence of Mc being double-valued is for stars
260 * that have a dredge-up. If Mt = Mc where Mc could be the value taken
261 * from CHeB or from the AGB we need to check the current stellar type.
262  if(mt.gt.mcbagb.or.(mt.ge.mc1.and.kw.gt.4))then
263  if(kw.eq.6)then
264  lambda = min(0.9d0,0.3d0+0.001d0*mass**5)
265  mc1 = (mt - lambda*mc1)/(1.d0 - lambda)
266  else
267  mc1 = mt
268  endif
269  if(mc1.le.gb(7))then
270  tn = tscls(10) - (1.d0/((gb(5)-1.d0)*gb(2)*gb(4)))*
271  & (mc1**(1.d0-gb(5)))
272  else
273  tn = tscls(11) - (1.d0/((gb(6)-1.d0)*gb(2)*gb(3)))*
274  & (mc1**(1.d0-gb(6)))
275  endif
276  else
277  if(mass.gt.zpars(3))then
278  mc1 = mcheif(mass,zpars(2),zpars(10))
279  if(mt.le.mc1)then
280  tn = tscls(2)
281  else
282  tn = tscls(2) + tscls(3)*((mt - mc1)/(mcbagb - mc1))
283  endif
284  elseif(mass.le.zpars(2))then
285  mc1 = mcgbf(lums(3),gb,lums(6))
286  mc2 = mcgbf(lums(4),gb,lums(6))
287  if(mt.le.mc1)then
288  tn = tscls(1)
289  elseif(mt.le.mc2)then
290  if(mt.le.gb(7))then
291  tn = tscls(4) - (1.d0/((gb(5)-1.d0)*gb(1)*gb(4)))*
292  & (mt**(1.d0-gb(5)))
293  else
294  tn = tscls(5) - (1.d0/((gb(6)-1.d0)*gb(1)*gb(3)))*
295  & (mt**(1.d0-gb(6)))
296  endif
297  else
298  tn = tscls(2) + tscls(3)*((mt - mc2)/(mcbagb - mc2))
299  endif
300  else
301  mc1 = mcheif(mass,zpars(2),zpars(9))
302  mc2 = mcheif(mass,zpars(2),zpars(10))
303  if(mt.le.mc1)then
304  tn = tscls(1)
305  elseif(mt.le.mc2)then
306  tn = tscls(1) + tgb*((mt - mc1)/(mc2 - mc1))
307  else
308  tn = tscls(2) + tscls(3)*((mt - mc2)/(mcbagb - mc2))
309  endif
310  endif
311  endif
312  endif
313  tn = min(tn,tscls(14))
314 *
315  goto 100
316 *
317  90 continue
318 *
319 * Calculate Helium star Main Sequence lifetime.
320 *
321  tm = themsf(mass)
322  tscls(1) = tm
323 *
324 * Zero- and terminal age Helium star main sequence luminosity
325 *
326  lums(1) = lzhef(mass)
327  am = max(0.d0,0.85d0-0.08d0*mass)
328  lums(2) = lums(1)*(1.d0+0.45d0+am)
329 *
330 * Set the Helium star GB parameters
331 *
332  gb(8) = 8.0d-05
333  gb(3) = 4.1d+04
334  gb(4) = 5.5d+04/(1.d0+0.4d0*mass**4)
335  gb(5) = 5.d0
336  gb(6) = 3.d0
337  gb(7) = (gb(3)/gb(4))**(1.d0/(gb(5)-gb(6)))
338 * Change in slope of giant L-Mc relation.
339  lums(6) = gb(4)*gb(7)**gb(5)
340 *
341 *** Set Helium star GB timescales
342 *
343  mc1 = mcgbf(lums(2),gb,lums(6))
344  tscls(4) = tm + (1.d0/((gb(5)-1.d0)*gb(8)*gb(4)))*
345  & mc1**(1.d0-gb(5))
346  tscls(6) = tscls(4) - (tscls(4) - tm)*((gb(7)/mc1)
347  & **(1.d0-gb(5)))
348  tscls(5) = tscls(6) + (1.d0/((gb(6)-1.d0)*gb(8)*gb(3)))*
349  & gb(7)**(1.d0-gb(6))
350 *
351 * Get an idea of when Mc = MIN(Mt,Mc,C,max) on the GB
352  mtc = min(mt,1.45d0*mt-0.31d0)
353  if(mtc.le.0.d0) mtc = mt
354  mcmax = min(mtc,max(mch,0.773d0*mass-0.35d0))
355  if(mcmax.le.gb(7))then
356  tscls(14) = tscls(4) - (1.d0/((gb(5)-1.d0)*gb(8)*gb(4)))*
357  & (mcmax**(1.d0-gb(5)))
358  else
359  tscls(14) = tscls(5) - (1.d0/((gb(6)-1.d0)*gb(8)*gb(3)))*
360  & (mcmax**(1.d0-gb(6)))
361  endif
362  tscls(14) = max(tscls(14),tm)
363  tn = tscls(14)
364 *
365  goto 100
366 *
367  95 continue
368  tm = 1.0d+10
369  96 continue
370  tn = 1.0d+10
371 *
372  100 continue
373  mass = mass0
374 *
375  return
376  end
377 ***