[yt-svn] commit/yt: 6 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Fri Mar 28 12:09:47 PDT 2014
6 new commits in yt:
https://bitbucket.org/yt_analysis/yt/commits/7a72b9113c4b/
Changeset: 7a72b9113c4b
Branch: yt-3.0
User: galtay
Date: 2014-03-26 00:27:19
Summary: Added a few things here.
1) OWLS ions are now available as particle fields
2) OWLS number density of ions are available as smoothed fields.
3) added routine to translate new ion naming convention (_p*_)
to roman numeral notation in slice and projection plots
Affected #: 6 files
diff -r 0536fb6516915620927dcffb4281021a46b73214 -r 7a72b9113c4bb3dda4a90d45886f407d7da9acdd yt/fields/species_fields.py
--- a/yt/fields/species_fields.py
+++ b/yt/fields/species_fields.py
@@ -59,34 +59,46 @@
* data[ftype,'density']
return _density
-def add_species_field_by_density(registry, ftype, species):
+def add_species_field_by_density(registry, ftype, species,
+ particle_type = False):
"""
This takes a field registry, a fluid type, and a species name and then
adds the other fluids based on that. This assumes that the field
"SPECIES_density" already exists and refers to mass density.
"""
registry.add_field((ftype, "%s_fraction" % species),
- function = _create_fraction_func(ftype, species),
- units = "")
+ function = _create_fraction_func(ftype, species),
+ particle_type = particle_type,
+ units = "")
+
registry.add_field((ftype, "%s_mass" % species),
- function = _create_mass_func(ftype, species),
- units = "g")
+ function = _create_mass_func(ftype, species),
+ particle_type = particle_type,
+ units = "g")
+
registry.add_field((ftype, "%s_number_density" % species),
- function = _create_number_density_func(ftype, species),
- units = "cm**-3")
+ function = _create_number_density_func(ftype, species),
+ particle_type = particle_type,
+ units = "cm**-3")
-def add_species_field_by_fraction(registry, ftype, species):
+def add_species_field_by_fraction(registry, ftype, species,
+ particle_type = False):
"""
This takes a field registry, a fluid type, and a species name and then
adds the other fluids based on that. This assumes that the field
"SPECIES_fraction" already exists and refers to mass fraction.
"""
registry.add_field((ftype, "%s_density" % species),
- function = _create_density_func(ftype, species),
- units = "g/cm**3")
+ function = _create_density_func(ftype, species),
+ particle_type = particle_type,
+ units = "g/cm**3")
+
registry.add_field((ftype, "%s_mass" % species),
- function = _create_mass_func(ftype, species),
- units = "g")
+ function = _create_mass_func(ftype, species),
+ particle_type = particle_type,
+ units = "g")
+
registry.add_field((ftype, "%s_number_density" % species),
- function = _create_number_density_func(ftype, species),
- units = "cm**-3")
+ function = _create_number_density_func(ftype, species),
+ particle_type = particle_type,
+ units = "cm**-3")
diff -r 0536fb6516915620927dcffb4281021a46b73214 -r 7a72b9113c4bb3dda4a90d45886f407d7da9acdd yt/frontends/sph/data_structures.py
--- a/yt/frontends/sph/data_structures.py
+++ b/yt/frontends/sph/data_structures.py
@@ -209,7 +209,10 @@
if "length" in unit_base:
length_unit = unit_base["length"]
elif "UnitLength_in_cm" in unit_base:
- length_unit = (unit_base["UnitLength_in_cm"], "cm")
+ if self.cosmological_simulation == 0:
+ length_unit = (unit_base["UnitLength_in_cm"], "cm")
+ else:
+ length_unit = (unit_base["UnitLength_in_cm"], "cmcm/h")
else:
raise RuntimeError
length_unit = _fix_unit_ordering(length_unit)
@@ -229,7 +232,10 @@
if "mass" in unit_base:
mass_unit = unit_base["mass"]
elif "UnitMass_in_g" in unit_base:
- mass_unit = (unit_base["UnitMass_in_g"], "g")
+ if self.cosmological_simulation == 0:
+ mass_unit = (unit_base["UnitMass_in_g"], "g")
+ else:
+ mass_unit = (unit_base["UnitMass_in_g"], "g/h")
else:
# Sane default
mass_unit = (1.0, "1e10*Msun/h")
@@ -286,6 +292,8 @@
_particle_mass_name = "Mass"
_field_info_class = OWLSFieldInfo
+
+
def _parse_parameter_file(self):
handle = h5py.File(self.parameter_filename, mode="r")
hvals = {}
diff -r 0536fb6516915620927dcffb4281021a46b73214 -r 7a72b9113c4bb3dda4a90d45886f407d7da9acdd yt/frontends/sph/fields.py
--- a/yt/frontends/sph/fields.py
+++ b/yt/frontends/sph/fields.py
@@ -14,17 +14,27 @@
# The full license is in the file COPYING.txt, distributed with this software.
#-----------------------------------------------------------------------------
+import os
import numpy as np
+import owls_ion_tables as oit
from yt.funcs import *
+
from yt.fields.field_info_container import \
FieldInfoContainer
+
from .definitions import \
gadget_ptypes, \
ghdf5_ptypes
-from yt.fields.species_fields import add_species_field_by_fraction
+from yt.config import ytcfg
+from yt.utilities.physical_constants import mh
+from yt.fields.species_fields import \
+ add_species_field_by_fraction, \
+ add_species_field_by_density
+from yt.fields.particle_fields import \
+ add_volume_weighted_smoothed_field
# Here are helper functions for things like vector fields and so on.
@@ -60,24 +70,27 @@
class OWLSFieldInfo(SPHFieldInfo):
- _species_fractions = ['H_fraction', 'He_fraction', 'C_fraction',
- 'N_fraction', 'O_fraction', 'Ne_fraction',
- 'Mg_fraction', 'Si_fraction', 'Fe_fraction']
+ _ions = ("c1", "c2", "c3", "c4", "c5", "c6",
+ "fe2", "fe17", "h1", "he1", "he2", "mg1", "mg2", "n2",
+ "n3", "n4", "n5", "n6", "n7", "ne8", "ne9", "ne10", "o1",
+ "o6", "o7", "o8", "si2", "si3", "si4", "si13")
+
+ _elements = ("H", "He", "C", "N", "O", "Ne", "Mg", "Si", "Fe")
# override
#--------------------------------------------------------------
def __init__(self, *args, **kwargs):
new_particle_fields = (
- ('Hydrogen', ('', ['H_fraction'], None)),
- ('Helium', ('', ['He_fraction'], None)),
- ('Carbon', ('', ['C_fraction'], None)),
- ('Nitrogen', ('', ['N_fraction'], None)),
- ('Oxygen', ('', ['O_fraction'], None)),
- ('Neon', ('', ['Ne_fraction'], None)),
- ('Magnesium', ('', ['Mg_fraction'], None)),
- ('Silicon', ('', ['Si_fraction'], None)),
- ('Iron', ('', ['Fe_fraction'], None))
+ ("Hydrogen", ("", ["H_fraction"], None)),
+ ("Helium", ("", ["He_fraction"], None)),
+ ("Carbon", ("", ["C_fraction"], None)),
+ ("Nitrogen", ("", ["N_fraction"], None)),
+ ("Oxygen", ("", ["O_fraction"], None)),
+ ("Neon", ("", ["Ne_fraction"], None)),
+ ("Magnesium", ("", ["Mg_fraction"], None)),
+ ("Silicon", ("", ["Si_fraction"], None)),
+ ("Iron", ("", ["Fe_fraction"], None))
)
self.known_particle_fields += new_particle_fields
@@ -85,9 +98,257 @@
super(OWLSFieldInfo,self).__init__( *args, **kwargs )
+
+
+
+
+ def setup_particle_fields(self, ptype):
+ """ additional particle fields that are not on disk """
+
+ super(OWLSFieldInfo,self).setup_particle_fields(ptype)
+
+ if ptype == "PartType0":
+
+ # this adds the particle element fields
+ # X_density, X_mass, and X_number_density
+ # where X is an element of self._elements.
+ # X_fraction are defined in file
+ #-----------------------------------------------
+ for s in self._elements:
+ add_species_field_by_fraction(self, ptype, s,
+ particle_type=True)
+
+ # this defines the ion density on particles
+ # X_density
+ #-----------------------------------------------
+ self.setup_ion_density_particle_fields()
+
+
+ # this adds the rest of the ion particle fields
+ # X_fraction, X_mass, X_number_density
+ #-----------------------------------------------
+ for ion in self._ions:
+
+ # construct yt name for ion
+ #---------------------------------------------------
+ if ion[0:2].isalpha():
+ symbol = ion[0:2].capitalize()
+ roman = int(ion[2:])
+ else:
+ symbol = ion[0:1].capitalize()
+ roman = int(ion[1:])
+
+ pstr = "_p" + str(roman-1)
+ yt_ion = symbol + pstr
+
+ add_species_field_by_density(self, ptype, yt_ion,
+ particle_type=True)
+
+
+
+ ptype = 'PartType0'
+ num_neighbors = 48
+ loaded = []
+
+ # add smoothed element fields
+ #-----------------------------------------------
+ for s in self._elements:
+ fname = s + '_number_density'
+ fn = add_volume_weighted_smoothed_field(
+ ptype, "particle_position", "particle_mass",
+ "smoothing_length", "density", fname, self,
+ num_neighbors)
+ loaded += fn
+
+ self.alias(("gas", fname), fn[0])
+
+ self._show_field_errors += loaded
+
+
+ # add smoothed ion fields
+ #-----------------------------------------------
+ for ion in self._ions:
+
+ if ion[0:2].isalpha():
+ symbol = ion[0:2].capitalize()
+ roman = int(ion[2:])
+ else:
+ symbol = ion[0:1].capitalize()
+ roman = int(ion[1:])
+
+ pstr = "_p" + str(roman-1)
+ yt_ion = symbol + pstr
+
+ fname = yt_ion + '_number_density'
+ fn = add_volume_weighted_smoothed_field(
+ ptype, "particle_position", "particle_mass",
+ "smoothing_length", "density", fname, self,
+ num_neighbors)
+ loaded += fn
+
+ self.alias(("gas", fname), fn[0])
+
+ self._show_field_errors += loaded
+
+
+ # find dependencies
+ #-----------------------------------------------
+ self.find_dependencies(loaded)
+
+
+ return
+
+
+
+
+
+
+ def setup_ion_density_particle_fields( self ):
+ """ Sets up particle fields for ion densities. """
+
+ # loop over all ions and make fields
+ #----------------------------------------------
+# for ion in ("h1",): # self._ions:
+ for ion in self._ions:
+
+ # construct yt name for ion
+ #---------------------------------------------------
+ if ion[0:2].isalpha():
+ symbol = ion[0:2].capitalize()
+ roman = int(ion[2:])
+ else:
+ symbol = ion[0:1].capitalize()
+ roman = int(ion[1:])
+
+ pstr = "_p" + str(roman-1)
+ yt_ion = symbol + pstr
+ ftype = "PartType0"
+
+ # add ion density field for particles
+ #---------------------------------------------------
+ fname = yt_ion + '_density'
+ dens_func = self._create_ion_density_func( ftype, ion )
+ self.add_field( (ftype, fname),
+ function = dens_func,
+ units="g/cm**3",
+ particle_type=True )
+ self._show_field_errors.append( (ftype,fname) )
+
+
+
+ def _create_ion_density_func( self, ftype, ion ):
+ """ returns a function that calculates the ion density of a particle.
+ """
+
+ def _ion_density(field, data):
+
+ # get element symbol from ion string. ion string will
+ # be a member of the tuple _ions (i.e. si13)
+ #--------------------------------------------------------
+ if ion[0:2].isalpha():
+ symbol = ion[0:2].capitalize()
+ else:
+ symbol = ion[0:1].capitalize()
+
+ # mass fraction for the element
+ #--------------------------------------------------------
+ m_frac = data[ftype, symbol+"_fraction"]
+
+ # get nH and T for lookup
+ #--------------------------------------------------------
+ log_nH = np.log10( data["PartType0", "H_number_density"] )
+ log_T = np.log10( data["PartType0", "Temperature"] )
+
+ # get name of owls_ion_file for given ion
+ #--------------------------------------------------------
+ owls_ion_path = self._get_owls_ion_data_dir()
+ fname = os.path.join( owls_ion_path, ion+".hdf5" )
+
+ # create ionization table for this redshift
+ #--------------------------------------------------------
+ itab = oit.IonTableOWLS( fname )
+ itab.set_iz( data.pf.current_redshift )
+
+ # find ion balance using log nH and log T
+ #--------------------------------------------------------
+ i_frac = itab.interp( log_nH, log_T )
+ return data[ftype,"Density"] * m_frac * i_frac
+
+ return _ion_density
+
+
+
+
+
+ # this function sets up the X_mass, X_density, X_fraction, and
+ # X_number_density fields where X is the name of an OWLS element.
+ #-------------------------------------------------------------
def setup_fluid_fields(self):
- # here species_name is "H", "He", etc
- for s in self._species_fractions:
- species_name = s.split('_')[0]
- add_species_field_by_fraction(self, "gas", species_name)
+
+ return
+
+
+ # these add the smoothed element fields
+ #-----------------------------------------------
+ for s in self._elements:
+ add_species_field_by_fraction(self, "gas", s)
+
+
+
+ # this should add the smoothed ion fields
+ #-----------------------------------------------
+ for ion in self._ions:
+
+ if ion[0:2].isalpha():
+ symbol = ion[0:2].capitalize()
+ roman = int(ion[2:])
+ else:
+ symbol = ion[0:1].capitalize()
+ roman = int(ion[1:])
+
+ pstr = "_p" + str(roman-1)
+ yt_ion = symbol + pstr
+
+ add_species_field_by_density(self, "gas", yt_ion)
+
+
+ # this function returns the owls_ion_data directory. if it doesn't
+ # exist it will download the data from http://yt-project.org/data
+ #-------------------------------------------------------------
+ def _get_owls_ion_data_dir(self):
+
+ txt = "Attempting to download ~ 30 Mb of owls ion data from %s to %s."
+ data_file = "owls_ion_data.tar.gz"
+ data_url = "http://yt-project.org/data"
+
+ # get test_data_dir from yt config (ytcgf)
+ #----------------------------------------------
+ tdir = ytcfg.get("yt","test_data_dir")
+
+ # set download destination to tdir or ./ if tdir isnt defined
+ #----------------------------------------------
+ if tdir == "/does/not/exist":
+ data_dir = "./"
+ else:
+ data_dir = tdir
+
+
+ # check for owls_ion_data directory in data_dir
+ # if not there download the tarball and untar it
+ #----------------------------------------------
+ owls_ion_path = os.path.join( data_dir, "owls_ion_data" )
+
+ if not os.path.exists(owls_ion_path):
+ mylog.info(txt % (data_url, data_dir))
+ fname = data_dir + "/" + data_file
+ fn = download_file(os.path.join(data_url, data_file), fname)
+
+ cmnd = "cd " + data_dir + "; " + "tar xf " + data_file
+ os.system(cmnd)
+
+
+ if not os.path.exists(owls_ion_path):
+ raise RuntimeError, "Failed to download owls ion data."
+
+ return owls_ion_path
diff -r 0536fb6516915620927dcffb4281021a46b73214 -r 7a72b9113c4bb3dda4a90d45886f407d7da9acdd yt/frontends/sph/io.py
--- a/yt/frontends/sph/io.py
+++ b/yt/frontends/sph/io.py
@@ -53,7 +53,7 @@
_vector_fields = ("Coordinates", "Velocity", "Velocities")
_known_ptypes = ghdf5_ptypes
_var_mass = None
- _element_fields = ('Hydrogen', 'Helium', 'Carbon', 'Nitrogen', 'Oxygen',
+ _element_names = ('Hydrogen', 'Helium', 'Carbon', 'Nitrogen', 'Oxygen',
'Neon', 'Magnesium', 'Silicon', 'Iron' )
@@ -110,7 +110,7 @@
ind = self._known_ptypes.index(ptype)
data[:] = self.pf["Massarr"][ind]
- elif field in self._element_fields:
+ elif field in self._element_names:
rfield = 'ElementAbundance/' + field
data = g[rfield][:][mask,...]
diff -r 0536fb6516915620927dcffb4281021a46b73214 -r 7a72b9113c4bb3dda4a90d45886f407d7da9acdd yt/frontends/sph/owls_ion_tables.py
--- /dev/null
+++ b/yt/frontends/sph/owls_ion_tables.py
@@ -0,0 +1,224 @@
+"""
+OWLS ion tables
+
+A module to handle the HM01 UV background spectra and ionization data from the
+OWLS photoionization equilibrium lookup tables.
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2013, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+import sys
+import h5py
+import numpy as np
+
+
+
+
+def h5rd( fname, path, dtype=None ):
+ """ Read Data. Return a dataset located at <path> in file <fname> as
+ a numpy array.
+ e.g. rd( fname, '/PartType0/Coordinates' ). """
+
+ data = None
+ with h5py.File( fname, 'r' ) as h5f:
+ ds = h5f[path]
+ if dtype == None:
+ dtype = ds.dtype
+ data = np.zeros( ds.shape, dtype=dtype )
+ data = ds.value
+ return data
+
+
+
+class IonTableSpectrum:
+
+ """ A class to handle the HM01 spectra in the OWLS ionization tables. """
+
+ def __init__(self, ion_file):
+
+ where = '/header/spectrum/gammahi'
+ self.GH1 = h5rd( ion_file, where ) # GH1[1/s]
+
+ where = '/header/spectrum/logenergy_ryd'
+ self.logryd = h5rd( ion_file, where ) # E[ryd]
+
+ where = '/header/spectrum/logflux'
+ self.logflux = h5rd( ion_file, where ) # J[ergs/s/Hz/Sr/cm^2]
+
+ where = '/header/spectrum/redshift'
+ self.z = h5rd( ion_file, where ) # z
+
+
+
+ def return_table_GH1_at_z(self,z):
+
+ # find redshift indices
+ #-----------------------------------------------------------------
+ i_zlo = np.argmin( np.abs( self.z - z ) )
+ if self.z[i_zlo] < z:
+ i_zhi = i_zlo + 1
+ else:
+ i_zhi = i_zlo
+ i_zlo = i_zlo - 1
+
+ z_frac = (z - self.z[i_zlo]) / (self.z[i_zhi] - self.z[i_zlo])
+
+ # find GH1 from table
+ #-----------------------------------------------------------------
+ logGH1_all = np.log10( self.GH1 )
+ dlog_GH1 = logGH1_all[i_zhi] - logGH1_all[i_zlo]
+
+ logGH1_table = logGH1_all[i_zlo] + z_frac * dlog_GH1
+ GH1_table = 10.0**logGH1_table
+
+ return GH1_table
+
+
+
+
+class IonTableOWLS:
+
+ """ A class to handle OWLS ionization tables. """
+
+ DELTA_nH = 0.25
+ DELTA_T = 0.1
+
+ def __init__(self, ion_file):
+
+ self.ion_file = ion_file
+
+ # ionbal is indexed like [nH, T, z]
+ # nH and T are log quantities
+ #---------------------------------------------------------------
+ self.nH = h5rd( ion_file, '/logd' ) # log nH [cm^-3]
+ self.T = h5rd( ion_file, '/logt' ) # log T [K]
+ self.z = h5rd( ion_file, '/redshift' ) # z
+
+ # read the ionization fractions
+ # linear values stored in file so take log here
+ # ionbal is the ionization balance (i.e. fraction)
+ #---------------------------------------------------------------
+ self.ionbal = h5rd( ion_file, '/ionbal' ).astype(np.float64)
+ self.ionbal_orig = self.ionbal.copy()
+
+ ipositive = np.where( self.ionbal > 0.0 )
+ izero = np.where( self.ionbal <= 0.0 )
+ self.ionbal[izero] = self.ionbal[ipositive].min()
+
+ self.ionbal = np.log10( self.ionbal )
+
+
+ # load in background spectrum
+ #---------------------------------------------------------------
+ self.spectrum = IonTableSpectrum( ion_file )
+
+ # calculate the spacing along each dimension
+ #---------------------------------------------------------------
+ self.dnH = self.nH[1:] - self.nH[0:-1]
+ self.dT = self.T[1:] - self.T[0:-1]
+ self.dz = self.z[1:] - self.z[0:-1]
+
+ self.order_str = '[log nH, log T, z]'
+
+
+
+
+
+ # sets iz and fz
+ #-----------------------------------------------------
+ def set_iz( self, z ):
+
+ if z <= self.z[0]:
+ self.iz = 0
+ self.fz = 0.0
+ elif z >= self.z[-1]:
+ self.iz = len(self.z) - 2
+ self.fz = 1.0
+ else:
+ for iz in range( len(self.z)-1 ):
+ if z < self.z[iz+1]:
+ self.iz = iz
+ self.fz = ( z - self.z[iz] ) / self.dz[iz]
+ break
+
+
+
+ # interpolate the table at a fixed redshift for the input
+ # values of nH and T ( input should be log ). A simple
+ # tri-linear interpolation is used.
+ #-----------------------------------------------------
+ def interp( self, nH, T ):
+
+ nH = np.array( nH )
+ T = np.array( T )
+
+ if nH.size != T.size:
+ print ' array size mismatch !!! '
+ sys.exit(1)
+
+ # field discovery will have nH.size == 1 and T.size == 1
+ # in that case we simply return 1.0
+
+ if nH.size == 1 and T.size == 1:
+ ionfrac = 1.0
+ return ionfrac
+
+
+ # find inH and fnH
+ #-----------------------------------------------------
+ inH = np.int32( ( nH - self.nH[0] ) / self.DELTA_nH )
+ fnH = ( nH - self.nH[inH] ) / self.dnH[inH]
+
+ indx = np.where( inH < 0 )[0]
+ if len(indx) > 0:
+ inH[indx] = 0
+ fnH[indx] = 0.0
+
+ indx = np.where( inH >= len(nH) )[0]
+ if len(indx) > 0:
+ inH[indx] = len(nH)-2
+ fnH[indx] = 1.0
+
+
+ # find iT and fT
+ #-----------------------------------------------------
+ iT = np.int32( ( T - self.T[0] ) / self.DELTA_T )
+ fT = ( T - self.T[iT] ) / self.dT[iT]
+
+ indx = np.where( iT < 0 )[0]
+ if len(indx) > 0:
+ iT[indx] = 0
+ fT[indx] = 0.0
+
+ indx = np.where( iT >= len(T) )[0]
+ if len(indx) > 0:
+ iT[indx] = len(T)-2
+ fT[indx] = 1.0
+
+
+ iz = self.iz
+ fz = self.fz
+
+ # calculate interpolated value
+ # use tri-linear interpolation on the log values
+ #-----------------------------------------------------
+
+ ionfrac = self.ionbal[inH, iT, iz ] * (1-fnH) * (1-fT) * (1-fz) + \
+ self.ionbal[inH+1, iT, iz ] * (fnH) * (1-fT) * (1-fz) + \
+ self.ionbal[inH, iT+1, iz ] * (1-fnH) * (fT) * (1-fz) + \
+ self.ionbal[inH, iT, iz+1] * (1-fnH) * (1-fT) * (fz) + \
+ self.ionbal[inH+1, iT, iz+1] * (fnH) * (1-fT) * (fz) + \
+ self.ionbal[inH, iT+1, iz+1] * (1-fnH) * (fT) * (fz) + \
+ self.ionbal[inH+1, iT+1, iz] * (fnH) * (fT) * (1-fz) + \
+ self.ionbal[inH+1, iT+1, iz+1] * (fnH) * (fT) * (fz)
+
+ return 10**ionfrac
diff -r 0536fb6516915620927dcffb4281021a46b73214 -r 7a72b9113c4bb3dda4a90d45886f407d7da9acdd yt/visualization/fixed_resolution.py
--- a/yt/visualization/fixed_resolution.py
+++ b/yt/visualization/fixed_resolution.py
@@ -25,6 +25,8 @@
import _MPL
import numpy as np
import weakref
+import re
+import string
class FixedResolutionBuffer(object):
r"""
@@ -147,6 +149,41 @@
if f not in exclude and f[0] not in self.data_source.pf.particle_types:
self[f]
+
+ def _is_ion( self, fname ):
+ p = re.compile("_p[0-9]+_")
+ result = False
+ if p.search( fname ) != None:
+ result = True
+ return result
+
+ def _ion_to_label( self, fname ):
+ pnum2rom = {
+ "0":"I", "1":"II", "2":"III", "3":"IV", "4":"V",
+ "5":"VI", "6":"VII", "7":"VIII", "8":"IX", "9":"X",
+ "10":"XI", "11":"XII", "12":"XIII", "13":"XIV", "14":"XV",
+ "15":"XVI", "16":"XVII", "17":"XVIII", "18":"XIX", "19":"XX"}
+
+ p = re.compile("_p[0-9]+_")
+ m = p.search( fname )
+ if m != None:
+ pstr = m.string[m.start()+1:m.end()-1]
+ segments = fname.split("_")
+ for i,s in enumerate(segments):
+ segments[i] = string.capitalize(s)
+ if s == pstr:
+ ipstr = i
+ element = segments[ipstr-1]
+ roman = pnum2rom[pstr[1:]]
+ label = element + '\/' + roman + '\/' + \
+ string.join( segments[ipstr+1:], '\/' )
+ else:
+ label = fname
+ return label
+
+
+
+
def _get_info(self, item):
info = {}
ftype, fname = field = self.data_source._determine_fields(item)[0]
@@ -172,8 +209,13 @@
info['label'] = finfo.display_name
if info['label'] is None:
- info['label'] = r'$\rm{'+fname+r'}$'
- info['label'] = r'$\rm{'+fname.replace('_','\/').title()+r'}$'
+ if self._is_ion( fname ):
+ fname = self._ion_to_label( fname )
+ info['label'] = r'$\rm{'+fname+r'}$'
+ info['label'] = r'$\rm{'+fname.replace('_','\/')+r'}$'
+ else:
+ info['label'] = r'$\rm{'+fname+r'}$'
+ info['label'] = r'$\rm{'+fname.replace('_','\/').title()+r'}$'
elif info['label'].find('$') == -1:
info['label'] = info['label'].replace(' ','\/')
info['label'] = r'$\rm{'+info['label']+r'}$'
https://bitbucket.org/yt_analysis/yt/commits/942d37c8d66a/
Changeset: 942d37c8d66a
Branch: yt-3.0
User: galtay
Date: 2014-03-26 08:19:46
Summary: added more OWLS ion smoothed fields and aliases (gas and stars)
Affected #: 1 file
diff -r 7a72b9113c4bb3dda4a90d45886f407d7da9acdd -r 942d37c8d66a40694be457ef94a5a399a369e24e yt/frontends/sph/fields.py
--- a/yt/frontends/sph/fields.py
+++ b/yt/frontends/sph/fields.py
@@ -77,8 +77,13 @@
_elements = ("H", "He", "C", "N", "O", "Ne", "Mg", "Si", "Fe")
- # override
- #--------------------------------------------------------------
+ _num_neighbors = 48
+
+ _add_elements = ("PartType0", "PartType4")
+
+ _add_ions = ("PartType0")
+
+
def __init__(self, *args, **kwargs):
new_particle_fields = (
@@ -99,30 +104,70 @@
+ def setup_particle_fields(self, ptype):
+ """ additional particle fields derived from those in snapshot.
+ we also need to add the smoothed fields here b/c setup_fluid_fields
+ is called before setup_particle_fields. """
+ smoothed_suffixes = \
+ { "PartType0":
+ ("_number_density", "_density", "_mass"),
+ "PartType4":
+ ("_number_density", "_density", "_mass", "_fraction") }
- def setup_particle_fields(self, ptype):
- """ additional particle fields that are not on disk """
+ loaded = []
- super(OWLSFieldInfo,self).setup_particle_fields(ptype)
+ # we add particle element fields for stars and gas
+ #-----------------------------------------------------
+ if ptype in self._add_elements:
- if ptype == "PartType0":
# this adds the particle element fields
# X_density, X_mass, and X_number_density
- # where X is an element of self._elements.
- # X_fraction are defined in file
+ # where X is an item of self._elements.
+ # X_fraction are defined in snapshot
#-----------------------------------------------
for s in self._elements:
add_species_field_by_fraction(self, ptype, s,
particle_type=True)
+ # this needs to be called after the call to
+ # add_species_field_by_fraction for some reason ...
+ # not sure why yet.
+ #-------------------------------------------------------
+ super(OWLSFieldInfo,self).setup_particle_fields(ptype)
+
+
+ # and now we add the smoothed versions
+ #-----------------------------------------------------
+ if ptype in self._add_elements:
+
+ for s in self._elements:
+ for sfx in smoothed_suffixes[ptype]:
+ fname = s + sfx
+ fn = add_volume_weighted_smoothed_field(
+ ptype, "particle_position", "particle_mass",
+ "smoothing_length", "density", fname, self,
+ self._num_neighbors)
+ loaded += fn
+
+ if ptype == "PartType0":
+ self.alias(("gas", fname), fn[0])
+ elif ptype == "PartType4":
+ self.alias(("star", fname), fn[0])
+
+
+ # we only add ion fields for gas. this takes some
+ # time as the ion abundances have to be interpolated
+ # from cloudy tables (optically thin)
+ #-----------------------------------------------------
+ if ptype in self._add_ions:
+
# this defines the ion density on particles
- # X_density
+ # X_density for all items in self._ions
#-----------------------------------------------
- self.setup_ion_density_particle_fields()
-
+ self.setup_gas_ion_density_particle_fields( ptype )
# this adds the rest of the ion particle fields
# X_fraction, X_mass, X_number_density
@@ -141,34 +186,18 @@
pstr = "_p" + str(roman-1)
yt_ion = symbol + pstr
+ # add particle field
+ #---------------------------------------------------
add_species_field_by_density(self, ptype, yt_ion,
particle_type=True)
-
- ptype = 'PartType0'
- num_neighbors = 48
- loaded = []
-
- # add smoothed element fields
- #-----------------------------------------------
- for s in self._elements:
- fname = s + '_number_density'
- fn = add_volume_weighted_smoothed_field(
- ptype, "particle_position", "particle_mass",
- "smoothing_length", "density", fname, self,
- num_neighbors)
- loaded += fn
-
- self.alias(("gas", fname), fn[0])
-
- self._show_field_errors += loaded
-
-
# add smoothed ion fields
#-----------------------------------------------
for ion in self._ions:
+ # construct yt name for ion
+ #---------------------------------------------------
if ion[0:2].isalpha():
symbol = ion[0:2].capitalize()
roman = int(ion[2:])
@@ -179,36 +208,27 @@
pstr = "_p" + str(roman-1)
yt_ion = symbol + pstr
- fname = yt_ion + '_number_density'
- fn = add_volume_weighted_smoothed_field(
- ptype, "particle_position", "particle_mass",
- "smoothing_length", "density", fname, self,
- num_neighbors)
- loaded += fn
+ for sfx in smoothed_suffixes[ptype]:
+ fname = yt_ion + sfx
+ fn = add_volume_weighted_smoothed_field(
+ ptype, "particle_position", "particle_mass",
+ "smoothing_length", "density", fname, self,
+ self._num_neighbors)
+ loaded += fn
- self.alias(("gas", fname), fn[0])
+ self.alias(("gas", fname), fn[0])
- self._show_field_errors += loaded
- # find dependencies
- #-----------------------------------------------
- self.find_dependencies(loaded)
+ self._show_field_errors += loaded
+ self.find_dependencies(loaded)
- return
-
-
-
-
-
-
- def setup_ion_density_particle_fields( self ):
- """ Sets up particle fields for ion densities. """
+ def setup_gas_ion_density_particle_fields( self, ptype ):
+ """ Sets up particle fields for gas ion densities. """
# loop over all ions and make fields
#----------------------------------------------
-# for ion in ("h1",): # self._ions:
for ion in self._ions:
# construct yt name for ion
@@ -222,7 +242,7 @@
pstr = "_p" + str(roman-1)
yt_ion = symbol + pstr
- ftype = "PartType0"
+ ftype = ptype
# add ion density field for particles
#---------------------------------------------------
@@ -289,29 +309,6 @@
return
- # these add the smoothed element fields
- #-----------------------------------------------
- for s in self._elements:
- add_species_field_by_fraction(self, "gas", s)
-
-
-
- # this should add the smoothed ion fields
- #-----------------------------------------------
- for ion in self._ions:
-
- if ion[0:2].isalpha():
- symbol = ion[0:2].capitalize()
- roman = int(ion[2:])
- else:
- symbol = ion[0:1].capitalize()
- roman = int(ion[1:])
-
- pstr = "_p" + str(roman-1)
- yt_ion = symbol + pstr
-
- add_species_field_by_density(self, "gas", yt_ion)
-
# this function returns the owls_ion_data directory. if it doesn't
# exist it will download the data from http://yt-project.org/data
https://bitbucket.org/yt_analysis/yt/commits/cc8e7be98b6c/
Changeset: cc8e7be98b6c
Branch: yt-3.0
User: galtay
Date: 2014-03-26 09:33:07
Summary: attempting to fix star aliases and smoothing
Affected #: 2 files
diff -r 942d37c8d66a40694be457ef94a5a399a369e24e -r cc8e7be98b6cbdcec40b51e2f39ac4fa11752d4a yt/fields/field_info_container.py
--- a/yt/fields/field_info_container.py
+++ b/yt/fields/field_info_container.py
@@ -66,7 +66,7 @@
def setup_fluid_fields(self):
pass
- def setup_particle_fields(self, ptype):
+ def setup_particle_fields(self, ptype, ftype='gas', num_neighbors=48 ):
for f, (units, aliases, dn) in sorted(self.known_particle_fields):
self.add_output_field((ptype, f),
units = units, particle_type = True, display_name = dn)
@@ -99,7 +99,9 @@
self.add_output_field(field,
units = self.pf.field_units.get(field, ""),
particle_type = True)
- self.setup_smoothed_fields(ptype)
+ self.setup_smoothed_fields(ptype,
+ num_neighbors=num_neighbors,
+ ftype=ftype)
def setup_smoothed_fields(self, ptype, num_neighbors = 64, ftype = "gas"):
# We can in principle compute this, but it is not yet implemented.
diff -r 942d37c8d66a40694be457ef94a5a399a369e24e -r cc8e7be98b6cbdcec40b51e2f39ac4fa11752d4a yt/frontends/sph/fields.py
--- a/yt/frontends/sph/fields.py
+++ b/yt/frontends/sph/fields.py
@@ -116,7 +116,6 @@
("_number_density", "_density", "_mass", "_fraction") }
- loaded = []
# we add particle element fields for stars and gas
#-----------------------------------------------------
@@ -136,13 +135,24 @@
# add_species_field_by_fraction for some reason ...
# not sure why yet.
#-------------------------------------------------------
- super(OWLSFieldInfo,self).setup_particle_fields(ptype)
+ if ptype == 'PartType0':
+ ftype='gas'
+ elif ptype == 'PartType1':
+ ftype='dm'
+ elif ptype == 'PartType4':
+ ftype='star'
+ elif ptype == 'all':
+ ftype='all'
+
+ super(OWLSFieldInfo,self).setup_particle_fields(
+ ptype, num_neighbors=self._num_neighbors, ftype=ftype)
# and now we add the smoothed versions
#-----------------------------------------------------
if ptype in self._add_elements:
+ loaded = []
for s in self._elements:
for sfx in smoothed_suffixes[ptype]:
fname = s + sfx
@@ -156,6 +166,8 @@
self.alias(("gas", fname), fn[0])
elif ptype == "PartType4":
self.alias(("star", fname), fn[0])
+ self._show_field_errors += loaded
+ self.find_dependencies(loaded)
# we only add ion fields for gas. this takes some
@@ -208,6 +220,7 @@
pstr = "_p" + str(roman-1)
yt_ion = symbol + pstr
+ loaded = []
for sfx in smoothed_suffixes[ptype]:
fname = yt_ion + sfx
fn = add_volume_weighted_smoothed_field(
@@ -218,11 +231,10 @@
self.alias(("gas", fname), fn[0])
+ self._show_field_errors += loaded
+ self.find_dependencies(loaded)
- self._show_field_errors += loaded
- self.find_dependencies(loaded)
-
def setup_gas_ion_density_particle_fields( self, ptype ):
""" Sets up particle fields for gas ion densities. """
https://bitbucket.org/yt_analysis/yt/commits/c92bbd6c541f/
Changeset: c92bbd6c541f
Branch: yt-3.0
User: galtay
Date: 2014-03-28 20:10:25
Summary: removed smoothed star fields and reerted default num neighbors to 64
Affected #: 2 files
diff -r cc8e7be98b6cbdcec40b51e2f39ac4fa11752d4a -r c92bbd6c541f5a366a5d86c085af6a68bc7e03ea yt/fields/field_info_container.py
--- a/yt/fields/field_info_container.py
+++ b/yt/fields/field_info_container.py
@@ -66,7 +66,7 @@
def setup_fluid_fields(self):
pass
- def setup_particle_fields(self, ptype, ftype='gas', num_neighbors=48 ):
+ def setup_particle_fields(self, ptype, ftype='gas', num_neighbors=64 ):
for f, (units, aliases, dn) in sorted(self.known_particle_fields):
self.add_output_field((ptype, f),
units = units, particle_type = True, display_name = dn)
diff -r cc8e7be98b6cbdcec40b51e2f39ac4fa11752d4a -r c92bbd6c541f5a366a5d86c085af6a68bc7e03ea yt/frontends/sph/fields.py
--- a/yt/frontends/sph/fields.py
+++ b/yt/frontends/sph/fields.py
@@ -109,11 +109,7 @@
we also need to add the smoothed fields here b/c setup_fluid_fields
is called before setup_particle_fields. """
- smoothed_suffixes = \
- { "PartType0":
- ("_number_density", "_density", "_mass"),
- "PartType4":
- ("_number_density", "_density", "_mass", "_fraction") }
+ smoothed_suffixes = ("_number_density", "_density", "_mass"),
@@ -148,13 +144,13 @@
ptype, num_neighbors=self._num_neighbors, ftype=ftype)
- # and now we add the smoothed versions
+ # and now we add the smoothed versions for PartType0
#-----------------------------------------------------
- if ptype in self._add_elements:
+ if ptype == 'PartType0':
loaded = []
for s in self._elements:
- for sfx in smoothed_suffixes[ptype]:
+ for sfx in smoothed_suffixes:
fname = s + sfx
fn = add_volume_weighted_smoothed_field(
ptype, "particle_position", "particle_mass",
@@ -162,19 +158,17 @@
self._num_neighbors)
loaded += fn
- if ptype == "PartType0":
- self.alias(("gas", fname), fn[0])
- elif ptype == "PartType4":
- self.alias(("star", fname), fn[0])
+ self.alias(("gas", fname), fn[0])
+
self._show_field_errors += loaded
self.find_dependencies(loaded)
- # we only add ion fields for gas. this takes some
- # time as the ion abundances have to be interpolated
- # from cloudy tables (optically thin)
- #-----------------------------------------------------
- if ptype in self._add_ions:
+ # we only add ion fields for gas. this takes some
+ # time as the ion abundances have to be interpolated
+ # from cloudy tables (optically thin)
+ #-----------------------------------------------------
+
# this defines the ion density on particles
# X_density for all items in self._ions
https://bitbucket.org/yt_analysis/yt/commits/1234d06b45b3/
Changeset: 1234d06b45b3
Branch: yt-3.0
User: galtay
Date: 2014-03-28 20:17:22
Summary: fixed typos associated with smoothed_suffixes
Affected #: 1 file
diff -r c92bbd6c541f5a366a5d86c085af6a68bc7e03ea -r 1234d06b45b315e6ab9142ac781e9db3a2816621 yt/frontends/sph/fields.py
--- a/yt/frontends/sph/fields.py
+++ b/yt/frontends/sph/fields.py
@@ -109,7 +109,7 @@
we also need to add the smoothed fields here b/c setup_fluid_fields
is called before setup_particle_fields. """
- smoothed_suffixes = ("_number_density", "_density", "_mass"),
+ smoothed_suffixes = ("_number_density", "_density", "_mass")
@@ -215,7 +215,7 @@
yt_ion = symbol + pstr
loaded = []
- for sfx in smoothed_suffixes[ptype]:
+ for sfx in smoothed_suffixes:
fname = yt_ion + sfx
fn = add_volume_weighted_smoothed_field(
ptype, "particle_position", "particle_mass",
https://bitbucket.org/yt_analysis/yt/commits/ff4b94238c78/
Changeset: ff4b94238c78
Branch: yt-3.0
User: MatthewTurk
Date: 2014-03-28 20:09:41
Summary: Merged in galtay/yt/yt-3.0 (pull request #765)
OWLS ions and ion labels on slice/projections
Affected #: 7 files
diff -r 333edf0db58a8316aab87eed13f9971d1eedd0db -r ff4b94238c78fb1f97e4d3cb535274b5fed23049 yt/fields/field_info_container.py
--- a/yt/fields/field_info_container.py
+++ b/yt/fields/field_info_container.py
@@ -66,7 +66,7 @@
def setup_fluid_fields(self):
pass
- def setup_particle_fields(self, ptype):
+ def setup_particle_fields(self, ptype, ftype='gas', num_neighbors=64 ):
for f, (units, aliases, dn) in sorted(self.known_particle_fields):
units = self.pf.field_units.get((ptype, f), units)
self.add_output_field((ptype, f),
@@ -100,7 +100,9 @@
self.add_output_field(field,
units = self.pf.field_units.get(field, ""),
particle_type = True)
- self.setup_smoothed_fields(ptype)
+ self.setup_smoothed_fields(ptype,
+ num_neighbors=num_neighbors,
+ ftype=ftype)
def setup_smoothed_fields(self, ptype, num_neighbors = 64, ftype = "gas"):
# We can in principle compute this, but it is not yet implemented.
diff -r 333edf0db58a8316aab87eed13f9971d1eedd0db -r ff4b94238c78fb1f97e4d3cb535274b5fed23049 yt/fields/species_fields.py
--- a/yt/fields/species_fields.py
+++ b/yt/fields/species_fields.py
@@ -59,34 +59,46 @@
* data[ftype,'density']
return _density
-def add_species_field_by_density(registry, ftype, species):
+def add_species_field_by_density(registry, ftype, species,
+ particle_type = False):
"""
This takes a field registry, a fluid type, and a species name and then
adds the other fluids based on that. This assumes that the field
"SPECIES_density" already exists and refers to mass density.
"""
registry.add_field((ftype, "%s_fraction" % species),
- function = _create_fraction_func(ftype, species),
- units = "")
+ function = _create_fraction_func(ftype, species),
+ particle_type = particle_type,
+ units = "")
+
registry.add_field((ftype, "%s_mass" % species),
- function = _create_mass_func(ftype, species),
- units = "g")
+ function = _create_mass_func(ftype, species),
+ particle_type = particle_type,
+ units = "g")
+
registry.add_field((ftype, "%s_number_density" % species),
- function = _create_number_density_func(ftype, species),
- units = "cm**-3")
+ function = _create_number_density_func(ftype, species),
+ particle_type = particle_type,
+ units = "cm**-3")
-def add_species_field_by_fraction(registry, ftype, species):
+def add_species_field_by_fraction(registry, ftype, species,
+ particle_type = False):
"""
This takes a field registry, a fluid type, and a species name and then
adds the other fluids based on that. This assumes that the field
"SPECIES_fraction" already exists and refers to mass fraction.
"""
registry.add_field((ftype, "%s_density" % species),
- function = _create_density_func(ftype, species),
- units = "g/cm**3")
+ function = _create_density_func(ftype, species),
+ particle_type = particle_type,
+ units = "g/cm**3")
+
registry.add_field((ftype, "%s_mass" % species),
- function = _create_mass_func(ftype, species),
- units = "g")
+ function = _create_mass_func(ftype, species),
+ particle_type = particle_type,
+ units = "g")
+
registry.add_field((ftype, "%s_number_density" % species),
- function = _create_number_density_func(ftype, species),
- units = "cm**-3")
+ function = _create_number_density_func(ftype, species),
+ particle_type = particle_type,
+ units = "cm**-3")
diff -r 333edf0db58a8316aab87eed13f9971d1eedd0db -r ff4b94238c78fb1f97e4d3cb535274b5fed23049 yt/frontends/sph/data_structures.py
--- a/yt/frontends/sph/data_structures.py
+++ b/yt/frontends/sph/data_structures.py
@@ -209,7 +209,10 @@
if "length" in unit_base:
length_unit = unit_base["length"]
elif "UnitLength_in_cm" in unit_base:
- length_unit = (unit_base["UnitLength_in_cm"], "cm")
+ if self.cosmological_simulation == 0:
+ length_unit = (unit_base["UnitLength_in_cm"], "cm")
+ else:
+ length_unit = (unit_base["UnitLength_in_cm"], "cmcm/h")
else:
raise RuntimeError
length_unit = _fix_unit_ordering(length_unit)
@@ -229,7 +232,10 @@
if "mass" in unit_base:
mass_unit = unit_base["mass"]
elif "UnitMass_in_g" in unit_base:
- mass_unit = (unit_base["UnitMass_in_g"], "g")
+ if self.cosmological_simulation == 0:
+ mass_unit = (unit_base["UnitMass_in_g"], "g")
+ else:
+ mass_unit = (unit_base["UnitMass_in_g"], "g/h")
else:
# Sane default
mass_unit = (1.0, "1e10*Msun/h")
@@ -286,6 +292,8 @@
_particle_mass_name = "Mass"
_field_info_class = OWLSFieldInfo
+
+
def _parse_parameter_file(self):
handle = h5py.File(self.parameter_filename, mode="r")
hvals = {}
diff -r 333edf0db58a8316aab87eed13f9971d1eedd0db -r ff4b94238c78fb1f97e4d3cb535274b5fed23049 yt/frontends/sph/fields.py
--- a/yt/frontends/sph/fields.py
+++ b/yt/frontends/sph/fields.py
@@ -14,17 +14,27 @@
# The full license is in the file COPYING.txt, distributed with this software.
#-----------------------------------------------------------------------------
+import os
import numpy as np
+import owls_ion_tables as oit
from yt.funcs import *
+
from yt.fields.field_info_container import \
FieldInfoContainer
+
from .definitions import \
gadget_ptypes, \
ghdf5_ptypes
-from yt.fields.species_fields import add_species_field_by_fraction
+from yt.config import ytcfg
+from yt.utilities.physical_constants import mh
+from yt.fields.species_fields import \
+ add_species_field_by_fraction, \
+ add_species_field_by_density
+from yt.fields.particle_fields import \
+ add_volume_weighted_smoothed_field
# Here are helper functions for things like vector fields and so on.
@@ -90,24 +100,32 @@
class OWLSFieldInfo(SPHFieldInfo):
- _species_fractions = ['H_fraction', 'He_fraction', 'C_fraction',
- 'N_fraction', 'O_fraction', 'Ne_fraction',
- 'Mg_fraction', 'Si_fraction', 'Fe_fraction']
+ _ions = ("c1", "c2", "c3", "c4", "c5", "c6",
+ "fe2", "fe17", "h1", "he1", "he2", "mg1", "mg2", "n2",
+ "n3", "n4", "n5", "n6", "n7", "ne8", "ne9", "ne10", "o1",
+ "o6", "o7", "o8", "si2", "si3", "si4", "si13")
- # override
- #--------------------------------------------------------------
+ _elements = ("H", "He", "C", "N", "O", "Ne", "Mg", "Si", "Fe")
+
+ _num_neighbors = 48
+
+ _add_elements = ("PartType0", "PartType4")
+
+ _add_ions = ("PartType0")
+
+
def __init__(self, *args, **kwargs):
new_particle_fields = (
- ('Hydrogen', ('', ['H_fraction'], None)),
- ('Helium', ('', ['He_fraction'], None)),
- ('Carbon', ('', ['C_fraction'], None)),
- ('Nitrogen', ('', ['N_fraction'], None)),
- ('Oxygen', ('', ['O_fraction'], None)),
- ('Neon', ('', ['Ne_fraction'], None)),
- ('Magnesium', ('', ['Mg_fraction'], None)),
- ('Silicon', ('', ['Si_fraction'], None)),
- ('Iron', ('', ['Fe_fraction'], None))
+ ("Hydrogen", ("", ["H_fraction"], None)),
+ ("Helium", ("", ["He_fraction"], None)),
+ ("Carbon", ("", ["C_fraction"], None)),
+ ("Nitrogen", ("", ["N_fraction"], None)),
+ ("Oxygen", ("", ["O_fraction"], None)),
+ ("Neon", ("", ["Ne_fraction"], None)),
+ ("Magnesium", ("", ["Mg_fraction"], None)),
+ ("Silicon", ("", ["Si_fraction"], None)),
+ ("Iron", ("", ["Fe_fraction"], None))
)
self.known_particle_fields += new_particle_fields
@@ -115,9 +133,255 @@
super(OWLSFieldInfo,self).__init__( *args, **kwargs )
+
+ def setup_particle_fields(self, ptype):
+ """ additional particle fields derived from those in snapshot.
+ we also need to add the smoothed fields here b/c setup_fluid_fields
+ is called before setup_particle_fields. """
+
+ smoothed_suffixes = ("_number_density", "_density", "_mass")
+
+
+
+ # we add particle element fields for stars and gas
+ #-----------------------------------------------------
+ if ptype in self._add_elements:
+
+
+ # this adds the particle element fields
+ # X_density, X_mass, and X_number_density
+ # where X is an item of self._elements.
+ # X_fraction are defined in snapshot
+ #-----------------------------------------------
+ for s in self._elements:
+ add_species_field_by_fraction(self, ptype, s,
+ particle_type=True)
+
+ # this needs to be called after the call to
+ # add_species_field_by_fraction for some reason ...
+ # not sure why yet.
+ #-------------------------------------------------------
+ if ptype == 'PartType0':
+ ftype='gas'
+ elif ptype == 'PartType1':
+ ftype='dm'
+ elif ptype == 'PartType4':
+ ftype='star'
+ elif ptype == 'all':
+ ftype='all'
+ super(OWLSFieldInfo,self).setup_particle_fields(
+ ptype, num_neighbors=self._num_neighbors, ftype=ftype)
+
+
+ # and now we add the smoothed versions for PartType0
+ #-----------------------------------------------------
+ if ptype == 'PartType0':
+
+ loaded = []
+ for s in self._elements:
+ for sfx in smoothed_suffixes:
+ fname = s + sfx
+ fn = add_volume_weighted_smoothed_field(
+ ptype, "particle_position", "particle_mass",
+ "smoothing_length", "density", fname, self,
+ self._num_neighbors)
+ loaded += fn
+
+ self.alias(("gas", fname), fn[0])
+
+ self._show_field_errors += loaded
+ self.find_dependencies(loaded)
+
+
+ # we only add ion fields for gas. this takes some
+ # time as the ion abundances have to be interpolated
+ # from cloudy tables (optically thin)
+ #-----------------------------------------------------
+
+
+ # this defines the ion density on particles
+ # X_density for all items in self._ions
+ #-----------------------------------------------
+ self.setup_gas_ion_density_particle_fields( ptype )
+
+ # this adds the rest of the ion particle fields
+ # X_fraction, X_mass, X_number_density
+ #-----------------------------------------------
+ for ion in self._ions:
+
+ # construct yt name for ion
+ #---------------------------------------------------
+ if ion[0:2].isalpha():
+ symbol = ion[0:2].capitalize()
+ roman = int(ion[2:])
+ else:
+ symbol = ion[0:1].capitalize()
+ roman = int(ion[1:])
+
+ pstr = "_p" + str(roman-1)
+ yt_ion = symbol + pstr
+
+ # add particle field
+ #---------------------------------------------------
+ add_species_field_by_density(self, ptype, yt_ion,
+ particle_type=True)
+
+
+ # add smoothed ion fields
+ #-----------------------------------------------
+ for ion in self._ions:
+
+ # construct yt name for ion
+ #---------------------------------------------------
+ if ion[0:2].isalpha():
+ symbol = ion[0:2].capitalize()
+ roman = int(ion[2:])
+ else:
+ symbol = ion[0:1].capitalize()
+ roman = int(ion[1:])
+
+ pstr = "_p" + str(roman-1)
+ yt_ion = symbol + pstr
+
+ loaded = []
+ for sfx in smoothed_suffixes:
+ fname = yt_ion + sfx
+ fn = add_volume_weighted_smoothed_field(
+ ptype, "particle_position", "particle_mass",
+ "smoothing_length", "density", fname, self,
+ self._num_neighbors)
+ loaded += fn
+
+ self.alias(("gas", fname), fn[0])
+
+ self._show_field_errors += loaded
+ self.find_dependencies(loaded)
+
+
+
+ def setup_gas_ion_density_particle_fields( self, ptype ):
+ """ Sets up particle fields for gas ion densities. """
+
+ # loop over all ions and make fields
+ #----------------------------------------------
+ for ion in self._ions:
+
+ # construct yt name for ion
+ #---------------------------------------------------
+ if ion[0:2].isalpha():
+ symbol = ion[0:2].capitalize()
+ roman = int(ion[2:])
+ else:
+ symbol = ion[0:1].capitalize()
+ roman = int(ion[1:])
+
+ pstr = "_p" + str(roman-1)
+ yt_ion = symbol + pstr
+ ftype = ptype
+
+ # add ion density field for particles
+ #---------------------------------------------------
+ fname = yt_ion + '_density'
+ dens_func = self._create_ion_density_func( ftype, ion )
+ self.add_field( (ftype, fname),
+ function = dens_func,
+ units="g/cm**3",
+ particle_type=True )
+ self._show_field_errors.append( (ftype,fname) )
+
+
+
+
+ def _create_ion_density_func( self, ftype, ion ):
+ """ returns a function that calculates the ion density of a particle.
+ """
+
+ def _ion_density(field, data):
+
+ # get element symbol from ion string. ion string will
+ # be a member of the tuple _ions (i.e. si13)
+ #--------------------------------------------------------
+ if ion[0:2].isalpha():
+ symbol = ion[0:2].capitalize()
+ else:
+ symbol = ion[0:1].capitalize()
+
+ # mass fraction for the element
+ #--------------------------------------------------------
+ m_frac = data[ftype, symbol+"_fraction"]
+
+ # get nH and T for lookup
+ #--------------------------------------------------------
+ log_nH = np.log10( data["PartType0", "H_number_density"] )
+ log_T = np.log10( data["PartType0", "Temperature"] )
+
+ # get name of owls_ion_file for given ion
+ #--------------------------------------------------------
+ owls_ion_path = self._get_owls_ion_data_dir()
+ fname = os.path.join( owls_ion_path, ion+".hdf5" )
+
+ # create ionization table for this redshift
+ #--------------------------------------------------------
+ itab = oit.IonTableOWLS( fname )
+ itab.set_iz( data.pf.current_redshift )
+
+ # find ion balance using log nH and log T
+ #--------------------------------------------------------
+ i_frac = itab.interp( log_nH, log_T )
+ return data[ftype,"Density"] * m_frac * i_frac
+
+ return _ion_density
+
+
+
+
+
+ # this function sets up the X_mass, X_density, X_fraction, and
+ # X_number_density fields where X is the name of an OWLS element.
+ #-------------------------------------------------------------
def setup_fluid_fields(self):
- # here species_name is "H", "He", etc
- for s in self._species_fractions:
- species_name = s.split('_')[0]
- add_species_field_by_fraction(self, "gas", species_name)
+
+ return
+
+
+
+ # this function returns the owls_ion_data directory. if it doesn't
+ # exist it will download the data from http://yt-project.org/data
+ #-------------------------------------------------------------
+ def _get_owls_ion_data_dir(self):
+
+ txt = "Attempting to download ~ 30 Mb of owls ion data from %s to %s."
+ data_file = "owls_ion_data.tar.gz"
+ data_url = "http://yt-project.org/data"
+
+ # get test_data_dir from yt config (ytcgf)
+ #----------------------------------------------
+ tdir = ytcfg.get("yt","test_data_dir")
+
+ # set download destination to tdir or ./ if tdir isnt defined
+ #----------------------------------------------
+ if tdir == "/does/not/exist":
+ data_dir = "./"
+ else:
+ data_dir = tdir
+
+
+ # check for owls_ion_data directory in data_dir
+ # if not there download the tarball and untar it
+ #----------------------------------------------
+ owls_ion_path = os.path.join( data_dir, "owls_ion_data" )
+
+ if not os.path.exists(owls_ion_path):
+ mylog.info(txt % (data_url, data_dir))
+ fname = data_dir + "/" + data_file
+ fn = download_file(os.path.join(data_url, data_file), fname)
+
+ cmnd = "cd " + data_dir + "; " + "tar xf " + data_file
+ os.system(cmnd)
+
+
+ if not os.path.exists(owls_ion_path):
+ raise RuntimeError, "Failed to download owls ion data."
+
+ return owls_ion_path
diff -r 333edf0db58a8316aab87eed13f9971d1eedd0db -r ff4b94238c78fb1f97e4d3cb535274b5fed23049 yt/frontends/sph/io.py
--- a/yt/frontends/sph/io.py
+++ b/yt/frontends/sph/io.py
@@ -53,7 +53,7 @@
_vector_fields = ("Coordinates", "Velocity", "Velocities")
_known_ptypes = ghdf5_ptypes
_var_mass = None
- _element_fields = ('Hydrogen', 'Helium', 'Carbon', 'Nitrogen', 'Oxygen',
+ _element_names = ('Hydrogen', 'Helium', 'Carbon', 'Nitrogen', 'Oxygen',
'Neon', 'Magnesium', 'Silicon', 'Iron' )
@@ -110,7 +110,7 @@
ind = self._known_ptypes.index(ptype)
data[:] = self.pf["Massarr"][ind]
- elif field in self._element_fields:
+ elif field in self._element_names:
rfield = 'ElementAbundance/' + field
data = g[rfield][:][mask,...]
diff -r 333edf0db58a8316aab87eed13f9971d1eedd0db -r ff4b94238c78fb1f97e4d3cb535274b5fed23049 yt/frontends/sph/owls_ion_tables.py
--- /dev/null
+++ b/yt/frontends/sph/owls_ion_tables.py
@@ -0,0 +1,224 @@
+"""
+OWLS ion tables
+
+A module to handle the HM01 UV background spectra and ionization data from the
+OWLS photoionization equilibrium lookup tables.
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2013, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+import sys
+import h5py
+import numpy as np
+
+
+
+
+def h5rd( fname, path, dtype=None ):
+ """ Read Data. Return a dataset located at <path> in file <fname> as
+ a numpy array.
+ e.g. rd( fname, '/PartType0/Coordinates' ). """
+
+ data = None
+ with h5py.File( fname, 'r' ) as h5f:
+ ds = h5f[path]
+ if dtype == None:
+ dtype = ds.dtype
+ data = np.zeros( ds.shape, dtype=dtype )
+ data = ds.value
+ return data
+
+
+
+class IonTableSpectrum:
+
+ """ A class to handle the HM01 spectra in the OWLS ionization tables. """
+
+ def __init__(self, ion_file):
+
+ where = '/header/spectrum/gammahi'
+ self.GH1 = h5rd( ion_file, where ) # GH1[1/s]
+
+ where = '/header/spectrum/logenergy_ryd'
+ self.logryd = h5rd( ion_file, where ) # E[ryd]
+
+ where = '/header/spectrum/logflux'
+ self.logflux = h5rd( ion_file, where ) # J[ergs/s/Hz/Sr/cm^2]
+
+ where = '/header/spectrum/redshift'
+ self.z = h5rd( ion_file, where ) # z
+
+
+
+ def return_table_GH1_at_z(self,z):
+
+ # find redshift indices
+ #-----------------------------------------------------------------
+ i_zlo = np.argmin( np.abs( self.z - z ) )
+ if self.z[i_zlo] < z:
+ i_zhi = i_zlo + 1
+ else:
+ i_zhi = i_zlo
+ i_zlo = i_zlo - 1
+
+ z_frac = (z - self.z[i_zlo]) / (self.z[i_zhi] - self.z[i_zlo])
+
+ # find GH1 from table
+ #-----------------------------------------------------------------
+ logGH1_all = np.log10( self.GH1 )
+ dlog_GH1 = logGH1_all[i_zhi] - logGH1_all[i_zlo]
+
+ logGH1_table = logGH1_all[i_zlo] + z_frac * dlog_GH1
+ GH1_table = 10.0**logGH1_table
+
+ return GH1_table
+
+
+
+
+class IonTableOWLS:
+
+ """ A class to handle OWLS ionization tables. """
+
+ DELTA_nH = 0.25
+ DELTA_T = 0.1
+
+ def __init__(self, ion_file):
+
+ self.ion_file = ion_file
+
+ # ionbal is indexed like [nH, T, z]
+ # nH and T are log quantities
+ #---------------------------------------------------------------
+ self.nH = h5rd( ion_file, '/logd' ) # log nH [cm^-3]
+ self.T = h5rd( ion_file, '/logt' ) # log T [K]
+ self.z = h5rd( ion_file, '/redshift' ) # z
+
+ # read the ionization fractions
+ # linear values stored in file so take log here
+ # ionbal is the ionization balance (i.e. fraction)
+ #---------------------------------------------------------------
+ self.ionbal = h5rd( ion_file, '/ionbal' ).astype(np.float64)
+ self.ionbal_orig = self.ionbal.copy()
+
+ ipositive = np.where( self.ionbal > 0.0 )
+ izero = np.where( self.ionbal <= 0.0 )
+ self.ionbal[izero] = self.ionbal[ipositive].min()
+
+ self.ionbal = np.log10( self.ionbal )
+
+
+ # load in background spectrum
+ #---------------------------------------------------------------
+ self.spectrum = IonTableSpectrum( ion_file )
+
+ # calculate the spacing along each dimension
+ #---------------------------------------------------------------
+ self.dnH = self.nH[1:] - self.nH[0:-1]
+ self.dT = self.T[1:] - self.T[0:-1]
+ self.dz = self.z[1:] - self.z[0:-1]
+
+ self.order_str = '[log nH, log T, z]'
+
+
+
+
+
+ # sets iz and fz
+ #-----------------------------------------------------
+ def set_iz( self, z ):
+
+ if z <= self.z[0]:
+ self.iz = 0
+ self.fz = 0.0
+ elif z >= self.z[-1]:
+ self.iz = len(self.z) - 2
+ self.fz = 1.0
+ else:
+ for iz in range( len(self.z)-1 ):
+ if z < self.z[iz+1]:
+ self.iz = iz
+ self.fz = ( z - self.z[iz] ) / self.dz[iz]
+ break
+
+
+
+ # interpolate the table at a fixed redshift for the input
+ # values of nH and T ( input should be log ). A simple
+ # tri-linear interpolation is used.
+ #-----------------------------------------------------
+ def interp( self, nH, T ):
+
+ nH = np.array( nH )
+ T = np.array( T )
+
+ if nH.size != T.size:
+ print ' array size mismatch !!! '
+ sys.exit(1)
+
+ # field discovery will have nH.size == 1 and T.size == 1
+ # in that case we simply return 1.0
+
+ if nH.size == 1 and T.size == 1:
+ ionfrac = 1.0
+ return ionfrac
+
+
+ # find inH and fnH
+ #-----------------------------------------------------
+ inH = np.int32( ( nH - self.nH[0] ) / self.DELTA_nH )
+ fnH = ( nH - self.nH[inH] ) / self.dnH[inH]
+
+ indx = np.where( inH < 0 )[0]
+ if len(indx) > 0:
+ inH[indx] = 0
+ fnH[indx] = 0.0
+
+ indx = np.where( inH >= len(nH) )[0]
+ if len(indx) > 0:
+ inH[indx] = len(nH)-2
+ fnH[indx] = 1.0
+
+
+ # find iT and fT
+ #-----------------------------------------------------
+ iT = np.int32( ( T - self.T[0] ) / self.DELTA_T )
+ fT = ( T - self.T[iT] ) / self.dT[iT]
+
+ indx = np.where( iT < 0 )[0]
+ if len(indx) > 0:
+ iT[indx] = 0
+ fT[indx] = 0.0
+
+ indx = np.where( iT >= len(T) )[0]
+ if len(indx) > 0:
+ iT[indx] = len(T)-2
+ fT[indx] = 1.0
+
+
+ iz = self.iz
+ fz = self.fz
+
+ # calculate interpolated value
+ # use tri-linear interpolation on the log values
+ #-----------------------------------------------------
+
+ ionfrac = self.ionbal[inH, iT, iz ] * (1-fnH) * (1-fT) * (1-fz) + \
+ self.ionbal[inH+1, iT, iz ] * (fnH) * (1-fT) * (1-fz) + \
+ self.ionbal[inH, iT+1, iz ] * (1-fnH) * (fT) * (1-fz) + \
+ self.ionbal[inH, iT, iz+1] * (1-fnH) * (1-fT) * (fz) + \
+ self.ionbal[inH+1, iT, iz+1] * (fnH) * (1-fT) * (fz) + \
+ self.ionbal[inH, iT+1, iz+1] * (1-fnH) * (fT) * (fz) + \
+ self.ionbal[inH+1, iT+1, iz] * (fnH) * (fT) * (1-fz) + \
+ self.ionbal[inH+1, iT+1, iz+1] * (fnH) * (fT) * (fz)
+
+ return 10**ionfrac
diff -r 333edf0db58a8316aab87eed13f9971d1eedd0db -r ff4b94238c78fb1f97e4d3cb535274b5fed23049 yt/visualization/fixed_resolution.py
--- a/yt/visualization/fixed_resolution.py
+++ b/yt/visualization/fixed_resolution.py
@@ -25,6 +25,8 @@
import _MPL
import numpy as np
import weakref
+import re
+import string
class FixedResolutionBuffer(object):
r"""
@@ -147,6 +149,41 @@
if f not in exclude and f[0] not in self.data_source.pf.particle_types:
self[f]
+
+ def _is_ion( self, fname ):
+ p = re.compile("_p[0-9]+_")
+ result = False
+ if p.search( fname ) != None:
+ result = True
+ return result
+
+ def _ion_to_label( self, fname ):
+ pnum2rom = {
+ "0":"I", "1":"II", "2":"III", "3":"IV", "4":"V",
+ "5":"VI", "6":"VII", "7":"VIII", "8":"IX", "9":"X",
+ "10":"XI", "11":"XII", "12":"XIII", "13":"XIV", "14":"XV",
+ "15":"XVI", "16":"XVII", "17":"XVIII", "18":"XIX", "19":"XX"}
+
+ p = re.compile("_p[0-9]+_")
+ m = p.search( fname )
+ if m != None:
+ pstr = m.string[m.start()+1:m.end()-1]
+ segments = fname.split("_")
+ for i,s in enumerate(segments):
+ segments[i] = string.capitalize(s)
+ if s == pstr:
+ ipstr = i
+ element = segments[ipstr-1]
+ roman = pnum2rom[pstr[1:]]
+ label = element + '\/' + roman + '\/' + \
+ string.join( segments[ipstr+1:], '\/' )
+ else:
+ label = fname
+ return label
+
+
+
+
def _get_info(self, item):
info = {}
ftype, fname = field = self.data_source._determine_fields(item)[0]
@@ -172,8 +209,13 @@
info['label'] = finfo.display_name
if info['label'] is None:
- info['label'] = r'$\rm{'+fname+r'}$'
- info['label'] = r'$\rm{'+fname.replace('_','\/').title()+r'}$'
+ if self._is_ion( fname ):
+ fname = self._ion_to_label( fname )
+ info['label'] = r'$\rm{'+fname+r'}$'
+ info['label'] = r'$\rm{'+fname.replace('_','\/')+r'}$'
+ else:
+ info['label'] = r'$\rm{'+fname+r'}$'
+ info['label'] = r'$\rm{'+fname.replace('_','\/').title()+r'}$'
elif info['label'].find('$') == -1:
info['label'] = info['label'].replace(' ','\/')
info['label'] = r'$\rm{'+info['label']+r'}$'
Repository URL: https://bitbucket.org/yt_analysis/yt/
--
This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.
More information about the yt-svn
mailing list