1 """Classes and routines to calculate enrichment and ionization histories.
2 """
3 import copy
4
5 import numpy
6 import scipy
7 import scipy.constants.codata as scc
8
9 import cosmolopy.reionization as cr
10 import cosmolopy.utils as utils
11 from cosmolopy.utils import Saveable
12
13 from cosmolopy import luminosityfunction
14 import padconvolve
15
17 """Make a normalized tophat function (nonzero from tmin to tmax).
18
19 t: array-like
20
21 The result is normalized so that the sum of the bins is equal to
22 norm (not the integral), i.e. there's no multiplication by the
23 width of the bins in the normalization.
24
25 Returns the kernel and the cdf.
26 """
27
28 kernel = numpy.ones(t.shape)
29 kernel[t>tmax] = 0.
30 kernel[t<tmin] = 0.
31 ksum = numpy.sum(kernel)
32 kernel = kernel/ksum
33 cdf = numpy.cumsum(kernel)
34 return kernel, cdf
35
37 data = numpy.loadtxt(filename, skiprows=1, delimiter=',',
38 usecols=range(5))
39 return data
40
41 -class LFIonHistory(Saveable):
42 """Calculate an ionization history based on luminosity function evolution.
43
44 LF parameters enter here, e.g. choice of LF, f_esc_gamma, M*-z slope, alpha
45
46 Calculate ionizing emissivity history (from LF, for example) Ndot(z).
47 Integrate ionizing emissivity to find x_tot(z).
48 Integrate recombinations to find x(z).
49 """
50
51 - def __init__(self, z, lf_params, cosmo,
52 lf_args={}, lf_maglim=None, MStar_z_slope=None, alpha=None,
53 f_esc_gamma=None, clump_fact=None, clump_fact_func=None):
54 """
55 z - array of redshift values
56 lf_params - dictionary with keys 'z', 'MStar', 'phiStar', 'alpha'
57 f_esc_gamma - ionizing photon escape fraction
58 cosmo - cosmological parameter dict (see cosmolopy.parameters)
59 lf_args - dict of keyword args for luminosityfunction.LFHistory
60 lf_maglim - magnitude limit for interating LF
61 MStar_z_slope - slope of high-z M* extrapolation
62 alpha - faint-end slope of the LF
63 """
64
65 self.z = z
66 self.lf_params = lf_params
67 self.cosmo = cosmo
68 self.lf_args = lf_args
69 self.lf_maglim = lf_maglim
70 self.MStar_z_slope = MStar_z_slope
71 self.alpha = alpha
72 self.f_esc_gamma = f_esc_gamma
73 self.clump_fact = clump_fact
74 self.clump_fact_func = clump_fact_func
75
76
77 if alpha is not None:
78 lf_params = copy.deepcopy(lf_params)
79 lf_params['alpha'] = numpy.array(lf_params['alpha'])
80 lf_params['alpha'][:] = alpha
81 self.lf_params = lf_params
82
83 if MStar_z_slope is not None:
84 lf_args = copy.deepcopy(lf_args)
85 if not 'extrap_args' in lf_args:
86 lf_args['extrap_args'] = {}
87 lf_args['extrap_args']['slopes'] = [0, MStar_z_slope]
88 lf_args['extrap_var'] = 'z'
89 self.lf_args = lf_args
90
91 args = copy.deepcopy(self.cosmo)
92 args.update(lf_args)
93 self.lfh = luminosityfunction.LFHistory(params=lf_params,**args)
94
95
96 self.NdotTot = self.lfh.iPhotonRateDensity_z(z, lf_maglim)
97
98
99 self.xTot = self.lfh.ionization(z, lf_maglim)
100
101 if f_esc_gamma is not None:
102 self.xEsc = self.xTot * f_esc_gamma
103 self.integrate_ion_recomb()
104 self.optical_depth()
105 self.optical_depth(use_recomb=False)
106
107 - def make_ion_func(self):
108 self.xEsc = self.xTot * self.f_esc_gamma
109 ion_func = utils.Extrapolate1d(self.z, self.xEsc,
110 bounds_behavior=[0.0,
111 'extrapolate']
112 )
113 return ion_func
114
116 if self.clump_fact is not None:
117 self.clump_fact_func = lambda z1: self.clump_fact
118 if self.clump_fact_func is None:
119 self.clump_fact_func = cr.clumping_factor_HB
120 return self.clump_fact_func
121
123 """Integrate the recombination rate.
124
125 Returns and stores results.
126 """
127 self.make_clump_fact_func()
128 ion_func = self.make_ion_func()
129
130 x, w, t = cr.integrate_ion_recomb(self.z,
131 ion_func=ion_func,
132 clump_fact_func =
133 self.clump_fact_func,
134 **self.cosmo)
135 self.xr = x
136 self.w = w
137 self.t = t
138 return x, w, t
139
140 - def optical_depth(self, use_recomb=True, store=True):
141 """Calculate Thompson electron scattering optical depth."""
142
143
144
145 tau_z0 = cr.optical_depth_instant(self.z[-1], **self.cosmo)
146
147
148 if use_recomb:
149 xPhys = self.xr.copy()
150 else:
151 xPhys = self.xEsc.copy()
152 xPhys[xPhys > 1.0] = 1.0
153 tau_later = (tau_z0 +
154 cr.integrate_optical_depth(xPhys[...,::-1],
155 xPhys[...,::-1],
156 self.z[::-1],
157 **self.cosmo))
158 tau = tau_later[...,::-1]
159 if store:
160 if use_recomb:
161 self.tau = tau
162 else:
163 self.tauEsc = tau
164 return tau
165
167 """Find the redshift at with xr reaches unity."""
168 z_ion_func = utils.Extrapolate1d(self.xr, self.z)
169 self.z_reion = z_ion_func(1.0)
170 return self.z_reion
171
172 -class EnrichmentHistory(Saveable):
173 """Calculate CIV histories based on ionization.
174 """
175
176 - def __init__(self, ionHist, f_x_Z, z_change=None, f_change=None,
177 delay_kernel=None,
178 ):
179 """Calculate instantaneos histories based on ionization.
180
181 Enrichment parameters enter here, e.g. f_x_Z, z_change, f_change.
182
183 Omega_Z is fraction of critical density contributed by metals
184 = nz * mz / rho_crit.
185
186 Z can just as well represent a specific metal element or
187 ionization species given the appropriatly decreased f_x_Z.
188
189 """
190 self.ionHist = ionHist
191 self.f_x_Z = f_x_Z
192 self.z_change = z_change
193 self.f_change=f_change
194 self.z = ionHist.z
195 self.calculate_Omega_Z()
196 self.delay_kernel = delay_kernel
197 if delay_kernel is not None:
198 self.convolve()
199
201 if self.z_change is not None:
202 print "f_x_Z changes from %.3g to %.3g (f = %.3g) at z<%.2f" % \
203 (self.f_x_Z,
204 self.f_x_Z * self.f_change,
205 self.f_change, self.z_change)
206 xfunc = scipy.interpolate.interp1d(self.ionHist.z[::-1],
207 self.ionHist.xTot[::-1])
208 x_change = xfunc(self.z_change)
209 xTot = self.ionHist.xTot
210
211 f_early = self.f_x_Z
212 f_late = self.f_x_Z * self.f_change
213 self.Omega_Z = ((self.z > self.z_change) * f_early * xTot +
214 (self.z <= self.z_change) *
215 (f_late * xTot + (f_early - f_late) * x_change))
216 else:
217 print "f_x_Z = %.3g" % (self.f_x_Z)
218 self.Omega_Z = self.f_x_Z * self.ionHist.xTot
219 return self.Omega_Z
220
221 - def convolve(self):
222 """Convolve delay curve with metal histories.
223 Convolve Omega_Z with delay curve to produce Omega_Z_convolved.
224 """
225 self.Omega_Z_convolved = padconvolve.padded_convolve(self.Omega_Z,
226 self.delay_kernel,
227 origin=0,
228 value=0.0)
229 return self.Omega_Z_convolved
230
231 - def binned_Omega_Z(self, zmin, zmax, useConvolved=True):
232 if useConvolved:
233 o = self.Omega_Z_convolved
234 else:
235 if hasattr(self, 'Omega_Z_unconvolved'):
236 o = self.Omega_Z_unconvolved
237 else:
238 o = self.Omega_Z
239 mask = numpy.logical_and(self.z >= zmin,
240 self.z < zmax)
241 return numpy.mean(o[mask])
242
243 - def fit_f_change(self, z_bin, Omega_Z_obs,
244 min_f_change=0.0,
245 max_f_change=None,
246 useConvolved=True):
247 """Fit f_change to an observed value of Omega_Z_obs in redshift bin
248 z_bin[0]--z_bin[1].
249 """
250 if not hasattr(self, 'original_f_change'):
251 self.original_f_change = self.f_change
252 def fitfunc(fc):
253 self.f_change = fc
254 self.calculate_Omega_Z()
255 if useConvolved:
256 self.convolve()
257 return (self.binned_Omega_Z(z_bin[0], z_bin[1], useConvolved) -
258 Omega_Z_obs)
259
260 if fitfunc(min_f_change) > Omega_Z_obs:
261 print "Warning: can't match Omega_Z with f_change > %.3g" \
262 % (min_f_change)
263 self.f_change = min_f_change
264 return min_f_change
265
266 if max_f_change is None:
267 max_f_change = 3. * self.original_f_change
268 ff = fitfunc(max_f_change)
269 while ff < 0:
270 print " max_f_change = %.3g too small (diff = %.3g)..." \
271 % (max_f_change, ff)
272 max_f_change = 3. * max_f_change
273 ff = fitfunc(max_f_change)
274 print " considering f_change between %.3g and %.3g" % (min_f_change,
275 max_f_change)
276 f_change_fit = scipy.optimize.minpack.bisection(fitfunc,
277 min_f_change,
278 max_f_change,
279 xtol=self.original_f_change*1e-3)
280 if useConvolved:
281 self.f_change = f_change_fit
282 else:
283 self.f_change_unconvolved = f_change_fit
284 self.Omega_Z_unconvolved = self.Omega_Z
285 return f_change_fit
286
287 - def fit_f_x_Z(self, z_bin, Omega_Z_obs,
288 min_f_x_Z=0.0,
289 max_f_x_Z=None,
290 useConvolved=True):
291 """Fit f_x_Z to an observed value of Omega_Z_obs in redshift bin
292 z_bin[0]--z_bin[1].
293 """
294 if not hasattr(self, 'original_f_x_Z'):
295 self.original_f_x_Z = self.f_x_Z
296 def fitfunc(fxc):
297 self.f_x_Z = fxc
298 self.calculate_Omega_Z()
299 if useConvolved:
300 self.convolve()
301 return (self.binned_Omega_Z(z_bin[0], z_bin[1], useConvolved) -
302 Omega_Z_obs)
303
304 if fitfunc(min_f_x_Z) > Omega_Z_obs:
305 print "Warning: can't match Omega_Z with f_x_Z > %.3g" \
306 % (min_f_x_Z)
307 self.f_x_Z = min_f_x_Z
308 return min_f_x_Z
309
310 if max_f_x_Z is None:
311 max_f_x_Z = 3. * self.original_f_x_Z
312 ff = fitfunc(max_f_x_Z)
313 while ff < 0:
314 print " max_f_x_Z = %.3g too small (diff = %.3g)..." \
315 % (max_f_x_Z, ff)
316 max_f_x_Z = 3. * max_f_x_Z
317 ff = fitfunc(max_f_x_Z)
318 print " considering f_x_Z between %.3g and %.3g" % (min_f_x_Z,
319 max_f_x_Z)
320 f_x_Z_fit = scipy.optimize.minpack.bisection(fitfunc,
321 min_f_x_Z,
322 max_f_x_Z,
323 xtol=self.original_f_x_Z*1e-3)
324 self.f_x_Z = f_x_Z_fit
325 return self.f_x_Z
326
328 - def __init__(self, ionHists, delay_kernels, f_x_Zs,
329 z_changes, f_changes,
330 lf_names=None, dl_names=None, en_names=None):
331 self.ionHists = ionHists
332 self.delay_kernels = delay_kernels
333 self.f_x_Zs = f_x_Zs
334 self.z_changes = z_changes
335 self.f_changes = f_changes
336 self.lf_names=lf_names
337 self.dl_names=dl_names
338 self.en_names=en_names
339
340 self.nLF = len(ionHists)
341 self.nDl = len(delay_kernels)
342 self.nEn = len(f_x_Zs)
343
344 EnHists = []
345 for iLF in range(self.nLF):
346 if lf_names is not None:
347 print "Using %s." % lf_names[iLF]
348 square = []
349 for iDl in range(self.nDl):
350 if lf_names is not None:
351 print "Using %s." % dl_names[iDl]
352 line = []
353 for iEn in range(self.nEn):
354 if en_names is not None:
355 print "Using %s." % en_names[iEn]
356 line.append(EnrichmentHistory(ionHists[iLF],
357 f_x_Zs[iEn],
358 z_changes[iEn],
359 f_changes[iEn],
360 delay_kernels[iDl]))
361 square.append(line)
362 EnHists.append(square)
363 self.EnHists = EnHists
364
365 - def fit_f_change(self, z_bin, Omega_Z_obs,
366 min_f_change=0.0, max_f_change=None):
367 f_changes = []
368 f_changes_unconvolved = []
369 for iLF in range(self.nLF):
370 square = []
371 line_unc = []
372 for iDl in range(self.nDl):
373 line = []
374 for iEn in range(self.nEn):
375 if iDl == 0:
376 line_unc.append(\
377 self.EnHists[iLF][iDl][iEn].fit_f_change(z_bin,
378 Omega_Z_obs,
379 min_f_change,
380 max_f_change,
381 useConvolved=False))
382 line.append(\
383 self.EnHists[iLF][iDl][iEn].fit_f_change(z_bin,
384 Omega_Z_obs,
385 min_f_change,
386 max_f_change,
387 useConvolved=True))
388
389 square.append(line)
390 f_changes.append(square)
391 f_changes_unconvolved.append(line_unc)
392 self.f_changes = numpy.array(f_changes)
393 self.f_changes_unconvolved = numpy.array(f_changes_unconvolved)
394 return f_changes, f_changes_unconvolved
395
396 - def fit_f_x_Z(self, z_bin, Omega_Z_obs,
397 min_f_x_Z=0.0, max_f_x_Z=None,
398 useConvolved=True):
399 f_x_Zs = []
400 for iLF in range(self.nLF):
401 square = []
402 for iDl in range(self.nDl):
403 line = []
404 for iEn in range(self.nEn):
405 line.append(\
406 self.EnHists[iLF][iDl][iEn].fit_f_x_Z(z_bin,
407 Omega_Z_obs,
408 min_f_x_Z,
409 max_f_x_Z,
410 useConvolved))
411 square.append(line)
412 f_x_Zs.append(square)
413 self.f_x_Zs = f_x_Zs
414 return f_x_Zs
415
416 - def convolve(self):
417 Omega_Zs_convolved = []
418 for iLF in range(self.nLF):
419 square = []
420 for iDl in range(self.nDl):
421 line = []
422 for iEn in range(self.nEn):
423 line.append(\
424 self.EnHists[iLF][iDl][iEn].convolve())
425 square.append(line)
426 Omega_Zs_convolved.append(square)
427 self.Omega_Zs_convolved = Omega_Zs_convolved
428 return self.Omega_Zs_convolved