Nbody6
gntage.f
Go to the documentation of this file.
1  SUBROUTINE gntage(mc,mt,kw,zpars,m0,aj)
2 *
3 * A routine to determine the age of a giant from its core mass and type.
4 * ----------------------------------------------------------------------
5 *
6 * Author : C. A. Tout
7 * Date : 24th September 1996
8 * Revised: 21st February 1997 to include core-helium-burning stars
9 *
10 * Rewritten: 2nd January 1998 by J. R. Hurley to be compatible with
11 * new evolution routines and to include new stellar types.
12 *
13  implicit none
14 *
15  integer kw
16  integer j,jmax
17  parameter(jmax=30)
18 *
19  real*8 mc,mt,m0,aj,tm,tn
20  real*8 tscls(20),lums(10),gb(10),zpars(20)
21  real*8 mmin,mmax,mmid,dm,f,fmid,dell,derl,lum
22  real*8 macc,lacc,tiny
23  parameter(macc=0.00001d0,lacc=0.0001d0,tiny=1.0d-14)
24  real*8 mcx,mcy
25 *
28 *
29 * This should only be entered with KW = 3, 4, 5, 6 or 9
30 *
31 * First we check that we don't have a CheB star
32 * with too small a core mass.
33  if(kw.eq.4)then
34 * Set the minimum CHeB core mass using M = Mflash
35  mcy = mcheif(zpars(2),zpars(2),zpars(10))
36  if(mc.le.mcy) kw = 3
37  endif
38 *
39 * Next we check that we don't have a GB star for M => Mfgb
40  if(kw.eq.3)then
41 * Set the maximum GB core mass using M = Mfgb
42  mcy = mcheif(zpars(3),zpars(2),zpars(9))
43  if(mc.ge.mcy)then
44  kw = 4
45  aj = 0.d0
46  endif
47  endif
48 *
49  if(kw.eq.6)then
50 *
51 * We try to start the star from the start of the SAGB by
52 * setting Mc = Mc,TP.
53 *
54  mcy = 0.44d0*2.25d0 + 0.448d0
55  if(mc.gt.mcy)then
56 * A type 6 with this sized core mass cannot exist as it should
57 * already have become a NS or BH as a type 5.
58 * We set it up so that it will.
59  mcx = (mc + 0.35d0)/0.773d0
60  elseif(mc.ge.0.8d0)then
61  mcx = (mc - 0.448d0)/0.44d0
62  else
63  mcx = mc
64  endif
65  m0 = mbagbf(mcx)
66  if(m0.lt.tiny)then
67 * Carbon core mass is less then the minimum for the start of SAGB.
68 * This must be the case of a low-mass C/O or O/Ne WD with only a
69 * very small envelope added or possibly the merger of a helium star
70 * with a main sequence star. We will set m0 = mt and then reset the
71 * core mass to allow for some helium to be added to the C/O core.
72  kw = 14
73  else
74  CALL star(kw,m0,mt,tm,tn,tscls,lums,gb,zpars)
75  aj = tscls(13)
76  endif
77  endif
78 *
79  if(kw.eq.5)then
80 *
81 * We fit a Helium core mass at the base of the AGB.
82 *
83  m0 = mbagbf(mc)
84  if(m0.lt.tiny)then
85 * Helium core mass is less then the BAGB minimum.
86  kw = 14
87  else
88  CALL star(kw,m0,mt,tm,tn,tscls,lums,gb,zpars)
89  aj = tscls(2) + tscls(3)
90  endif
91  endif
92 *
93 *
94  if(kw.eq.4)then
95 *
96 * The supplied age is actually the fractional age, fage, of CHeB lifetime
97 * that has been completed, ie. 0 <= aj <= 1.
98 *
99 * Get the minimum, fage=1, and maximum, fage=0, allowable masses
100  mcy = mcagbf(zpars(2))
101  if(mc.ge.mcy)then
102  mmin = mbagbf(mc)
103  else
104  mmin = zpars(2)
105  endif
106  mmax = mheif(mc,zpars(2),zpars(10))
107  if(aj.lt.tiny)then
108  m0 = mmax
109  goto 20
110  elseif(aj.ge.1.d0)then
111  m0 = mmin
112  goto 20
113  endif
114 * Use the bisection method to find m0
115  fmid = (1.d0-aj)*mcheif(mmax,zpars(2),zpars(10)) +
116  & aj*mcagbf(mmax) - mc
117  f = (1.d0-aj)*mcheif(mmin,zpars(2),zpars(10)) +
118  & aj*mcagbf(mmin) - mc
119  if(f*fmid.ge.0.d0)then
120 * This will probably occur if mc is just greater than the minimum
121 * allowed mass for a CHeB star and fage > 0.
122  kw = 3
123  goto 90
124  endif
125  m0 = mmin
126  dm = mmax - mmin
127  do 10 j = 1,jmax
128  dm = 0.5d0*dm
129  mmid = m0 + dm
130  fmid = (1.d0-aj)*mcheif(mmid,zpars(2),zpars(10)) +
131  & aj*mcagbf(mmid) - mc
132  if(fmid.lt.0.d0) m0 = mmid
133  if(abs(dm).lt.macc.or.abs(fmid).lt.tiny) goto 20
134  10 continue
135  20 continue
136 *
137  CALL star(kw,m0,mt,tm,tn,tscls,lums,gb,zpars)
138  aj = tscls(2) + aj*tscls(3)
139 *
140  endif
141 *
142  90 continue
143 *
144  if(kw.eq.3)then
145 *
146 * First we double check that we don't have a GB star for M => Mfgb
147  mcy = mcheif(zpars(3),zpars(2),zpars(9))
148 * Next we find an m0 so as to place the star at the BGB
149  mcx = mcheif(zpars(2),zpars(2),zpars(9))
150  if(mc.gt.mcx)then
151  m0 = mheif(mc,zpars(2),zpars(9))
152  else
153 * Use Newton-Raphson to find m0 from Lbgb
154  m0 = zpars(2)
155  CALL star(kw,m0,mt,tm,tn,tscls,lums,gb,zpars)
156  lum = lmcgbf(mc,gb)
157  j = 0
158  30 continue
159  dell = lbgbf(m0) - lum
160  if(abs(dell/lum).le.lacc) goto 40
161  derl = lbgbdf(m0)
162  m0 = m0 - dell/derl
163  j = j + 1
164  goto 30
165  40 continue
166  endif
167  CALL star(kw,m0,mt,tm,tn,tscls,lums,gb,zpars)
168  aj = tscls(1) + 1.0d-06*(tscls(2) - tscls(1))
169 *
170  endif
171 *
172  if(kw.eq.8.or.kw.eq.9)then
173 *
174 * We make a post-MS naked helium star.
175 * To make things easier we put the star at the TMS point
176 * so it actually begins as type 8.
177 *
178  kw = 8
179  mmin = mc
180  CALL star(kw,mmin,mc,tm,tn,tscls,lums,gb,zpars)
181  mcx = mcgbf(lums(2),gb,lums(6))
182  f = mcx - mc
183  mmax = mt
184  do 50 j = 1,jmax
185  CALL star(kw,mmax,mc,tm,tn,tscls,lums,gb,zpars)
186  mcy = mcgbf(lums(2),gb,lums(6))
187  if(mcy.gt.mc) goto 60
188  mmax = 2.d0*mmax
189  50 continue
190  60 continue
191  fmid = mcy - mc
192 * Use the bisection method to find m0
193  m0 = mmin
194  dm = mmax - mmin
195  do 70 j = 1,jmax
196  dm = 0.5d0*dm
197  mmid = m0 + dm
198  CALL star(kw,mmid,mc,tm,tn,tscls,lums,gb,zpars)
199  mcy = mcgbf(lums(2),gb,lums(6))
200  fmid = mcy - mc
201  if(fmid.lt.0.d0) m0 = mmid
202  if(abs(dm).lt.macc.or.abs(fmid).lt.tiny) goto 80
203  70 continue
204  80 continue
205 *
206  CALL star(kw,m0,mt,tm,tn,tscls,lums,gb,zpars)
207  aj = tm + 1.0d-10*tm
208 *
209  endif
210 *
211  if(kw.eq.14)then
212 *
213  kw = 4
214  m0 = mt
215  mcy = mcagbf(m0)
216  aj = mc/mcy
217  CALL star(kw,m0,mt,tm,tn,tscls,lums,gb,zpars)
218  if(m0.le.zpars(2))then
219  mcx = mcgbf(lums(4),gb,lums(6))
220  else
221  mcx = mcheif(m0,zpars(2),zpars(10))
222  end if
223  mc = mcx + (mcy - mcx)*aj
224  aj = tscls(2) + aj*tscls(3)
225  endif
226 *
227  RETURN
228  END