1 """Implements stellar initial mass functions from Romano et al. 2005.
2
3 Full Citation: Romano D., Chiappini C., Matteucci F., Tosi M., 2005,
4 A\&A, 430, 491 (2005A&A...430..491R)
5
6 http://adsabs.harvard.edu/abs/2005A&A...430..491R
7 """
8
9 import math
10
11 import numpy
12 import scipy
13 import scipy.integrate
14
15 import cosmolopy.utils as utils
16 from cosmolopy.utils import Normalize
17
18
19
21 """Convert an IMF by mass to an IMF by number.
22 """
23 def imf_number_function(mass):
24 return imf_mass_function(mass)/mass
25 imf_number_function.__name__ = \
26 imf_mass_function.__name__.replace('imf_mass_',
27 'imf_number_')
28 imf_number_function.__dict__.update(imf_mass_function.__dict__)
29 imf_number_function.__doc__ = imf_mass_function.__doc__
30 return imf_number_function
31
33 """Return a function consisting of seperate powerlaws.
34
35
36 F(x) = a_i x^{p_i} for l_i < x < l_{i+1}
37
38 Parameters
39 ----------
40
41 limits: array (length n+1)
42 boundaries of the specified powerlaws. Must be one greater in
43 length than coefficents and powers. Specify -numpy.infty for
44 the first limit or numpy.infty for the last limit for
45 unbounded powerlaws.
46
47 coefficients: array (length n)
48 values of the coefficient a_i
49
50 powers: array (length n)
51 values of the powerlaw indices p_i
52
53 The returned function takes a single, one-dimensional array of
54 values on which to operate.
55
56 """
57
58 limits = numpy.atleast_1d(limits)
59 coefficients = numpy.atleast_1d(coefficients)
60 powers = numpy.atleast_1d(powers)
61 for array in [limits, coefficients, powers]:
62 if array.ndim > 1:
63 raise ValueError("arguments must be a 1D arrays or scalars.")
64
65 if not len(coefficients) == len(powers):
66 raise ValueError("coefficients and powers must be the same length.")
67
68 if not len(limits) == len(powers)+1:
69 raise ValueError("limits must be one longer than powers.")
70
71 limits = limits.reshape((-1,1))
72 coefficients = coefficients.reshape((-1,1))
73 powers = powers.reshape((-1,1))
74
75 def pf_func(x):
76 x = numpy.atleast_1d(x)
77 if x.ndim > 1:
78 raise ValueError("argument must be a 1D array or scalar.")
79 y = numpy.sum((coefficients * x**powers) *
80 (x > limits[0:-1]) * (x < limits[1:]),
81 axis=0)
82 y[x < limits[0]] = None
83 y[x > limits[-1]] = None
84 return y
85 return pf_func
86
87
88
89
90 coeff_Scalo86 = numpy.array([0.19, 0.24])
91 pow_Scalo86 = numpy.array([-1.35, -1.70])
92 lims_Scalo86 = numpy.array([0, 2., numpy.infty])
93 _imf_mass_Scalo86 = broken_powerlaw_function(lims_Scalo86,
94 coeff_Scalo86,
95 pow_Scalo86)
96
98 """Scalo 1986 stellar Initial Mass function by mass via Romano et al. (2005).
99
100 Parameters
101 ----------
102
103 mass: 1d array
104 mass in Solar masses.
105
106 """
107 return _imf_mass_Scalo86(mass)
108
109
110 coeff_Scalo98 = numpy.array([0.39, 0.39, 0.16])
111 pow_Scalo98 = numpy.array([-0.2, -1.7, -1.3])
112 lims_Scalo98 = numpy.array([0, 1., 10., numpy.infty])
113 _imf_mass_Scalo98 = broken_powerlaw_function(lims_Scalo98,
114 coeff_Scalo98,
115 pow_Scalo98)
116
118 """Scalo 1998 stellar Initial Mass function by mass via Romano et al. (2005).
119
120 Parameters
121 ----------
122
123 mass: 1d array
124 mass in Solar masses.
125
126 """
127 return _imf_mass_Scalo98(mass)
128
129
130
131 coeff_Tinsley = numpy.array([0.21, 0.26, 2.6])
132 pow_Tinsley = numpy.array([-1., -1.3, -2.3])
133 lims_Tinsley = numpy.array([0, 2., 10., numpy.infty])
134 _imf_mass_Tinsley = broken_powerlaw_function(lims_Tinsley,
135 coeff_Tinsley,
136 pow_Tinsley)
137
139 """Tinsley stellar Initial Mass function by mass via Romano et al. (2005).
140
141 Parameters
142 ----------
143
144 mass: 1d array
145 mass in Solar masses.
146
147 """
148 return _imf_mass_Tinsley(mass)
149
150
151
152 coeff_Kroupa = numpy.array([0.58, 0.31, 0.31])
153 pow_Kroupa = numpy.array([-0.3, -1.2, -1.7])
154 lims_Kroupa = numpy.array([0, 0.5, 1., numpy.infty])
155 _imf_mass_Kroupa = broken_powerlaw_function(lims_Kroupa,
156 coeff_Kroupa,
157 pow_Kroupa)
158
160 """Kroupa stellar Initial Mass function by mass via Romano et al. (2005).
161
162 Parameters
163 ----------
164
165 mass: 1d array
166 mass in Solar masses.
167
168 """
169 return _imf_mass_Kroupa(mass)
170
171
173 """Salpeter stellar Initial Mass function by mass via Romano et al. (2005).
174
175 """
176 return 0.17 * mass**-1.35
177
178 mc = 0.079
179 logmc = math.log10(mc)
180 sigma = 0.69
181 twosigsq = 2. * sigma**2
182
183 fa = 0.85
184 fb = 0.24
185
187 return fa * numpy.exp(-1 * ((numpy.log10(mass) -
188 logmc)**2. /
189 twosigsq))
191 return fb * mass**-1.3
193 """Chabrier stellar Initial Mass function by mass via Romano et al. (2005).
194
195 With slope -1.3 above 1 MSun.
196
197 """
198 if numpy.isscalar(mass):
199 if mass < 1.:
200 return _imf_lowmass_Chabrier(mass)
201 else:
202 return _imf_highmass_Chabrier(mass)
203 mass = numpy.asarray(mass)
204 imf = numpy.empty(mass.shape)
205 lowmask = mass < 1.
206 imf[lowmask] = _imf_lowmass_Chabrier(mass[lowmask])
207 imf[mass >= 1.] = _imf_highmass_Chabrier(mass[mass >= 1.])
208
209 return imf
210
211
212 imf_number_Chabrier = by_number(imf_mass_Chabrier)
213
214
215
217 """Plot and test the normalization of imf_function."""
218 label = imf_function.__name__.replace('imf_', '')
219 label = label.replace('_', ' ')
220 print label
221 dm = 0.005
222 mmin = 0.1
223 mmax = 100.
224 mass = numpy.arange(mmin, mmax + 1.1 * dm, dm)
225 imf = imf_function(mass)
226
227
228 liimf = utils.logquad(imf_function, mmin, mmax)
229 print liimf[0], ", ",
230 print 1.-liimf[0]
231 assert abs(liimf[0] - 1) < 1e-10
232
233
234 cimf = scipy.integrate.cumtrapz(imf * dm)
235 print cimf[-1], ", ",
236 print 1.-cimf[-1]
237 assert abs(cimf[-1] - 1) < 1e-3
238
239 pylab.subplot(121)
240 l = pylab.plot(mass, imf, label=label)
241 pylab.plot(mass[1:], cimf, l[0].get_c())
242 pylab.gca().set_yscale('log')
243 pylab.legend(loc='best')
244
245 pylab.subplot(122)
246 pylab.plot(mass, imf, l[0].get_c(), label=label)
247 pylab.plot(mass[1:], cimf, l[0].get_c())
248 pylab.gca().set_yscale('log')
249 pylab.gca().set_xscale('log')
250
252
253 import sys
254 module = sys.modules[__name__]
255 functionlist = [getattr(module, name) for name in dir(module) if name.startswith('imf_mass_')]
256
257 import pylab
258 import scipy
259 import scipy.integrate
260 pylab.figure(figsize=(11,6))
261
262
263 for function in functionlist:
264 test_plot_norm_imf(Normalize(0.1,100.)(function))
265
266
267 pylab.figure(figsize=(11,6))
268 for function in functionlist:
269 test_plot_norm_imf(Normalize(0.1,100.)(by_number(function)))
270
271 if __name__ == '__main__':
272 import matplotlib.pyplot as pylab
273 import os
274 import sys
275
276
277 if len(sys.argv)==1:
278 print "Run with a filename argument to produce image files, e.g.:"
279 print " python initialmassfunction.py imf.png"
280 print " python initialmassfunction.py imf.eps"
281 if len(sys.argv) > 1:
282 filename = sys.argv[1]
283 prefix, extension = os.path.splitext(filename)
284 else:
285 filename = None
286
287
288
289 test_plot_norm_all()
290
291
292 if filename is None:
293 pylab.show()
294 else:
295 from matplotlib import _pylab_helpers
296 for manager in _pylab_helpers.Gcf.get_all_fig_managers():
297 fig = manager.canvas.figure
298 if len(fig.get_label()) > 0:
299 label = fig.get_label()
300 else:
301 label = '_Fig' + str(manager.num)
302 newfilename = prefix + '_' + label + extension
303 fig.savefig(newfilename, bbox_inches="tight")
304