This notebook will help you to access the chemical abundance tracks that you see in paper 1 figure 17 but for all elements. You can use it to compare it to your own model/data.



In [1]:

    
%pylab inline









    



Populating the interactive namespace from numpy and matplotlib

We load the results from the best parameter Chempy run:



In [2]:

    
# Single zone models for sun, arcturus and cas with default and alternative yield set
sun_def_ab = np.load('data/paper_1_abundance_tracks/single_sun_default_abundances.npy')
sun_def_cube = np.load('data/paper_1_abundance_tracks/single_sun_default_cube.npy')
sun_alt_ab = np.load('data/paper_1_abundance_tracks/single_sun_alternative_abundances.npy')
sun_alt_cube = np.load('data/paper_1_abundance_tracks/single_sun_alternative_cube.npy')
arc_def_ab = np.load('data/paper_1_abundance_tracks/single_arc_default_abundances.npy')
arc_def_cube = np.load('data/paper_1_abundance_tracks/single_arc_default_cube.npy')
arc_alt_ab = np.load('data/paper_1_abundance_tracks/single_arc_alternative_abundances.npy')
arc_alt_cube = np.load('data/paper_1_abundance_tracks/single_arc_alternative_cube.npy')
cas_def_ab = np.load('data/paper_1_abundance_tracks/single_cas_default_abundances.npy')
cas_def_cube = np.load('data/paper_1_abundance_tracks/single_cas_default_cube.npy')
cas_alt_ab = np.load('data/paper_1_abundance_tracks/single_cas_alternative_abundances.npy')
cas_alt_cube = np.load('data/paper_1_abundance_tracks/single_cas_alternative_cube.npy')
# multizone models for sun, arcturus and cas with default and alternative yield set
mult_sun_def_ab = np.load('data/paper_1_abundance_tracks/sun_default_abundances.npy')
mult_sun_def_cube = np.load('data/paper_1_abundance_tracks/sun_default_cube.npy')
mult_sun_alt_ab = np.load('data/paper_1_abundance_tracks/sun_alternative_abundances.npy')
mult_sun_alt_cube = np.load('data/paper_1_abundance_tracks/sun_alternative_cube.npy')
mult_arc_def_ab = np.load('data/paper_1_abundance_tracks/arc_default_abundances.npy')
mult_arc_def_cube = np.load('data/paper_1_abundance_tracks/arc_default_cube.npy')
mult_arc_alt_ab = np.load('data/paper_1_abundance_tracks/arc_alternative_abundances.npy')
mult_arc_alt_cube = np.load('data/paper_1_abundance_tracks/arc_alternative_cube.npy')
mult_cas_def_ab = np.load('data/paper_1_abundance_tracks/cas_default_abundances.npy')
mult_cas_def_cube = np.load('data/paper_1_abundance_tracks/cas_default_cube.npy')
mult_cas_alt_ab = np.load('data/paper_1_abundance_tracks/cas_alternative_abundances.npy')
mult_cas_alt_cube = np.load('data/paper_1_abundance_tracks/cas_alternative_cube.npy')

These are the available elements (not all are trustworthy)



In [3]:

    
print(sun_def_ab.dtype.names)









    



('H', 'He', 'Li', 'Be', 'B', 'C', 'N', 'O', 'F', 'Ne', 'Na', 'Mg', 'Al', 'Si', 'P', 'S', 'Cl', 'Ar', 'K', 'Ca', 'Sc', 'Ti', 'V', 'Cr', 'Mn', 'Fe', 'Co', 'Ni', 'Cu', 'Zn', 'Ga', 'Ge', 'As', 'Se', 'Br', 'Kr', 'Rb', 'Sr', 'Y', 'Zr', 'Nb', 'Mo', 'weights')

These are the corresponding time-steps in Gyrs, 13.5 being present-day.



In [4]:

    
print(sun_def_cube['time'])









    



[  0.    0.1   0.2   0.3   0.4   0.5   0.6   0.7   0.8   0.9   1.    1.1
   1.2   1.3   1.4   1.5   1.6   1.7   1.8   1.9   2.    2.1   2.2   2.3
   2.4   2.5   2.6   2.7   2.8   2.9   3.    3.1   3.2   3.3   3.4   3.5
   3.6   3.7   3.8   3.9   4.    4.1   4.2   4.3   4.4   4.5   4.6   4.7
   4.8   4.9   5.    5.1   5.2   5.3   5.4   5.5   5.6   5.7   5.8   5.9
   6.    6.1   6.2   6.3   6.4   6.5   6.6   6.7   6.8   6.9   7.    7.1
   7.2   7.3   7.4   7.5   7.6   7.7   7.8   7.9   8.    8.1   8.2   8.3
   8.4   8.5   8.6   8.7   8.8   8.9   9.    9.1   9.2   9.3   9.4   9.5
   9.6   9.7   9.8   9.9  10.   10.1  10.2  10.3  10.4  10.5  10.6  10.7
  10.8  10.9  11.   11.1  11.2  11.3  11.4  11.5  11.6  11.7  11.8  11.9
  12.   12.1  12.2  12.3  12.4  12.5  12.6  12.7  12.8  12.9  13.   13.1
  13.2  13.3  13.4  13.5]

Here is how you can plot them:



In [5]:

    
plt.plot(sun_def_ab['Fe'][1:],sun_def_ab['Mg'][1:]-sun_def_ab['Fe'][1:], label = 'single sun')
plt.plot(mult_sun_def_ab['Fe'][1:],mult_sun_def_ab['Mg'][1:]-mult_sun_def_ab['Fe'][1:], label = 'multi sun')
plt.plot(arc_def_ab['Fe'][1:],arc_def_ab['Mg'][1:]-arc_def_ab['Fe'][1:], label = 'single arc')
plt.plot(mult_arc_def_ab['Fe'][1:],mult_arc_def_ab['Mg'][1:]-mult_arc_def_ab['Fe'][1:], label = 'multi arc')
plt.plot(cas_def_ab['Fe'][1:],cas_def_ab['Mg'][1:]-cas_def_ab['Fe'][1:], label = 'single cas')
plt.plot(mult_cas_def_ab['Fe'][1:],mult_cas_def_ab['Mg'][1:]-mult_cas_def_ab['Fe'][1:], label = 'multi cas')
plt.plot(sun_alt_ab['Fe'][1:],sun_alt_ab['Mg'][1:]-sun_alt_ab['Fe'][1:],linestyle = '--', label = 'single sun alternative')
plt.plot(mult_sun_alt_ab['Fe'][1:],mult_sun_alt_ab['Mg'][1:]-mult_sun_alt_ab['Fe'][1:],linestyle = '--', label = 'multi sun alternative')
plt.plot(arc_alt_ab['Fe'][1:],arc_alt_ab['Mg'][1:]-arc_alt_ab['Fe'][1:],linestyle = '--', label = 'single arc alternative')
plt.plot(mult_arc_alt_ab['Fe'][1:],mult_arc_alt_ab['Mg'][1:]-mult_arc_alt_ab['Fe'][1:],linestyle = '--', label = 'multi arc alternative')
plt.plot(cas_alt_ab['Fe'][1:],cas_alt_ab['Mg'][1:]-cas_alt_ab['Fe'][1:],linestyle = '--', label = 'single cas alternative')
plt.plot(mult_cas_alt_ab['Fe'][1:],mult_cas_alt_ab['Mg'][1:]-mult_cas_alt_ab['Fe'][1:],linestyle = '--', label = 'multi cas alternative')
plt.xlabel('[Fe/H]')
plt.ylabel('[Mg/Fe]')
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))









    Out[5]:





<matplotlib.legend.Legend at 0x7f84c2b80940>



In [6]:

    
plt.plot(sun_def_cube['time'][1:],sun_def_ab['Mg'][1:]-sun_def_ab['Fe'][1:])
plt.xlabel('time in Gyr')
plt.ylabel('[Mg/Fe]')









    Out[6]:





Text(0,0.5,'[Mg/Fe]')

You can do this for all the available elements. I would see the results for [Fe/H] < 1 as rough extrapolations. Beware that you have to weight in the SFR and age-distribution of your tracer stars if you want to compare e.g. to metallicity distribution function. See https://github.com/jan-rybizki/Chempy/blob/master/tutorials/5-Chempy_function.ipynb the section called "A note on chemical evolution tracks and 'by eye' fit".



In [7]:

    
print(sun_def_cube.dtype.names)









    



('sfr', 'infall', 'time', 'feedback', 'mass_in_remnants', 'stars', 'gas', 'Z', 'alpha', 'sn1a', 'sn2', 'pn', 'bh', 'hn', 'Al', 'Ar', 'As', 'B', 'Be', 'Br', 'C', 'Ca', 'Cl', 'Co', 'Cr', 'Cu', 'F', 'Fe', 'Ga', 'Ge', 'H', 'He', 'K', 'Kr', 'Li', 'Mg', 'Mn', 'Mo', 'N', 'Na', 'Nb', 'Ne', 'Ni', 'O', 'P', 'Rb', 'S', 'Sc', 'Se', 'Si', 'Sr', 'Ti', 'V', 'Y', 'Zn', 'Zr')

You can of course also compare to all other values saved in the cube e.g. SN rates or sfr/gas/stellar mass



In [8]:

    
plt.plot(sun_def_cube['time'],sun_def_cube['sfr'], label = "SFR")
plt.plot(sun_def_cube['time'],sun_def_cube['infall'], label = "infall")
plt.plot(sun_def_cube['time'],sun_def_cube['stars'], label = "stars")
plt.plot(sun_def_cube['time'],sun_def_cube['gas'], label = "gas")
plt.plot(sun_def_cube['time'],sun_def_cube['sn1a'], label = "sn1a")
plt.plot(sun_def_cube['time'],sun_def_cube['sn2'], label = "sn2")
plt.plot(sun_def_cube['time'],sun_def_cube['mass_in_remnants'], label = "mass in remnants")
plt.yscale('log')
plt.legend()









    Out[8]:





<matplotlib.legend.Legend at 0x7f84c2a39da0>