enrichpy.yields

1 """Encapsulates stellar chemical yield data, currently only carbon 12. 2 3 General scheme: 4 5 For each data source, a Yield_Data subclass handles reading in the 6 data and converting it into numpy array of yields, stellar masses, and 7 stellar metallicites. 8 9 The Interpolate_Yields class take a Yield_Data subclass instance and 10 handles the creation of an interpolation function that returns the 11 yield as a function of stellar mass and metallicity. 12 13 """ 14 import re 15 import optparse 16 import os 17 from pkg_resources import resource_stream 18 19 import scipy 20 import scipy.interpolate as si 21 import scipy.interpolate.fitpack2 as sif2 22 import numpy 23 24 import initialmassfunction 25

26 -def ejection_from_yield(mYield, mIni, mRem, xIni):

27 """Convert net yield of an element to total ejected mass. 28 29 mYield, mIni, and mRem should be in the same units (which will also 30 be the units of the returned ejected mass. 31 32 xIni is the initial mass fraction of the element. 33 34 See GBM equation 11. 35 36 Gavilan M., Buell J.F., Molla M. Astron. Astrophys. 2005, 432, 861 37 (2005A&A...432..861G) 38 """ 39 return mYield + ((mIni - mRem) * xIni)

40

41 -def yield_from_ejection(mEj, mIni, mRem, xIni):

42 """Convert total ejected mass to net yield of an element. 43 44 mEj, mIni, and mRem should be in the same units (which will also 45 be the units of the returned yield. 46 47 xIni is the initial mass fraction of the element. 48 49 See GBM equation 11. 50 51 Gavilan M., Buell J.F., Molla M. Astron. Astrophys. 2005, 432, 861 52 (2005A&A...432..861G) 53 """ 54 return mEj - ((mIni - mRem) * xIni)

55

56 -class Yield_Data:

57 """Base class for yield data.""" 58 59 ### Subclasses should define the following, where applicable: ### 60 #shortname = 61 #filename = 62 #metal_col = # stellar inition metallicity 63 #mass_col = # stellar initial mass 64 #rem_col = # stellar remnant mass 65 #yield_col = # net mass yield of element 66 #ejection_col = # total mass of element ejected 67 ### 68

69 - def __init__(self, interp_args={}):

70 self.interp_args=interp_args 71 self.loadtxt_args={} 72 self.metal_mass_frac_C = 0.178 73 self.load_data()

74

75 - def load_data(self):

76 self.read_table() 77 self.select_data() 78 self.regrid_data()

79

80 - def read_table(self):

81 """Read data in self.filename. 82 83 Uses the pkg_resources.resource_stream function to access the 84 data file. 85 """ 86 stream = resource_stream(__name__, self.filename) 87 self.data = numpy.loadtxt(stream, **self.loadtxt_args) 88 stream.close()

89

90 - def select_data(self):

91 if hasattr(self, 'mass_col'): 92 self.mass_star = self.data[:,self.mass_col] 93 if hasattr(self, 'rem_col'): 94 self.mass_rem = self.data[:,self.rem_col] 95 if hasattr(self, 'metal_col'): 96 self.metal_star = self.data[:,self.metal_col] 97 if hasattr(self, 'yield_col'): 98 self.mass_C = self.data[:,self.yield_col] 99 if hasattr(self, 'ejection_col'): 100 self.mass_C_total = self.data[:,self.ejection_col] 101 102 # If we have net yield and remnant mass, but no total ejection mass: 103 if (hasattr(self, 'yield_col') and hasattr(self, 'rem_col') and 104 not hasattr(self, 'ejection_col')): 105 self.mass_C_total = ejection_from_yield(self.mass_C, 106 self.mass_star, 107 self.mass_rem, 108 self.metal_mass_frac_C * 109 self.metal_star) 110 # If we total ejection and remnant mass, but no net yield mass: 111 if (hasattr(self, 'ejection_col') and hasattr(self, 'rem_col') and 112 not hasattr(self, 'yield_col')): 113 self.mass_C = yield_from_ejection(self.mass_C_total, 114 self.mass_star, 115 self.mass_rem, 116 self.metal_mass_frac_C * 117 self.metal_star)

118 119

120 - def regrid_data(self):

121 """Subclasses can define this if they want.""" 122 pass

123

124 - def __call__(self, ejections=False):

125 if ejections: 126 return self.mass_star, self.metal_star, self.mass_C_total 127 else: 128 return self.mass_star, self.metal_star, self.mass_C

129

130 - def interpolate_function(self, ejections=False):

131 return Interpolate_Yields(*self(ejections=ejections), 132 **self.interp_args)

133 #return Interpolate_Yields(*self(ejections=ejections)) 134

135 -class Data_vdHG(Yield_Data):

136 """van den Hoek and Groenewegen (1997): 0.8 -- 8 Msun; Z = 0.001 -- 0.04 137 138 van den Hoek, L.B., \& Groenewegen, M.A.T.\ 1997, \aaps, 123, 305 139 (1997A&AS..123..305V) 140 141 [http://vizier.cfa.harvard.edu/viz-bin/VizieR?-source=J/A+AS/123/305] 142 143 Notes 144 ----- 145 146 metal_mass_frac_C is infered to be 0.224 from the table below: 147 148 Table1: Initial element abundances adopted 149 -------------------------------------------------------------- 150 Element. Z=0.001 0.004 0.008 0.02 0.04 151 -------------------------------------------------------------- 152 H 0.756 0.744 0.728 0.68 0.62 153 4He 0.243 0.252 0.264 0.30 0.34 154 12C 2.24E-4 9.73E-4 1.79E-3 4.47E-3 9.73E-3 155 13C 0.04E-4 0.16E-4 0.29E-4 0.72E-4 1.56E-4 156 14N 0.70E-4 2.47E-4 5.59E-4 1.40E-3 2.47E-3 157 16O 5.31E-4 2.11E-3 4.24E-3 1.06E-2 2.11E-2 158 159 """ 160 161 shortname = 'vdH&G1997' 162 filename = 'yield_data/1997A&AS..123..305V_tab2-22_total.dat' 163 metal_col=0 164 mass_col=1 165 yield_col=4 166 mass_Ej_col = 10 167

168 - def __init__(self, interp_args={}):

169 self.interp_args=interp_args 170 self.loadtxt_args = {'usecols' : range(1,13)} 171 self.metal_mass_frac_C = 0.224 172 self.load_data()

173

174 - def select_data(self):

175 self.mass_star = self.data[:,self.mass_col] 176 self.mass_star_ejected = self.data[:,self.mass_Ej_col] 177 self.mass_rem = self.mass_star - self.mass_star_ejected 178 self.metal_star = self.data[:,self.metal_col] 179 self.mass_C = self.data[:,self.yield_col] 180 181 # If we have net yield and remnant mass, but no total ejection mass: 182 self.mass_C_total = ejection_from_yield(self.mass_C, 183 self.mass_star, 184 self.mass_rem, 185 self.metal_mass_frac_C * 186 self.metal_star)

187 188

189 -class Data_H(Yield_Data):

190 """Herwig 2004: 2-6 Msun; Z=0.0001 191 192 Herwig, F.\ 2004, \apjs, 155, 651 (2004ApJS..155..651H) 193 194 2-6 Msun; Z=0.0001 195 [http://vizier.cfa.harvard.edu/viz-bin/VizieR?-source=J/ApJS/155/651] 196 197 Intermediate-mass stellar evolution tracks from the main 198 sequence to the tip of the AGB for five initial masses (2-6 199 Msolar) and metallicity Z=0.0001 have been computed. ...yields 200 are presented. 201 202 """ 203 204 shortname = 'H2004' 205 filename = 'yield_data/2004ApJS..155..651H_table5.dat'

206 - def _H_iso_converter(*args):

207 """Strip element names from isotope in ChLi table.""" 208 isotope = args[-1] 209 alphanum = re.compile('(\D*)(\d*)') # non-digit followed by digit 210 element, number = alphanum.match(isotope).groups() 211 if number == '': 212 number=1 213 return int(number)

214

215 - def __init__(self, interp_args={}):

216 self.interp_args=interp_args 217 self.loadtxt_args={'converters' : {1: self._H_iso_converter}} 218 self.metal_mass_frac_C = 0.178 219 self.load_data()

220

221 - def select_data(self):

222 #self.mass_star = self.data[:,self.mass_col] 223 #self.metal_star = self.data[:,self.metal_col] 224 self.mass_C = self.data[2,2:] 225 self.mass_star = numpy.array([2.,3.,4.,5.,6.]) 226 self.metal_star = numpy.ones(self.mass_star.shape) * 0.0001

227

228 -class Data_GBM(Yield_Data):

229 """Gavilan, Buell, & Molla (2005): 0.8 -- 100 Msun; log(Z/Zsun)=-0.2--+0.2 230 231 Gavilan M., Buell J.F., Molla M. Astron. Astrophys. 2005, 432, 861 232 (2005A&A...432..861G) 233 234 log(Z/Zsun)=-0.2, -0.1, 0.0, +0.1 and +0.2 235 236 [http://cdsarc.u-strasbg.fr/viz-bin/ftp-index?J/A%2bA/432/861] 237 238 """ 239 240 shortname = 'GB&M2005' 241 filename = 'yield_data/2005A&A...432..861G/table2.dat' 242 metal_col=0 243 mass_col=1 244 rem_col=2 245 ejection_col=5 246

247 - def regrid_data(self):

248 """ Interpolate data onto a common mass grid. 249 """ 250 mass_C = self.mass_C 251 mass_C_total = self.mass_C_total 252 mass_star = self.mass_star 253 metal_star = self.metal_star 254 255 ### Set up a mass scale with reasonable spacing. ### 256 257 # Get all unique mass values in the data set. 258 mass_scale = numpy.sort(numpy.unique(mass_star)) 259 260 # Now step through the values, discarding ones that are too close 261 # together. 262 263 last_mass = mass_scale[0] 264 new_mass_scale = [last_mass] 265 numavg = 0 266 for m in list(mass_scale[1:]): 267 if m - last_mass < (0.04 * m): 268 # If the next value is too close, we average the positions. 269 numavg = numavg + 1 270 new_mass_scale[-1] = (numavg * new_mass_scale[-1] + m)/(numavg+1) 271 else: 272 numavg = 0 273 new_mass_scale.append(m) 274 last_mass = m 275 mass_scale = numpy.array(new_mass_scale) 276 277 ZVals = numpy.unique(metal_star) 278 279 new_len = len(ZVals) * len(mass_scale) 280 new_mass_C = numpy.zeros(new_len) 281 new_mass_C_total = numpy.zeros(new_len) 282 repeated_inds = range(len(mass_scale)) * len(ZVals) 283 new_mass_star = mass_scale[repeated_inds] 284 new_metal_star = numpy.zeros(new_len) 285 286 doplot = False 287 if doplot: 288 print new_mass_star 289 import pylab 290 291 for Z, Zind in zip(list(ZVals), range(len(ZVals))): 292 mask = metal_star == Z 293 ifunc = si.interp1d(numpy.log10(mass_star[mask]), 294 mass_C[mask], 295 kind=3) 296 ifunc_total = si.interp1d(numpy.log10(mass_star[mask]), 297 mass_C_total[mask], 298 kind=3) 299 i0 = Zind * len(mass_scale) 300 i1 = (Zind+1) * len(mass_scale) 301 302 new_mass_C[i0:i1] = ifunc(numpy.log10(mass_scale)) 303 new_mass_C_total[i0:i1] = ifunc_total(numpy.log10(mass_scale)) 304 new_metal_star[i0:i1] = Z 305 306 if doplot: 307 print mass_scale.shape 308 print new_mass_C[i0:i1].shape 309 pylab.figure() 310 pylab.scatter(mass_star[mask], mass_C[mask], marker='x') 311 pylab.plot(new_mass_star[i0:i1], new_mass_C[i0:i1], '-+') 312 pylab.title(str(Z)) 313 314 self.orig_mass_star = self.mass_star 315 self.orig_metal_star = self.metal_star 316 self.orig_mass_C = self.mass_C 317 self.orig_mass_C_total = self.mass_C_total 318 319 self.mass_star = new_mass_star 320 self.metal_star = new_metal_star 321 self.mass_C = new_mass_C 322 self.mass_C_total = new_mass_C_total

323

324 -class Exponentiate:

325 """Decorator to exponentiate some arguments before calling the function. 326 """

327 - def __init__(self, xexp=True, yexp=True):

328 self.xexp = xexp 329 self.yexp = yexp

330

331 - def __call__(self, function):

332 if (not self.xexp) and (not self.yexp): 333 return function 334 335 if self.xexp and (not self.yexp): 336 newfunction = lambda x, y: function(10**x, y) 337 elif (not self.xexp) and self.yexp: 338 newfunction = lambda x, y: function(x, 10**y) 339 elif self.xexp and self.yexp: 340 newfunction = lambda x, y: function(10**x, 10**y) 341 342 newfunction.__name__ = function.__name__ 343 newfunction.__dict__.update(function.__dict__) 344 newfunction.__doc__ = function.__doc__ 345 return newfunction

346

347 -class Interpolate_Yields:

348 """The ejected carbon mass is interpolated on a logarithmic 349 metalicity scale and mass scale. 350 351 Essentially a convenience class to encapsulate my current choices 352 about the interpolation. 353 354 """ 355

356 - def __init__(self, mass_star, metal_star, mass_C, weight_function=None, 357 mass_scale='linear', metal_scale='linear', 358 kx=1, 359 ky=1, 360 swapxy = False):

361 self.swapxy = swapxy 362 self.mass_scale = mass_scale 363 self.metal_scale = metal_scale 364 365 # Find unique m and Z values. 366 mVals = numpy.unique(mass_star) 367 mVals.sort() 368 369 ZVals = numpy.unique(metal_star) 370 ZVals.sort() 371 372 self.x = mass_star 373 if not mass_scale=='linear': 374 self.x = numpy.log10(x) 375 376 self.y = metal_star 377 if not metal_scale=='linear': 378 self.y = numpy.log10() 379 380 self.xknots = numpy.unique(self.x) 381 self.xknots.sort() 382 self.yknots = numpy.unique(self.y) 383 self.yknots.sort() 384 385 self.mVals = mVals 386 self.ZVals = ZVals 387 388 self.mass_C = mass_C 389 self.z = self.mass_C 390 if weight_function is not None: 391 self.z *= weight_function(self.x) 392 393 ### bicubic spline interpolation ### 394 395 if self.swapxy: 396 ## Note that x and y are swapped due to an error I encountered 397 ## trying to use kx > ky. 398 self._interpfunc = sif2.LSQBivariateSpline(self.y, self.x, 399 self.z, 400 self.yknots, 401 self.xknots, 402 kx=ky, ky=kx) 403 self._swap_mCfunc = Exponentiate(not metal_scale=='linear', 404 not mass_scale=='linear')\ 405 (self._interpfunc) 406 self._mCfunc = lambda x, y: self._swap_mCfunc(y,x).transpose() 407 else: 408 self._interpfunc = sif2.LSQBivariateSpline(self.x, self.y, 409 self.mass_C, 410 self.xknots, 411 self.yknots, 412 kx=kx, ky=ky) 413 self._mCfunc = Exponentiate(not mass_scale=='linear', 414 not metal_scale=='linear')\ 415 (self._interpfunc)

416

417 - def __call__(self, m, Z):

418 ## Note that x and y are swapped due to an error I encountered 419 ## trying to use kx > ky. 420 return self._mCfunc(m, Z)

421

422 -def linear_coefficients(x,y):

423 slope = (y[1:]-y[:-1]) / (x[1:] - x[:-1]) 424 intercept = y[:-1] - slope * x[:-1] 425 return slope, intercept

426

427 -def print_linear_info(slope, intercept):

428 print " mass_C = %.3g m_star + %3.g Msun" % (slope, intercept) 429 print " mass_C = 0 when m = %.3g Msun" % \ 430 (-1. * intercept / slope)

431

432 -def extrapolate_function(yield_function, metallicity, return_coeffs=False):

433 """Linearly extrapolate yield data to high and low masses. 434 435 Returns two functions: the high mass extrapolation and the low 436 mass extrapolation. 437 """ 438 439 # High mass extrapolation 440 highM = numpy.sort(yield_function.mVals)[-2:] 441 high_mass_C = yield_function(highM, metallicity).flatten() 442 highSlope, highIntercept = linear_coefficients(highM, high_mass_C) 443 print "Extrapolating from points at m_star = %.3g and %.3g Msun:" % \ 444 (highM[0], highM[1]) 445 print_linear_info(highSlope, highIntercept) 446 exfuncHigh = lambda mass: mass * highSlope + highIntercept 447 448 # Low mass extrapolation 449 lowM = numpy.sort(yield_function.mVals)[:2] 450 low_mass_C = yield_function(lowM, metallicity).flatten() 451 lowSlope, lowIntercept = linear_coefficients(lowM, low_mass_C) 452 print "Extrapolating from points at m_star = %.3g and %.3g Msun:" % \ 453 (lowM[0], lowM[1]) 454 print_linear_info(lowSlope, lowIntercept) 455 exfuncLow = lambda mass: mass * lowSlope + lowIntercept 456 457 if return_coeffs: 458 return (exfuncHigh, exfuncLow, 459 highSlope, highIntercept, lowSlope, lowIntercept) 460 else: 461 return exfuncHigh, exfuncLow

462

463 -class Data_ChLi(Yield_Data):

464 """Chieffi & Limongi (2004) yields: 13--35 Msun; Z = 0, 1e-6 -- 0.02 465 466 Chieffi, A., \& Limongi, M.\ 2004, \apj, 608, 405 467 (2004ApJ...608..405C) 468 469 [http://vizier.cfa.harvard.edu/viz-bin/Cat?J/ApJ/608/405] 470 471 """ 472 shortname = 'Ch&Li2004' 473 filename = 'yield_data/2004ApJ...608..405C_yields.dat'

474 - def _ChLi_iso_converter(*args):

475 """Strip element names from isotope in ChLi table.""" 476 isotope = args[-1] 477 alphanum = re.compile('(\d*)(\D*)') # digit followed by non-digit 478 number, element = alphanum.match(isotope).groups() 479 return int(number)

480

481 - def __init__(self, interp_args={}):

482 self.interp_args=interp_args 483 self.loadtxt_args={'converters' : {2: self._ChLi_iso_converter}} 484 self.metal_mass_frac_C = 0.178 485 self.load_data()

486

487 - def select_data(self):

488 """Return mass, metallicity, and C yield. 489 """ 490 data = self.data 491 isonumber = data[:,2] 492 # Select 12C. 493 cmask = isonumber == 12 494 495 # Set of mass values for the columns 4-9. 496 mass_star = numpy.array([13., 15., 20., 25., 30., 35.]) 497 498 # Metallicity values. 499 metal_star = data[cmask][:,0:1] 500 501 # Yield data 502 mass_C = data[cmask][:,3:] 503 504 # Form compatible grid of mass_star and metal_star values. 505 ones = numpy.ones(mass_C.shape) 506 mass_star = ones * mass_star 507 metal_star = ones * metal_star 508 self.mass_star = mass_star.flatten() 509 self.metal_star = metal_star.flatten() 510 self.mass_C = mass_C.flatten() 511 return self.mass_star, self.metal_star, self.mass_C

512

513 -def test_plot_data_interp(data, nm=1e3, nZ=1e3, **args):

514 515 mCfunc = data.interpolate_function() 516 517 mass_star, metal_star, mass_C = data() 518 519 m = numpy.linspace(0.75 * min(mass_star), 1.05 * max(mass_star), nm) 520 from collections import deque 521 marks = deque(['.', 'o', 'v','+','x','*']) 522 colors = deque(['r', 'y', 'g','c','b','k', 'm']) 523 524 import matplotlib.pyplot as pylab 525 pylab.figure() 526 for Z in numpy.unique(metal_star): 527 mark = marks[0] 528 col = colors[0] 529 marks.rotate() 530 colors.rotate() 531 datamask = metal_star==Z 532 #print mass_star[datamask], mass_C[datamask] 533 pylab.plot(mass_star[datamask], mass_C[datamask], 534 marker=mark, label=str(Z), c=col, ls='') 535 pylab.plot(m, mCfunc(m, Z), c=col, ls=':') 536 537 pylab.legend(loc='best') 538 539 Z = numpy.logspace(numpy.min(numpy.log10(metal_star[metal_star>0])) - 1, 540 numpy.max(numpy.log10(metal_star)) + 0.5, nZ) 541 if numpy.any(metal_star <= 0): 542 Z = numpy.hstack(([0], Z)) 543 544 logZ = numpy.log10(Z) 545 pylab.figure() 546 for M in numpy.unique(mass_star): 547 mark = marks[0] 548 col = colors[0] 549 marks.rotate() 550 colors.rotate() 551 datamask = numpy.logical_and(mass_star==M, metal_star>0) 552 pylab.plot(numpy.log10(metal_star[datamask]), mass_C[datamask], 553 marker=mark, label=str(M), c=col, ls='') 554 pylab.plot(logZ[Z>0], mCfunc(M, Z).flatten()[Z>0], c=col, ls=':') 555 556 pylab.legend(loc='best') 557 558 559 pylab.figure() 560 for M in numpy.unique(mass_star): 561 mark = marks[0] 562 col = colors[0] 563 marks.rotate() 564 colors.rotate() 565 datamask = mass_star==M 566 pylab.plot(metal_star[datamask], mass_C[datamask], 567 marker=mark, label=str(M), c=col, ls='') 568 pylab.plot(Z, mCfunc(M, Z).flatten(), c=col, ls=':') 569 #print mCfunc(M, 10**logZ) 570 571 #pylab.xlim(Z[0], Z[-1]) 572 pylab.legend(loc='best')

573 574

575 -def data_overlap(data_list):

576 names = ['m*', 'Z*', 'mC'] 577 mins = [numpy.inf, numpy.inf, numpy.inf] 578 maxa = [None, None, None] 579 maxmin = [None, None, None] 580 minmax = [None, None, None] 581 for data in data_list: 582 print data.shortname 583 datarrays = data() 584 for i, dat in enumerate(datarrays): 585 imin = numpy.nanmin(dat) 586 imax = numpy.nanmax(dat) 587 print(("%5.3g <= " + names[i] + " <= %5.3g ") % (imin, imax)), 588 mins[i] = min(imin, mins[i]) 589 maxa[i] = max(imax, maxa[i]) 590 maxmin[i] = max(imin, maxmin[i]) 591 minmax[i] = min(imax, maxmin[i]) 592 print 593 print mins 594 print maxa 595 print maxmin 596 print minmax 597 return mins, maxa, maxmin, minmax

598

599 -def test_plot_interp(data_list, Z_list, nm=1e3, 600 maxmass=None, ejections=False, title=None, 601 logx=True, logy=True):

602 import matplotlib.pyplot as pylab 603 fig = pylab.figure(figsize=(9,6)) 604 if title is not None: 605 fig.set_label(title) 606 ax = fig.add_subplot(111) 607 608 from collections import deque 609 styles = deque(['-', '--', ':', '-.']) 610 611 mins, maxa, maxmin, minmax = data_overlap(data_list) 612 613 for data in data_list: 614 print data.shortname 615 style = styles[0] 616 styles.rotate(-1) 617 618 if ejections: 619 if not hasattr(data, 'mass_C_total'): 620 print "No total ejected masses for " + data.shortname 621 continue 622 mCfunc_ej = data.interpolate_function(ejections=True) 623 mass_star, metal_star, mass_C_ej = data(ejections=True) 624 625 mCfunc = data.interpolate_function() 626 mass_star, metal_star, mass_C = data() 627 628 #m = numpy.linspace(0.75 * min(mass_star), 1.05 * max(mass_star), nm) 629 m = numpy.linspace(min(mass_star), max(mass_star), nm) 630 colors = deque(['k','b', 'c', 'g','r']) 631 for Z in Z_list: 632 color = colors[0] 633 colors.rotate(-1) 634 if (Z < min(metal_star)) or (Z > max(metal_star)): 635 print "skipping ", 636 print Z, 637 print data.shortname 638 639 continue 640 Zlab = "Z = %.3g" % (Z) 641 pylab.plot(m, mCfunc(m, Z), ls=style, c=color, 642 label=data.shortname + " " + Zlab) 643 if ejections: 644 color = colors[0] 645 colors.rotate(-1) 646 pylab.plot(m, mCfunc_ej(m, Z), ls=style, c=color, 647 label=data.shortname + " ej " + Zlab) 648 if maxmass is not None: 649 pylab.xlim(xmax=maxmass) 650 #pylab.legend(loc='best') 651 labels=[] 652 lines=[] 653 for line in ax.lines: 654 label = line.get_label() 655 if label[0] == '_': 656 continue 657 labels.append(label) 658 lines.append(line) 659 if logx: 660 pylab.xscale('log') 661 if logy: 662 pylab.yscale('log') 663 fig.legend(lines, labels, 'lower right') 664 fig.subplots_adjust(left=0.07, bottom=0.07, top=0.95, right=0.62)

665 666

667 -class Data_CaLa(Yield_Data):

668 """Campbell & Lattanzio 2008 0.85 -- 3 Msun; Z = 0, & [Fe/H] = -6.5 -- -3.0 669 670 Campbell, S.~W., \& Lattanzio, J.~C.\ 2008, \aap, 490, 769 671 (2008A&A...490..769C) 672 673 [http://vizier.u-strasbg.fr/viz-bin/VizieR?-source=J/A+A/490/769] 674 675 """ 676 shortname = 'Ca&La2008' 677 basefilename = 'yield_data/2008A&A...490..769C_table' 678

679 - def __init__(self, interp_args={}):

680 self.interp_args={'swapxy':True} 681 self.interp_args.update(interp_args) 682 self.loadtxt_args={} 683 self.metal_mass_frac_C = 0.178 684 self.load_data()

685 686 # def _CaLa_iso_converter(self, isotope): 687 # return 0.0 688 # def _CaLa_iso_converter(self, isotope): 689 # """Strip element names from isotope in CaLa table.""" 690 # alphanum = re.compile('(\D*)(\d*)') # non-digit followed by digit 691 # element, number = alphanum.match(isotope).groups() 692 # print isotope 693 # print element 694 # print number 695 # return int(number) 696

697 - def read_table(self):

698 self.tables = [] 699 for i in range(1,6): 700 cols = [1,2,3,4,5,6] 701 if i==4: 702 cols = [1,2,3,4,5] 703 stream = resource_stream(__name__, 704 self.basefilename + str(i) + '.dat') 705 table = numpy.loadtxt(stream, 706 #converters={0: self._CaLa_iso_converter}, 707 usecols=cols 708 ) 709 stream.close() 710 if i==4: 711 table = numpy.hstack((table, 712 numpy.nan * numpy.ones((len(table),1)))) 713 self.tables.append(table) 714 #print "Table " + str(i) 715 #print table.shape 716 #import newio 717 #table = newio.loadtxt(self.basefilename + '6.dat') 718 stream = resource_stream(__name__, self.basefilename + '6.dat') 719 table = numpy.loadtxt(stream) 720 stream.close() 721 self.tables.append(table) 722 self.data = numpy.vstack(self.tables[0:5]) 723 self.remnant_data = self.tables[5] 724 return self.tables

725

726 - def select_data(self):

727 """Return mass, metallicity, and C yield. 728 """ 729 data = self.data 730 remnant_data = self.remnant_data 731 732 # Set of mass values for the columns 2-5. 733 mass_star = numpy.array([0.85, 1.0, 2.0, 3.0]) 734 735 # Metallicity values. 736 Z_solar = 0.02 737 738 metal_star = Z_solar * numpy.array([[0, 739 10**-6.5, 740 10**-5.45, 741 10**-4.0, 742 10**-3.0]]).transpose() 743 744 # Select 12C. 745 isonumber = data[:,0] 746 cmask = isonumber == 12 747 748 # Fractional mass values. Columns are: 749 # A; Initial; 0.85 Msun; 1.0 Msun; 2.0 Msun; 3.0 Msun 750 init_frac_mass_C = data[cmask][:,1:2] 751 frac_mass_C = data[cmask][:,2:] 752 753 # Table 6 contains Remnant Masses: 754 # Met 0.85 1.0 2.0 3.0 755 mass_Re = remnant_data[:, 1:] 756 # Ejected mass = init mass - remnant mass 757 mass_Ej = mass_star - mass_Re 758 759 # initial mass of carbon 760 init_mass_C = init_frac_mass_C * mass_star 761 762 # total ejected mass of carbon 763 mass_C_tot = frac_mass_C * mass_Ej 764 765 # net mass of carbon produced 766 mass_C_net = mass_C_tot - init_mass_C 767 768 # Form compatible grid of mass_star and metal_star values. 769 ones = numpy.ones(mass_C_tot.shape) 770 mass_star = ones * mass_star 771 metal_star = ones * metal_star 772 773 # Convert to 1D arrays and mask out missing data: 774 data_mask = numpy.isfinite(mass_Ej.flatten()) 775 self.mass_star = mass_star.flatten()[data_mask] 776 self.metal_star = metal_star.flatten()[data_mask] 777 self.mass_Ej = mass_Ej.flatten()[data_mask] 778 self.mass_C_total = mass_C_tot.flatten()[data_mask] 779 self.mass_C = mass_C_net.flatten()[data_mask] 780 781 return self.mass_star, self.metal_star, self.mass_C

782 783

784 -def test_all():

785 imf = initialmassfunction.imf_number_Chabrier 786 args = {'weight_function' : imf} 787 weighted_list = [Data_GBM(args), 788 Data_ChLi(args), 789 Data_CaLa(args), 790 Data_vdHG(args), 791 Data_H(args)] 792 unweighted_list = [Data_GBM(), 793 Data_ChLi(), 794 Data_CaLa(), 795 Data_vdHG(), 796 Data_H()] 797 metal_list = [0.02, 1e-3, 1e-4, 1e-5, 0.0] 798 for ej in [True, False]: 799 if ej: 800 pt = 'Ej' 801 else: 802 pt = 'Yi' 803 test_plot_interp(unweighted_list, metal_list, 804 ejections=ej, title='unwFull'+pt) 805 # test_plot_interp(unweighted_list, metal_list, 806 # maxmass=8.0, 807 # ejections=ej, title='unwZoom'+pt) 808 809 test_plot_interp(weighted_list, metal_list, 810 ejections=ej, title='wFull'+pt)

811 # test_plot_interp(weighted_list, metal_list, 812 # maxmass=8.0, 813 # ejections=ej, title='wZoom'+pt) 814 815 816 # test_plot_data_interp(Data_GBM()) 817 # test_plot_data_interp(Data_ChLi()) 818 # test_plot_data_interp(Data_CaLa()) 819 820 if __name__ == '__main__': 821 822 usage = """ """ 823 824 ### Argument parsing. ### 825 parser = optparse.OptionParser(usage) 826 parser.add_option("-f", "--file", action='store_true', dest="filename", 827 default=None) 828 (options, args) = parser.parse_args() 829 if (len(args) > 0) and (options.filename is None): 830 options.filename = args[0] 831 if options.filename is None: 832 print "No filename given." 833 print usage 834 else: 835 prefix, extension = os.path.splitext(options.filename) 836 837 ### Main code ### 838 839 import pylab 840 test_all() 841 842 ### Plot output code. ### 843 if options.filename is None: 844 pylab.show() 845 else: 846 from matplotlib import _pylab_helpers 847 for manager in _pylab_helpers.Gcf.get_all_fig_managers(): 848 fig = manager.canvas.figure 849 if len(fig.get_label()) > 0: 850 label = fig.get_label() 851 else: 852 label = '_Fig' + str(manager.num) 853 newfilename = prefix + '_' + label + extension 854 fig.savefig(newfilename, dpi=75)#, bbox_inches="tight") 855

Source Code for Module enrichpy.yields