In [ ]:
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import sygma as s
import os
print (s.__file__)
%matplotlib inline
In [ ]:
%run read_wiersma_figures
In [ ]:
## Parameter for standard best fit
iniZ_st=0.0126
imf_type_st='chabrier'
imf_bdys_st=[0.1,100] #based on paper p. 18 (and tests on p.19)
imf_yields_range_st=[0.8,100] #0.8 based on calculations in paper p.19, paper uses 40Msun but then messses up plots, 100Msun good agreement.
sn1a_rate_st='exp' #based on paper .p. 23
exp_dtd_st=2e9 #based on paper p. 23
direct_norm_1a_st=0.02 #0.016 #better fit #0.02, #based on paper p. 23 (a parameter, see also p22)
dt_st=7e6
transitionmass_st=8 #6.8 #chosen to match the AGB and massive star intersection.
yield_interp_st='wiersma' # means we correct for initial abundance as in wiersma (independent if yields are net yields or not.)
netyields_on_st=True #tells the code we actual use tabula.ted net yields.
wiersmamod_st=True
table_st='yield_tables/agb_and_massive_stars_portinari98_marigo01_net_yields.txt'
exclude_masses_st=[] #needs to be empty else 60Msun exlcuded!, cannot justify NOT to use 6,7Msun
sn1a_table_st='yield_tables/sn1a_t03.txt'
iniabu_table_st='yield_tables/iniabu/iniab_solar_Wiersma_iso.ppn'
mgal_st = 1e4 #made to sygma.py defualt, #reasonable choice, high resolution of mass of a particle ejected in a hydro simulation.
hardsetZ = 0.02
In [ ]:
s0=s.sygma(iolevel=0,iniZ=iniZ_st,imf_type=imf_type_st,
imf_bdys=imf_bdys_st,
imf_yields_range=imf_yields_range_st,
sn1a_rate=sn1a_rate_st,
exp_dtd=exp_dtd_st,
direct_norm_1a=direct_norm_1a_st,
dt=dt_st,
transitionmass=transitionmass_st,
yield_interp=yield_interp_st,
netyields_on=netyields_on_st,
wiersmamod=wiersmamod_st,
table=table_st, \
exclude_masses=exclude_masses_st,
sn1a_table=sn1a_table_st,iniabu_table=iniabu_table_st,\
hardsetZ=hardsetZ)
In [ ]:
#%matplotlib nbagg
fign = 1
markevery=5
plt.close(2)
props = dict(boxstyle='square', facecolor='w', alpha=1)
f, axarr = plt.subplots(3, 2)
elements=['Total','C','N','O','Si','Fe']
for specie in elements:
specieW=specie
if specie=='Total':
plt.sca(axarr[0, 0])
if specie=='C':
plt.sca(axarr[0, 1])
elif specie=='N':
plt.sca(axarr[1, 0])
elif specie=='O':
plt.sca(axarr[1, 1])
elif specie=='Si':
plt.sca(axarr[2, 0])
elif specie=='Fe':
plt.sca(axarr[2, 1])
if not 'G-1' in elements and specie=='Total':
specieW='tot'
s0.plot_totmasses(fig=fign,source='agb',norm='current',color='b',
label='AGB, Z=1e-4',shape='--',marker='x',log=False,markevery=markevery)
s0.plot_totmasses(fig=1,source='massive',norm='current',color='b',
label='Massive, Z=1e-4',shape='-',marker='x',log=False,markevery=markevery)
s0.plot_totmasses(fig=fign,source='sn1a',norm='current',color='b',
label='SN1a, Z=1e-4',shape=':',marker='x',log=False,markevery=markevery)
[x,y]=getW(fig=2,specie=specieW,source='agb',Z=0.02)
plt.plot(x,y,color='r',linestyle='--')
[x,y]=getW(fig=2,specie=specieW,source='massive',Z=0.02)
plt.plot(x,y,color='r')
else:
if specie=='G-1':
specieW='tot'
s0.plot_mass(fig=fign,specie=specie,source='agb',norm='current',color='b',
label=specie+' (AGB), Z=2e-2',shape='--',marker='x',markevery=markevery)
s0.plot_mass(fig=fign,specie=specie,source='massive',norm='current',color='b',
label=specie+' (massive), Z=2e-2',shape='-',marker='x',markevery=markevery)
s0.plot_mass(fig=fign,specie=specie,source='sn1a',norm='current',color='b',
label=specie+' (SN1a), Z=2e-2',shape=':',marker='x',markevery=markevery)
[x,y]=getW(fig=2,specie=specieW,source='agb',Z=0.02)
plt.plot(x,y,color='r',linestyle='--')
[x,y]=getW(fig=2,specie=specieW,source='massive',Z=0.02)
plt.plot(x,y,color='r',linestyle='-')
if not specie in ['N','C','O','G-1']:
[x,y]=getW(fig=2,specie=specieW,source='sn1a',Z=0.02)
plt.plot(x,y,color='r',linestyle=':')
ax=plt.gca()
ax.text(0.1, 0.85, specie, transform=ax.transAxes, fontsize=16,verticalalignment='top', bbox=props)
ax.minorticks_on()
plt.title('')
plt.legend().set_visible(False)
if (specie=='N'):
ax.set_ylabel('Fraction of accumulated ejecta',size=16)
elif not ( ((specie=='Total') or ( specie=='N')) or ( specie=='Si')):
#ax.get_yaxis().set_ticklabels([]) #set_visible(False)
plt.ylabel('')
else:
plt.ylabel('')
if specie=='Total':
plt.plot([],[],marker='',linestyle='--',color='k',label='AGB',linewidth=3)
plt.plot([],[],marker='',linestyle='-',color='k',label='Massive',linewidth=3)
plt.plot([],[],marker='',linestyle=':',color='k',label='SNIa',linewidth=3)
l2=ax.lines[-3:]
leg=plt.legend(l2,('AGB','Massive','SNIa'),loc=6,prop={'size': 14})
leg.set_visible(True)
leg.get_frame().set_edgecolor('k')
leg.get_frame().set_alpha(1.)
if specie in ['Total','C','N','O']:
ax.get_xaxis().set_ticklabels([])
plt.xlabel('')
else:
plt.xlabel('SSP age [yr]',size=16)
plt.ylim(0,1)
#general settings
#plt.subplots_adjust(hspace=0.1, bottom=0.125,wspace=0.05,left=0.05,right=0.95)
plt.subplots_adjust(hspace=0.1, bottom=0.125,wspace=0.15,left=0.1,right=0.9)
fig=plt.gcf()
fig.set_size_inches(14,10,forward=True) #figsize=(6, 4)
#plt.tight_layout()
matplotlib.rcParams.update({'font.size': 16})
In [ ]:
s1=s.sygma(iolevel=0,iniZ=iniZ_st,imf_type=imf_type_st,
imf_bdys=imf_bdys_st,
imf_yields_range=imf_yields_range_st,
sn1a_rate=sn1a_rate_st,
exp_dtd=exp_dtd_st,
direct_norm_1a=direct_norm_1a_st,
dt=dt_st,
transitionmass=7.5,
yield_interp=yield_interp_st,
netyields_on=netyields_on_st,
wiersmamod=wiersmamod_st,
table=table_st, \
exclude_masses=exclude_masses_st,
sn1a_table=sn1a_table_st,iniabu_table=iniabu_table_st,\
hardsetZ=hardsetZ)
s2=s.sygma(iolevel=0,iniZ=iniZ_st,imf_type=imf_type_st,
imf_bdys=imf_bdys_st,
imf_yields_range=imf_yields_range_st,
sn1a_rate=sn1a_rate_st,
exp_dtd=exp_dtd_st,
direct_norm_1a=direct_norm_1a_st,
dt=dt_st,
transitionmass=8.5,
yield_interp=yield_interp_st,
netyields_on=netyields_on_st,
wiersmamod=wiersmamod_st,
table=table_st, \
exclude_masses=exclude_masses_st,
sn1a_table=sn1a_table_st,iniabu_table=iniabu_table_st,\
hardsetZ=hardsetZ)
In [ ]:
fign = fign + 1
figW=2
fig=plt.figure(fign)
ax=plt.gca()
Z=0.02
specie='tot'
[age,y_agb]=getW(fig=figW,specie=specie,source='agb',Z=Z)
l1,=ax.plot(age,y_agb,color='r',linestyle='--',linewidth=4,label='AGB, W09')
[age,y_massive]=getW(fig=figW,specie=specie,source='massive',Z=Z)
l2,=ax.plot(age,y_massive,color='r',linestyle='-',linewidth=4,label='Massive, W09')
#ax.add_artist(leg)
leg=ax.legend(loc=1,fontsize=17)
leg.get_frame().set_edgecolor('k')
leg.get_frame().set_alpha(1.)
s2.plot_totmasses(fig=fign,source='agb',color='g',shape='--',norm='current',fsize=[10,6],marker='s')
s2.plot_totmasses(fig=fign,source='massive',color='g',shape='-',norm='current',fsize=[10,6],marker='o')
s0.plot_totmasses(fig=fign,source='agb',color='k',shape='--',norm='current',fsize=[10,6],marker='^')
s0.plot_totmasses(fig=fign,source='massive',color='k',shape='-',norm='current',fsize=[10,6],marker='>')
s1.plot_totmasses(fig=fign,source='agb',color='b',shape='--',norm='current',fsize=[10,6],marker='+')
s1.plot_totmasses(fig=fign,source='massive',color='b',shape='-',norm='current',fsize=[10,6],marker='x')
plt.yscale('linear')
plt.ylabel('IMF-weighted fraction of the ejecta')
plt.title('')
ax.minorticks_on()
l2=ax.lines[2:]
ax.add_artist(leg)
leg=ax.legend(l2,('AGB, '+str(s2.transitionmass)+'M$_{\odot}$','Massive, '+str(s2.transitionmass)+'M$_{\odot}$',
'AGB, '+str(s0.transitionmass)+'M$_{\odot}$','Massive, '+str(s0.transitionmass)+'M$_{\odot}$',
'AGB, '+str(s1.transitionmass)+'M$_{\odot}$','Massive, '+str(s1.transitionmass)+'M$_{\odot}$')
,loc=6,prop={'size': 13})
leg.get_frame().set_edgecolor('k')
leg.get_frame().set_alpha(1.)
ax.set_ylim(0,1)
fig.set_size_inches(8,5,forward=True)
plt.tight_layout()
plt.xlim(2e6,13e9)
plt.ylabel('Fraction of accumulated ejecta', fontsize=19)
plt.xlabel('SSP age [yr]', fontsize=19)
fig=plt.gcf();
leg=plt.legend(loc=6,fontsize=17)
leg.get_frame().set_edgecolor('k')
leg.get_frame().set_alpha(1.)
plt.tight_layout()
fig.set_size_inches(8,6,forward=True)
plt.subplots_adjust(hspace=0.1, bottom=0.225, wspace=0.05,left=0.15,right=0.9)
matplotlib.rcParams.update({'font.size': 19})
In [ ]: