enrichpy.delay

1 """Class Delay represents the delay distribution for metals emerging from stars. 2 3 """ 4 import math 5 import sys 6 import cPickle 7 import os 8 9 import scipy 10 import scipy.integrate 11 import scipy.optimize 12 import numpy 13 14 import cosmolopy.utils as utils 15 import yields 16 import initialmassfunction 17 import lifetime 18

19 -class Delay:

20 """Represents the delay distribution for metals emerging from stars. 21 22 Takes the IMF, stellar metal yield function, and stellar lifetime 23 function and calculates various the distribution of ejected metal 24 mass with respect to stellar mass and stellar lifetime. 25 """ 26

27 - def __init__(self, metallicity, 28 imfFunc = initialmassfunction.imf_number_Chabrier, 29 yieldClass = yields.Data_GBM, 30 ejections=False, 31 lifetimeFunc = lifetime.main_sequence_life_K97, 32 lowMass = None, 33 highMass = None, 34 minTime = None, 35 maxTime = None, 36 ):

37 """Initialize the distribution. 38 39 Parameters 40 ---------- 41 42 metallicity: scalar 43 metallicity value for the yield calculation. 44 45 imfFunc: function 46 should return IMF by number given stellar mass in Msun. 47 48 yieldClass: 49 class encapsulating the yield data (see yields.py) 50 51 lifetimeFunc: function 52 should return lifetime of a star given mass in Msun. 53 54 lowMass, highMass: float 55 define the interval over which to normalize the distributions. 56 Default to the range of masses in yieldClass().mVals 57 58 """ 59 60 self.metallicity = metallicity 61 self.yieldClass = yieldClass 62 self.lifetimeFunc = lifetimeFunc 63 self.ejections = ejections 64 65 ## Stellar mass as a function of lifetime. ## 66 self._m_star_func = \ 67 lifetime.mass_from_main_sequence_life_function(self.lifetimeFunc) 68 69 ### Initialize some functions that require interpolation grids. ### 70 71 ## Yield function. 72 self._yields = yieldClass() 73 self._yieldFunc_2d = \ 74 self._yields.interpolate_function(ejections=ejections) 75 76 ## Normalization constant (determined later) to turn yieldFunc 77 ## * imf into the fraction of carbon produced by stars between 78 ## M and M+dm. 79 self._fCnorm = 1.0 80 81 ## Convert minTime and maxTime to lowMass and highMass. ## 82 if minTime is not None: 83 highMass = self.m_star(minTime) 84 if maxTime is not None: 85 lowMass = self.m_star(maxTime) 86 87 ## Choose low and high masses using the range of yield data.## 88 self.lowMass = lowMass 89 if lowMass is None: 90 self.lowMass = numpy.min(self._yieldFunc_2d.mVals) 91 self.highMass = highMass 92 if highMass is None: 93 self.highMass = numpy.max(self._yieldFunc_2d.mVals) 94 95 ## Calculate the equivalent range of stellar lifetimes. ## 96 self.maxTime = self.t_life(self.lowMass) 97 self.minTime = self.t_life(self.highMass) 98 99 print "Normalizing delay function from %.3g -- %.3g Msun or " %\ 100 (self.lowMass, self.highMass) 101 print " from %.3g -- %.3g years." %\ 102 (self.maxTime, self.minTime) 103 104 ## Normalize the IMF. ## 105 self.imfFunc = utils.Normalize(self.lowMass, self.highMass)(imfFunc) 106 107 ## Normalize the CCDF to unity at lowMass. ## 108 # Use an array of masses because the integration works better 109 # when split into pieces. 110 marray = [] 111 if (lowMass is None) and (highMass is None): 112 # Use masses at which yield is defined (rather than 113 # interpolated). 114 marray = self._yieldFunc_2d.mVals 115 else: 116 # Or just use an array of values covering the given mass range. 117 marray = numpy.linspace(self.lowMass, self.highMass, 1000.) 118 ccdf = self.metal_mass_CCDF(marray) 119 self._fCnorm = ccdf[0] 120 print " delay curve normalization = %.3g" % self._fCnorm 121 122 mass_limit_error = self.mass_limit_error() 123 print "Estimated error from exclusion of m_star > %.3g Msun = %.3g" % \ 124 (self.highMass, mass_limit_error) 125 self.printstats() 126 sys.stdout.flush()

127 128 # def dumb_pickle_filter(self, prefix='temp_pic_test'): 129 # """Filter out attributes that can't be pickled. 130 # """ 131 # picname = prefix + '.pkl' 132 # picfile = open(picname, 'w') 133 # for k, v in self.__dict__.items(): 134 # # Avoid self references (and thereby infinite recurion). 135 # if v is self: 136 # delattr(self, k) 137 # continue 138 # # Remove any attributes that can't be pickled. 139 # try: 140 # cPickle.dump(v, picfile) 141 # except TypeError: 142 # print "Can't pickle: ", 143 # print k, ": ", 144 # print type(v), " " 145 # delattr(self, k) 146 # picfile.close() 147 # os.remove(picname) 148 149 # def __getstate__(self): 150 # """Prepare a state of pickling.""" 151 # self.dumb_pickle_filter() 152 # return self.__dict__ 153

154 - def mass_limit_error(self, m_max=None, imf_index=-2.3):

155 """Estimate the error incurred by excluding stars above m_max. 156 157 m_max defaults to self.highMass. 158 159 imf_index is the powerlaw index of the imf at high masses. 160 """ 161 if m_max is None: 162 m_max = self.highMass 163 164 (exfuncHigh, exfuncLow, 165 highSlope, highIntercept, lowSlope, lowIntercept) = \ 166 yields.extrapolate_function(self._yieldFunc_2d, 167 self.metallicity, 168 return_coeffs=True) 169 170 phi0 = self.imf(m_max) * (m_max**2.3) 171 phi0 *= self._fCnorm 172 plustwo = imf_index + 2. 173 plusone = imf_index + 1. 174 175 ## The error is simply the integral of the linear yield times 176 ## the powelaw IMF from m_max to infinity. 177 error = ((highSlope * phi0 * m_max**(plustwo) / (-1.* plustwo)) + 178 (highIntercept * phi0 * m_max**(plusone) / (-1.* plusone))) 179 return error

180

181 - def metal_mass_dist_postinterp(self, m):

182 """Like metal_mass_dist, but with interpolation after IMF 183 multiplication. 184 """ 185 if not hasattr(self, '_fCfunc_2d'): 186 self._fCfunc_2d = \ 187 self.yieldClass(interp_args={'weight_function':self.imf}, 188 ).interpolate_function(ejections=self.ejections) 189 return (self._fCfunc_2d(m, self.metallicity).flatten() / self._fCnorm)

190

191 - def metal_mass_dist(self, m):

192 """Fraction of total carbon yield from stars of mass m, 193 dfC/dm, i.e. the distribution of carbon production as a 194 function of stellar mass. 195 196 This is weighted by the stellar initial mass function (IMF), 197 so that metal_mass_dist(m) * dm is the fraction of carbon 198 emitted by stars with masses between m and m+dm. 199 200 Note that this is normalized over the interval from 201 self.lowMass to self.highMass. 202 """ 203 dist = (self.imf(m) * 204 self._yieldFunc_2d(m, self.metallicity).flatten() / 205 self._fCnorm) 206 #dist[m<self.lowMass] = None 207 #dist[m>self.highMass] = None 208 return dist

209

210 - def yield_mass(self, m):

211 """Carbon yield from a star of mass m. 212 213 This is the total mass of carbon emitted by a single star of 214 mass m over its lifetime. 215 216 """ 217 return self._yieldFunc_2d(m, self.metallicity).flatten()

218

219 - def metal_mass_CCDF(self, m, highMass=None, method='romberg', **kwargs):

220 """Complementary cumulative distribution of metal yield by mass, 221 i.e. the fraction of the total metal output produced by stars 222 above mass m. 223 224 This is the integral of metal_mass_dist from m to highMass (or 225 self.highMass). 226 227 Note that this is normalized over the interval from 228 self.lowMass to self.highMass, so values may be negative above 229 highMass and greater than one below lowMass. 230 231 kwargs get passed to utils.ccumulate, which may passes them to 232 utils.integrate_piecewise, which may pass them to the 233 integration routine. 234 235 See also: metal_mass_dist. 236 237 """ 238 if highMass is None: 239 highMass = self.highMass 240 ccdf = utils.ccumulate(self.metal_mass_dist, m, max=highMass, 241 method=method, 242 **kwargs) 243 ccdf[m<self.lowMass] = numpy.nanmax(ccdf) 244 ccdf[m>self.highMass] = 0.0 245 return ccdf

246 247

248 - def metal_mass_dist_binned(self, m, ccdf=None, **kwargs):

249 """Calculate the values of a binned enrichment function. 250 251 Each returned value bdf[i] is the integral of the mass dist 252 function from m[i-1] to m[i]. bdf[0] is the integral from 253 self.lowMass to m[0]. 254 255 Also returned are delta_m values (the width of each bin in 256 units of mass), so the average derivative can be calculated, 257 if desired. 258 259 Note that this is normalized over the interval from 260 self.lowMass to self.highMass. 261 262 See also: metal_mass_CCDF 263 """ 264 265 if ccdf is None: 266 # Cumulative delay function. 267 cdf = 1. - numpy.nan_to_num(self.metal_mass_CCDF(m, **kwargs)) 268 else: 269 cdf = 1. - numpy.nan_to_num(ccdf) 270 271 # Break up into bins. 272 bdf = numpy.empty(cdf.shape) 273 bdf[0] = cdf[0] 274 bdf[1:] = cdf[1:] - cdf[:-1] 275 276 # Width of the bins in time. 277 dm = numpy.empty(cdf.shape) 278 dm[0] = m[0] - self.lowMass 279 dm[1:] = m[1:] - m[:-1] 280 return bdf, dm

281

282 - def metal_time_CDF(self, t, ccdf=None, return_mass=False, **kwargs):

283 """The delay CDF of metal yield as a function of stellar age. 284 285 I.e. the fraction of the total metal yield produced by stars 286 with a lifetime of t or less. 287 288 Note that this is normalized over the interval from 289 self.lowMass to self.highMass, so values may be negative above 290 the lifetime corresponding to highMass and greater than one 291 below the lowMass lifetime. 292 293 Calculated using metal_mass_CCDF and m_star. 294 295 kwargs get passed to metal_mass_CCDF. 296 """ 297 m = self.m_star(t) 298 if ccdf is None: 299 cdf = self.metal_mass_CCDF(m, **kwargs) 300 else: 301 cdf = ccdf 302 303 if return_mass: 304 return cdf, m 305 return cdf

306

307 - def metal_time_dist_binned(self, t, ccdf=None, **kwargs):

308 """Calculate the values of a binned delay function. 309 310 Each returned value bdf[i] is the integral of the delay 311 function from t[i-1] to t[i]. bdf[0] is the integral from 312 self.minTime to t[0]. 313 314 Also returned are delta_t values (the width of each bin in 315 units of time), so the average derivative can be calculated, 316 if desired. 317 318 Note that this is normalized over the interval from 319 self.lowMass to self.highMass. 320 321 See also: metal_time_CDF 322 """ 323 324 if ccdf is None: 325 # Cumulative delay function. 326 cdf = numpy.nan_to_num(self.metal_time_CDF(t, **kwargs)) 327 else: 328 cdf = ccdf 329 # Break up into bins. 330 bdf = numpy.empty(cdf.shape) 331 bdf[0] = cdf[0] 332 bdf[1:] = cdf[1:] - cdf[:-1] 333 334 # Width of the bins in time. 335 dt = numpy.empty(cdf.shape) 336 dt[0] = t[0] - self.minTime 337 dt[1:] = t[1:] - t[:-1] 338 return bdf, dt

339

340 - def metal_time_dist_unbinned(self, t, **kwargs):

341 """Delay function inferred from metal_time_dist_binned. 342 343 Simply returns the bin values divided by the bin size. Note 344 that the first bin is from self.minTime to t[0]. 345 346 Notes 347 ----- 348 349 This is the most accurate way to get the distribution of metal 350 emission delay times. A more analytical treatment would be to 351 convert directly from the continuous metal-mass distribution, 352 but this involves the derivative of the stellar lifetime as a 353 funtion of mass, which can be numerically problematic. 354 355 Note that this is normalized over the interval from 356 self.lowMass to self.highMass. 357 """ 358 359 bdf, dt = self.metal_time_dist_binned(t, **kwargs) 360 df = bdf/dt 361 return df, dt

362

363 - def imf(self, mass):

364 """The initial mass function. 365 366 Note that the normalization may include stars that do not 367 contribute to the yield of metals. 368 """ 369 return self.imfFunc(mass)

370

371 - def m_star(self, t_life):

372 """Stellar mass as a function of lifetime. """ 373 return self._m_star_func(t_life)

374

375 - def t_life(self, m_star):

376 """Stellar lifetime as a function of mass. """ 377 return self.lifetimeFunc(m_star)

378

379 - def stats(self, meantol=1e-3):

380 """Calculate median and mean delay.""" 381 382 npieces = 1./meantol 383 tarray = numpy.logspace(numpy.log10(self.minTime), 384 numpy.log10(self.maxTime), 385 npieces) 386 marray = self.m_star(tarray) 387 388 ccdf = self.metal_mass_CCDF(marray) 389 cdf = 1. - ccdf 390 391 # Find median: 392 medindex = numpy.max(numpy.where(ccdf<0.5)) 393 # Ambiguity of the exact index: 394 mediana = tarray[medindex] 395 medianb = tarray[medindex+1] 396 medianMassa = marray[medindex] 397 medianMassb = marray[medindex+1] 398 399 ## Calculate mean time. ### 400 dist, dt = self.metal_time_dist_binned(tarray, ccdf=ccdf) 401 tarrayb = numpy.empty(tarray.shape) 402 tarrayb[0] = 0 403 tarrayb[1:] = tarray[:-1] 404 meana = numpy.sum(dist * tarray) / numpy.sum(dist) 405 meanb = numpy.sum(dist * tarrayb) / numpy.sum(dist) 406 407 ## Calculate mean mass. ### 408 marray = marray[::-1] 409 ccdf = ccdf[::-1] 410 dist, dm = self.metal_mass_dist_binned(marray, ccdf=ccdf) 411 marrayb = numpy.empty(marray.shape) 412 marrayb[0] = 0 413 marrayb[1:] = marray[:-1] 414 meanMassa = numpy.sum(dist * marray) / numpy.sum(dist) 415 meanMassb = numpy.sum(dist * marrayb) / numpy.sum(dist) 416 return (mediana, (medianb-mediana)/min(mediana,medianb), 417 medianMassb, (medianMassa-medianMassb)/min(medianMassb,mediana), 418 meana, 419 (meanb-meana)/min(meana,meanb), 420 meanMassa, 421 (meanMassb-meanMassa)/min(meanMassa,meanMassb))

422

423 - def printstats(self):

424 median, medianerr, medianMass, medianMasserr, mean, meanerr, meanMass, meanMasserr = self.stats() 425 print "Median and mean delays are %.4g and %.4g Gyr" % (median/1e9, 426 mean/1e9) 427 print " (frac. med err = %.3g)" % medianerr 428 print " (frac. mean err = %.3g)" % meanerr 429 430 print " Mass from mean delay is %.3g Msun" % self.m_star(mean) 431 print "Median and mean masses are %.4g and %.4g Msun" % (medianMass, 432 meanMass) 433 print " (frac. med err = %.3g)" % medianMasserr 434 print " (frac. mean err = %.3g)" % meanMasserr 435 print " Age from mean mass is %.3g Gyr" % (self.t_life(meanMass)/1e9)

436 437 438

439 -class MockDelay(Delay):

440 _fCnorm = numpy.nan

441 - def __init__(self, binmasses, fractions, 442 lifetimeFunc = lifetime.main_sequence_life_K97, 443 maxTime=None, 444 **args):

445 """ Mimics the Delay class. 446 """ 447 448 self.lifetimeFunc = lifetimeFunc 449 ## Stellar mass as a function of lifetime. ## 450 self._m_star_func = \ 451 lifetime.mass_from_main_sequence_life_function(self.lifetimeFunc) 452 453 454 if maxTime is not None: 455 lowMass = self.m_star(maxTime) 456 print "Resetting low mass limit to %.3g (from %.3g)." % \ 457 (lowMass, binmasses[0]) 458 binmasses[0] = lowMass 459 460 if not sum(fractions) == 1: 461 norm = 1./sum(fractions) 462 print "Normalizing fractions with a factor of %.3g" % norm 463 fractions = fractions * norm 464 465 assert len(fractions) == len(binmasses) - 1,\ 466 "len(fractions) must be one less than len(splitmasses)." 467 assert numpy.all(numpy.sort(binmasses) == binmasses),\ 468 "binmasses must be ordered from low to high." 469 470 self.binmasses = binmasses 471 self.fractions = fractions 472 473 ## Choose low and high masses using the range of yield data.## 474 self.lowMass = numpy.min(self.binmasses) 475 self.highMass = numpy.max(self.binmasses) 476 477 self.minTime = self.t_life(self.highMass) 478 self.maxTime = self.t_life(self.lowMass) 479 self.printstats() 480 sys.stdout.flush()

481

482 - def metal_mass_dist(self, m):

483 # Cumulative delay function. 484 ccdf = numpy.nan_to_num(self.metal_mass_CCDF(m)) 485 cdf = numpy.nanmax(ccdf) - ccdf 486 487 # Break up into bins. 488 bdf = numpy.empty(cdf.shape) 489 bdf[0] = cdf[0] 490 bdf[1:] = cdf[1:] - cdf[:-1] 491 492 # Width of the bins in time. 493 dm = numpy.empty(cdf.shape) 494 dm[0] = m[0] - self.lowMass 495 dm[1:] = m[1:] - m[:-1] 496 497 df = bdf/dm 498 return df

499

500 - def metal_mass_CCDF(self, mass):

501 ccdf = numpy.zeros(mass.shape) 502 503 # Step through the fraction indices in reverse 504 for ibin in range(len(self.fractions)-1, -1, -1): 505 # Add the current fraction to all lower bins. 506 ccdf[mass < self.binmasses[ibin]] += self.fractions[ibin] 507 508 # And add the correct portion of the current fraction to 509 # points in the current bin. 510 binmask = numpy.logical_and(mass >= self.binmasses[ibin], 511 mass < self.binmasses[ibin+1]) 512 binwidth = self.binmasses[ibin+1] - self.binmasses[ibin] 513 ccdf[binmask] += (self.fractions[ibin] * 514 (self.binmasses[ibin+1] - mass[binmask]) / 515 binwidth) 516 return ccdf

517 - def metal_mass_dist_postinterp(self, m):

518 return self.metal_mass_dist(m)

519

520 - def imf(self, mass):

521 return numpy.nan * mass

522

523 - def yield_mass(self, m):

524 return numpy.nan * m

525

526 -def test_plot_binned_delay_function():

527 ddist = Delay(0.02) 528 t = numpy.arange(0., 5.e10, 1.33e6) 529 bdf, dt = ddist.metal_time_dist_binned(t) 530 cbdf = numpy.cumsum(bdf) 531 print "Sum of binned function = ", numpy.sum(bdf) 532 pylab.subplot(221) 533 pylab.plot(t, t * bdf/dt,'.', alpha=0.2) 534 pylab.subplot(222) 535 pylab.plot(t, cbdf,'.', alpha=0.2)

536

537 -def test_plot_delay_function(dlogt=1e-3):

538 Z = 0.02 539 ftfunc = delay_function(0.02) 540 t = 10**numpy.arange(6.0,10.2,dlogt) 541 dt = t[1:] - t[:-1] 542 ft = ftfunc(t[:-1], dt) 543 544 ft2 = ft 545 dlnt = numpy.log(t[1:]) - numpy.log(t[:-1]) 546 ft2[numpy.logical_not(numpy.isfinite(ft))] = 0 547 ftcum = scipy.integrate.cumtrapz(ft2 * dt) 548 549 lifetimeFunc = lifetime.main_sequence_life_K97 550 m_star_func = \ 551 lifetime.mass_from_main_sequence_life_function(lifetimeFunc, 552 logmmin=10.) 553 m = m_star_func(t[:-1]) 554 dmdt = (m - m_star_func(t[:-1] + dt))/dt 555 556 # print 'ft=',ft 557 Z = 0.02 558 559 import pylab 560 #pylab.figure() 561 pylab.subplot(221) 562 pylab.plot(t[:-1], t[:-1]*ft, alpha=0.7) 563 ax = pylab.gca() 564 ax.set_xscale('log') 565 pylab.xlim(t.min(), t.max(), alpha=0.7) 566 567 #pylab.twinx() 568 pylab.subplot(222) 569 pylab.plot(t[1:-1], ftcum, alpha=0.7) 570 if dlogt > 1e-2: 571 pylab.plot(t[1:-1], ftcum, '.', alpha=0.7) 572 pylab.plot(t[1:-1], 1.0 - ftcum, alpha=0.7) 573 pylab.axhline(y=0.5) 574 ax = pylab.gca() 575 ax.set_xscale('log') 576 pylab.xlim(t.min(), t.max(), alpha=0.7) 577 #pylab.gca().set_yscale('log') 578 579 pylab.subplot(223) 580 pylab.plot(t[:-1], m, alpha=0.7) 581 pylab.gca().set_xscale('log') 582 pylab.xlim(t.min(), t.max(), alpha=0.7) 583 pylab.gca().set_yscale('log') 584 pylab.axhline(y=100.0) 585 #pylab.twinx() 586 587 pylab.subplot(224) 588 pylab.plot(t[:-1], dmdt, alpha=0.7) 589 pylab.gca().set_xscale('log') 590 pylab.gca().set_yscale('log') 591 pylab.xlim(t.min(), t.max(), alpha=0.7)

592

593 -def test_plot_converge_delay_function():

594 test_plot_delay_function(1.e-4) 595 test_plot_delay_function(1.e-3) 596 test_plot_delay_function(1.e-2)

597 #test_plot_binned_delay_function() 598

599 -def plot_delay_function(metallicity, 600 imfFunc = initialmassfunction.imf_number_Chabrier, 601 yieldClass = yields.Data_GBM, 602 lifetimeFunc = lifetime.main_sequence_life_K97):

603 #lifetimeFunc = lifetime.main_sequence_life_MM89): 604 """ Plot delay function, yield, IMF, and lifetime. 605 606 """ 607 608 dm = 0.01 609 mass = numpy.arange(0, 110., dm) 610 611 time = 10**numpy.arange(6., 10.5, 0.01) 612 dt = time[1:] - time[:-1] 613 time = time[1:] 614 615 m_star_func = \ 616 lifetime.mass_from_main_sequence_life_function(lifetimeFunc) 617 mass_from_time = m_star_func(time) 618 619 imf = imfFunc(mass) 620 621 mCfunc_weighted = yieldClass(weight_function=imfFunc) 622 mCfunc_noweight = yieldClass() 623 624 mC_weighted = mCfunc_weighted(mass, metallicity) 625 mC_noweight = mCfunc_noweight(mass, metallicity) 626 627 tLife = lifetimeFunc(mass) 628 dmdt = (mass[1:] - mass[:-1])/(tLife[1:] - tLife[:-1]) 629 630 ddist = Delay(0.02, 631 imfFunc = imfFunc, 632 yieldClass = yieldClass, 633 lifetimeFunc = lifetimeFunc) 634 fdelay = ddist.metal_time_dist_unbinned(time) 635 fdelay[numpy.isnan(fdelay)] = 0.0 636 cumfdelay = numpy.cumsum(fdelay * dt) 637 638 print mass_from_time 639 print "Min time = %.3g" % min(time[numpy.isfinite(cumfdelay)]) 640 print "Max time = %.3g " % max(time[numpy.isfinite(cumfdelay)]) 641 print "Max mass = %.5g" % numpy.nanmax(mass_from_time[numpy.isfinite(cumfdelay)]) 642 print "Min mass = %.5g" % numpy.nanmin(mass_from_time[numpy.isfinite(cumfdelay)]) 643 print "Min cdf = %.5g " % min(cumfdelay[numpy.isfinite(cumfdelay)]) 644 print "Max cdf = %.5g" % max(cumfdelay[numpy.isfinite(cumfdelay)]) 645 646 delay_fig = pylab.figure(figsize=(11,6)) 647 delay_fig.set_label('vsMass') 648 649 pylab.subplot(221) 650 pylab.plot(mass, imf) 651 pylab.xlabel('mass (Msun)') 652 pylab.ylabel('IMF') 653 pylab.gca().set_yscale('log') 654 655 pylab.subplot(222) 656 pylab.plot(mass, mC_noweight, label='unweighted', ls='--') 657 pylab.xlabel('mass (Msun)') 658 pylab.ylabel(r'$m_C$ unweighted (--)') 659 pylab.gca().set_yscale('log') 660 #pylab.legend(loc='best') 661 662 pylab.twinx() 663 pylab.plot(mass, mC_weighted, label='IMF weighted', ls=':') 664 pylab.ylabel(r'$m_C$ weighted (:)') 665 pylab.gca().set_yscale('log') 666 #pylab.legend(loc='best') 667 668 pylab.subplot(223) 669 pylab.plot(mass, tLife) 670 pylab.xlabel('mass (Msun)') 671 pylab.ylabel('tLife (years)') 672 pylab.gca().set_yscale('log') 673 674 pylab.twinx() 675 pylab.plot(mass[1:], -1. * dmdt, label='-dmdt', ls=':') 676 pylab.gca().set_yscale('log') 677 pylab.ylabel(r'dmdt(:)') 678 679 pylab.subplot(224) 680 pylab.plot(mass_from_time, fdelay) 681 pylab.gca().set_yscale('log') 682 pylab.xlabel('mass (Msun)') 683 pylab.ylabel('fdelay') 684 685 pylab.twinx() 686 pylab.plot(mass_from_time, cumfdelay) 687 pylab.ylabel('cumfdelay') 688 pylab.gca().set_yscale('log') 689 pylab.xlim(numpy.min(mass), numpy.max(mass)) 690 691 pylab.subplots_adjust(wspace=0.5, hspace=0.29) 692 693 delay_fig_time = pylab.figure(figsize=(11,6)) 694 delay_fig_time.set_label('vsTime') 695 696 pylab.subplot(221) 697 pylab.plot(tLife, imf) 698 pylab.xlabel('tLife (yr)') 699 pylab.ylabel('IMF') 700 pylab.gca().set_yscale('log') 701 pylab.gca().set_xscale('log') 702 703 pylab.twinx() 704 pylab.plot(tLife[1:], -1. * tLife[1:] * imf[1:] * dmdt, ls=':') 705 pylab.ylabel('IMF * t * -dm/dt') 706 pylab.gca().set_yscale('log') 707 708 pylab.subplot(222) 709 pylab.plot(tLife, mC_noweight, label='unweighted', ls='--') 710 pylab.xlabel('tLife (yr)') 711 pylab.ylabel(r'$m_C$ unweighted (--)') 712 pylab.gca().set_yscale('log') 713 pylab.gca().set_xscale('log') 714 #pylab.legend(loc='best') 715 716 pylab.twinx() 717 pylab.plot(tLife, mC_weighted, label='IMF weighted', ls=':') 718 pylab.ylabel(r'$m_C$ weighted (:)') 719 pylab.gca().set_yscale('log') 720 #pylab.legend(loc='best') 721 722 pylab.subplot(223) 723 pylab.plot(tLife, mass) 724 pylab.xlabel('tLife (yr)') 725 pylab.ylabel('mass (Msun)') 726 pylab.gca().set_yscale('log') 727 pylab.gca().set_xscale('log') 728 729 pylab.twinx() 730 pylab.plot(tLife[1:], -1. * dmdt, label='-dmdt', ls=':') 731 pylab.gca().set_yscale('log') 732 pylab.ylabel(r'-dmdt(:)') 733 734 pylab.subplot(224) 735 pylab.plot(time, time * fdelay) 736 #pylab.gca().set_yscale('log') 737 pylab.gca().set_xscale('log') 738 pylab.xlabel('tLife (years)') 739 pylab.ylabel('t * fdelay') 740 pylab.subplots_adjust(wspace=0.5, hspace=0.29) 741 742 pylab.twinx() 743 pylab.plot(time, cumfdelay) 744 pylab.ylabel('cumfdelay')

745 #pylab.gca().set_yscale('log') 746

747 -def test_Delay(plots=True, skip_asserts=False, delayfunc=None):

748 import numpy.testing.utils as ntest 749 750 maxTime = delayfunc.maxTime 751 minTime = delayfunc.minTime 752 753 # Set up mass and time coordinates. 754 npoints = 1000 755 #dmass = 0.01 756 #mass = numpy.arange(delayfunc.lowMass, delayfunc.highMass + dmass, dmass) 757 mass, dmass = numpy.linspace(delayfunc.lowMass, delayfunc.highMass, npoints, 758 retstep=True) 759 760 #dlogt = 0.01 761 #time = 10**numpy.arange(numpy.log10(minTime), numpy.log10(maxTime) + dlogt, 762 # dlogt) 763 #time = numpy.logspace(numpy.log10(minTime), numpy.log10(maxTime), npoints) 764 time = numpy.linspace(0, maxTime, npoints) 765 766 mass_from_time = delayfunc.m_star(time) 767 time_from_mass = delayfunc.t_life(mass) 768 769 # Calculate delay-function-related quantities. 770 metal_mass_dist = delayfunc.metal_mass_dist(mass) 771 metal_mass_dist_postinterp = delayfunc.metal_mass_dist_postinterp(mass) 772 metal_mass_CCDF = delayfunc.metal_mass_CCDF(mass) 773 metal_time_CDF = delayfunc.metal_time_CDF(time) 774 metal_time_CDF_time_from_mass = delayfunc.metal_time_CDF(time_from_mass) 775 metal_time_dist_binned, dtime_bin = delayfunc.metal_time_dist_binned(time) 776 metal_time_dist_unbinned, dtime_ubin = \ 777 delayfunc.metal_time_dist_unbinned(time) 778 779 imf = delayfunc.imf(mass) 780 yield_mass = delayfunc.yield_mass(mass) 781 782 ### Test some relationships between the quantities. ### 783 import scipy.integrate 784 785 # Do a simple integral over the PDF: 786 test_mass_CCDF = scipy.integrate.cumtrapz((dmass * 787 metal_mass_dist)[::-1])[::-1] 788 diff_mass_CCDF = test_mass_CCDF - metal_mass_CCDF[:-1] 789 790 sum_metal_time_dist = \ 791 numpy.cumsum(numpy.nan_to_num(metal_time_dist_binned)) 792 793 test_time_CDF = numpy.cumsum(metal_time_dist_unbinned * dtime_ubin) 794 795 test_metal_mass_dist = imf * yield_mass / delayfunc._fCnorm 796 797 if not skip_asserts: 798 799 ### The CDF should be the integral of the metal dist. ### 800 801 # The CDF should be the integral of the metal dist: 802 ntest.assert_almost_equal(metal_mass_CCDF[:-1], test_mass_CCDF, 2) 803 804 # The metal time CDF and metal mass CCDF should be 805 # equal. (This could be 10 digits rather than four, except for 806 # a small hicup in the time from mass conversion.)# 807 ntest.assert_almost_equal(metal_time_CDF_time_from_mass, 808 metal_mass_CCDF, 4) 809 810 ### The CCDF and CDF should be close one at lowMass. ### 811 marray = numpy.linspace(delayfunc.lowMass, delayfunc.highMass, 100) 812 ccdf = delayfunc.metal_mass_CCDF(marray) 813 cdf = delayfunc.metal_time_CDF(delayfunc.t_life(marray)) 814 ntest.assert_approx_equal(ccdf[0], 1., 4) 815 ntest.assert_approx_equal(cdf[0], 1., 4) 816 817 # The sum of the binned dist should be one over the norm 818 # range. 819 ntest.assert_approx_equal(sum_metal_time_dist[-1], 1., 4) 820 821 ### The sum of the binned dist should be equal to the CDF ### 822 ntest.assert_almost_equal(\ 823 numpy.cumsum(numpy.nan_to_num(metal_time_dist_binned)), 824 numpy.nan_to_num(metal_time_CDF)) 825 826 ### The integral of the time dist should be equal to the CDF ### 827 ntest.assert_almost_equal(numpy.nan_to_num(metal_time_CDF), 828 test_time_CDF, 10) 829 830 ### The product of IMF and yield should be proportional to the metal 831 ### mass distribution. ### 832 ntest.assert_almost_equal(metal_mass_dist, test_metal_mass_dist, 10) 833 834 ### The metal mass dist should be similar if we interpolate before 835 ### or after the IMF multiplication. ### 836 ntest.assert_almost_equal(metal_mass_dist, 837 metal_mass_dist_postinterp, 2) 838 839 ### Make some plots ### 840 841 if plots: 842 pylab.figure().set_label('metal_mass_CCDF') 843 844 pylab.subplot(221) 845 pylab.plot(mass[:-1], test_mass_CCDF, label="test") 846 pylab.plot(mass, metal_mass_CCDF, label="CCDF") 847 pylab.legend(loc='best') 848 pylab.axvline(x=delayfunc.lowMass) 849 pylab.axvline(x=delayfunc.highMass) 850 pylab.twinx() 851 pylab.plot(mass, metal_mass_dist, ':', label="PDF") 852 853 pylab.subplot(222) 854 pylab.plot(mass[:-1], diff_mass_CCDF, label="test diff") 855 pylab.legend(loc='best') 856 pylab.twinx() 857 pylab.plot(mass, metal_mass_dist, ":", label="PDF") 858 859 pylab.subplot(223) 860 pylab.plot(time_from_mass, metal_mass_CCDF, label="CCDF") 861 pylab.xscale('log') 862 pylab.legend(loc='best') 863 864 pylab.subplot(224) 865 pylab.plot(mass, metal_mass_dist, ':', label="PDF") 866 pylab.axvline(x=delayfunc.lowMass) 867 pylab.axvline(x=delayfunc.highMass) 868 pylab.legend(loc='best') 869 870 pylab.figure().set_label('metal_mass_PDF') 871 pylab.subplot(211) 872 pylab.plot(mass, metal_mass_dist, ':', label="PDF") 873 pylab.plot(mass, metal_mass_dist_postinterp, '-.', label="postinterp") 874 pylab.plot(mass, test_metal_mass_dist, '--', label="~ IMF * yield") 875 pylab.axvline(x=delayfunc.lowMass) 876 pylab.axvline(x=delayfunc.highMass) 877 pylab.legend(loc='best') 878 pylab.subplot(212) 879 pylab.plot(mass, metal_mass_dist - 880 test_metal_mass_dist, ':', label="PDF - c*(IMF * yield)") 881 pylab.plot(mass, metal_mass_dist - 882 metal_mass_dist_postinterp, '-', label="PDF - postinterp") 883 pylab.axvline(x=delayfunc.lowMass) 884 pylab.axvline(x=delayfunc.highMass) 885 pylab.legend(loc='best') 886 887 pylab.figure().set_label('mass_time') 888 pylab.subplot(211) 889 pylab.plot(time, mass_from_time, label="mass from time") 890 pylab.plot(time_from_mass, mass, label="time from mass") 891 pylab.xscale('log') 892 pylab.subplot(212) 893 pylab.plot(mass_from_time, time, label="mass from time") 894 pylab.plot(mass, time_from_mass, label="time from mass") 895 pylab.yscale('log') 896 pylab.legend(loc='best') 897 898 pylab.figure().set_label('metal_time_CDF') 899 900 pylab.subplot(221) 901 pylab.plot(time, metal_time_CDF, label="CDF") 902 pylab.plot(time, sum_metal_time_dist, label="sum of binned") 903 pylab.plot(time, test_time_CDF, label="int of unbinned") 904 pylab.legend(loc='best') 905 pylab.xscale('log') 906 pylab.axvline(x=minTime) 907 pylab.axvline(x=maxTime) 908 909 pylab.subplot(222) 910 pylab.plot(mass_from_time, metal_time_CDF, label="CDF") 911 pylab.plot(mass_from_time, sum_metal_time_dist, label="sum of binned") 912 pylab.plot(mass_from_time, test_time_CDF, label="int of unbinned") 913 pylab.axvline(x=delayfunc.lowMass) 914 pylab.axvline(x=delayfunc.highMass) 915 pylab.legend(loc='best') 916 917 pylab.subplot(223) 918 pylab.plot(time, sum_metal_time_dist - metal_time_CDF, label="sum-CDF") 919 pylab.plot(time, test_time_CDF - metal_time_CDF, label="int-CDF") 920 pylab.legend(loc='best') 921 pylab.xscale('log') 922 pylab.axvline(x=minTime) 923 pylab.axvline(x=maxTime) 924 925 pylab.subplot(224) 926 pylab.plot(mass_from_time, sum_metal_time_dist - metal_time_CDF, 927 label="sum-CDF") 928 pylab.plot(mass_from_time, 929 test_time_CDF - metal_time_CDF, label="int-CDF") 930 pylab.legend(loc='best') 931 pylab.axvline(x=delayfunc.lowMass) 932 pylab.axvline(x=delayfunc.highMass) 933 934 935 pylab.figure().set_label('metal_time_PDF_bin') 936 pylab.plot(time, metal_time_dist_binned, label="binned") 937 pylab.xscale('log') 938 pylab.axvline(x=minTime) 939 pylab.axvline(x=maxTime) 940 pylab.legend(loc='best') 941 942 pylab.figure().set_label('metal_time_PDF_unbin') 943 pylab.plot(time, metal_time_dist_unbinned, label="unbinned") 944 pylab.xscale('log') 945 pylab.legend(loc='best') 946 947 pylab.figure().set_label('CDF_CCDF') 948 pylab.subplot(211) 949 pylab.plot(mass, metal_mass_CCDF, label='mass CCDF') 950 pylab.plot(mass, metal_time_CDF_time_from_mass, label='time CDF') 951 pylab.legend(loc='best') 952 pylab.subplot(212) 953 pylab.plot(mass, metal_time_CDF_time_from_mass - metal_mass_CCDF, ':', 954 label='CDF - CCDF') 955 pylab.legend(loc='best')

956 957

958 -def compare_dist_times():

959 import time 960 timelist = [] 961 names = [] 962 963 names.append('init') 964 timelist.append(time.time()) 965 delayfunc = Delay(0.02) 966 967 names.append('linspace') 968 timelist.append(time.time()) 969 m = numpy.linspace(delayfunc.lowMass, delayfunc.highMass, 1e3) 970 971 names.append('postinterp') 972 timelist.append(time.time()) 973 delayfunc.metal_mass_dist_postinterp(m) 974 975 names.append('dist') 976 timelist.append(time.time()) 977 delayfunc.metal_mass_dist(m) 978 979 names.append('ccumulate dist') 980 timelist.append(time.time()) 981 utils.ccumulate(delayfunc.metal_mass_dist, m, max=delayfunc.highMass) 982 983 names.append('ccumulate dist romberg') 984 timelist.append(time.time()) 985 utils.ccumulate(delayfunc.metal_mass_dist, m, max=delayfunc.highMass, 986 method='romberg') 987 988 names.append('ccumulate postintepr') 989 timelist.append(time.time()) 990 utils.ccumulate(delayfunc.metal_mass_dist_postinterp, m, 991 max=delayfunc.highMass) 992 993 names.append('ccumulate postintepr romberg') 994 timelist.append(time.time()) 995 utils.ccumulate(delayfunc.metal_mass_dist_postinterp, m, 996 max=delayfunc.highMass, 997 method='romberg') 998 999 timelist.append(time.time()) 1000 timearray = numpy.asarray(timelist) 1001 1002 time_diffs = 1e3 * (timearray[1:] - timearray[:-1]) 1003 1004 for diff, name in zip(time_diffs, names): 1005 print "%s: %.4g ms" % (name, diff)

1006

1007 -def run_profs():

1008 import cProfile 1009 cProfile.run('delayfunc = Delay(0.02)', 'Delay.prof') 1010 cProfile.run('m = numpy.linspace(delayfunc.lowMass, delayfunc.highMass, 1e3)') 1011 1012 cProfile.run("""utils.ccumulate(delayfunc.metal_mass_dist, m, max=delayfunc.highMass)""", 'cumulate_dist.prof') 1013 1014 cProfile.run("""utils.ccumulate(delayfunc.metal_mass_dist_postinterp, m, max=delayfunc.highMass)""", 'cumulate_dist_postinterp.prof')

1015 1016 1017 if __name__ == '__main__': 1018 import matplotlib.pyplot as pylab 1019 import os 1020 import sys 1021 1022 ### Argument parsing. ### 1023 if len(sys.argv)==1: 1024 print "Run with a filename argument to produce image files, e.g.:" 1025 print " python delay.py delay.png" 1026 print " python delay.py delay.eps" 1027 if len(sys.argv) > 1: 1028 filename = sys.argv[1] 1029 prefix, extension = os.path.splitext(filename) 1030 else: 1031 filename = None 1032 1033 ### Main code area. ### 1034 1035 #run_profs() 1036 #compare_dist_times() 1037 1038 test_Delay(skip_asserts=False, 1039 delayfunc = Delay(0.02, maxTime=2e9, ejections=False)) 1040 test_Delay(skip_asserts=False, 1041 delayfunc = Delay(0.02, maxTime=2e9, ejections=True)) 1042 1043 ##test_Delay(skip_asserts=False, delayfunc = Delay(0.02)) 1044 1045 #test_Delay(skip_asserts=True, 1046 # delayfunc = Delay(0.0001, maxTime=2e9, 1047 # yieldClass = yields.Data_H)) 1048 1049 #bins = numpy.array([1.4, 5, 100.]) 1050 #fracs = numpy.array([0.5, 0.5]) 1051 #test_Delay(skip_asserts=True, delayfunc = MockDelay(bins, fracs)) 1052 1053 1054 #test_plot_converge_delay_function() 1055 #plot_delay_function(0.02) 1056 #test_plot_binned_delay_function() 1057 1058 ### Plot output code. ### 1059 if filename is None: 1060 pylab.show() 1061 else: 1062 from matplotlib import _pylab_helpers 1063 for manager in _pylab_helpers.Gcf.get_all_fig_managers(): 1064 fig = manager.canvas.figure 1065 if len(fig.get_label()) > 0: 1066 label = fig.get_label() 1067 else: 1068 label = '_Fig' + str(manager.num) 1069 newfilename = prefix + '_' + label + extension 1070 fig.savefig(newfilename, bbox_inches="tight") 1071

Source Code for Module enrichpy.delay