-
Notifications
You must be signed in to change notification settings - Fork 24
/
gntage.f
299 lines (299 loc) · 8.8 KB
/
gntage.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
***
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
***