In [1]:
%matplotlib inline
import deepdish as dd
import h5py
import numpy as np
from galaxy_analysis.plot.plot_styles import *
import matplotlib.pyplot as plt
In [2]:
#run = 'NSNS1'
#run = 'HNE1'
#run = 'AGB1'
run = 'AGB1'
run2 = 'SNE1'
infile = '/home/aemerick/work/enzo_runs/mixing_experiment/' + run + '/' + run + '_gas_abundances.h5'
infile = '/home/aemerick/work/enzo_runs/mixing_experiment/' + run + '/abundances/gas_abundances.h5'
h5data = h5py.File(infile,'r')
zbins = dd.io.load(infile,'/DD0500/CNM/mass_fraction/bins')
data_list = np.sort([x for x in h5data.keys() if 'DD' in x])
infile2 = '/home/aemerick/work/enzo_runs/mixing_experiment/' + run2 + '/' + run2 + '_gas_abundances.h5'
infile2 = '/home/aemerick/work/enzo_runs/mixing_experiment/' + run2 + '/abundances/gas_abundances.h5'
h5data2 = h5py.File(infile2,'r')
data_list2 = np.sort([x for x in h5data2.keys() if 'DD' in x])
In [3]:
run_names = ['AGB1','NSNS1','SNE1','HNE1']
h5data_dict = {}
data_list_dict = {}
infile_dict = {}
for run in run_names:
#infile_dict[run] = '/home/aemerick/work/enzo_runs/mixing_experiment/' + run + '/' + run + '_gas_abundances.h5'
infile_dict[run] = '/home/aemerick/work/enzo_runs/mixing_experiment/' + run + '/abundances/gas_abundances.h5'
h5data_dict[run] = h5py.File(infile_dict[run],'r')
data_list_dict[run] = np.sort([x for x in h5data_dict[run].keys() if 'DD' in x])
abundances_to_average = {}
abundances_to_average['AGB1'] = ['C','N','O','Na','Mg','Si','S','Ca','Mn','Fe','Ni','As','Sr','Y','Ba','La','Ce','Pr','Nd']
abundances_to_average['NSNS1'] = ['Na','Si','S','Ca','Mn','Ni','As']
abundances_to_average['SNE1'] = ['Na','Si','S','Ca','Mn','Ni','As']
abundances_to_average['HNE1'] = ['Na']
data_list = data_list_dict['AGB1']
infile = infile_dict['AGB1']
h5data = h5data_dict['AGB1']
zbins = dd.io.load(infile,'/DD0500/CNM/mass_fraction/bins')
In [4]:
run_label_dict = {"AGB1":"II_E46", "NSNS1" :"II_E49", "NSNS3" : "II_E50", "SNE1" : "II_E51", "HNE1":"II_E52",
"AGB1_400" : "III_E46", "NSNS1_400" : "II_E49", "SNE1_400" : "III_E51",
"AGB1_151" : "I_E46", "NSNS1_151" : "I_E50", "SNE1_151" : "I_E51"}
In [5]:
h5data['DD0672'].keys()
Out[5]:
In [6]:
def _old_get_mass_fraction(zdex=None, mf = None, data_names = data_list, fname = infile,
SN_element = 'O', AGB_element = 'Ba', bins = zbins,
norm_stat = 'median'):
result = {'AGB' : np.zeros(np.size(data_names)),
'SN' : np.zeros(np.size(data_names)),
'time' : np.zeros(np.size(data_names))}
log_z = np.log10( 0.5 * (bins[1:] + bins[:-1]))
for i,d in enumerate(data_names):
#
# load mass fraction
#
dname = str(d)
SN = dd.io.load(fname, '/' + dname + '/CNM/mass_fraction/' + SN_element + '_Fraction')
AGB = dd.io.load(fname, '/' + dname + '/CNM/mass_fraction/' + AGB_element + '_Fraction')
# np.interp(mf, np.cumsum(SN['hist']), log_z - SN['median'])
if (norm_stat is None) or (norm_stat == ""):
SN_norm = np.zeros(np.size(log_z))
AGB_norm = np.zeros(np.size(log_z))
else:
SN_norm = np.log10(SN[norm_stat])
AGB_norm = np.log10(AGB[norm_stat])
# find the 1 dex limit
if (not (zdex is None)):
result['SN'][i]= 1.0 - np.interp(zdex, log_z - SN_norm, np.cumsum(SN['hist']))
result['AGB'][i]= 1.0 - np.interp(zdex, log_z - AGB_norm, np.cumsum(AGB['hist']))
else:
result['SN'][i] = np.interp(1.0-mf, np.cumsum( SN['hist']), log_z - SN_norm)
result['AGB'][i] = np.interp(1.0-mf, np.cumsum(AGB['hist']), log_z - AGB_norm)
result['time'][i] = dd.io.load(fname, '/'+dname+'/general/Time')
# slightly better will be to just track the time evolution of 0.10 and 0.9
# mass fractions (how far from mediam)
result['time'] = result['time'] - result['time'][0]
return result
def get_mass_fraction(zdex=None, mf = None, data_names = data_list, fname = infile,
bins = zbins, elements = None,
norm_stat = 'median', phase = 'CNM'):
if elements is None:
elements = dd.io.load(fname, '/metal_species')
elif isinstance(elements, basestring):
elements = [elements]
result = {}
for e in elements:
e = e.decode('utf-8')
result[e] = np.zeros(np.size(data_names))
result['time'] = np.zeros(np.size(data_names))
log_z = np.log10( 0.5 * (bins[1:] + bins[:-1]))
for i,d in enumerate(data_names):
#
# load mass fraction
#
dname = str(d)
for e in elements:
e = e.decode('utf-8')
mass_data = dd.io.load(fname, '/' + dname + '/' + phase + '/mass_fraction/' + e + '_Fraction')
# np.interp(mf, np.cumsum(SN['hist']), log_z - SN['median'])
if (norm_stat is None) or (norm_stat == ""):
norm = np.zeros(np.size(log_z))
else:
norm = np.log10(mass_data[norm_stat])
# find the 1 dex limit
if (not (zdex is None)):
result[e][i]= 1.0 - np.interp(zdex, log_z - SN_norm, np.cumsum(mass_data['hist']))
else:
result[e][i] = np.interp(1.0-mf, np.cumsum( mass_data['hist']), log_z - norm)
result['time'][i] = dd.io.load(fname, '/'+dname+'/general/Time')
# slightly better will be to just track the time evolution of 0.10 and 0.9
# mass fractions (how far from mediam)
result['time'] = result['time'] - result['time'][0]
return result
def get_stat(statname, data_names = data_list, fname = infile,
bins = zbins, elements = None, logval = False, average_all= False, average_type = None, phase = 'CNM'):
if elements is None:
elements = dd.io.load(fname, '/metal_species')
elif isinstance(elements, basestring):
elements = [elements]
result = {}
for e in elements:
e = e.decode('utf-8')
if statname != 'hist':
result[e] = np.zeros(np.size(data_names))
else:
result[e] = [None]*np.size(data_names)
result['time'] = np.zeros(np.size(data_names))
log_z = np.log10( 0.5 * (bins[1:] + bins[:-1]))
for i,d in enumerate(data_names):
#
# load mass fraction
#
dname = str(d)
for e in elements:
e = e.decode('utf-8')
#try:
mass_data = dd.io.load(fname, '/' + dname + '/' + phase + '/mass_fraction/' + e + '_Fraction')
#except:
# print("Failing for " + dname, e)
# continue
# np.interp(mf, np.cumsum(SN['hist']), log_z - SN['median'])
# find the 1 dex limit
result[e][i] = mass_data[statname]
result['time'][i] = dd.io.load(fname, '/'+dname+'/general/Time')
# slightly better will be to just track the time evolution of 0.10 and 0.9
# mass fractions (how far from mediam)
if logval:
for e in elements:
e = e.decode('utf-8')
result[e] = np.log10(result[e])
result['time'] = result['time'] - result['time'][0]
if average_all:# and statname != 'hist':
if statname != 'hist':
result['average'] = np.zeros(np.size(result['time']))
else:
result['average'] = np.zeros(np.shape(result[e]))
count = 0
for k in result.keys():
if average_type is None:
if k == 'time' or k =='average' or np.isnan(result[k][0]):
continue
else:
if not (k in abundances_to_average[average_type]):
continue
if statname == 'hist':
#temp = result[k]*1.0
#temp[ np.isinf(temp) ] = 0.0
temp = 10.0**(result[k])
result['average'] = result['average'] + temp
result['average'] = result['average'] + result[k]
count = count + 1
result['average'] = result['average']/(1.0 * count)
if statname == 'hist':
# result['average'][result['average']==0.0] = -np.inf
result['average'] = np.log10(result['average'])
return result
def list_stats(data_name, fname = infile):
if not (isinstance(data_name,basestring)):
if isinstance(data_name[0],basestring):
dname = data_name[0]
else:
raise ValueError
else:
dname = data_name
return (dd.io.load(fname, '/' + str(dname) +'/CNM/mass_fraction/O_Fraction')).keys()
In [ ]:
In [8]:
SN = dd.io.load('gas_abundances.h5', '/DD0500/CNM/mass_fraction/O_Fraction')
zval = np.log10(0.5*(zbins[1:]+zbins[:-1])) - np.log10(SN['median'])
xval = np.cumsum(SN['hist'])
print np.interp(1.0, zval, xval)
d = 'DD0640'
dd.io.load('gas_abundances.h5', '/'+d+'/general/Time')
In [55]:
result = old_get_mass_fraction(zdex=1.0)
result2 = old_get_mass_fraction(zdex=0.5)
result3 = old_get_mass_fraction(zdex=0.25)
In [57]:
fig, ax = plt.subplots()
fig.set_size_inches(6,6)
ax.plot(result['time'], result['SN'], lw = 3, color = 'black', label = '1.00 dex')
ax.plot(result['time'], result['AGB'], lw = 3, color = 'black', ls = '--')
ax.plot(result2['time'], result2['SN'], lw = 3, color = 'C0', label = '0.50 dex')
ax.plot(result2['time'], result2['AGB'], lw = 3, color = 'C0', ls = '--')
ax.plot(result3['time'], result3['SN'], lw = 3, color = 'C1', label = '0.25 dex')
ax.plot(result3['time'], result3['AGB'], lw = 3, color = 'C1', ls = '--')
ax.legend(loc='best')
plt.tight_layout()
ax.set_xlim(result['time'][0], result['time'][-1])
ax.set_xlabel('Time (Myr)')
ax.set_ylabel('Mass Fraction with log(Z) - log(Zmed) > X dex')
ax.semilogy()
Out[57]:
In [79]:
result = get_mass_fraction(mf=0.1)
result2 = get_mass_fraction(mf=0.9)
result3 = get_mass_fraction(mf=0.25)
result4 = get_mass_fraction(mf=0.75)
result5 = get_mass_fraction(mf=0.01)
result6 = get_mass_fraction(mf=0.99)
In [80]:
fig, ax = plt.subplots()
fig.set_size_inches(6,6)
ax.plot(result['time'], result['SN'], lw = 3, color = 'black', label = '10%')
#ax.plot(result['time'], result['AGB'], lw = 3, color = 'black', ls = '--')
ax.plot(result2['time'], result2['SN'], lw = 3, color = 'C0', label = '90%')
#ax.plot(result2['time'], result2['AGB'], lw = 3, color = 'C0', ls = '--')
ax.plot(result3['time'], result3['SN'], lw = 3, color = 'C1', label = '25%')
#ax.plot(result['time'], result['AGB'], lw = 3, color = 'black', ls = '--')
ax.plot(result4['time'], result4['SN'], lw = 3, color = 'C2', label = '75%')
#ax.plot(result2['time'], result2['AGB'], lw = 3, color = 'C0', ls = '--')
ax.plot(result5['time'], result5['SN'], lw = 3, color = 'C3', label = '1%')
ax.plot(result6['time'], result6['SN'], lw = 3, color = 'C4', label = '99%')
ax.legend(loc='best')
plt.tight_layout()
ax.set_xlim(result['time'][0], result['time'][-1])
ax.set_xlabel('Time (Myr)')
ax.set_ylabel('log(Z) - log(Z_med) [dex]')
ax.set_ylim(-1,1)
Out[80]:
In [77]:
fig, ax = plt.subplots()
fig.set_size_inches(6,6)
#ax.plot(result['time'], result['SN'], lw = 3, color = 'black', label = '10%')
ax.plot(result['time'], result['AGB'], lw = 3, color = 'black', ls = '--')
#ax.plot(result2['time'], result2['SN'], lw = 3, color = 'C0', label = '90%')
ax.plot(result2['time'], result2['AGB'], lw = 3, color = 'C0', ls = '--')
#ax.plot(result3['time'], result3['SN'], lw = 3, color = 'C1', label = '25%')
ax.plot(result3['time'], result3['AGB'], lw = 3, color = 'C1', ls = '--')
#ax.plot(result4['time'], result4['SN'], lw = 3, color = 'C2', label = '75%')
ax.plot(result4['time'], result4['AGB'], lw = 3, color = 'C2', ls = '--')
ax.legend(loc='best')
plt.tight_layout()
ax.set_xlim(result['time'][0], result['time'][-1])
ax.set_xlabel('Time (Myr)')
ax.set_ylabel('log(Z) - log(Z_med) [dex]')
ax.set_ylim(-1,1)
Out[77]:
In [18]:
list_stats(data_list)
Out[18]:
In [70]:
median = get_stat('median', logval = True, average_all = True)
In [61]:
mean = get_stat('mean', logval = True, average_all = True)
Q1 = get_stat('Q1', logval = True, average_all = True)
Q3 = get_stat('Q3', logval = True, average_all = True)
d9 = get_stat('decile_9',logval=True,average_all=True)
d1 = get_stat('decile_1',logval=True,average_all=True)
In [7]:
def _get_plot_values(phase, run, stat, species = 'average'):
if stat in all_data[phase][run].keys():
y = all_data[phase][run][stat][species]
elif stat == 'mean-median':
y = all_data[phase][run]['mean'][species] - all_data[phase][run]['median'][species]
elif stat == 'max-min':
y = all_data[phase][run]['max'][species] - all_data[phase][run]['min'][species]
elif stat == 'IQR':
y = all_data[phase][run]['Q3'][species] - all_data[phase][run]['Q1'][species]
elif stat == 'IDR':
y = all_data[phase][run]['decile_9'][species] - all_data[phase][run]['decile_1'][species]
else:
raise ValueError
time = all_data[phase][run]['mean']['time']
return time, y
In [8]:
all_data = {}
phases = ['Disk','CNM','WNM','WIM','HIM','halo']
phases = ['CNM']
for phase in phases:
all_data[phase] = {}
print(phase)
for run in run_names:
i = 0
all_data[phase][run] = {}
for stat in ['hist','mean','median','Q1','Q3','decile_9','decile_1','min','max']:
all_data[phase][run][stat] = get_stat(stat, logval = True,
average_all = True, average_type = run,
data_names = data_list_dict[run],
fname = infile_dict[run], phase = phase)
In [36]:
event_data = {}
for run in run_names + ['HNE2','NSNS3','AGB1_400','NSNS1_400','SNE1_400']:
event_data[run] = np.genfromtxt('/home/aemerick/work/enzo_runs/mixing_experiment/'+run+'/event_data_table.dat',
#+\
#run + '_event_data_table.dat',
names = True, dtype=None)
In [95]:
event_data['SNE1']
Out[95]:
In [109]:
#
# ty and plot a PDF or two
#
fig,ax = plt.subplots(4,1)
fig.set_size_inches(6,24)
def plotable_hist(x):
x2 = np.ones(np.size(x)+1)
x2[:-1] = x
x2[-1] = x[-2]
return x2
run = 'AGB1'
e = 'C'
agb_hist = all_data['CNM'][run]['hist']
median = all_data['CNM'][run]['median']
mean = all_data['CNM'][run]['mean']
times = all_data['CNM'][run]['hist']['time']
plot_times = [0.5, 10.2, 50.5, 75.5]
plot_i = [0]*np.size(plot_times)
for i,t in enumerate(plot_times):
plot_i[i] = int(np.argmin( np.abs(times - t) ))
ci = 0
for i,index in enumerate(plot_i):
ax[i].set_xlabel('log(Tracer Abundance) [Arbitrary]')
ax[i].set_ylabel('PDF')
ax[i].set_xlim(-21,-5)
ax[i].set_ylim(-4,0)
for run in ['AGB1','SNE1','NSNS1','HNE1']:
color = "C%0i"%(ci)
if run == 'AGB1':
e = 'C'
else:
e = 'Na'
for i,index in enumerate(plot_i):
agb_hist = all_data['CNM'][run]['hist']
median = all_data['CNM'][run]['median']
mean = all_data['CNM'][run]['mean']
times = all_data['CNM'][run]['hist']['time']
ax[i].step(np.log10(zbins), plotable_hist(agb_hist[e][index]), where='post', lw = 3,
color = color, label = run) #label = '%0.2f Myr'%times[index],
#for i,index in enumerate(plot_i):
ax[i].plot( [median[e][index]]*2, ax[i].get_ylim(), ls = '-', lw = 3, color = color)
ax[i].plot( [mean[e][index]]*2, ax[i].get_ylim(), ls = '--', lw = 3, color = color)
ci = ci + 1
ax[0].legend(loc='best')
plt.minorticks_on()
plt.tight_layout()
#print(zbins)
In [44]:
In [6]:
run = 'AGB1'
select = (event_data[run]['z'] == 0.0)
select = (event_data[run]['z'] > -10000)
tracers = event_data[run]['Element'][select]
r = event_data[run]['r_cyl'][select]
tracers = tracers[np.argsort(r)]
r = r[np.argsort(r)]
fig, ax = plt.subplots()
fig.set_size_inches(8,8)
radial_colors = {"0" : 'black', "100" : "C0", "300" : "C1" , "600" : "C2"}
plot_time = 10.0
yvals = np.zeros(np.size(r))
for i,t in enumerate(tracers):
data = np.genfromtxt('/home/aemerick/work/enzo_runs/mixing_experiment/' + run +\
'/' + t + '_enrichment_evolution.dat', names = True)
func = lambda x, phase : np.interp(x, data['time'], data[phase])
yvals[i] = func(plot_time,'CGM')
ax.plot(data['time'],data['CGM'], color = radial_colors["%0.0f"%(r[i])], label = "%.0f"%(r[i]))
#ax.scatter( r, yvals, label = run + " %.0f"%(plot_time), s = 20)
ax.set_xlabel('r (pc)')
ax.set_ylabel(r'f$_{\rm ej}$')
ax.set_ylim(0,1.0)
ax.set_xlim(0.0, 150.0)
ax.legend(loc='best',ncol=2)
plt.minorticks_on()
plt.tight_layout()
In [49]:
run = 'AGB1'
#select = (event_data[run]['z'] == 0.0)
fig, axes = plt.subplots(1,3)
fig.set_size_inches(18,6)
#radial_colors = {"0" : 'black', "100" : "C0", "300" : "C1" , "600" : "C2"}
for ax,run in zip(axes,['AGB1','NSNS1','SNE1']):
select = event_data[run]['z'] == 0.0
tracers = event_data[run]['Element'][select]
n = event_data[run]['n_m'][select]
for plot_time in [5,10.0, 50.0, 100.0]:
yvals = np.zeros(np.size(n))
for i,t in enumerate(tracers):
data = np.genfromtxt('/home/aemerick/work/enzo_runs/mixing_experiment/' + run +\
'/' + t + '_enrichment_evolution.dat', names = True)
func = lambda x, phase : np.interp(x, data['time'], data[phase])
yvals[i] = func(plot_time,'CGM')
#ax.plot(data['time'],data['CGM'], color = radial_colors["%0.0f"%(r[i])], label = "%.0f"%(r[i]))
ax.scatter( n, yvals, label = run + " %.0f"%(plot_time), s = 20)
ax.set_xlabel(r'<n> (cm$^{-3}$)')
ax.set_ylabel(r'f$_{\rm ej}$')
ax.set_ylim(0,1.0)
ax.set_xlim(1.0E-5,1.0E2)
ax.semilogx()
ax.legend(loc='best',ncol=1)
plt.minorticks_on()
plt.tight_layout()
In [48]:
run = 'AGB1'
#select = (event_data[run]['z'] == 0.0)
fig, axes = plt.subplots(1,3)
fig.set_size_inches(18,6)
#radial_colors = {"0" : 'black', "100" : "C0", "300" : "C1" , "600" : "C2"}
for ax,run in zip(axes,['AGB1_400','NSNS1_400','SNE1_400']):
select = event_data[run]['z'] == 0.0
tracers = event_data[run]['Element'][select]
n = event_data[run]['n_m'][select]
for plot_time in [5,10.0, 50.0, 100.0]:
yvals = np.zeros(np.size(n))
for i,t in enumerate(tracers):
data = np.genfromtxt('/home/aemerick/work/enzo_runs/mixing_experiment/' + run +\
'/' + t + '_enrichment_evolution.dat', names = True)
func = lambda x, phase : np.interp(x, data['time'], data[phase])
yvals[i] = func(plot_time,'CGM')
#ax.plot(data['time'],data['CGM'], color = radial_colors["%0.0f"%(r[i])], label = "%.0f"%(r[i]))
ax.scatter( n, yvals, label = run + " %.0f"%(plot_time), s = 20)
ax.set_xlabel(r'<n> (cm$^{-3}$)')
ax.set_ylabel(r'f$_{\rm ej}$')
ax.set_ylim(0,1.0)
ax.set_xlim(1.0E-5,1.0E2)
ax.semilogx()
ax.legend(loc='best',ncol=1)
plt.minorticks_on()
plt.tight_layout()
In [7]:
run = 'AGB1'
#select = (event_data[run]['z'] == 0.0)
fig, axes = plt.subplots(1,3)
fig.set_size_inches(18,6)
#radial_colors = {"0" : 'black', "100" : "C0", "300" : "C1" , "600" : "C2"}
k = 0
for ax,run in zip(axes,['AGB1_400','NSNS1_400','SNE1_400']):
select = event_data[run]['z'] == 0.0
tracers = event_data[run]['Element'][select]
r = event_data[run]['r_cyl'][select]
for j,plot_time in enumerate([ 10.0,100]):
yvals = np.zeros(np.size(r))
for i,t in enumerate(tracers):
data = np.genfromtxt('/home/aemerick/work/enzo_runs/mixing_experiment/' + run +\
'/' + t + '_enrichment_evolution.dat', names = True)
func = lambda x, phase : np.interp(x, data['time'], data[phase])
yvals[i] = func(plot_time,'CGM')
#ax.plot(data['time'],data['CGM'], color = radial_colors["%0.0f"%(r[i])], label = "%.0f"%(r[i]))
l = 0
for rval in np.unique(r):
x = rval
y = np.average( yvals[ r == rval])
std = np.std(yvals[r == rval])
if k == 0 and l == 0:
ax.scatter( x, y, label = " %.0f"%(plot_time), s = 20, color = "C%i"%(j))
else:
ax.scatter( x, y, s = 20, color = "C%i"%(j))
ax.errorbar( x, y, std, color = 'C%i'%(j))
l = l + 1
ax.set_xlabel(r'r (pc)')
ax.set_ylabel(r'f$_{\rm ej}$')
ax.set_ylim(0,1.0)
ax.set_xlim(-1, 610.0)
ax.legend(loc='best',ncol=1)
k = k + 1
plt.minorticks_on()
plt.tight_layout()
In [18]:
run = 'AGB1'
#select = (event_data[run]['z'] == 0.0)
runs_to_plot = ['AGB1','NSNS1','NSNS3','SNE1']
fig, all_axes = plt.subplots(3,len(runs_to_plot), sharex=True,sharey=False)
fig.set_size_inches(6*len(runs_to_plot),6*3)
fig.subplots_adjust(wspace=0,hspace=0)
#radial_colors = {"0" : 'black', "100" : "C0", "300" : "C1" , "600" : "C2"}
m = 0
for phase,axes in zip(['CGM','CNM','mean-median'],all_axes):
k = 0
for ax,run in zip(axes,runs_to_plot):
select = event_data[run]['z'] >=-9999 #== 0.0
tracers = event_data[run]['Element'][select]
r = event_data[run]['r_cyl'][select]
if m == 0:
xy = (400.0, 0.85)
ax.annotate(run_label_dict[run], xy, xytext=xy, textcoords = 'data')
for j,plot_time in enumerate([ 20.0,100]):
yvals = np.zeros(np.size(r))
for i,t in enumerate(tracers):
if phase == 'mean-median':
if run in all_data['CNM'].keys():
x,y = _get_plot_values('CNM', run, 'mean-median', t)
func = lambda xval : np.interp(xval, x, y)
yvals[i] = func(plot_time)
else:
continue
else:
data = np.genfromtxt('/home/aemerick/work/enzo_runs/mixing_experiment/' + run +\
'/' + t + '_enrichment_evolution.dat', names = True)
func = lambda x, pname : np.interp(x, data['time'], data[pname])
if phase in ['CNM','WNM','WIM','HIM']:
yvals[i] = func(plot_time,phase) / func(plot_time,'Disk')
elif phase == 'CGM':
yvals[i] = func(plot_time, phase)
#ax.plot(data['time'],data['CGM'], color = radial_colors["%0.0f"%(r[i])], label = "%.0f"%(r[i]))
l = 0
for rval in np.unique(r):
x = rval + (j*2 - 1)*7.5
select = r == rval
y = np.average( yvals[select])
std = np.std(yvals[select])
if k == 0 and l == 0:
ax.scatter( x, y, label = " %.0f Myr"%(plot_time), s = 20, color = "C%i"%(j))
else:
ax.scatter( x, y, s = 20, color = "C%i"%(j))
#print run, np.max(yvals[select])
ax.errorbar(x, y,
yerr = [[y- np.min(yvals[select])], [np.max(yvals[select]) - y]], color = 'C%i'%(j),
lw = 3)
l = l + 1
if m == 1:
ax.set_xlabel(r'r (pc)')
if m == 0 and (k == 0):
ax.set_ylabel(r'f$_{\rm ' + phase + '}$')
elif m == 1 and (k == 0):
ax.set_ylabel(r'f$_{\rm ' + phase + r'}$ / f$_{\rm Disk}$')
elif m == 2 and (k == 0):
ax.set_ylabel(r'log(mean) - log(median) [dex]')
if m == 2:
ax.set_ylim(0,14)
else:
ax.set_ylim(0,1.02)
#ax.set_ylim(1.0E-3,2.0)
#ax.semilogy()
ax.set_xlim(-10, 630.0)
k = k + 1
plt.minorticks_on()
#plt.tight_layout()
ax.minorticks_on()
m = m + 1
all_axes[(0,0)].legend(loc='lower right',ncol=1)
fig.savefig("II_radial_dependence.png")
In [90]:
#select = (event_data[run]['z'] == 0.0)
fig, all_axes = plt.subplots(2,3, sharex=True,sharey=True)
fig.set_size_inches(18,12)
fig.subplots_adjust(wspace=0,hspace=0)
#radial_colors = {"0" : 'black', "100" : "C0", "300" : "C1" , "600" : "C2"}
m = 0
for phase,axes in zip(['CGM','CNM'],all_axes):
k = 0
for ax,run in zip(axes,['AGB1_400','NSNS1_400','SNE1_400']):
select = event_data[run]['z'] >= -9999 #== 0.0
tracers = event_data[run]['Element'][select]
r = event_data[run]['r_cyl'][select]
if m == 0:
xy = (400.0, 0.85)
ax.annotate(run_label_dict[run], xy, xytext=xy, textcoords = 'data')
for j,plot_time in enumerate([ 10.0,100]):
yvals = np.zeros(np.size(r))
for i,t in enumerate(tracers):
data = np.genfromtxt('/home/aemerick/work/enzo_runs/mixing_experiment/' + run +\
'/' + t + '_enrichment_evolution.dat', names = True)
func = lambda x, pname : np.interp(x, data['time'], data[pname])
if phase in ['CNM','WNM','WIM','HIM']:
yvals[i] = func(plot_time,phase) / func(plot_time,'Disk')
elif phase == 'CGM':
yvals[i] = func(plot_time, phase)
#ax.plot(data['time'],data['CGM'], color = radial_colors["%0.0f"%(r[i])], label = "%.0f"%(r[i]))
l = 0
for rval in np.unique(r):
x = rval + (j*2 - 1)*7.5
select = r == rval
y = np.average( yvals[select])
std = np.std(yvals[select])
if k == 0 and l == 0:
ax.scatter( x, y, label = " %.0f Myr"%(plot_time), s = 20, color = "C%i"%(j))
else:
ax.scatter( x, y, s = 20, color = "C%i"%(j))
#print run, np.max(yvals[select])
ax.errorbar(x, y,
yerr = [[y- np.min(yvals[select])], [np.max(yvals[select]) - y]], color = 'C%i'%(j),
lw = 3)
l = l + 1
if m == 1:
ax.set_xlabel(r'r (pc)')
if m == 0 and (k == 0):
ax.set_ylabel(r'f$_{\rm ' + phase + '}$')
elif m == 1 and (k == 0):
ax.set_ylabel(r'f$_{\rm ' + phase + r'}$ / f$_{\rm Disk}$')
ax.set_ylim(0,1.02)
#ax.set_ylim(1.0E-3,2.0)
#ax.semilogy()
ax.set_xlim(-1, 610.0)
k = k + 1
plt.minorticks_on()
#plt.tight_layout()
ax.minorticks_on()
m = m + 1
all_axes[(0,0)].legend(loc='lower right',ncol=1)
all_axes[(0,0)].legend(loc='lower right',ncol=1)
all_axes[(0,0)].legend(loc='lower right',ncol=1)
fig.savefig("III_radial_dependence.png")
In [16]:
#select = (event_data[run]['z'] == 0.0)
runs_to_plot = ['AGB1','NSNS1','NSNS3','SNE1']
fig, all_axes = plt.subplots(3,len(runs_to_plot), sharex=True,sharey=False)
fig.set_size_inches(6*len(runs_to_plot),6 * 3)
fig.subplots_adjust(wspace=0,hspace=0)
#radial_colors = {"0" : 'black', "100" : "C0", "300" : "C1" , "600" : "C2"}
m = 0
for phase,axes in zip(['CGM','CNM','mean-median'],all_axes):
k = 0
for ax,run in zip(axes,runs_to_plot):
select = event_data[run]['z'] > -99999 #== 0.0
tracers = event_data[run]['Element'][select]
n = np.log10(event_data[run]['n_m'][select])
if m == 0:
xy = (0, 0.925)
ax.annotate(run_label_dict[run], xy, xytext=xy, textcoords = 'data')
for j,plot_time in enumerate([ 10.0,100]):
yvals = np.zeros(np.size(n))
for i,t in enumerate(tracers):
if phase == 'mean-median':
if run in all_data['CNM'].keys():
x,y = _get_plot_values('CNM', run, 'mean-median', t)
func = lambda xval : np.interp(xval, x, y)
yvals[i] = func(plot_time)
else:
continue
else:
data = np.genfromtxt('/home/aemerick/work/enzo_runs/mixing_experiment/' + run +\
'/' + t + '_enrichment_evolution.dat', names = True)
func = lambda x, pname : np.interp(x, data['time'], data[pname])
if phase in ['CNM','WNM','WIM','HIM']:
yvals[i] = func(plot_time,phase) / func(plot_time,'Disk')
elif phase == 'CGM':
yvals[i] = func(plot_time, phase)
#ax.plot(data['time'],data['CGM'], color = radial_colors["%0.0f"%(r[i])], label = "%.0f"%(r[i]))
l = 0
num_bins = 14
min_bin = -5.0
max_bin = 2.0
db = (max_bin - min_bin) / (1.0*num_bins)
for ibin in np.arange(num_bins):
x = min_bin + db * ibin
select = (n >= x) * (n < x + db)
if np.size(yvals[select]) == 0:
continue
y = np.average( yvals[select])
x = x + (j*2 - 1)*0.1*db
if l == 0:
ax.scatter( x+0.5*db, y, label = " %.0f Myr"%(plot_time), s = 20, color = "C%i"%(j))
else:
ax.scatter( x+0.5*db, y, s = 20, color = "C%i"%(j))
#print run, np.max(yvals[select])
ax.errorbar(x+0.5*db, y,
yerr = [[y- np.min(yvals[select])], [np.max(yvals[select]) - y]], color = 'C%i'%(j),
lw = 3)
l = l + 1
if m == 1:
ax.set_xlabel(r'<n> (cm$^{-3}$)')
#ax.set_ylim(1.0E-3,2.0)
#ax.semilogy()
ax.set_xlim(-5,2)
if m < 2:
ax.set_ylim(0,1.02)
if m == 2:
ax.set_ylim(0,14)
if m == 0 and (k == 0):
ax.set_ylabel(r'f$_{\rm ' + phase + '}$')
elif m == 1 and (k == 0):
ax.set_ylabel(r'f$_{\rm ' + phase + r'}$ / f$_{\rm Disk}$')
elif m == 2 and (k == 0):
ax.set_ylabel(r' log(mean) - log(median) [dex]')
k = k + 1
plt.minorticks_on()
#plt.tight_layout()
ax.minorticks_on()
m = m + 1
all_axes[(0,3)].legend(loc='lower right',ncol=1)
fig.savefig("II_n_dependence.png")
In [101]:
In [92]:
#select = (event_data[run]['z'] == 0.0)
runs_to_plot = ['AGB1_400','NSNS1_400','SNE1_400']
fig, all_axes = plt.subplots(2,len(runs_to_plot), sharex=True,sharey=True)
fig.set_size_inches(6*len(runs_to_plot),12)
fig.subplots_adjust(wspace=0,hspace=0)
#radial_colors = {"0" : 'black', "100" : "C0", "300" : "C1" , "600" : "C2"}
m = 0
for phase,axes in zip(['CGM','CNM'],all_axes):
k = 0
for ax,run in zip(axes,runs_to_plot):
select = event_data[run]['z'] > -99999 #== 0.0
tracers = event_data[run]['Element'][select]
n = np.log10(event_data[run]['n_m'][select])
if m == 0:
xy = (0, 0.925)
ax.annotate(run_label_dict[run], xy, xytext=xy, textcoords = 'data')
for j,plot_time in enumerate([ 10.0,100]):
yvals = np.zeros(np.size(n))
for i,t in enumerate(tracers):
data = np.genfromtxt('/home/aemerick/work/enzo_runs/mixing_experiment/' + run +\
'/' + t + '_enrichment_evolution.dat', names = True)
func = lambda x, pname : np.interp(x, data['time'], data[pname])
if phase in ['CNM','WNM','WIM','HIM']:
yvals[i] = func(plot_time,phase) / func(plot_time,'Disk')
elif phase == 'CGM':
yvals[i] = func(plot_time, phase)
#ax.plot(data['time'],data['CGM'], color = radial_colors["%0.0f"%(r[i])], label = "%.0f"%(r[i]))
l = 0
num_bins = 14
min_bin = -5.0
max_bin = 2.0
db = (max_bin - min_bin) / (1.0*num_bins)
for ibin in np.arange(num_bins):
x = min_bin + db * ibin
select = (n >= x) * (n < x + db)
if np.size(yvals[select]) == 0:
continue
y = np.average( yvals[select])
x = x + (j*2 - 1)*0.1*db
if k == 0 and l == 0:
ax.scatter( x+0.5*db, y, label = " %.0f Myr"%(plot_time), s = 20, color = "C%i"%(j))
else:
ax.scatter( x+0.5*db, y, s = 20, color = "C%i"%(j))
#print run, np.max(yvals[select])
ax.errorbar(x+0.5*db, y,
yerr = [[y- np.min(yvals[select])], [np.max(yvals[select]) - y]], color = 'C%i'%(j),
lw = 3)
l = l + 1
if m == 1:
ax.set_xlabel(r'<n> (cm$^{-3}$)')
if m == 0 and (k == 0):
ax.set_ylabel(r'f$_{\rm ' + phase + '}$')
elif m == 1 and (k == 0):
ax.set_ylabel(r'f$_{\rm ' + phase + r'}$ / f$_{\rm Disk}$')
ax.set_ylim(0,1.02)
#ax.set_ylim(1.0E-3,2.0)
#ax.semilogy()
ax.set_xlim(-5,2)
k = k + 1
plt.minorticks_on()
#plt.tight_layout()
ax.minorticks_on()
m = m + 1
all_axes[(0,0)].legend(loc='lower right',ncol=1)
fig.savefig("III_n_dependence.png")
In [20]:
run = 'AGB1'
#select = (event_data[run]['z'] == 0.0)
fig, axes = plt.subplots(1,3)
fig.set_size_inches(18,6)
#radial_colors = {"0" : 'black', "100" : "C0", "300" : "C1" , "600" : "C2"}
m = 0
for ax,phase in zip(axes,['CGM','CNM','mean-median']):
for k,run in enumerate(['AGB1','NSNS1','NSNS3','SNE1','HNE1','HNE2']):
select = event_data[run]['z'] == 0.0
tracers = event_data[run]['Element'][select]
E = np.log10(event_data[run]['E51'][select])
#print E
for j,plot_time in enumerate([20.0, 100.0]):
yvals = np.zeros(np.size(E))
for i,t in enumerate(tracers):
if phase == 'mean-median':
if run in all_data['CNM'].keys():
x,y = _get_plot_values('CNM', run, 'mean-median', t)
func = lambda xval : np.interp(xval, x, y)
yvals[i] = func(plot_time)
else:
continue
else:
data = np.genfromtxt('/home/aemerick/work/enzo_runs/mixing_experiment/' + run +\
'/' + t + '_enrichment_evolution.dat', names = True)
func = lambda x, phase : np.interp(x, data['time'], data[phase])
if phase == 'CGM':
yvals[i] = func(plot_time,phase)
else:
yvals[i] = func(plot_time,phase) / func(plot_time,'Disk')
#ax.plot(data['time'],data['CGM'], color = radial_colors["%0.0f"%(r[i])], label = "%.0f"%(r[i]))
l = 0
min_bin = -6.25
max_bin = 2.25
db = 0.5
num_bins = (max_bin - min_bin) / db
for ibin in np.arange(num_bins):
x = min_bin + db * ibin
select = (E >= x) * (E < x + db)
if np.size(yvals[select]) == 0:
continue
y = np.average( yvals[select])
x = x + (j*2 - 1)*0.15*db
if k == 0 and l == 0:
ax.scatter( x+0.5*db, y, label = " %.0f Myr"%(plot_time), s = 20, color = "C%i"%(j))
else:
ax.scatter( x+0.5*db, y, s = 20, color = "C%i"%(j))
#print run, np.max(yvals[select])
ax.errorbar(x+0.5*db, y,
yerr = [[y- np.min(yvals[select])], [np.max(yvals[select]) - y]], color = 'C%i'%(j),
lw = 3)
l = l + 1
ax.set_xlabel(r'log(E$_{\rm ej}$ [erg])')
if m < 2:
ax.set_ylim(0,1.0)
elif m == 2:
ax.set_ylim(0,14.0)
ax.set_xlim(-6, 2)
ax.minorticks_on()
m = m + 1
plt.tight_layout()
#ax.semilogx()
axes[0].legend(loc='lower right',ncol=1)
axes[0].set_ylabel(r'f$_{\rm CGM}$')
axes[1].set_ylabel(r'f$_{\rm CNM}$ / f$_{\rm Disk}$')
axes[2].set_ylabel(r'log(mean)-log(median) [dex]')
fig.savefig("II_Eej_CNM_avg.png")
In [23]:
#
#
#
# Now I want to plot f_ej
#
#
# gather which tracers are at the desired r or z
#
# pick any r, only midplane though
#
run = 'AGB1'
select = event_data[run]['z'] == 0.0
tracers = event_data[run]['Element'][select]
r = event_data[run]['r_cyl'][select]
tracers = tracers[np.argsort(r)]
r = r[np.argsort(r)]
print tracers
# plot evolution?
fig, ax = plt.subplots(1,3)
fig.set_size_inches(24,8)
radial_ls = ['-','-.','--',':']
radial_colors = {"0" : 'black', "100" : "C0", "300" : "C1" , "600" : "C2"}
for i,e in enumerate(tracers):
color = radial_colors["%.0f"%(r[i])]
if r[i] > 0:
ls = radial_ls[i%4]
else:
ls = '-'
#viridis( i / (1.0*np.size(tracers)))
stat = 'mean-median'
x,y = _get_plot_values('CNM', run, stat, e)
ax[0].plot(x,y,lw = 3,label = "%.1f"%(r[i]), color = color)
ax[0].set_ylabel(stat)
stat = 'mean'
x,y = _get_plot_values('CNM', run, stat, e)
ax[1].plot(x,y,lw=3, color = color, label = "%.0f"%(r[i]))
ax[1].set_ylabel(stat)
stat = 'median'
x,y = _get_plot_values('CNM', run, stat, e)
ax[2].plot(x,y,lw=3,color = color)
ax[2].set_ylabel(stat)
for a in ax:
a.set_xlabel('Time (Myr)')
a.set_xlim(0,150.0)
ax[0].set_ylim(0,14)
#ax[1].set_ylim()
ax[1].legend(loc='best', ncol = 2)
#ax
plt.minorticks_on()
plt.tight_layout()
In [24]:
#
#
#
# Now I want to plot f_ej
#
#
# gather which tracers are at the desired r or z
#
# pick any r, only midplane though
#
run = 'AGB1'
select = event_data[run]['z'] == 0.0
tracers = event_data[run]['Element'][select]
r = event_data[run]['r_cyl'][select]
tracers = tracers[np.argsort(r)]
r = r[np.argsort(r)]
print tracers
# plot evolution?
fig, ax = plt.subplots(1,3)
fig.set_size_inches(24,8)
radial_ls = ['-','-.','--',':']
radial_colors = {"0" : 'black', "100" : "C0", "300" : "C1" , "600" : "C2"}
for i,e in enumerate(tracers):
color = radial_colors["%.0f"%(r[i])]
if r[i] > 0:
ls = radial_ls[i%4]
else:
ls = '-'
#viridis( i / (1.0*np.size(tracers)))
stat = 'mean-median'
phase = 'Disk'
x,y = _get_plot_values(phase, run, stat, e)
ax[0].plot(x,y,lw = 3,label = "%.1f"%(r[i]), color = color)
ax[0].set_ylabel(stat)
stat = 'mean'
x,y = _get_plot_values(phase, run, stat, e)
ax[1].plot(x,y,lw=3, color = color, label = "%.0f"%(r[i]))
ax[1].set_ylabel(stat)
stat = 'median'
x,y = _get_plot_values(phase, run, stat, e)
ax[2].plot(x,y,lw=3,color = color)
ax[2].set_ylabel(stat)
for a in ax:
a.set_xlabel('Time (Myr)')
a.set_xlim(0,150.0)
ax[0].set_ylim(0,14)
#ax[1].set_ylim()
ax[1].legend(loc='best', ncol = 2)
#ax
plt.minorticks_on()
plt.tight_layout()
In [ ]:
print all_data['CNM']['AGB1'].keys()
In [19]:
fig, ax = plt.subplots(1,2)
fig.set_size_inches(16,8)
lines = {'AGB1' : '-', 'NSNS1' : '--', 'SNE1' : '-.' , 'HNE1' : ':'}
labels = {'AGB1' : 'II_E46', 'NSNS1' : 'II_E49' , 'SNE1' : 'II_E51', 'HNE1' : 'II_E52'}
for run in run_names:
x,y = _get_plot_values('Disk', run, 'mean-median', 'average')
ax[0].plot(x, y, color = 'black', lw = 3, ls = lines[run], label = labels[run])
for run in run_names:
x,y = _get_plot_values('CNM', run, 'mean-median', 'average')
ax[1].plot(x, y, color = color_dict['CNM'], lw = 3, ls = lines[run], label = labels[run])
for a in ax:
a.set_xlim(0,150.0)
a.set_ylim(0,14)
a.set_ylabel(r'log(Mean) - log(Median) (dex)')
a.set_xlabel(r'Time (Myr)')
ax[0].legend(ncol=2)
plt.minorticks_on()
plt.tight_layout()
ax[0].annotate("ISM", xy = (10,1), xycoords = 'data')
ax[1].annotate("CNM", xy = (10,1), xycoords = 'data')
fig.savefig("ISM_CNM_average_mean-median.png")
In [ ]:
fig, ax = plt.subplots(1,2)
fig.set_size_inches(16,8)
lines = {'AGB1_400' : '-', 'NSNS1_400' : '--', 'SNE1_400' : '-.' }
labels = {'AGB1_400' : 'III_E46', 'NSNS1_400' : 'III_E49' , 'SNE1' : 'III_E51'}
for run in run_names:
x,y = _get_plot_values('Disk', run, 'mean-median', 'average')
ax[0].plot(x, y, color = 'black', lw = 3, ls = lines[run], label = labels[run])
for run in run_names:
x,y = _get_plot_values('CNM', run, 'mean-median', 'average')
ax[1].plot(x, y, color = color_dict['CNM'], lw = 3, ls = lines[run], label = labels[run])
for a in ax:
a.set_xlim(0,150.0)
a.set_ylim(0,14)
a.set_ylabel(r'log(Mean) - log(Median) (dex)')
a.set_xlabel(r'Time (Myr)')
ax[0].legend(ncol=2)
plt.minorticks_on()
plt.tight_layout()
ax[0].annotate("ISM", xy = (10,1), xycoords = 'data')
ax[1].annotate("CNM", xy = (10,1), xycoords = 'data')
fig.savefig("ISM_CNM_400_average_mean-median.png")
In [14]:
fig, all_ax = plt.subplots(6,4)
fig.set_size_inches(24,36)
phases = ['Disk','CNM','WNM','WIM','HIM','halo']
#
# Note mean reflects ejection fraction and initial mixing mass
#
colors = {'AGB1':'C0','NSNS1':'C1','SNE1':'C2','HNE1':'C3'}
minx = 1000
maxx = -1
for ax,phase in zip(all_ax,phases):
for run in run_names:
x,y = _get_plot_values(phase, run, 'IQR', 'average')
ax[0].plot(x, y, color = colors[run], lw = 3, ls = '-', label = run)
minx = np.min( [minx,np.min(x)])
maxx = np.max( [maxx,np.max(x)])
for run in run_names:
x,y = _get_plot_values(phase, run, 'IDR', 'average')
ax[1].plot(x, y, color = colors[run], lw = 3, ls = '-', label = run)
for run in run_names:
x,y = _get_plot_values(phase, run, 'mean-median', 'average')
ax[2].plot(x, y, color = colors[run], lw = 3, ls = '-', label = run)
for run in run_names:
x,y = _get_plot_values(phase, run, 'max-min', 'average')
ax[3].plot(x, y, color = colors[run], lw = 3, ls = '-', label = run)
for a in ax:
a.set_xlabel(r'Time (Myr)')
a.set_xlim(2,maxx)
ax[0].set_ylabel(phase + 'IQR')
ax[1].set_ylabel('IDR')
ax[2].set_ylabel('Mean - Median')
ax[3].set_ylabel('Max - Min')
ax[0].legend(loc= 'best')
ax[0].set_ylim(0,14)
ax[1].set_ylim(0,14)
ax[2].set_ylim(0,14)
ax[3].set_ylim(10,30)
plt.minorticks_on()
plt.tight_layout()
In [12]:
fig, all_ax = plt.subplots(6,3)
fig.set_size_inches(18,36)
phases = ['Disk','CNM','WNM','WIM','HIM','halo']
#
# Note mean reflects ejection fraction and initial mixing mass
#
colors = {'AGB1':'C0','NSNS1':'C1','SNE1':'C2','HNE1':'C3'}
minx = 1000
maxx = -1
for ax,phase in zip(all_ax,phases):
for run in run_names:
x,y = _get_plot_values(phase, run, 'mean-median', 'average')
ax[0].plot(x, y, color = colors[run], lw = 3, ls = '-', label = run)
minx = np.min( [minx,np.min(x)])
maxx = np.max( [maxx,np.max(x)])
for run in run_names:
x,y = _get_plot_values(phase, run, 'mean', 'average')
ax[1].plot(x, y, color = colors[run], lw = 3, ls = '-', label = run)
for run in run_names:
x,y = _get_plot_values(phase, run, 'median', 'average')
ax[2].plot(x, y, color = colors[run], lw = 3, ls = '-', label = run)
for a in ax:
a.set_xlabel(r'Time (Myr)')
a.set_xlim(1,maxx)
ax[0].set_ylabel(phase + ' Mean - Median')
ax[1].set_ylabel('Mean')
ax[2].set_ylabel('Median')
ax[0].legend(loc= 'best')
ax[0].set_ylim(0,14)
ax[1].set_ylim(-15,-5)
ax[2].set_ylim(-15,-5)
plt.minorticks_on()
plt.tight_layout()
Interpretation so far from the above plots:
1) The mean and median values themselves are affected by the initial mass of the total ejecta, but the difference is NOT. The IQR and IDR are also somewhat sensitive to this in the early evolution (i.e. during the very first enrichment moments) but not afterwards. So looking at the differences should be independent of exactly how much enrichment occurs.
2) The initial very large spikes in IQR and IDR are very sensitive to the very very high metal fractions reaction by these events in the injection sites. They are able to clear out much of the ambient ISM leaving behind a very low density, metal enriched gas. This does not occur in the other events, which have more tame initial spreads.
3) The trends for mean - median vs. IQR and ID are different. This is likely because of the above. The median is less sensitive to the spikes?
4) Mean values in the ISM tend to converge in each phase, indicating that the total amount of metal in a given phase is set only by the total amount ejected from the galaxy and the time it takes for gas to change phases. The evolution of the spreads is driven entirely from changes in the distribution from mixing / homogenization.
5) For mean - median, spreads are rank ordered by inverse injection energy for all but the very initial phase of the evolution. This is again, not true for IQR and IDR. This could maybe becuase of differences in how the metals in each event are strewn about in the initial phases, but not sure. Need to think about this a bit more maybe....
NEW:
Actually it might be that IQR and IDR are bad statistics for this case. Take a case with an initial enrichment of small mass of ISM to some metal fraction Z (delta function or maybe a narrow gaussian at that Z), the PDF will be the sum of the ~delta function at zero abundance and ~delta at Z. But if M_Z <<< M_ISM, then median will be in the middle of the ~zero abundance delta, and IQR and IDR will be very narrow around it. In this case I'd argue that IQR and IDR do a bad job of capturing a measure of homogeneity. But Mean - median still will work. Maybe just look at min - max?
Maybe compute volume fraction of the ISM that contains > X % of the metals in the ISM? where X is like 0.1 or 0.01?
In [11]:
fig, all_ax = plt.subplots(6,3)
fig.set_size_inches(18,36)
phases = ['Disk','CNM','WNM','WIM','HIM','halo']
#
# Note mean reflects ejection fraction and initial mixing mass
#
colors = {'AGB1':'C0','NSNS1':'C1','SNE1':'C2','HNE1':'C3'}
minx = 1000
maxx = -1
for ax,phase in zip(all_ax,phases):
for run in run_names:
x,y = _get_plot_values(phase, run, 'mean', 'average')
ax[0].plot(x, y, color = colors[run], lw = 3, ls = '-', label = run)
minx = np.min( [minx,np.min(x)])
maxx = np.max( [maxx,np.max(x)])
for run in run_names:
x,y = _get_plot_values(phase, run, 'median', 'average')
ax[1].plot(x, y, color = colors[run], lw = 3, ls = '-', label = run)
for run in run_names:
x,y = _get_plot_values(phase, run, 'mean-median', 'average')
ax[2].plot(x, y, color = colors[run], lw = 3, ls = '-', label = run)
for a in ax:
a.set_xlabel(r'Time (Myr)')
a.set_xlim(5,maxx)
ax[0].set_ylabel(phase + 'IQR')
ax[1].set_ylabel('IDR')
ax[2].set_ylabel('Mean - Median')
ax[0].legend(loc= 'best')
plt.tight_layout()
In [71]:
fig, ax = plt.subplots()
fig.set_size_inches(8,8)
ax.plot(mean['time'], Q3['Na'], lw = 3, color = 'C0', ls = '--', label = 'SNE IQR')
ax.plot(mean['time'], Q1['Na'], lw = 3, color = 'C0', ls = '--', label = 'SNE IQR')
ax.plot(mean['time'], mean['Na'], lw = 3, color = 'C0', ls = '-', label = 'SNE IQR')
ax.plot(mean['time'], median['Na'], lw = 3, color = 'black', ls = '-', label = 'SNE IQR')
#ax.plot(mean_2['time'], Q3_2['average'], lw = 3, color = 'C2', ls = '--', label = 'HNE IQR')
#ax.plot(mean_2['time'], Q1_2['average'], lw = 3, color = 'C2', ls = '--', label = 'HNE IQR')
#ax.plot(mean_2['time'], mean_2['average'], lw = 3, color = 'C2', ls = '-', label = 'HNE IQR')
ax.set_xlim(mean['time'][0], mean['time'][-1])
plt.tight_layout()
ax.set_xlabel(r'Time (Myr)')
ax.set_ylabel(r'Stat')
Out[71]:
In [ ]:
In [ ]: