In [3]:
%matplotlib inline
import yt
import numpy as np
import matplotlib.pyplot as plt
from yt.visualization.base_plot_types import get_multi_plot
import matplotlib.colorbar as cb
from matplotlib.colors import LogNorm
from galaxy_analysis.yt_fields import field_generators as fg
from galaxy_analysis.utilities import convert_abundances as ca
from galaxy_analysis.utilities import utilities as galutil
from galaxy_analysis.plot.plot_styles import *
from matplotlib.ticker import (MultipleLocator, FormatStrFormatter,
AutoMinorLocator)
def plot_settings():
fsize = 21
rc('text',usetex=False)
rc('font',size=fsize)
return
plot_settings()
rc('text',usetex=False)
rc('font',size=fsize)
In [4]:
workdir = "/home/aemerick/work/enzo_runs/cosmo_testing/feedback/grid/"
In [11]:
dsname = 'DD0045'
dsname = dsname + '/' + dsname
ds = yt.load(workdir + dsname)
fg.generate_derived_fields(ds)
ds = yt.load(workdir + dsname)
fg.generate_particle_filters(ds)
data = ds.all_data()
In [12]:
field_units = {'Density' : 'g/cm**3',
'x_velocity' : 'km/s',
'y_velocity' : 'km/s','z_velocity':'km/s','vx' : 'km/s',
"Temperature" : 'K'}
print("%20s - %30s: %5s %5s %5s"%("Type 0","Field","Sum","Min","Max"))
for field in ds.field_list:
vals = data[field]
if np.size(vals) > 0:
if field[1] in field_units:
vals = vals.to(field_units[field[1]])
elif 'kph' in field[1]:
vals = vals.to('1/s')
elif '_Density' in field[1]:
try:
vals = vals.to('g/cm**3')
except:
vals = vals
print("%20s - %30s: %5.3E %5.3E %5.3E"%(field[0],field[1], np.sum(vals), np.min(vals), np.max(vals)))
else:
print("%20s - %30s: None None None"%(field[0],field[1]))
In [13]:
pt = data['particle_type']
print(np.unique(pt))
In [14]:
#
# Check to make sure no particles that are SNIa
# get double assigned SNIa types
#
num_doubled = 0
nsnia = 0
for i in np.arange(np.size(pt)):
num_neg = 0
if (data['snia_metal_fraction'][i] < 0):
num_neg += 1
if (data['snia_sch_metal_fraction'][i] < 0):
num_neg += 1
if (data['snia_sds_metal_fraction'][i] < 0):
num_neg += 1
if (data['snia_hers_metal_fraction'][i] < 0):
num_neg += 1
if num_neg == 1:
nsnia += 1
elif num_neg > 1:
num_doubled += 1
print("doubled up on particle",i,num_neg)
print(num_doubled)
In [15]:
#
# Check out the WD particles
#
def check_SNIa_properties(data):
pt = data['particle_type']
m = data['particle_mass'].to('Msun')
WD_select = pt == 12
complete_SNIa = (pt==12)*(m==0.0)
n_part = np.size(pt)
ideal = {'DDS' : 0.05126958193824957,'SDS': 0.005242912915506227, 'sCh' :0.09073471211157359,
'HeRS' : 0.0035997607307840106, 'total' : 0.15084696769611342}
print("Number of particles ", n_part)
print("Number of WD particles ", np.size(pt[WD_select]))
print("Number of particles that will (or have) gone SNIa ", nsnia)
print("Number of particles that have gone SNIa ", np.size(pt[complete_SNIa]))
sntypes = ['DDS','sCh','SDS','HeRS']
#
# Count the types
#
#
num_snia = {}
data_field = {'DDS' : 'snia_metal_fraction', 'sCh' : 'snia_sch_metal_fraction',
'SDS' : 'snia_sds_metal_fraction', 'HeRS' : 'snia_hers_metal_fraction'}
for sntype in sntypes:
select = data[ data_field[sntype] ] < 0
num_snia[sntype] = np.size(pt[select])
print("Fraction of SNIa Types")
total_num = 0
for sntype in sntypes:
total_num += num_snia[sntype]
print("%10s: Number = %4i - Fraction = %.5f - Target Fraction %.5f"%(sntype, num_snia[sntype], num_snia[sntype]/(1.0*n_part), ideal[sntype]))
print("%10s: Number = %4i - Fraction = %.5f - Target Fraction %.5f"%("Total", total_num, total_num/(1.0*n_part), ideal['total']))
check_SNIa_properties(data)
In [16]:
#
# Check lifetimes
#
lt = data['dynamical_time'].to('Myr')
select = lt < 15.0E3 # Myr
print(np.size(lt[select]))
print(np.min(lt[select]), np.max(lt[select]))
print(np.median(lt[select]))
print(np.average(lt[select]))
In [ ]:
In [ ]: