Test of POPIII star input
Test of SSP with POPIII yields. Focus are basic GCE features. You can find the documentation here. Note that all POPIII stars are massive stars. Hence we cannot extend the IMF down to for example 1Msun.
In [2]:
%pylab nbagg
import sygma as s
reload(s)
s.__file__
#from imp import *
#s=load_source('sygma','/home/nugrid/nugrid/SYGMA/SYGMA_online/SYGMA_dev/sygma.py')
from scipy.integrate import quad
from scipy.interpolate import UnivariateSpline
import matplotlib.pyplot as plt
import numpy as np
$\odot$ Evolution of ISM fine
$\odot$ Sources of massive and AGB stars distinguished
$\odot$ Test of final mass of ISM for different IMF boundaries
$\odot$ Test of Salpeter, Chabrier, Kroupa IMF by checking the evolution of ISM mass (incl. alphaimf)
$\odot$ Test if SNIa on/off works
$\odot$ Test of the three SNIa implementations, the evolution of SN1a contributions
$\odot$ Test of parameter tend, dt and special_timesteps
$\odot$ Test of parmeter mgal
$\odot$ Test of netyields_on
TODO: test non-linear yield fitting (hard set in code right now, no input parameter provided)
The IMF allows to calculate the number of stars $N_{12}$ in the mass interval [m1,m2] with
(I) $N_{12}$ = kN $\int {m1}^{m2} m^{-2.35} dm$
Where kN is the normalization constant. It can be derived from the total amount of mass of the system $M{tot}$ since the total mass $M_{12}$ in the mass interval above can be estimated with
(II) $M_{12}$ = kN $\int {m1}^{m2} m^{-1.35} dm$
With a total mass interval of [1,30] and $M_{tot}=1e11$ the $k_N$ can be derived:
$1e11 = k_N/0.35 * (1^{-0.35} - 30^{-0.35})$
In [60]:
k_N=1e11*0.35/ (10**-0.35 - 30**-0.35) #(I)
N_tot=k_N/1.35 * (10**-1.35 - 30**-1.35) #(II)
Yield_tot=0.1*N_tot
Includes stars from 10Msun to 30Msun (upper end consistent with higher Z).
At ~5e6 M30 star starts to contribute. Need to resolve steps between masses and hence chhose small constant time interval for s2 run. There are
In [61]:
reload(s)
s1=s.sygma(iolevel=0,mgal=1e11,dt=1e7,tend=1.3e10,special_timesteps=-1,imf_type='salpeter',
imf_yields_range_pop3=[10,30],imf_bdys_pop3=[10,30],
pop3_table='yield_tables/popIII_h1.txt',table='yield_tables/isotope_yield_table_h1.txt',
sn1a_on=False, sn1a_table='yield_tables/sn1a_h1.txt',
iniabu_table='yield_tables/iniabu/iniab1.0E-04GN93_alpha_h1.ppn')
s2=s.sygma(iolevel=0,mgal=1e11,dt=1e5,tend=1e7,special_timesteps=-1,imf_type='salpeter',
imf_yields_range_pop3=[10,30],imf_bdys_pop3=[10,30],
pop3_table='yield_tables/popIII_h1.txt',table='yield_tables/isotope_yield_table_h1.txt',
sn1a_on=False, sn1a_table='yield_tables/sn1a_h1.txt',
iniabu_table='yield_tables/iniabu/iniab1.0E-04GN93_alpha_h1.ppn')
Yield_tot_sim=s1.history.ism_iso_yield[-1][0]
In [62]:
print Yield_tot_sim
print Yield_tot
print 'ratio should be 1 : ',Yield_tot_sim/Yield_tot
No production of any source, except massive stars expected. Massive stars produce all H.
In [63]:
print s1.history.ism_iso_yield_agb[-1][0]
print 'should be 1: ',s1.history.ism_iso_yield_massive[-1][0]/Yield_tot
print 'No SNIa contribution:',s1.history.ism_iso_yield_1a[-1][0]
As expected massive stars contribute to 'All'.
In [64]:
s1.plot_totmasses(fig=1,source='all',markevery=2,marker='^')
s1.plot_totmasses(fig=1,source='agb')
s1.plot_totmasses(fig=1,source='massive',marker='x',markevery=3)
s1.plot_totmasses(fig=1,source='sn1a',marker='D')
s1.plot_mass(fig=1,specie='H-1',marker='+',markevery=1)
#plt.legend(loc=7,fontsize=14)
#mpld3.display()
#plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.xlim(1e6,1e9)
Out[64]:
Higher resolution:
In [65]:
s2.plot_totmasses(fig=2,source='all',markevery=4,marker='^')
s2.plot_totmasses(fig=2,source='massive',marker='x',markevery=6)
s2.plot_mass(fig=2,specie='H-1',marker='+',markevery=10)
#plt.legend(loc=7,fontsize=14)
#mpld3.display()
#plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.xlim(6e6,3e7)
Out[65]:
In [66]:
import read_yields as ry
y=ry.read_nugrid_yields('yield_tables/popIII_h1.txt')
zm_lifetime_grid=s1.zm_lifetime_grid_current
idx_z = (np.abs(zm_lifetime_grid[0]-0.0001)).argmin() #Z=0
grid_masses=zm_lifetime_grid[1][::-1]
grid_lifetimes=zm_lifetime_grid[2][idx_z][::-1]
plt.figure(981)
plt.plot(grid_masses,grid_lifetimes,label='spline fit grid points (SYGMA)',marker='x')
m=[]
ages=[]
for k in range(len(y.table_mz)):
m_ini=float(y.table_mz[k].split(',')[0].split('=')[1])
if m_ini>=30:
continue
m.append(m_ini)
ages.append(y.age[k])
plt.plot(np.array(m),np.log10(np.array(ages)),marker='+',markersize=20,label='input yield grid',linestyle='None')
plt.xlabel('Mini/Msun')
plt.ylabel('log lifetime');plt.legend(prop={'size':14})
Out[66]:
In [85]:
s2.plot_totmasses(fig=561,marker='x',label='totmass',markevery=2)
s2.plot_mass(fig=561,specie='H',label='H, sim',color='k',shape='-',marker='o',markevery=3)
import read_yields as ry
y=ry.read_nugrid_yields('yield_tables/popIII_h1.txt')
m=[]
ages=[]
for k in range(len(y.table_mz)):
m_ini=float(y.table_mz[k].split(',')[0].split('=')[1])
if m_ini>=30:
continue
m.append(m_ini)
ages.append(y.age[k])
#print m[-1],ages[-1]
def yields(m,k_N):
return ( k_N/1.35 * (m**-1.35 - 30.**-1.35) ) * 0.1
yields1=[]
for m1 in m:
yields1.append(yields(m1,k_N))
plt.figure(561)
plt.plot(ages,yields1,marker='+',linestyle='',markersize=15,label='H, semi')
plt.legend(loc=4)
plt.xlim(5e6,3e7)
Out[85]:
In [68]:
k_N=1e11*0.35/ (15**-0.35 - 30**-0.35)
N_tot=k_N/1.35 * (15**-1.35 - 30**-1.35)
Yield_tot=0.1*N_tot
In [69]:
# imf_bdys_pop3=[15,30]
In [70]:
##reload(chem_evol)
#dreload(s)
import sygma as s
reload(s)
s1=s.sygma(iolevel=1,mgal=1e11,iniZ=0,dt=1e7,tend=1.3e10,imf_type='salpeter',
imf_bdys=[10,30],imf_bdys_pop3=[15,30],pop3_table='yield_tables/popIII_h1.txt',
table='yield_tables/isotope_yield_table_h1.txt',sn1a_on=False, sn1a_table='yield_tables/sn1a_h1.txt',
iniabu_table='yield_tables/iniabu/iniab1.0E-04GN93_alpha_h1.ppn')
Yield_tot_sim=s1.history.ism_iso_yield[-1][0]
In [71]:
print 'Sould be 1:' ,Yield_tot_sim/Yield_tot
In [72]:
k_N=1e11*0.35/ (10**-0.35 - 15**-0.35)
N_tot=k_N/1.35 * (10**-1.35 - 15**-1.35)
Yield_tot=0.1*N_tot
In [73]:
s1=s.sygma(iolevel=1,mgal=1e11,dt=1e7,tend=1.3e10,imf_type='salpeter',imf_yields_range_pop3=[10,30],imf_bdys_pop3=[10,15],pop3_table='yield_tables/popIII_h1.txt',table='yield_tables/isotope_yield_table_h1.txt',sn1a_on=False, sn1a_table='yield_tables/sn1a_h1.txt', iniabu_table='yield_tables/iniabu/iniab1.0E-04GN93_alpha_h1.ppn')
Yield_tot_sim=s1.history.ism_iso_yield[-1][0]
In [74]:
k_N=1e11*0.35/ (5**-0.35 - 100**-0.35) # IMF range
N_tot=k_N/1.35 * (10**-1.35 - 30**-1.35) # yield range
Yield_tot=0.1*N_tot
In [75]:
s1=s.sygma(iolevel=1,mgal=1e11,dt=1e7,tend=1.3e10,imf_type='salpeter',imf_bdys_pop3=[5,100],
pop3_table='yield_tables/popIII_h1.txt',table='yield_tables/isotope_yield_table_h1.txt',
sn1a_on=False, sn1a_table='yield_tables/sn1a_h1.txt',
iniabu_table='yield_tables/iniabu/iniab1.0E-04GN93_alpha_h1.ppn')
#imf_yields_range_pop3=[10,30],
Yield_tot_sim=s1.history.ism_iso_yield[-1][0]
The code should ignore s1.imf_bdys and use s1.imf_bdys_pop3
In [76]:
print s1.imf_bdys,s1.imf_bdys_pop3,s1.imf_yields_range_pop3
In [77]:
print Yield_tot_sim
print Yield_tot
print 'Sould be 1:' ,Yield_tot_sim/Yield_tot
In [78]:
alphaimf = 1.5 #Set test alphaimf
In [79]:
k_N=1e11*(alphaimf-2)/ (-10**-(alphaimf-2) + 30**-(alphaimf-2))
N_tot=k_N/(alphaimf-1) * (-10**-(alphaimf-1) + 30**-(alphaimf-1))
Yield_tot=0.1*N_tot
In [80]:
s1=s.sygma(iolevel=0,mgal=1e11,dt=1e7,tend=1.3e10,imf_type='alphaimf',alphaimf=1.5,imf_bdys_pop3=[10,30],pop3_table='yield_tables/popIII_h1.txt',table='yield_tables/isotope_yield_table_h1.txt',sn1a_on=False, sn1a_table='yield_tables/sn1a_h1.txt', iniabu_table='yield_tables/iniabu/iniab1.0E-04GN93_alpha_h1.ppn')
Yield_tot_sim=s1.history.ism_iso_yield[-1][0]
In [81]:
print 'Should be 1 :',Yield_tot/Yield_tot_sim
In [82]:
reload(s)
s1=s.sygma(iolevel=0,mgal=1e11,dt=1e7,tend=1.3e10,imf_type='alphaimf',imf_bdys_pop3=[10,30],pop3_table='yield_tables/popIII_h1.txt',table='yield_tables/isotope_yield_table_h1.txt',sn1a_on=False, sn1a_table='yield_tables/sn1a_h1.txt', iniabu_table='yield_tables/iniabu/iniab1.0E-04GN93_alpha_h1.ppn')
s2=s.sygma(iolevel=0,mgal=1e11,dt=1e7,tend=1.3e10,imf_type='alphaimf',imf_bdys_pop3=[10,30],pop3_table='yield_tables/popIII_h1.txt',table='yield_tables/isotope_yield_table_h1.txt',sn1a_on=True, sn1a_table='yield_tables/sn1a_h1.txt', iniabu_table='yield_tables/iniabu/iniab1.0E-04GN93_alpha_h1.ppn')
In [83]:
print 'Should be 0:',(s1.history.ism_elem_yield_1a[0]),(s1.history.ism_elem_yield_1a[-1])
print (s1.history.ism_elem_yield[0]),(s1.history.ism_elem_yield[-1])
print 'Should be 0:',(s2.history.ism_elem_yield_1a[0]),(s2.history.ism_elem_yield_1a[-1])
print (s2.history.ism_elem_yield[0]),(s2.history.ism_elem_yield[-1])
print (s1.history.ism_elem_yield[-1][0] + s2.history.ism_elem_yield_1a[-1][0])/s2.history.ism_elem_yield[-1][0]
#s2.plot_mass(specie='H-1',source='sn1a') #plot s1 data (without sn) cannot be plotted -> error, maybe change plot function?
In [84]:
s0=s.sygma(iolevel=0,imf_bdys=[0.01,100],imf_yields_range=[0.02,99],imf_type='chabrier',transitionmass=6,sfr='input',iniZ=0.0,\
dt=1e7,tend=1.3e10, mgal=1e1,sn1a_on=True,sn1a_rate='exp',exp_dtd=2e9,exclude_masses=[100,6,7],netyields_on=True,pop3_table='yield_tables/popIII_h1.txt')
In [94]:
s1=s.sygma(iolevel=1,mgal=1e11,dt=1e7,tend=1.3e10,imf_type='salpeter',imf_bdys_pop3=[10,30],imf_yields_range_pop3=[20,30],iniZ=-1,pop3_table='yield_tables/popIII_h1.txt',table='yield_tables/isotope_yield_table_h1.txt',sn1a_on=False, sn1a_table='yield_tables/sn1a_h1.txt', iniabu_table='yield_tables/iniabu/iniab1.0E-04GN93_alpha_h1.ppn')
#Yield_tot_sim=s1.history.ism_iso_yield[-1][0]
In [330]:
k_N=1e11*0.35/ (10**-0.35 - 30**-0.35)
N_tot=k_N/1.35 * (20**-1.35 - 30**-1.35)
Yield_tot=0.1*N_tot
In [331]:
s1.imf_bdys,s1.imf_bdys_pop3,s1.imf_yields_range_pop3
Out[331]:
In [332]:
print 'Sould be 1:' ,Yield_tot_sim/Yield_tot
In [ ]:
yield range [20,30]