| Home | Trees | Indices | Help |
|---|
|
|
1 #!/usr/bin/env python3
2 # *-* Coding: UTF-8 *-*
3
4 from __future__ import division, absolute_import
5
6 __author__ = "Eduardo dos Santos Pereira"
7 __email__ = "pereira.somoza@gmail.com"
8 __credits__ = ["Eduardo dos Santos Pereira"]
9 __license__ = "GPLV3"
10 __version__ = "1.0.1"
11 __maintainer__ = "Eduardo dos Santos Pereira"
12 __status__ = "Stable"
13
14
15 """Cosmic Star Formation Rate
16
17 This file is part of pystar.
18 copyright : Eduardo dos Santos Pereira
19
20 pystar is free software: you can redistribute it and/or modify
21 it under the terms of the GNU General Public License as published by
22 the Free Software Foundation, either version 3 of the License.
23 pystar is distributed in the hope that it will be useful,
24 but WITHOUT ANY WARRANTY; without even the implied warranty of
25 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
26 GNU General Public License for more details.
27
28 You should have received a copy of the GNU General Public License
29 along with Foobar. If not, see <http://www.gnu.org/licenses/>.
30
31 """
32
33 from numpy import log10, sqrt, array
34 from numpy.numarray import zeros, Float64
35 from .structures import structures
36 import scipy.interpolate as spint
37 from scipy.integrate import romberg
38 from .run_kut4 import rk4_int
39
40 import sys
41 pyversion = sys.version_info
42 if pyversion[0] >= 3:
43 from . import filedict
44 else:
45 print("Importing filedict for python2.7")
46 from . import filedict_old as filedict
47
48
50 """The Cosmic Star Formation rate
51 The model used to develop this class was presented by the first time
52 in the article of Pereira and Miranda (2010) - (MNRAS, 401, 1924, 2010).
53
54 The cosmologic background model is passed as a instance parameter:
55 cosmology
56
57 Keyword arguments:
58 tau -- (default - 2.5) time scale, in Gyr, of the CSFR.
59
60 eimf -- (default 1.35) exponent of the Initial Mass Function
61
62 nsch -- (default 1) the normalization factor in the CSFR model
63
64 imfType -- (default S - Salpeter) the Initial Mass Function Type.
65 Possible values:
66 S: Salpeter
67 K: Kroupa
68
69 lmin -- (default 6.0) log10 of the minal mass of the dark halo
70 where it is possible to have star formation.
71
72 zmax -- (defaul 20.0) - the maximum redshift to be considered
73
74 omegam -- (default 0.24) - The dark matter parameter
75
76 omegab -- (default 0.04) - The barionic parameter
77
78 omegal -- (default 0.73) - The dark energy parameter
79
80 h -- (default 0.7) - The h of the Hubble constant (H = h * 100)
81
82 massFunctionType -- (default \"ST\") The type of mass
83 function of dark matter halos used. Possibles values:
84 \"ST\" for Seth and Thormen mass function.
85 \"TK\" for Tinker et al. mass function.
86 """
87
88 - def __init__(self, cosmology,
89 tau=2.29, eimf=1.35, nsch=1, lmin=6.0, zmax=20.0,
90 imfType="S", **kwargs):
91
92 structures.__init__(self, cosmology, **kwargs)
93
94 cacheFile = self._cacheFIle + "_CSFR_" + str(tau)\
95 + "_" + str(eimf) + "_" + str(nsch) + "_" + imfType
96
97 self._cache_dictS = filedict.FileDict(filename=cacheFile + ".cache")
98
99 tau = tau * 1.0e9
100 self._cc = self._tck_ab[1]
101
102 #Cosmic Star Formation Rate normalization
103 self.__esnor = 1.0
104 self.__nsch = nsch
105 self.__tau = tau
106
107 self.imfType = imfType
108
109 self.__aminf1 = 2.5e1
110 self.__amsup1 = 1.4e2
111 self.__amin = 1.0e-1
112
113 self.__imfDict = {"S": self.__imfSalpeter,
114 "K": self.__imfKroupa
115 }
116
117 self.__anorm1 = None
118 self.__eimf = eimf
119 self.__eimf0 = eimf - 1.0
120
121 try:
122
123 self.__astar = self._cache_dictS['astar']
124 self.__csfr = self._cache_dictS['csfr']
125 self.__rho_gas = self._cache_dictS['rho_gas']
126 print("Data CSFR in Cache")
127
128 except:
129 self.__csfr, self.__rho_gas, self.__astar = self.__sfr()
130
131 tck_sf = spint.splrep(self.__astar, self.__csfr)
132 self.__cs = tck_sf[1]
133
134 tck_sg = spint.splrep(self.__astar, self.__rho_gas)
135 self.__cs2 = tck_sg[1]
136
138 """
139 Return a list with keys and functions of IMF's
140 """
141 dic = []
142 for key, value in list(self.__imfDict.items()):
143 dic.append([key, value])
144 return dic
145
152
154 """Return the interpolated, by cubi-spline, barionic accretion
155 rate into structures
156 """
157 j = -1
158
159 while(1):
160 j = j + 1
161 if(a < self._ascale[j + 1] and
162 a >= self._ascale[j] and
163 j <= (len(self._ascale) - 3)
164 ):
165 resp = (self._cc[j] / 6.0) * \
166 (((a - self._ascale[j + 1]) ** 3.0) /
167 (self._ascale[j] - self._ascale[j + 1])
168 - (a - self._ascale[j + 1]) *
169 (self._ascale[j] - self._ascale[j + 1])) \
170 - (self._cc[j + 1] / 6.0) * \
171 (((a - self._ascale[j]) ** 3.0) /
172 (self._ascale[j] - self._ascale[j + 1])
173 - (a - self._ascale[j]) *
174 (self._ascale[j] - self._ascale[j + 1])) \
175 + (self._abt2[j] * (a - self._ascale[j + 1])
176 - self._abt2[j + 1] * (a - self._ascale[j]))\
177 / (self._ascale[j] - self._ascale[j + 1])
178 return resp
179
180 elif(j >= (len(self._ascale) - 2) and j <= len(self._ascale)):
181 return self._abt2[j]
182 elif(a < self._ascale[0]):
183 raise NameError("Error in the spline function")
184 break
185
187 """Return the numerical function to be integrated to calculate
188 the mass density of barions into structures.
189 """
190
191 z = 1.0 / a - 1.0
192
193 if(z < 0.0):
194 z = 0.0
195
196 tage = self._cosmology.age(z)
197
198 age01 = 4.0 * log10(tage) - 2.704e+01
199 age02 = (3.6 - sqrt(age01)) / 2.0
200
201 mi_1 = 1.0e+01 ** age02
202 yr = self.__massEjected(mi_1)
203
204 if(self.__nsch == 1.0):
205 sexp = (1.0 - yr) / self.__tau
206 else:
207 sexp = (1.0 - yr) / self.__tau /\
208 self._cosmology.getRobr0() ** (self.__nsch - 1.0)
209
210 F = zeros(1, type=Float64)
211 F[0] = (-sexp * (rho_g[0]) ** self.__nsch
212 + self.__esnor * self.__spn(a) /
213 self._cosmology.getRobr0()) \
214 * self._cosmology.dt_dz(z) / a ** 2.0
215 return F
216
218 """Return the Cosmic Star Formation Rate
219 from the barionic gas into structures
220 """
221 if(self.__nsch == 1):
222 roes = rg / self.__tau
223 else:
224 roes = (rg ** self.__nsch) / self.__tau\
225 / self._cosmology.getRobr0() ** (self.__nsch - 1.0)
226 return roes
227
229 """Return the Cosmic Star Formation rate, density of barionic
230 gas into structures
231 """
232
233 #Normalization of the cosmic star formation rate
234 rho_g0 = array([1.0e-9])
235 a0 = self._ascale[0]
236 nf = len(self._ascale) - 1
237 af = self._ascale[nf]
238 step = (af - a0) / 100.0
239
240 A, R_g = rk4_int(self.__fcn, a0, rho_g0, af, step)
241
242 ng = len(A) - 1
243
244 if(self.__nsch == 1):
245 roes = R_g[ng] / self.__tau
246 else:
247 roes = (R_g[ng] ** self.__nsch) / self.__tau /\
248 self._cosmology.getRobr0() ** (self.__nsch - 1.0)
249
250 self.__esnor = 1.62593696e-2 / roes
251
252 rho_g0 = array([1.0e-9])
253 a0 = 1.0 / (self._zmax + 1.0)
254 nf = len(self._ascale) - 1
255 af = self._ascale[nf]
256 step = (af - a0) / 5000.
257 A, R_g = rk4_int(self.__fcn, a0, rho_g0, af, step)
258
259 rho_s = self.__csfr_gas(R_g)
260
261 self._cache_dictS['astar'] = A
262 self._cache_dictS['csfr'] = rho_s
263 self._cache_dictS['rho_gas'] = R_g
264
265 return rho_s, R_g, A
266
268 if(self.__anorm1 is None):
269 self.__imfSalpeter(10)
270 amexp1 = (1.0 / m_min) ** self.__eimf0
271 amexp2 = (1.0 / self.__amsup1) ** self.__eimf0
272 amexp3 = (1.0 / 8.0) ** self.__eimf0
273 amexp4 = (1.0 / self.__aminf1) ** self.__eimf0
274 amexp5 = (1.0 / m_min) ** self.__eimf
275 amexp6 = (1.0 / 8.0) ** self.__eimf
276 amexp7 = (1.0 / 10.0) ** self.__eimf
277 amexp8 = (1.0 / self.__aminf1) ** self.__eimf
278 amexp9 = (1.0 / self.__amsup1) ** self.__eimf
279
280 yrem1 = (amexp1 - amexp2) / self.__eimf0
281 yrem2 = 1.156e-01 * (amexp1 - amexp3) / self.__eimf0
282 yrem3 = 1.3e+01 * (amexp4 - amexp2) / self.__eimf0 / 2.4e+01
283 yrem4 = 4.551e-01 * (amexp5 - amexp6) / self.__eimf
284 yrem5 = 1.35e+00 * (amexp6 - amexp7) / self.__eimf
285 yrem6 = 1.40e+00 * (amexp7 - amexp8) / self.__eimf
286 yrem7 = 6.5e+01 * (amexp8 - amexp9) / self.__eimf / 6.0
287
288 yr = self.__anorm1 * (yrem1 - yrem2 - yrem3
289 - yrem4 - yrem5 - yrem6 + yrem7)
290
291 return yr
292
294 """
295 Return the mass integration of the mass ejected by the collapse of the
296 star
297 """
298 if(self.imfType == "S"):
299 return self.__massEjectedSalpeter(m_min)
300
301 else:
302 mEject = (romberg(self.__mPhi, m_min, self.__amsup1, tol=1.48e-04)
303 - romberg(self.__mrPhi, m_min, self.__amsup1, tol=1.48e-04)
304 )
305 return mEject
306
308 return m * self.phi(m)
309
312
314 """
315 Return the remnant mass of the object after the colapse of the star
316 with mass m
317 """
318
319 if(m > 0 and m < 1):
320 return 0
321 elif(m >= 1 and m <= 8):
322 return 0.1156 * m + 0.4551
323 elif(m > 8 and m <= 10):
324 return 1.35
325 elif(m > 10 and m < 25):
326 return 1.4
327 elif(m >= 25 and m <= 145):
328 return (13.0 / 24.0) * (m - 20)
329 else:
330 raise NameError("Error: Out of the mass range...")
331
333 if(self.__anorm1 is None):
334 alpha0 = 0.3
335 alpha1 = 1.3
336 alpha2 = 2.3
337 alpha3 = alpha2
338
339 k0 = 1
340 k1 = k0 * 0.08
341 k2 = k1 * 0.5
342 k3 = k2
343 A = [k0 * (self.__amsup1 ** (1.0 - alpha0)
344 - self.__amin ** (1.0 - alpha0)) / (1.0 - alpha0),
345 k1 * (self.__amsup1 ** (1.0 - alpha1)
346 - self.__amin ** (1.0 - alpha1)) / (1.0 - alpha1),
347 k2 * (self.__amsup1 ** (1.0 - alpha2)
348 - self.__amin ** (1.0 - alpha2)) / (1.0 - alpha2),
349 k3 * (self.__amsup1 ** (1.0 - alpha3)
350 - self.__amin ** (1.0 - alpha3)) / (1.0 - alpha3)
351 ]
352
353 self.__anorm1 = 1 / sum(A)
354
355 if(m > 0.08 and m <= 0.5):
356 return self.__anorm1 * m ** (-1.3)
357 elif(m > 0.5):
358 return self.__anorm1 * m ** (-2.3)
359 else:
360 raise NameError("Mass out of the range")
361
363
364 if(self.__anorm1 is None):
365 self.__amin = 1.0e-1
366 self.__anorm1 = self.__eimf0 / (1.0 / self.__amin ** self.__eimf0
367 - 1.0 / self.__amsup1 ** self.__eimf0)
368
369 return self.__anorm1 * m ** (-(1.0 + self.__eimf))
370
372 """Return the Initial Mass Function
373 """
374 try:
375 return self.__imfDict[self.imfType](m)
376 except:
377 raise NameError("No Defined Initial Mass Function")
378
380 """Return the Cosmic Star Formation rate as a function of z
381 """
382
383 a = 1.0 / (1.0 + z)
384
385 j = -1
386
387 while(1):
388 j = j + 1
389 if(a < self.__astar[j + 1] and
390 a >= self.__astar[j] and
391 j <= (len(self.__astar) - 3)):
392 resp = (self.__cs[j] / 6.0) * \
393 (((a - self.__astar[j + 1]) ** 3.0) /
394 (self.__astar[j] - self.__astar[j + 1])
395 - (a - self.__astar[j + 1]) *
396 (self.__astar[j] - self.__astar[j + 1])) \
397 - (self.__cs[j + 1] / 6.0) * \
398 (((a - self.__astar[j]) ** 3.0) /
399 (self.__astar[j] - self.__astar[j + 1])
400 - (a - self.__astar[j]) *
401 (self.__astar[j] - self.__astar[j + 1])) \
402 + (self.__csfr[j] * (a - self.__astar[j + 1])
403 - self.__csfr[j + 1] * (a - self.__astar[j])) / \
404 (self.__astar[j] - self.__astar[j + 1])
405 return resp[0]
406
407 elif(j >= (len(self.__astar) - 2) and j <= len(self.__astar)):
408 return self.__csfr[j]
409 elif(a < self.__astar[0]):
410 raise NameError("Error in spline csfr")
411 break
412
414 """Return the barionic gas density into structures
415 """
416 a = 1.0 / (1.0 + z)
417
418 j = -1
419
420 while(1):
421 j = j + 1
422 if(a < self.__astar[j + 1] and
423 a >= self.__astar[j] and
424 j <= (len(self.__astar) - 3)):
425 resp = (self.__cs2[j] / 6.0) * ((
426 (a - self.__astar[j + 1]) ** 3.0) /
427 (self.__astar[j] - self.__astar[j + 1])
428 - (a - self.__astar[j + 1]) *
429 (self.__astar[j] - self.__astar[j + 1])) \
430 - (self.__cs2[j + 1] / 6.0) * \
431 (((a - self.__astar[j]) ** 3.0) /
432 (self.__astar[j] - self.__astar[j + 1])
433 - (a - self.__astar[j]) *
434 (self.__astar[j] - self.__astar[j + 1])) \
435 + (self.__rho_gas[j] * (a - self.__astar[j + 1])
436 - self.__rho_gas[j + 1] * (a - self.__astar[j])) \
437 / (self.__astar[j] - self.__astar[j + 1])
438 return resp
439
440 elif(j >= (len(self.__astar) - 2) and j <= len(self.__astar)):
441 return self.__rho_gas[j]
442 elif(a < self.__astar[0]):
443 raise NameError("Error spline gas density")
444 break
445
| Home | Trees | Indices | Help |
|---|
| Generated by Epydoc 3.0.1 on Mon Aug 11 14:43:29 2014 | http://epydoc.sourceforge.net |