In [1]:

    

%matplotlib inline
from galaxy_analysis.plot.plot_styles import *
import matplotlib.pyplot as plt
from galaxy_analysis.analysis import Galaxy
import numpy as np
import yt
from galaxy_analysis.utilities import convert_abundances


    

In [2]:

    

#gal = Galaxy('DD0559', wdir = '/home/emerick/work/enzo_runs/pleiades/starIC/run11_30km_3pc/otrad_ion-nosn/')

#gal = Galaxy('DD0752', wdir = '/home/emerick/work/enzo_runs/pleiades/starIC/run11_30km/final_sndriving/') #_3pc/otrad_ion-nosn/')

gal = Galaxy('DD1053', wdir = '/home/aemerick/work/enzo_runs/leo_p/fiducial/sn_H2atten_H2sh')


    

    




Parsing Hierarchy : 100%|██████████| 1109/1109 [00:00<00:00, 15397.72it/s]
/home/aemerick/.local/lib/python3.7/site-packages/yt/fields/local_fields.py:46: UserWarning: Because 'sampling_type' not specified, yt will assume a cell 'sampling_type'
  warnings.warn("Because 'sampling_type' not specified, yt will "






    




tracer species present:  ['C', 'N', 'O', 'Na', 'Mg', 'Si', 'S', 'Ca', 'Mn', 'Fe', 'Ni', 'As', 'Sr', 'Y', 'Ba']
16 mass fields defined
16 mass fraction fields defined
15 number density fields defined
141 abundance ratio fields defined
141 particle abundance ratio fields defined
5 additional helper fields defined






    




Parsing Hierarchy : 100%|██████████| 1109/1109 [00:00<00:00, 8686.68it/s]
/home/aemerick/code/onezone/onezone/data_tables.py:98: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
  if  c == flag or id == flag:

In [3]:

    

def is_alive(data, t):
    
    lifetime = data[('io','particle_model_lifetime')].convert_to_units('Myr')
    t0       = data['creation_time'].convert_to_units('Myr')
    alive    = np.array([True] * np.size(t0))
    
    not_born    = t <  t0
    dead        = t > (t0 + lifetime)
    alive[not_born] = False
    alive[dead]     = False
    
    return alive


    

In [4]:

    

gal.ds.field_list


    

    
Out[4]:





[('all', 'birth_mass'),
 ('all', 'creation_time'),
 ('all', 'dynamical_time'),
 ('all', 'metallicity_fraction'),
 ('all', 'particle_As_fraction'),
 ('all', 'particle_Ba_fraction'),
 ('all', 'particle_C_fraction'),
 ('all', 'particle_Ca_fraction'),
 ('all', 'particle_Fe_fraction'),
 ('all', 'particle_H_fraction'),
 ('all', 'particle_He_fraction'),
 ('all', 'particle_Mg_fraction'),
 ('all', 'particle_Mn_fraction'),
 ('all', 'particle_N_fraction'),
 ('all', 'particle_Na_fraction'),
 ('all', 'particle_Ni_fraction'),
 ('all', 'particle_O_fraction'),
 ('all', 'particle_S_fraction'),
 ('all', 'particle_Si_fraction'),
 ('all', 'particle_Sr_fraction'),
 ('all', 'particle_Y_fraction'),
 ('all', 'particle_index'),
 ('all', 'particle_mass'),
 ('all', 'particle_position_x'),
 ('all', 'particle_position_y'),
 ('all', 'particle_position_z'),
 ('all', 'particle_type'),
 ('all', 'particle_velocity_x'),
 ('all', 'particle_velocity_y'),
 ('all', 'particle_velocity_z'),
 ('all', 'rad_table_T_pos'),
 ('all', 'rad_table_Z_pos'),
 ('all', 'rad_table_g_pos'),
 ('all', 'se_table_M_pos'),
 ('all', 'se_table_Z_pos'),
 ('all', 'sn_mass_ejected'),
 ('all', 'wind_mass_ejected'),
 ('all', 'yield_table_M_pos'),
 ('all', 'yield_table_Z_pos'),
 ('enzo', 'As_Density'),
 ('enzo', 'Ba_Density'),
 ('enzo', 'C_Density'),
 ('enzo', 'Ca_Density'),
 ('enzo', 'Dark_Matter_Density'),
 ('enzo', 'Density'),
 ('enzo', 'Electron_Density'),
 ('enzo', 'Fe_Density'),
 ('enzo', 'GasEnergy'),
 ('enzo', 'GravPotential'),
 ('enzo', 'H2II_Density'),
 ('enzo', 'H2I_Density'),
 ('enzo', 'H2I_kdiss'),
 ('enzo', 'HII_Density'),
 ('enzo', 'HI_Density'),
 ('enzo', 'HI_kph'),
 ('enzo', 'HM_Density'),
 ('enzo', 'HeIII_Density'),
 ('enzo', 'HeII_Density'),
 ('enzo', 'HeII_kph'),
 ('enzo', 'HeI_Density'),
 ('enzo', 'HeI_kph'),
 ('enzo', 'Metal_Density'),
 ('enzo', 'Mg_Density'),
 ('enzo', 'Mn_Density'),
 ('enzo', 'N_Density'),
 ('enzo', 'Na_Density'),
 ('enzo', 'Ni_Density'),
 ('enzo', 'OTLW_kdissH2I'),
 ('enzo', 'O_Density'),
 ('enzo', 'Pe_heating_rate'),
 ('enzo', 'PhotoGamma'),
 ('enzo', 'RadAccel1'),
 ('enzo', 'RadAccel2'),
 ('enzo', 'RadAccel3'),
 ('enzo', 'Ray_Segments'),
 ('enzo', 'S_Density'),
 ('enzo', 'Si_Density'),
 ('enzo', 'Sr_Density'),
 ('enzo', 'Temperature'),
 ('enzo', 'TotalEnergy'),
 ('enzo', 'Y_Density'),
 ('enzo', 'x-velocity'),
 ('enzo', 'y-velocity'),
 ('enzo', 'z-velocity'),
 ('io', 'birth_mass'),
 ('io', 'creation_time'),
 ('io', 'dynamical_time'),
 ('io', 'metallicity_fraction'),
 ('io', 'particle_As_fraction'),
 ('io', 'particle_Ba_fraction'),
 ('io', 'particle_C_fraction'),
 ('io', 'particle_Ca_fraction'),
 ('io', 'particle_Fe_fraction'),
 ('io', 'particle_H_fraction'),
 ('io', 'particle_He_fraction'),
 ('io', 'particle_Mg_fraction'),
 ('io', 'particle_Mn_fraction'),
 ('io', 'particle_N_fraction'),
 ('io', 'particle_Na_fraction'),
 ('io', 'particle_Ni_fraction'),
 ('io', 'particle_O_fraction'),
 ('io', 'particle_S_fraction'),
 ('io', 'particle_Si_fraction'),
 ('io', 'particle_Sr_fraction'),
 ('io', 'particle_Y_fraction'),
 ('io', 'particle_index'),
 ('io', 'particle_mass'),
 ('io', 'particle_position_x'),
 ('io', 'particle_position_y'),
 ('io', 'particle_position_z'),
 ('io', 'particle_type'),
 ('io', 'particle_velocity_x'),
 ('io', 'particle_velocity_y'),
 ('io', 'particle_velocity_z'),
 ('io', 'rad_table_T_pos'),
 ('io', 'rad_table_Z_pos'),
 ('io', 'rad_table_g_pos'),
 ('io', 'se_table_M_pos'),
 ('io', 'se_table_Z_pos'),
 ('io', 'sn_mass_ejected'),
 ('io', 'wind_mass_ejected'),
 ('io', 'yield_table_M_pos'),
 ('io', 'yield_table_Z_pos')]

In [31]:

    

t = 1250.0 * yt.units.Myr
select = is_alive(gal.df, t)

elements = ['N','O','Sr','Mg','Mn','Fe','Ba']

nrow = 5
ncol = 7
fig,ax = plt.subplots(nrow, ncol)
fig.set_size_inches(ncol*5,nrow*5)

axi = 0
axj = 0
bins = np.arange(-30,0, 0.1)
norm = 1.0

for denom in [None, 'H','Fe','Mg','N']:
    
    for e in elements:
        index = (axi,axj)
        
        if denom is None:
            field = 'particle_' + e + '_fraction'
            xdata = np.log10(gal.df[('all',field)])
        else:
            if e == denom:
                xdata = np.zeros(np.size(xdata))
            else:
                field = 'particle_' + e + '_over_' + denom
                xdata = gal.df[('io',field)]
            
    
        hist,bins = np.histogram(xdata[select], bins = bins)
        #dx   = 10.0**(bins[1:]) - 10.0**(bins[:-1])
        norm = 1.0*np.sum(hist)
        plot_histogram(ax[index], bins, hist / (norm), lw = 3, color = 'black', ls = '-')

        #hist, bins = np.histogram(xdata, bins = bins)
        #dx   = 10.0**(bins[1:]) - 10.0**(bins[:-1])
        #norm = 1.0*np.sum(hist)*dx
        #plot_histogram(ax[index], bins, hist / (norm), lw = 3, color = 'red', ls = '--')

        ax[index].set_ylim(0,0.2)
        #ax[index].semilogy()
        
        if denom is None:
            ax[index].set_xlabel(e + ' Mass Fraction')
            ax[index].set_xlim(-20,0)
        else:
            ax[index].set_xlabel('[X/' + denom +']')
            if denom == 'H':
                ax[index].set_xlim(-8,0)
            else:
                ax[index].set_xlim(-4,4)
            
        #ax[index].set_ylabel('Fraction of Stars')
    
        axj = axj + 1
        if axj >= ncol:
            axj = 0
            axi = axi + 1
            

plt.tight_layout()
fig.savefig('MDFs.png')


    

    




/home/aemerick/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:35: RuntimeWarning: invalid value encountered in true_divide






    

In [6]:

    

t = 500.0 * yt.units.Myr
select = is_alive(gal.df, t)

elements = ['N','O','Sr','Mg','Mn','Fe','Ba']

nrow = 2
ncol = 4
fig,ax = plt.subplots(nrow, ncol)
fig.set_size_inches(ncol*4,nrow*4)

axi = 0
axj = 0

denom_element = 'Mg'
xdata = gal.df['creation_time'].convert_to_units('Myr')
xdata = xdata - np.min(xdata)

bins = np.arange(-30,0, 0.1)
for e in elements:
    index = (axi,axj)
    if e == denom_element:
        continue
        
    ydata = gal.df[('io','particle_' + e + '_over_' + denom_element)]
    
    ax[index].scatter(xdata[select], ydata[select])
    

    #ax[index].set_xlim(-20,0)
    #ax[index].set_ylim(0,1)
    #ax[index].semilogy()
    ax[index].set_ylabel(e + ' over ' + denom_element)
    ax[index].set_xlabel('Formation Time')
    ax[index].set_xlim(0,500)

    
    axj = axj + 1
    if axj >= ncol:
        axj = 0
        axi = axi + 1
        
plt.tight_layout()


    

In [8]:

    

t = 500.0 * yt.units.Myr
select = is_alive(gal.df, t)

elements = ['N','O','Sr','Mg','Mn','Fe','Ba']

nrow = 2
ncol = 4
fig,ax = plt.subplots(nrow, ncol)
fig.set_size_inches(ncol*4,nrow*4)

axi = 0
axj = 0

denom_element = 'Mg'
xdata = gal.df['creation_time'].convert_to_units('Myr')
xdata = xdata - np.min(xdata)

bins = np.arange(-30,0, 0.1)
for e in elements:
    index = (axi,axj)
    if e == denom_element:
        continue
        
    ydata = gal.df[('io','particle_' + e + '_over_' + denom_element)]
    
    ydata = convert_abundances.renormalize(ydata, e, denom_element)
    
    ax[index].scatter(xdata[select], ydata[select])

    #ax[index].set_xlim(-20,0)
    #ax[index].set_ylim(0,1)
    #ax[index].semilogy()
    ax[index].set_ylabel(e + ' over ' + denom_element)
    ax[index].set_xlabel('Formation Time')
    ax[index].set_xlim(0,500)

    ax[index].set_ylim(np.average(ydata[select]) - 2, np.average(ydata[select]) + 2)

    axj = axj + 1
    if axj >= ncol:
        axj = 0
        axi = axi + 1
        
plt.tight_layout()


    

In [32]:

    

t = 1250.0 * yt.units.Myr
select = is_alive(gal.df, t)

fsize = 17
rc('text',usetex=False)
rc('font',size=fsize)

elements = ['N','O','Mg','Ca','Mn','Fe','Ba']

nrow = 4
ncol = 7
fig,ax = plt.subplots(nrow, ncol, sharex=True)
fig.set_size_inches(ncol*4,nrow*4)


   
axi = 0
axj = 0
    


bins = np.arange(-30,0, 0.1)
for denom_element in ['O','Fe','N','H']:
    xdata = gal.df['creation_time'].convert_to_units('Myr')
    xdata = xdata - np.min(xdata) 
 
    axj = 0
    for e in elements:
        index = (axi,axj)
        
        if e == denom_element:
            ydata = np.ones(np.size(xdata))*-100
        else:
            ydata = gal.df[('io','particle_' + e + '_over_' + denom_element)]
        
     
        ax[index].scatter(xdata[select], ydata[select], color = 'black', s = 20, alpha = 0.5)
    

    #ax[index].set_xlim(-20,0)
        if denom_element == 'H':
            ax[index].set_ylim(-6, 0)
            xy = (400, -0.5)
            ax[index].annotate(e, xy, xy)
        else:
            ax[index].set_ylim(-3,3)
    #ax[index].semilogy()
            xy = (400,2.5)
            ax[index].annotate(e, xy, xy)
        ax[index].set_xlim(0,1200)

        if axj == 0:
            ax[index].set_ylabel('[X/' + denom_element + ']')
        if axi == nrow - 1:
            ax[index].set_xlabel('Formation Time')
    
        axj = axj + 1
        if axj >= ncol:
            axj = 0
            axi = axi + 1
        
fig.subplots_adjust(hspace = 0, wspace = 0)

plt.setp([a.get_xticklabels() for a in fig.axes[:-ncol]], visible = False)

axi = 0
axj = 0
for a in fig.axes:

    if axj != 0:
        plt.setp(a.get_yticklabels(), visible = False)
    
    axj = axj + 1
    if axj >= ncol:
        axj = 0
        axi = axi + 1
#plt.setp([a.get_yticklabels() for a in fig.axes[:-ncol]], visible = False)

plt.minorticks_on()

#plt.tight_layout()
fig.savefig('stellar_abundances_time.png')


    

In [20]:

    

t = 1200.0 * yt.units.Myr
select = is_alive(gal.df, t)

elements = ['N','O','Sr','Mg','Mn','Fe','Ba']

nrow = 4
ncol = 7
fig,ax = plt.subplots(nrow, ncol, sharex=True)
fig.set_size_inches(ncol*4,nrow*4)

cdata = gal.df['creation_time'].convert_to_units('Myr')
cdata = cdata - np.min(cdata)
   
axi = 0
axj = 0
    
x1 = 'O'
x2 = 'Mg'

bins = np.arange(-30,0, 0.1)
for denom_element in ['O','Mg','N','H']:
#    xdata = gal.df['creation_time'].convert_to_units('Myr')
#    xdata = xdata - np.min(xdata) 
    xdata = gal.df[('io','particle_'+ x1 + '_over_' + x2)]

    
    axj = 0
    for e in elements:
        index = (axi,axj)
        
        if e == denom_element:
            ydata = np.ones(np.size(xdata))*-100
        else:
            ydata = gal.df[('io','particle_' + e + '_over_' + denom_element)]
        
     
        ax[index].scatter(xdata[select], ydata[select], s = 20, c=cdata[select])# color = 'black')
    

    #ax[index].set_xlim(-20,0)
        if denom_element == 'H':
            ax[index].set_ylim(-6, 0)
            xy = (-0.5, -0.5)
            ax[index].annotate(e, xy, xy)
        else:
            ax[index].set_ylim(-3,3)
    #ax[index].semilogy()
            xy = (-0.5,2.5)
            ax[index].annotate(e, xy, xy)
            
        if x1 == 'O' and x2 == 'H':
            ax[index].set_xlim(-6.5, 0.5)
        if x1 == 'Mg' and x2 == 'H':
            ax[index].set_xlim(-6.5, 0.5)

        if axj == 0:
            ax[index].set_ylabel('[X/' + denom_element + ']')
        if axi == nrow - 1:
            ax[index].set_xlabel('[' + x1 + '/' + x2 + ']')
    
        axj = axj + 1
        if axj >= ncol:
            axj = 0
            axi = axi + 1
        
fig.subplots_adjust(hspace = 0, wspace = 0)

plt.setp([a.get_xticklabels() for a in fig.axes[:-ncol]], visible = False)

axi = 0
axj = 0
for a in fig.axes:

    if axj != 0:
        plt.setp(a.get_yticklabels(), visible = False)
    
    axj = axj + 1
    if axj >= ncol:
        axj = 0
        axi = axi + 1
#plt.setp([a.get_yticklabels() for a in fig.axes[:-ncol]], visible = False)

plt.minorticks_on()

#plt.tight_layout()
fig.savefig('stellar_abundances_' + x1 + '_over_' + x2 + '.png')


    

In [ ]:

    

t = 1200.0 * yt.units.Myr
select = is_alive(gal.df, t)

elements = ['N','O','Ca','Mg','Mn','Fe','Ba']

nrow = 2
ncol = 4
fig,ax = plt.subplots(nrow, ncol)
fig.set_size_inches(ncol*5,nrow*5)

axi = 0
axj = 0
cdata = gal.df['creation_time'].convert_to_units('Myr')
cdata = cdata - np.min(cdata)

denom_element = 'Fe'
xdata = gal.df['creation_time'].convert_to_units('Myr')
#xdata = xdata - np.min(xdata)

e1 = 'Fe'
e2 = 'H'
xdata = gal.df[('io','particle_' + e1 + '_over_' + e2)]
#xdata = convert_abundances.renormalize(xdata, e1,e2)

bins = np.arange(-30,0, 0.1)
for e in elements:
    index = (axi,axj)
    if e == denom_element:
        continue
        
    ydata = gal.df[('io','particle_' + e + '_over_' + denom_element)]
    #ydata = convert_abundances.renormalize(ydata, e, denom_element)
    
    ax[index].scatter(xdata[select], ydata[select], c = cdata[select])
    

    #ax[index].set_xlim(-20,0)
    #ax[index].set_ylim(0,1)
    #ax[index].semilogy()
    ax[index].set_ylabel(e + ' over ' + denom_element)
    ax[index].set_xlabel(e1 + ' over ' + e2)
    ax[index].set_xlim(-8, -1)

    ax[index].set_ylim(np.average(ydata[select]) - 2, np.average(ydata[select]) + 2)
    
    axj = axj + 1
    if axj >= ncol:
        axj = 0
        axi = axi + 1
    plt.minorticks_on()
    plt.tight_layout()


    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [33]:

    

axj = 0
        axi = axi + 1
plt.minorticks_on()
plt.tight_layout()


    

In [ ]:

    

t = 1200.0 * yt.units.Myr
select = is_alive(gal.df, t)

elements = ['N','O','Sr','Mg','Mn','Fe','Ba']

nrow = 2
ncol = 4
fig,ax = plt.subplots(nrow, ncol)
fig.set_size_inches(ncol*5,nrow*5)


    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    

axi = 0
axj = 0

denom_element = 'Mg'
cdata = gal.df['creation_time'].convert_to_units('Myr')


    

In [ ]:

    




    

In [ ]:

    




    

In [26]:

    

cdata = cdata - np.min(cdata)

xdata = gal.df[('io','particle_' + denom_element + '_over_H')]

bins = np.arange(-30,0, 0.1)
for e in elements:
    index = (axi,axj)
    if e == denom_element:
        continue
        
    ydata = gal.df[('io','particle_' + e + '_over_' + denom_element)]
    
  
    ax[index].scatter(xdata[select], ydata[select], c = cdata[select])
    

    #ax[index].set_xlim(-20,0)
    #ax[index].set_ylim(0,1)
    #ax[index].semilogy()
    ax[index].set_ylabel(e + ' over ' + denom_element)
    ax[index].set_xlabel(denom_element + ' over H')
    ax[index].set_xlim(-5,0)

    
    axj = axj + 1
    if axj >= ncol:
        axj = 0
        axi = axi + 1
plt.minorticks_on()
        
plt.tight_layout()


    

In [ ]: