In [19]:
%matplotlib inline
from galaxy_analysis.plot.plot_styles import *
import matplotlib.pyplot as plt
import numpy as np
import onezone.data_tables as DT
from galaxy_analysis.utilities import convert_abundances as convert
import onezone.imf as imf
import h5py
def plot_settings():
fsize = 21
rc('text',usetex=False)
rc('font',size=fsize)
return
In [2]:
CCSN = DT.StellarYieldsTable('SNII')
wind = DT.StellarYieldsTable('wind')
fsize = 17
rc('text',usetex=False)
rc('font',size=fsize)
In [3]:
CCSN.x['mass'].size
Out[3]:
In [4]:
def make_manual_table(table, group, filename, datadir = './'):
#
# Load from the h5 file
#
hf = h5py.File(datadir + filename, 'r')
grp = hf[group]
table.yield_type = group
table.data_dir = datadir
table.x['mass'] = np.unique( grp['M'][()] )
table.x['metallicity'] = np.unique( grp['Z'][()] )
table.nbins['mass'] = np.size(table.x['mass'])
table.nbins['metallicity'] = np.size(table.x['metallicity'])
table._array_size = int( np.prod( list(table.nbins.values())))
for i,be in enumerate(grp['yield_names']):
e = be.decode('utf-8')
table.y[e] = grp['yields'][:,:,i]
table._available_yields = table.y.keys()
hf.close()
return
In [11]:
ddir = "/home/aemerick/code/galaxy_analysis/individualstar_data/construct_tables/"
#
# load LC18 yields
#
LC18 = {}
for ytype in ['SN','Wind','SN+Wind']:
LC18[ytype] = {}
for vr in ['R_0','R_150','R_300','mixture']:
fname = 'LC18_' + vr + '-' + 'LC18_winds_' + vr + '-FRUITY-PopIII.h5'
LC18[ytype][vr] = DT.StellarYieldsTable(ytype,fname,manual_table=True,data_dir=ddir)
make_manual_table(LC18[ytype][vr], ytype, fname, datadir = ddir)
In [ ]:
In [12]:
def interpolate_many(vals_list, fields, table):
return np.array([table.interpolate( [x,y], fields)[0] for x,y in vals_list])
def get_abundance_ratio(vals_list, element1, element2, yield_type):
#mass1 = interpolate_many(vals_list, element1, table)
#mass2 = interpolate_many(vals_list, element2, table)
aratios = np.array( [convert.get_yield_ratio(element1,element2,x,y,yield_type) for x,y in vals_list] )
return aratios
In [13]:
#
#
# LC18
#
#
def make_grid_points(table):
x,y = np.meshgrid(table.x['mass'], table.x['metallicity'])
points = [(0,0)]*np.size(table.x['mass'])*np.size(table.x['metallicity'])
i = 0
for yi in table.x['metallicity']:
for xi in table.x['mass']:
points[i] = (xi,yi)
i = i + 1
#z = get_abundance_ratio( points, 'O', 'Fe', 'SNII')
return x, y, points
In [30]:
plot_settings()
def plot_ratios_vs_mass(table, outname=None,title=False, denominators=['Fe']):
numerators = ['C','N','O','Mg','Si','Ca','Mn','Sr','Ba']
denominators = ['Fe']
ncol = 3
nrow = 3
fig,ax = plt.subplots(nrow, ncol, sharex=True,sharey=True)
fig.set_size_inches(ncol*5,nrow*5)
fig.subplots_adjust(hspace = 0, wspace = 0)
axi = 0
axj = 0
if 'SN' in outname:
xlim = (0,30)
if "Wind" in outname:
xlim = (0,125)
x,y,points = make_grid_points(table)
bins = np.arange(-30,0, 0.1)
for denom_element in denominators:
for num_element in numerators:
index = (axi,axj)
if num_element == denom_element:
axj = axj + 1
continue
#z = get_abundance_ratio( points, num_element, denom_element, table.yield_type)
z = table.interpolate_yield_ratio(num_element, denom_element, points)
newshape = (table.x['metallicity'].size, table.x['mass'].size)
ax[index].scatter(x, z.reshape(newshape), c = y, s = 50)#, alpha = 0.5)
for i in np.arange(table.x['metallicity'].size):
ax[index].plot(x[i],z.reshape(newshape)[i], color = viridis(i/(5*1.0)),
label = "[Z] = %0.2f"%(np.log10(y.reshape(newshape)[i][0]/0.0134)))
ax[index].set_xlim(xlim)
ax[index].plot(xlim,[0,0],lw=2,color='black',ls='--')
ax[index].set_ylim(-3.2,3.2)
ax[index].plot([25,25],ax[index].get_ylim(),lw=1,color='black',ls=':')
#ax[index].set_ylabel()
if axj == 0:
ax[index].set_ylabel("[X/" + denom_element + "]")
xy = (0.7,0.1)
ax[index].annotate("["+num_element + "/" + denom_element +"]",
xy,xy, xycoords = 'axes fraction')
axj = axj + 1
if axj >= ncol:
axj = 0
axi = axi + 1
plt.minorticks_on()
ax[(2,2)].legend(loc='upper right')
for i in np.arange(ncol):
ax[(nrow-1,i)].set_xlabel('Mass (Msun)')
if title:
xy = (0.15,0.9)
ax[(0,1)].annotate( outname.split('.p')[0], xy,xy,xycoords='axes fraction')
if not (outname is None):
fig.savefig(outname,bbox_inches="tight")
return
for t in ['SN+Wind','SN','Wind']:
plot_ratios_vs_mass(LC18[t]['R_0'], outname = 'LC18_'+t+'_R0.png', title=True)
plot_ratios_vs_mass(LC18[t]['R_150'], outname = 'LC18_'+t+'_R150.png',title=True)
plot_ratios_vs_mass(LC18[t]['R_300'], outname = 'LC18_'+t+'_R300.png',title=True)
plot_ratios_vs_mass(LC18[t]['mixture'], outname = 'LC18_'+t+'_mixture.png',title=True)
#plot_ratios_vs_mass(LC18['SN']['R_300'], outname = 'LC18_SN_R300.png')
In [36]:
plot_settings()
def plot_velocity_comparison(outname=None,title=False, ytype='SN+Wind'):
numerators = ['C','N','Mg','O','Mg','Ba','Na']
denominators = ['N','O','Ca','Fe','Fe','Fe','Fe']
ncol = 4
nrow = len(numerators)
fig,ax = plt.subplots(nrow, ncol, sharex=True,sharey=True)
fig.set_size_inches(ncol*5,nrow*5)
fig.subplots_adjust(hspace = 0, wspace = 0)
axi = 0
axj = 0
if ytype == 'SN':
xlim = (0,30)
if ytype == "Wind":
xlim = (0,125)
if ytype == 'SN+Wind':
xlim = (0,125)
bins = np.arange(-30,0, 0.1)
for num_element,denom_element in zip(numerators,denominators):
#z = get_abundance_ratio( points, num_element, denom_element, table.yield_type)
for tname in LC18[ytype].keys():
index = (axi,axj)
x,y,points = make_grid_points(LC18[ytype][tname])
z = LC18[ytype][tname].interpolate_yield_ratio(num_element, denom_element, points)
newshape = (LC18[ytype][tname].x['metallicity'].size, LC18[ytype][tname].x['mass'].size)
ax[index].scatter(x, z.reshape(newshape), c = y, s = 50)#, alpha = 0.5)
for i in np.arange(LC18[ytype][tname].x['metallicity'].size):
ax[index].plot(x[i],z.reshape(newshape)[i], color = viridis(i/(5*1.0)),
label = "[Z] = %0.2f"%(np.log10(y.reshape(newshape)[i][0]/0.0134)))
if axj == 0:
ax[index].set_ylabel("["+num_element + "/" + denom_element +"]")
xy = (0.1,0.9)
ax[index].annotate(tname,
xy,xy, xycoords = 'axes fraction')
ax[index].set_xlim(xlim)
ax[index].plot(xlim,[0,0],lw=2,color='black',ls='--')
ax[index].set_ylim(-3.2,3.2)
ax[index].plot([25,25],ax[index].get_ylim(),lw=1,color='black',ls=':')
#ax[index].set_ylabel()
axj = axj + 1
if axj >= ncol:
axj = 0
axi = axi + 1
plt.minorticks_on()
ax[(nrow-1,0)].legend(loc='upper right')
for i in np.arange(ncol):
ax[(nrow-1,i)].set_xlabel('Mass (Msun)')
#if title:
# fig.suptitle(outname.split('.p')[0])
if not (outname is None):
fig.savefig(outname,bbox_inches="tight")
return
plot_velocity_comparison(outname = 'LC18_SN_subset_velocity_comparison.png',title=True,ytype='SN')
plot_velocity_comparison(outname = 'LC18_Wind_subset_velocity_comparison.png',title=True,ytype='Wind')
plot_velocity_comparison(outname = 'LC18_SN+Wind_subset_velocity_comparison.png',title=True,ytype='SN+Wind')
#plot_ratios_vs_mass(LC18['SN']['R_300'], outname = 'LC18_SN_R300.png')
In [35]:
x,y = np.meshgrid(CCSN.x['mass'], CCSN.x['metallicity'])
points = [(0,0)]*np.size(CCSN.x['mass'])*np.size(CCSN.x['metallicity'])
i = 0
for yi in CCSN.x['metallicity']:
for xi in CCSN.x['mass']:
points[i] = (xi,yi)
i = i + 1
z = get_abundance_ratio( points, 'O', 'Fe', 'SNII')
In [89]:
numerators = ['O', 'Mg', 'Na', 'Ca','Ba']
denominators = ['Fe','Mg','O']
ncol = len(numerators)
nrow = len(denominators)
fig,ax = plt.subplots(nrow, ncol, sharex=True,)
fig.set_size_inches(ncol*6,nrow*6)
fig.subplots_adjust(hspace = 0)# wspace = 0)
axi = 0
axj = 0
bins = np.arange(-30,0, 0.1)
for denom_element in denominators:
for num_element in numerators:
index = (axi,axj)
if num_element == denom_element:
axj = axj + 1
continue
z = get_abundance_ratio( points, num_element, denom_element, 'SNII')
ax[index].scatter(x, z.reshape(5,12), c = y, s = 50)#, alpha = 0.5)
for i in np.arange(5):
ax[index].plot(x[i],z.reshape(5,12)[i], color = viridis(i/(5*1.0)),
label = "log$_{10}$(Z) = %0.2f"%(np.log10(y.reshape(5,12)[i][0])))
ax[index].set_xlim(0,26)
ax[index].set_ylabel("["+num_element + "/" + denom_element +"]")
axj = axj + 1
if axj >= ncol:
axj = 0
axi = axi + 1
plt.minorticks_on()
ax[(0,0)].legend(loc='best')
for i in np.arange(np.size(ncol)):
ax[(nrow-1,i)].set_xlabel('Mass (Msun)')
In [86]:
numerators = ['O', 'Mg', 'Sr','Ba']
denominators = ['Fe','Mg','Ba']
ncol = len(numerators)
nrow = len(denominators)
fig,ax = plt.subplots(nrow, ncol, sharex=True)
fig.set_size_inches(ncol*6,nrow*6)
fig.subplots_adjust(hspace = 0)# wspace = 0)
axi = 0
axj = 0
bins = np.arange(-30,0, 0.1)
for denom_element in denominators:
for num_element in numerators:
index = (axi,axj)
if num_element == denom_element:
axj = axj + 1
continue
z = get_abundance_ratio( points, num_element, denom_element, 'wind')
ax[index].scatter(x, z.reshape(5,12), c = y, s = 50)#, alpha = 0.5)
for i in np.arange(5):
ax[index].plot(x[i],z.reshape(5,12)[i], color = viridis(i/(5*1.0)),
label = "log$_{10}$(Z) = %0.2f"%(np.log10(y.reshape(5,12)[i][0])))
ax[index].set_xlim(0,26)
ax[index].set_ylabel("["+num_element + "/" + denom_element +"]")
ylim = ax[index].get_ylim()
ax[index].plot([2.0,2.0],ylim,color='black',lw=3,ls='--')
ax[index].set_ylim(ylim)
axj = axj + 1
if axj >= ncol:
axj = 0
axi = axi + 1
plt.minorticks_on()
ax[(0,0)].legend(loc='best')
for i in np.arange(np.size(ncol)):
ax[(nrow-1,i)].set_xlabel('Mass (Msun)')
In [ ]:
masses = CCSN.masses
points = zip(masses, 1.0E-4)
get_abundance_ratio( [[10.0,0.001], [15.0,0.001]], 'O','Fe','SNII')
In [ ]:
mass_values = np.linspace(1.0,100.0,128.0)
z_values = np.ones(128)*1.0E-4
#np.linspace(1.0E-4, 0.01,128.0)
mass_centers = 0.5*(mass_values[1:] + mass_values[:-1])
z_centers = 0.5*(z_values[1:] + z_values[:-1])
plot_mmesh, plot_zmesh = np.meshgrid(mass_values, z_values)
mmesh, zmesh = np.meshgrid( mass_centers, z_centers)
mflat, zflat = mmesh.flatten(), zmesh.flatten()
vals = zip(mflat, zflat)
SNII_aratios = get_abundance_ratio( vals, 'Ba','Fe','SNII').reshape(127,127)
wind_aratios = get_abundance_ratio( vals, 'Ba','Fe','wind').reshape(127,127)
In [ ]:
fig, ax = plt.subplots()
fig.set_size_inches(8,8)
p = ax.pcolormesh(plot_mmesh, plot_zmesh, SNII_aratios)
ax.set_xlim(1.0,100.0)
ax.set_ylim(0.0001,0.01)
plt.colorbar(p)
ax.colorbar
In [ ]:
fig, ax = plt.subplots()
fig.set_size_inches(8,8)
SNII_aratios[0][(mass_centers < 8.0) + (mass_centers>25.0)] = np.nan
ax.plot(mass_centers, SNII_aratios[0], label = 'SNII')
ax.set_xlim(1.0,25.0)
select = mass_centers > 8.0
ax.plot(mass_centers[select], wind_aratios[0][select], label = 'Winds')
select = mass_centers < 8.0
ax.plot(mass_centers[select], wind_aratios[0][select], label = 'AGB Winds')
ax.legend(loc='best')
In [ ]:
onez
In [5]:
In [17]:
def imf_weighted_yields(IMF, yield_table = wind):
"""
IMF is some funtion that can be evaluated to give an arbitrarily
normalized IMF weighted at some mass "m"
"""
z_values = 1.0*yield_table.x['metallicity']
y_names = yield_table.y_names()
m_values = 1.0* yield_table.x['mass']
m_values[-1] = 0.999*m_values[-1]
z_values[-1] = 0.999*z_values[-1]
print np.size(m_values),np.size(z_values)
mmesh, zmesh = np.meshgrid(m_values, z_values)
masses, z = mmesh.flatten(), zmesh.flatten()
vals = zip(masses,z)
result = {}
for x in y_names:
temp = interpolate_many(vals, x, yield_table) * IMF(masses)
result[x] = temp.reshape(np.size(z_values),np.size(m_values))
print x, result[x]
return result
In [14]:
salpeter = imf.salpeter(M_min = 1.0, M_max = 100.0)
interpolate_many([(1.0,0.0001),(22.5,0.01802)], 'O', CCSN)
Out[14]:
In [18]:
salpeter = imf.salpeter(M_min = 1.0, M_max = 100.0)
f_imf = lambda x : salpeter.imf(x)
result = imf_weighted_yields(f_imf, yield_table = CCSN)
In [19]:
result['O']
Out[19]:
In [ ]: