[yt-svn] commit/yt: 10 new changesets

commits-noreply at bitbucket.org commits-noreply at bitbucket.org
Thu Mar 20 04:45:10 PDT 2014


10 new commits in yt:

https://bitbucket.org/yt_analysis/yt/commits/a7257198c122/
Changeset:   a7257198c122
Branch:      yt-3.0
User:        MatthewTurk
Date:        2014-03-14 01:58:38
Summary:     Adding periodic table module
Affected #:  3 files

diff -r cc4e15d9f4d95b558690cd94e4f74b30bed6a02b -r a7257198c1227679c148c8a6aea54e96db8ea8bd yt/fields/species_fields.py
--- a/yt/fields/species_fields.py
+++ b/yt/fields/species_fields.py
@@ -31,6 +31,7 @@
     standard_particle_fields
 from yt.utilities.physical_constants import \
     mh, \
-    mass_sun_cgs
+    mass_sun_cgs, \
+    periodic_table
 from yt.funcs import *
 

diff -r cc4e15d9f4d95b558690cd94e4f74b30bed6a02b -r a7257198c1227679c148c8a6aea54e96db8ea8bd yt/utilities/periodic_table.py
--- /dev/null
+++ b/yt/utilities/periodic_table.py
@@ -0,0 +1,175 @@
+"""
+A simple periodic table.
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# 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 numpy as np
+import numbers
+import types
+
+_elements = (
+    (1, 1.0079400000, "Hydrogen", "H"),
+    (2, 4.0026020000, "Helium", "He"),
+    (3, 6.9410000000, "Lithium", "Li"),
+    (4, 9.0121820000, "Beryllium", "Be"),
+    (5, 10.8110000000, "Boron", "B"),
+    (6, 12.0107000000, "Carbon", "C"),
+    (7, 14.0067000000, "Nitrogen", "N"),
+    (8, 15.9994000000, "Oxygen", "O"),
+    (9, 18.9994000000, "Fluorine", "F"),
+    (10, 20.1797000000, "Neon", "Ne"),
+    (11, 22.9897692800, "Sodium", "Na"),
+    (12, 24.3050000000, "Magnesium", "Mg"),
+    (13, 26.9815386000, "Aluminium", "Al"),
+    (14, 28.0855000000, "Silicon", "Si"),
+    (15, 30.9737620000, "Phosphorus", "P"),
+    (16, 32.0650000000, "Sulphur", "S"),
+    (17, 35.4530000000, "Chlorine", "Cl"),
+    (18, 39.9480000000, "Argon", "Ar"),
+    (19, 39.0983000000, "Potassium", "K"),
+    (20, 40.0780000000, "Calcium", "Ca"),
+    (21, 44.9559120000, "Scandium", "Sc"),
+    (22, 47.8670000000, "Titanium", "Ti"),
+    (23, 50.9415000000, "Vanadium", "V"),
+    (24, 51.9961000000, "Chromium", "Cr"),
+    (25, 54.9380450000, "Manganese", "Mn"),
+    (26, 55.8450000000, "Iron", "Fe"),
+    (27, 58.9331950000, "Cobalt", "Co"),
+    (28, 58.6934000000, "Nickel", "Ni"),
+    (29, 63.5460000000, "Copper", "Cu"),
+    (30, 65.3800000000, "Zinc", "Zn"),
+    (31, 69.7230000000, "Gallium", "Ga"),
+    (32, 72.6400000000, "Germanium", "Ge"),
+    (33, 74.9216000000, "Arsenic", "As"),
+    (34, 78.9600000000, "Selenium", "Se"),
+    (35, 79.9040000000, "Bromine", "Br"),
+    (36, 83.7980000000, "Krypton", "Kr"),
+    (37, 85.4678000000, "Rubidium", "Rb"),
+    (38, 87.6200000000, "Strontium", "Sr"),
+    (39, 88.9058500000, "Yttrium", "Y"),
+    (40, 91.2240000000, "Zirkonium", "Zr"),
+    (41, 92.9063800000, "Niobium", "Nb"),
+    (42, 95.9600000000, "Molybdaenum", "Mo"),
+    (43, 98.0000000000, "Technetium", "Tc"),
+    (44, 101.0700000000, "Ruthenium", "Ru"),
+    (45, 102.9055000000, "Rhodium", "Rh"),
+    (46, 106.4200000000, "Palladium", "Pd"),
+    (47, 107.8682000000, "Silver", "Ag"),
+    (48, 112.4110000000, "Cadmium", "Cd"),
+    (49, 114.8180000000, "Indium", "In"),
+    (50, 118.7100000000, "Tin", "Sn"),
+    (51, 121.7600000000, "Antimony", "Sb"),
+    (52, 127.6000000000, "Tellurium", "Te"),
+    (53, 126.9044700000, "Iodine", "I"),
+    (54, 131.2930000000, "Xenon", "Xe"),
+    (55, 132.9054519000, "Cesium", "Cs"),
+    (56, 137.3270000000, "Barium", "Ba"),
+    (57, 138.9054700000, "Lanthanum", "La"),
+    (58, 140.1160000000, "Cerium", "Ce"),
+    (59, 140.9076500000, "Praseodymium", "Pr"),
+    (60, 144.2420000000, "Neodymium", "Nd"),
+    (61, 145.0000000000, "Promethium", "Pm"),
+    (62, 150.3600000000, "Samarium", "Sm"),
+    (63, 151.9640000000, "Europium", "Eu"),
+    (64, 157.2500000000, "Gadolinium", "Gd"),
+    (65, 158.9253500000, "Terbium", "Tb"),
+    (66, 162.5001000000, "Dysprosium", "Dy"),
+    (67, 164.9303200000, "Holmium", "Ho"),
+    (68, 167.2590000000, "Erbium", "Er"),
+    (69, 168.9342100000, "Thulium", "Tm"),
+    (70, 173.0540000000, "Ytterbium", "Yb"),
+    (71, 174.9668000000, "Lutetium", "Lu"),
+    (72, 178.4900000000, "Hafnium", "Hf"),
+    (73, 180.9478800000, "Tantalum", "Ta"),
+    (74, 183.8400000000, "Tungsten", "W"),
+    (75, 186.2070000000, "Rhenium", "Re"),
+    (76, 190.2300000000, "Osmium", "Os"),
+    (77, 192.2170000000, "Iridium", "Ir"),
+    (78, 192.0840000000, "Platinum", "Pt"),
+    (79, 196.9665690000, "Gold", "Au"),
+    (80, 200.5900000000, "Hydrargyrum", "Hg"),
+    (81, 204.3833000000, "Thallium", "Tl"),
+    (82, 207.2000000000, "Lead", "Pb"),
+    (83, 208.9804010000, "Bismuth", "Bi"),
+    (84, 210.0000000000, "Polonium", "Po"),
+    (85, 210.0000000000, "Astatine", "At"),
+    (86, 220.0000000000, "Radon", "Rn"),
+    (87, 223.0000000000, "Francium", "Fr"),
+    (88, 226.0000000000, "Radium", "Ra"),
+    (89, 227.0000000000, "Actinium", "Ac"),
+    (90, 232.0380600000, "Thorium", "Th"),
+    (91, 231.0358800000, "Protactinium", "Pa"),
+    (92, 238.0289100000, "Uranium", "U"),
+    (93, 237.0000000000, "Neptunium", "Np"),
+    (94, 244.0000000000, "Plutonium", "Pu"),
+    (95, 243.0000000000, "Americium", "Am"),
+    (96, 247.0000000000, "Curium", "Cm"),
+    (97, 247.0000000000, "Berkelium", "Bk"),
+    (98, 251.0000000000, "Californium", "Cf"),
+    (99, 252.0000000000, "Einsteinium", "Es"),
+    (100, 257.0000000000, "Fermium", "Fm"),
+    (101, 258.0000000000, "Mendelevium", "Md"),
+    (102, 259.0000000000, "Nobelium", "No"),
+    (103, 262.0000000000, "Lawrencium", "Lr"),
+    (104, 261.0000000000, "Rutherfordium", "Rf"),
+    (105, 262.0000000000, "Dubnium", "Db"),
+    (106, 266.0000000000, "Seaborgium", "Sg"),
+    (107, 264.0000000000, "Bohrium", "Bh"),
+    (108, 277.0000000000, "Hassium", "Hs"),
+    (109, 268.0000000000, "Meitnerium", "Mt"),
+    (110, 271.0000000000, "Ununnilium", "Ds"),
+    (111, 272.0000000000, "Unununium", "Rg"),
+    (112, 285.0000000000, "Ununbium", "Uub"),
+    (113, 284.0000000000, "Ununtrium", "Uut"),
+    (114, 289.0000000000, "Ununquadium", "Uuq"),
+    (115, 288.0000000000, "Ununpentium", "Uup"),
+    (116, 292.0000000000, "Ununhexium", "Uuh"),
+    (118, 294.0000000000, "Ununoctium", "Uuo"),
+)
+
+class Element:
+    def __init__(self, num, weight, name, symbol):
+        self.num = num
+        self.weight = weight
+        self.name = name
+        self.symbol = symbol
+
+    def __repr__(self):
+        return "Element: %s (%s)" % (self.symbol, self.name)
+
+class PeriodicTable:
+    def __init__(self):
+        self.elements_by_number = {}
+        self.elements_by_name = {}
+        self.elements_by_symbol = {}
+        for num, weight, name, symbol in _elements:
+            e = Element(num, weight, name, symbol)
+            self.elements_by_number[num] = e
+            self.elements_by_name[name] = e
+            self.elements_by_symbol[symbol] = e
+
+    def __getitem__(self, key):
+        if isinstance(key, (np.number, numbers.Number)):
+            d = self.elements_by_number
+        elif isinstance(key, types.StringTypes):
+            if len(key) <= 2:
+                d = self.elements_by_symbol
+            elif len(key) == 3 and key[0] == "U":
+                d = self.elements_by_symbol
+            else:
+                d = self.elements_by_name
+        else:
+            raise KeyError(key)
+        return d[key]
+
+periodic_table = PeriodicTable()

diff -r cc4e15d9f4d95b558690cd94e4f74b30bed6a02b -r a7257198c1227679c148c8a6aea54e96db8ea8bd yt/utilities/tests/test_periodic_table.py
--- /dev/null
+++ b/yt/utilities/tests/test_periodic_table.py
@@ -0,0 +1,14 @@
+from yt.testing import *
+from yt.utilities.periodic_table import _elements, periodic_table
+
+def test_element_accuracy():
+    for num, w, name, sym in _elements:
+        e0 = periodic_table[num]
+        e1 = periodic_table[name]
+        e2 = periodic_table[sym]
+        yield assert_equal, id(e0), id(e1)
+        yield assert_equal, id(e0), id(e2)
+        yield assert_equal, e0.num, num
+        yield assert_equal, e0.weight, w
+        yield assert_equal, e0.name, name
+        yield assert_equal, e0.symbol, sym


https://bitbucket.org/yt_analysis/yt/commits/dbbc2643cba8/
Changeset:   dbbc2643cba8
Branch:      yt-3.0
User:        MatthewTurk
Date:        2014-03-14 02:36:49
Summary:     Adding chemical formula parser
Affected #:  2 files

diff -r a7257198c1227679c148c8a6aea54e96db8ea8bd -r dbbc2643cba8902d32b146945a4ada361c7e862d yt/utilities/chemical_formulas.py
--- /dev/null
+++ b/yt/utilities/chemical_formulas.py
@@ -0,0 +1,44 @@
+"""
+A simple periodic table.
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# 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 numpy as np
+import string
+import re
+from .periodic_table import periodic_table
+
+class ChemicalFormula:
+    def __init__(self, formula_string):
+        # See YTEP-0003 for information on the format.
+        self.formula_string = formula_string
+        self.elements = []
+        if "_" in self.formula_string:
+            molecule, ionization = self.formula_string.split("_")
+            if ionization[0] == "p":
+                charge = int(ionization[1:])
+            elif ionization[0] == "m":
+                charge = -int(ionization[1:])
+            else:
+                raise NotImplementedError
+        else:
+            molecule = self.formula_string
+            charge = 0
+        self.charge = charge
+        for element, count in re.findall(r'([A-Z][a-z]*)(\d*)', molecule):
+            if count == '': count = 1
+            self.elements.append((periodic_table[element], int(count)))
+        self.weight = sum(n * e.weight for e, n in self.elements)
+
+    def __repr__(self):
+        return self.formula_string

diff -r a7257198c1227679c148c8a6aea54e96db8ea8bd -r dbbc2643cba8902d32b146945a4ada361c7e862d yt/utilities/tests/test_chemical_formulas.py
--- /dev/null
+++ b/yt/utilities/tests/test_chemical_formulas.py
@@ -0,0 +1,23 @@
+from yt.testing import *
+from yt.utilities.chemical_formulas import ChemicalFormula
+from yt.utilities.periodic_table import periodic_table
+
+_molecules = (
+
+    ("H2O_p1", (("H", 2), ("O", 1)), 1),
+    ("H2O_m1", (("H", 2), ("O", 1)), -1),
+    ("H2O", (("H", 2), ("O", 1)), 0),
+    ("H2SO4", (("H", 2), ("S", 1), ("O", 4)), 0),
+    # Now a harder one
+    ("UuoMtUuq3", (("Uuo", 1), ("Mt", 1), ("Uuq", 3)), 0)
+)
+
+def test_formulas():
+    for formula, components, charge in _molecules:
+        f = ChemicalFormula(formula)
+        w = sum( n * periodic_table[e].weight for e, n in components)
+        yield assert_equal, f.charge, charge
+        yield assert_equal, f.weight, w
+        for (n, c1), (e, c2) in zip(components, f.elements):
+            yield assert_equal, n, e.symbol
+            yield assert_equal, c1, c2


https://bitbucket.org/yt_analysis/yt/commits/a60cb6982818/
Changeset:   a60cb6982818
Branch:      yt-3.0
User:        MatthewTurk
Date:        2014-03-14 17:43:28
Summary:     Updating species fields for Enzo.
Affected #:  2 files

diff -r dbbc2643cba8902d32b146945a4ada361c7e862d -r a60cb6982818ddd2cd82949c4fd4a792cbf44a68 yt/fields/species_fields.py
--- a/yt/fields/species_fields.py
+++ b/yt/fields/species_fields.py
@@ -15,23 +15,56 @@
 
 import numpy as np
 
-from yt.fields.field_info_container import \
-    FieldInfoContainer, \
-    NullFunc, \
-    TranslationFunc, \
-    FieldInfo, \
-    ValidateParameter, \
-    ValidateDataField, \
-    ValidateProperty, \
-    ValidateSpatial, \
-    ValidateGridType
-from yt.fields.particle_fields import \
-    particle_deposition_functions, \
-    particle_vector_functions, \
-    standard_particle_fields
 from yt.utilities.physical_constants import \
     mh, \
     mass_sun_cgs, \
-    periodic_table
+    amu_cgs
 from yt.funcs import *
+from yt.utilities.chemical_formulas import \
+    ChemicalFormula
 
+# See YTEP-0003 for details, but we want to ensure these fields are all
+# populated:
+#
+#   * _mass
+#   * _density
+#   * _fraction
+#   * _number_density
+#
+
+def _create_fraction_func(ftype, species):
+    def _frac(field, data):
+        return data[ftype, "%s_density" % species] \
+             / data[ftype, "density"]
+    return _frac
+
+def _create_mass_func(ftype, species):
+    def _mass(field, data):
+        return data[ftype, "%s_density" % species] \
+             * data["index", "cell_volume"]
+    return _mass
+
+def _create_number_density_func(ftype, species):
+    formula = ChemicalFormula(species)
+    weight = formula.weight # This is in AMU
+    weight *= amu_cgs
+    def _number_density(field, data):
+        return data[ftype, "%s_density" % species] \
+             / amu_cgs
+    return _number_density
+
+def add_species_field_by_density(registry, ftype, species):
+    """
+    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 = "")
+    registry.add_field((ftype, "%s_mass" % species),
+                        function = _create_mass_func(ftype, species),
+                        units = "g")
+    registry.add_field((ftype, "%s_number_density" % species),
+                        function = _create_number_density_func(ftype, species),
+                        units = "cm**-3")

diff -r dbbc2643cba8902d32b146945a4ada361c7e862d -r a60cb6982818ddd2cd82949c4fd4a792cbf44a68 yt/frontends/enzo/fields.py
--- a/yt/frontends/enzo/fields.py
+++ b/yt/frontends/enzo/fields.py
@@ -20,6 +20,8 @@
     FieldInfoContainer
 from yt.units.yt_array import \
     YTArray
+from yt.fields.species_fields import \
+    add_species_field_by_density
 
 from yt.utilities.physical_constants import \
     mh, \
@@ -46,6 +48,19 @@
                 ("HDI", 3.0),
     ])
 
+known_species_names = {
+    'HI'    : 'H',
+    'HII'   : 'H_p1',
+    'HeI'   : 'He',
+    'HeII'  : 'He_p1',
+    'HeIII' : 'He_p2',
+    'H2I'   : 'H2',
+    'H2II'  : 'H2_p1',
+    'HM'    : 'H_m1',
+    'DI'    : 'D',
+    'DII'   : 'D_p1',
+    'HD'    : 'HD'
+}
 
 class EnzoFieldInfo(FieldInfoContainer):
     known_other_fields = (
@@ -115,32 +130,16 @@
         self.add_output_field(("enzo", "%s_Density" % species),
                            take_log=True,
                            units="code_mass/code_length**3")
-        self.alias(("gas", "%s_density" % species),
+        yt_name = known_species_names[species]
+        self.alias(("gas", "%s_density" % yt_name),
                    ("enzo", "%s_Density" % species))
-        def _species_mass(field, data):
-            return data["gas", "%s_density" % species] \
-                 * data["cell_volume"]
-        self.add_field(("gas", "%s_mass" % species),
-                           function=_species_mass,
-                           units = "g")
-        def _species_fraction(field, data):
-            return data["gas", "%s_density" % species] \
-                 / data["gas","density"]
-        self.add_field(("gas", "%s_fraction" % species),
-                           function=_species_fraction,
-                           units = "")
-        def _species_number_density(field, data):
-            return data["gas", "%s_density" % species] \
-                / known_species_masses[species]
-        self.add_field(("gas", "%s_number_density" % species),
-                           function=_species_number_density,
-                           units = "1/cm**3")
+        add_species_field_by_density(self, "gas", yt_name)
 
     def setup_species_fields(self):
         species_names = [fn.rsplit("_Density")[0] for ft, fn in 
                          self.field_list if fn.endswith("_Density")]
         species_names = [sp for sp in species_names
-                         if sp in known_species_masses]
+                         if sp in known_species_names]
         for sp in species_names:
             self.add_species_field(sp)
         def _number_density(_sp_list, masses):


https://bitbucket.org/yt_analysis/yt/commits/8363a2ae3ea2/
Changeset:   8363a2ae3ea2
Branch:      yt-3.0
User:        MatthewTurk
Date:        2014-03-14 17:45:27
Summary:     Updating description for chemical formula module.
Affected #:  1 file

diff -r a60cb6982818ddd2cd82949c4fd4a792cbf44a68 -r 8363a2ae3ea2aaa95ca2922d103d24e976d3fb1b yt/utilities/chemical_formulas.py
--- a/yt/utilities/chemical_formulas.py
+++ b/yt/utilities/chemical_formulas.py
@@ -1,5 +1,5 @@
 """
-A simple periodic table.
+Very basic chemical formula parser.
 
 
 


https://bitbucket.org/yt_analysis/yt/commits/64b10e309640/
Changeset:   64b10e309640
Branch:      yt-3.0
User:        MatthewTurk
Date:        2014-03-14 18:07:46
Summary:     This implements a potentially-controversial change of adding electrons as De
and adding both D and De into the "periodic table".
Affected #:  2 files

diff -r 8363a2ae3ea2aaa95ca2922d103d24e976d3fb1b -r 64b10e30964093603f90b79c83171a0b6fcbb66d yt/frontends/enzo/fields.py
--- a/yt/frontends/enzo/fields.py
+++ b/yt/frontends/enzo/fields.py
@@ -22,9 +22,8 @@
     YTArray
 from yt.fields.species_fields import \
     add_species_field_by_density
-
 from yt.utilities.physical_constants import \
-    mh, \
+    mh, me, mp, \
     mass_sun_cgs
 
 b_units = "code_magnetic"
@@ -36,7 +35,7 @@
   (sp, mh * v) for sp, v in [
                 ("HI", 1.0),
                 ("HII", 1.0),
-                ("Electron", 1.0),
+                ("De", 1.0),
                 ("HeI", 4.0),
                 ("HeII", 4.0),
                 ("HeIII", 4.0),
@@ -49,17 +48,18 @@
     ])
 
 known_species_names = {
-    'HI'    : 'H',
-    'HII'   : 'H_p1',
-    'HeI'   : 'He',
-    'HeII'  : 'He_p1',
-    'HeIII' : 'He_p2',
-    'H2I'   : 'H2',
-    'H2II'  : 'H2_p1',
-    'HM'    : 'H_m1',
-    'DI'    : 'D',
-    'DII'   : 'D_p1',
-    'HD'    : 'HD'
+    'HI'      : 'H',
+    'HII'     : 'H_p1',
+    'HeI'     : 'He',
+    'HeII'    : 'He_p1',
+    'HeIII'   : 'He_p2',
+    'H2I'     : 'H2',
+    'H2II'    : 'H2_p1',
+    'HM'      : 'H_m1',
+    'DI'      : 'D',
+    'DII'     : 'D_p1',
+    'HD'      : 'HD',
+    'Electron': 'De'
 }
 
 class EnzoFieldInfo(FieldInfoContainer):
@@ -85,6 +85,8 @@
         ("PhotoGamma", (ra_units, ["photo_gamma"], None)),
         ("Density", (rho_units, ["density"], None)),
         ("Metal_Density", (rho_units, ["metal_density"], None)),
+        # Note: we do not alias Electron_Density to anything
+        ("Electron_Density", (rho_units, [], None)),
     )
 
     known_particle_fields = (
@@ -140,20 +142,14 @@
                          self.field_list if fn.endswith("_Density")]
         species_names = [sp for sp in species_names
                          if sp in known_species_names]
+        def _electron_density(field, data):
+            return data["Electron_Density"] * (me/mp)
+        self.add_field(("gas", "De_density"),
+                       function = _electron_density,
+                       units = "g/cm**3")
         for sp in species_names:
             self.add_species_field(sp)
-        def _number_density(_sp_list, masses):
-            def _num_dens_func(field, data):
-                num = data.pf.arr(np.zeros_like(data["density"], np.float64),
-                                  "1/cm**3")
-                for sp in _sp_list:
-                    num += data["%s_density" % sp] / masses[sp]
-                return num
-            return _num_dens_func
-        func = _number_density(species_names, known_species_masses)
-        self.add_field(("gas", "number_density"),
-                           function = func,
-                           units = "1 / cm**3")
+
 
     def setup_fluid_fields(self):
         # Now we conditionally load a few other things.

diff -r 8363a2ae3ea2aaa95ca2922d103d24e976d3fb1b -r 64b10e30964093603f90b79c83171a0b6fcbb66d yt/utilities/periodic_table.py
--- a/yt/utilities/periodic_table.py
+++ b/yt/utilities/periodic_table.py
@@ -135,6 +135,9 @@
     (115, 288.0000000000, "Ununpentium", "Uup"),
     (116, 292.0000000000, "Ununhexium", "Uuh"),
     (118, 294.0000000000, "Ununoctium", "Uuo"),
+    # Now some special cases that are *not* elements
+    (-1, 2.014102, "Deuterium", "D"),
+    (-1, 0.00054858, "Electron", "De"),
 )
 
 class Element:


https://bitbucket.org/yt_analysis/yt/commits/5c40f4d50013/
Changeset:   5c40f4d50013
Branch:      yt-3.0
User:        MatthewTurk
Date:        2014-03-16 14:02:51
Summary:     Resolving conflict.
Affected #:  1 file

diff -r 64b10e30964093603f90b79c83171a0b6fcbb66d -r 5c40f4d50013ec830e3c37296196b2754fee685c yt/frontends/enzo/fields.py
--- a/yt/frontends/enzo/fields.py
+++ b/yt/frontends/enzo/fields.py
@@ -85,6 +85,7 @@
         ("PhotoGamma", (ra_units, ["photo_gamma"], None)),
         ("Density", (rho_units, ["density"], None)),
         ("Metal_Density", (rho_units, ["metal_density"], None)),
+        ("SN_Colour", (rho_units, [], None)),
         # Note: we do not alias Electron_Density to anything
         ("Electron_Density", (rho_units, [], None)),
     )


https://bitbucket.org/yt_analysis/yt/commits/732be20f7ed1/
Changeset:   732be20f7ed1
Branch:      yt-3.0
User:        MatthewTurk
Date:        2014-03-16 14:04:33
Summary:     Merging
Affected #:  17 files

diff -r 5c40f4d50013ec830e3c37296196b2754fee685c -r 732be20f7ed12269689d3e5040c7e9e2acbce121 yt/analysis_modules/halo_analysis/halo_callbacks.py
--- a/yt/analysis_modules/halo_analysis/halo_callbacks.py
+++ b/yt/analysis_modules/halo_analysis/halo_callbacks.py
@@ -427,8 +427,12 @@
             return            
     else:
         # take first instance of downward intersection with critical value
-        index = np.where((vod[:-1] >= critical_overdensity) &
-                         (vod[1:] < critical_overdensity))[0][0]
+        intersections = (vod[:-1] >= critical_overdensity) & \
+            (vod[1:] < critical_overdensity)
+        if not intersections.any():
+            halo.quantities.update(vquantities)
+            return            
+        index = np.where(intersections)[0][0]
 
     for field in fields:
         v_prof = profile_data[field][dfilter].to_ndarray()

diff -r 5c40f4d50013ec830e3c37296196b2754fee685c -r 732be20f7ed12269689d3e5040c7e9e2acbce121 yt/analysis_modules/halo_finding/rockstar/rockstar.py
--- a/yt/analysis_modules/halo_finding/rockstar/rockstar.py
+++ b/yt/analysis_modules/halo_finding/rockstar/rockstar.py
@@ -255,7 +255,7 @@
         if particle_mass is None:
             pmass_min, pmass_max = dd.quantities.extrema(
                 (ptype, "particle_mass"), non_zero = True)
-            if pmass_min != pmass_max:
+            if np.abs(pmass_max - pmass_min) / pmass_max > 0.01:
                 raise YTRockstarMultiMassNotSupported(pmass_min, pmass_max,
                     ptype)
             particle_mass = pmass_min

diff -r 5c40f4d50013ec830e3c37296196b2754fee685c -r 732be20f7ed12269689d3e5040c7e9e2acbce121 yt/data_objects/data_containers.py
--- a/yt/data_objects/data_containers.py
+++ b/yt/data_objects/data_containers.py
@@ -381,11 +381,13 @@
         else:
             self.index.save_object(self, name)
 
-    def to_glue(self, fields, label="yt"):
+    def to_glue(self, fields, label="yt", data_collection=None):
         """
         Takes specific *fields* in the container and exports them to
         Glue (http://www.glueviz.org) for interactive
-        analysis. Optionally add a *label*.  
+        analysis. Optionally add a *label*. If you are already within
+        the Glue environment, you can pass a *data_collection* object,
+        otherwise Glue will be started.
         """
         from glue.core import DataCollection, Data
         from glue.core.coordinates import coordinates_from_header
@@ -394,11 +396,14 @@
         gdata = Data(label=label)
         for component_name in fields:
             gdata.add_component(self[component_name], component_name)
-        dc = DataCollection([gdata])
 
-        app = GlueApplication(dc)
-        app.start()
-
+        if data_collection is None:
+            dc = DataCollection([gdata])
+            app = GlueApplication(dc)
+            app.start()
+        else:
+            data_collection.append(gdata)
+        
     def __reduce__(self):
         args = tuple([self.pf._hash(), self._type_name] +
                      [getattr(self, n) for n in self._con_args] +

diff -r 5c40f4d50013ec830e3c37296196b2754fee685c -r 732be20f7ed12269689d3e5040c7e9e2acbce121 yt/data_objects/static_output.py
--- a/yt/data_objects/static_output.py
+++ b/yt/data_objects/static_output.py
@@ -94,6 +94,7 @@
     particle_unions = None
     known_filters = None
     _index_class = None
+    field_units = None
 
     class __metaclass__(type):
         def __init__(cls, name, b, d):
@@ -133,6 +134,7 @@
         self.parameters = {}
         self.known_filters = self.known_filters or {}
         self.particle_unions = self.particle_unions or {}
+        self.field_units = self.field_units or {}
 
         # path stuff
         self.parameter_filename = str(filename)
@@ -345,6 +347,12 @@
         f = self.particle_fields_by_type
         fields = set_intersection([f[s] for s in union
                                    if s in self.particle_types_raw])
+        for field in fields:
+            units = set([])
+            for s in union:
+                units.add(self.field_units.get((s, field), ""))
+            if len(units) == 1:
+                self.field_units[union.name, field] = list(units)[0]
         self.particle_types += (union.name,)
         self.particle_unions[union.name] = union
         fields = [ (union.name, field) for field in fields]

diff -r 5c40f4d50013ec830e3c37296196b2754fee685c -r 732be20f7ed12269689d3e5040c7e9e2acbce121 yt/fields/field_info_container.py
--- a/yt/fields/field_info_container.py
+++ b/yt/fields/field_info_container.py
@@ -90,7 +90,8 @@
                 raise RuntimeError
             if field[0] not in self.pf.particle_types:
                 continue
-            self.add_output_field(field, units = "",
+            self.add_output_field(field, 
+                                  units = self.pf.field_units.get(field, ""),
                                   particle_type = True)
 
     def setup_fluid_aliases(self):

diff -r 5c40f4d50013ec830e3c37296196b2754fee685c -r 732be20f7ed12269689d3e5040c7e9e2acbce121 yt/fields/geometric_fields.py
--- a/yt/fields/geometric_fields.py
+++ b/yt/fields/geometric_fields.py
@@ -140,11 +140,11 @@
     def _cylindrical_r(field, data):
         center = data.get_field_parameter("center")
         normal = data.get_field_parameter("normal")
-        coords = data.pf.arr(obtain_rvec(data), "code_length")
+        coords = obtain_rvec(data)
         coords[0,...] -= center[0]
         coords[1,...] -= center[1]
         coords[2,...] -= center[2]
-        return get_cyl_r(coords, normal).in_cgs()
+        return data.pf.arr(get_cyl_r(coords, normal), "code_length").in_cgs()
 
     registry.add_field(("index", "cylindrical_r"),
              function=_cylindrical_r,

diff -r 5c40f4d50013ec830e3c37296196b2754fee685c -r 732be20f7ed12269689d3e5040c7e9e2acbce121 yt/frontends/halo_catalogs/halo_catalog/io.py
--- a/yt/frontends/halo_catalogs/halo_catalog/io.py
+++ b/yt/frontends/halo_catalogs/halo_catalog/io.py
@@ -82,6 +82,7 @@
                     data_file.file_id, pcount)
         ind = 0
         with h5py.File(data_file.filename, "r") as f:
+            if not f.keys(): return None
             pos = np.empty((pcount, 3), dtype="float64")
             pos = data_file.pf.arr(pos, "code_length")
             dx = np.finfo(f['particle_position_x'].dtype).eps
@@ -113,4 +114,6 @@
     def _identify_fields(self, data_file):
         with h5py.File(data_file.filename, "r") as f:
             fields = [("halos", field) for field in f]
-        return fields
+            units = dict([(("halos", field), 
+                           f[field].attrs["units"]) for field in f])
+        return fields, units

diff -r 5c40f4d50013ec830e3c37296196b2754fee685c -r 732be20f7ed12269689d3e5040c7e9e2acbce121 yt/frontends/halo_catalogs/rockstar/io.py
--- a/yt/frontends/halo_catalogs/rockstar/io.py
+++ b/yt/frontends/halo_catalogs/rockstar/io.py
@@ -123,4 +123,4 @@
     def _identify_fields(self, data_file):
         fields = [("halos", f) for f in halo_dt.fields if
                   "padding" not in f]
-        return fields
+        return fields, {}

diff -r 5c40f4d50013ec830e3c37296196b2754fee685c -r 732be20f7ed12269689d3e5040c7e9e2acbce121 yt/frontends/sph/io.py
--- a/yt/frontends/sph/io.py
+++ b/yt/frontends/sph/io.py
@@ -162,7 +162,7 @@
                 # We'll append it anyway.
                 fields.append((ptype, mname))
         f.close()
-        return fields
+        return fields, {}
 
 class IOHandlerGadgetHDF5(IOHandlerOWLS):
     _dataset_type = "gadget_hdf5"
@@ -345,7 +345,7 @@
                     elif req != ptype:
                         continue
                 field_list.append((ptype, field))
-        return field_list
+        return field_list, {}
 
 class IOHandlerTipsyBinary(BaseIOHandler):
     _dataset_type = "tipsy"
@@ -527,7 +527,7 @@
         return self._field_list
 
     def _identify_fields(self, data_file):
-        return self._field_list
+        return self._field_list, {}
 
     def _calculate_particle_offsets(self, data_file):
         field_offsets = {}
@@ -567,7 +567,7 @@
         f = []
         for ftype, fname in self.pf.parameters["field_list"]:
             f.append((str(ftype), str(fname)))
-        return f
+        return f, {}
 
     def _read_particle_coords(self, chunks, ptf):
         chunks = list(chunks)

diff -r 5c40f4d50013ec830e3c37296196b2754fee685c -r 732be20f7ed12269689d3e5040c7e9e2acbce121 yt/frontends/stream/io.py
--- a/yt/frontends/stream/io.py
+++ b/yt/frontends/stream/io.py
@@ -161,7 +161,7 @@
         return {'io': npart}
 
     def _identify_fields(self, data_file):
-        return self.fields[data_file.filename].keys()
+        return self.fields[data_file.filename].keys(), {}
 
 class IOHandlerStreamHexahedral(BaseIOHandler):
     _dataset_type = "stream_hexahedral"

diff -r 5c40f4d50013ec830e3c37296196b2754fee685c -r 732be20f7ed12269689d3e5040c7e9e2acbce121 yt/geometry/particle_geometry_handler.py
--- a/yt/geometry/particle_geometry_handler.py
+++ b/yt/geometry/particle_geometry_handler.py
@@ -117,8 +117,10 @@
     def _detect_output_fields(self):
         # TODO: Add additional fields
         pfl = []
+        units = {}
         for dom in self.data_files:
-            fl = self.io._identify_fields(dom)
+            fl, _units = self.io._identify_fields(dom)
+            units.update(_units)
             dom._calculate_offsets(fl)
             for f in fl:
                 if f not in pfl: pfl.append(f)
@@ -127,6 +129,7 @@
         pf.particle_types = tuple(set(pt for pt, pf in pfl))
         # This is an attribute that means these particle types *actually*
         # exist.  As in, they are real, in the dataset.
+        pf.field_units.update(units)
         pf.particle_types_raw = pf.particle_types
 
     def _identify_base_chunk(self, dobj):

diff -r 5c40f4d50013ec830e3c37296196b2754fee685c -r 732be20f7ed12269689d3e5040c7e9e2acbce121 yt/utilities/exceptions.py
--- a/yt/utilities/exceptions.py
+++ b/yt/utilities/exceptions.py
@@ -349,12 +349,6 @@
             self.ma)
         return v
 
-class YTFITSHeaderNotUnderstood(YTException):
-    def __str__(self):
-        return "This FITS header is not recognizable in its current form.\n" + \
-                "If you would like to force loading, specify: \n" + \
-                "ignore_unit_names = True"
-
 class YTEmptyProfileData(Exception):
     pass
 

diff -r 5c40f4d50013ec830e3c37296196b2754fee685c -r 732be20f7ed12269689d3e5040c7e9e2acbce121 yt/utilities/fits_image.py
--- a/yt/utilities/fits_image.py
+++ b/yt/utilities/fits_image.py
@@ -24,7 +24,7 @@
 class FITSImageBuffer(HDUList):
 
     def __init__(self, data, fields=None, units="cm",
-                 center=None, scale=None):
+                 center=None, scale=None, wcs=None):
         r""" Initialize a FITSImageBuffer object.
 
         FITSImageBuffer contains a list of FITS ImageHDU instances, and optionally includes
@@ -51,6 +51,8 @@
             Pixel scale in unit *units*. Will be ignored if *data* is
             a FixedResolutionBuffer or a YTCoveringGrid. Must be
             specified otherwise, or if *units* is "deg".
+        wcs : `astropy.wcs.WCS` instance, optional
+            Supply an AstroPy WCS instance to override automatic WCS creation.
 
         Examples
         --------
@@ -128,38 +130,43 @@
 
         proj_type = ["linear"]*self.dimensionality
 
-        if isinstance(img_data, FixedResolutionBuffer) and units != "deg":
-            # FRBs are a special case where we have coordinate
-            # information, so we take advantage of this and
-            # construct the WCS object
-            dx = (img_data.bounds[1]-img_data.bounds[0]).in_units(units)/self.nx
-            dy = (img_data.bounds[3]-img_data.bounds[2]).in_units(units)/self.ny
-            xctr = 0.5*(img_data.bounds[1]+img_data.bounds[0]).in_units(units)
-            yctr = 0.5*(img_data.bounds[3]+img_data.bounds[2]).in_units(units)
-            center = [xctr, yctr]
-        elif isinstance(img_data, YTCoveringGridBase):
-            dx, dy, dz = img_data.dds.in_units(units)
-            center = 0.5*(img_data.left_edge+img_data.right_edge).in_units(units)
-        elif units == "deg" and self.dimensionality == 2:
-            dx = -scale[0]
-            dy = scale[1]
-            proj_type = ["RA---TAN","DEC--TAN"]
+        if wcs is None:
+            if isinstance(img_data, FixedResolutionBuffer) and units != "deg":
+                # FRBs are a special case where we have coordinate
+                # information, so we take advantage of this and
+                # construct the WCS object
+                dx = (img_data.bounds[1]-img_data.bounds[0]).in_units(units)/self.nx
+                dy = (img_data.bounds[3]-img_data.bounds[2]).in_units(units)/self.ny
+                xctr = 0.5*(img_data.bounds[1]+img_data.bounds[0]).in_units(units)
+                yctr = 0.5*(img_data.bounds[3]+img_data.bounds[2]).in_units(units)
+                center = [xctr, yctr]
+            elif isinstance(img_data, YTCoveringGridBase):
+                dx, dy, dz = img_data.dds.in_units(units)
+                center = 0.5*(img_data.left_edge+img_data.right_edge).in_units(units)
+            elif units == "deg" and self.dimensionality == 2:
+                dx = -scale[0]
+                dy = scale[1]
+                proj_type = ["RA---TAN","DEC--TAN"]
+            else:
+                dx = scale[0]
+                dy = scale[1]
+                if self.dimensionality == 3: dz = scale[2]
+            
+            w.wcs.crval = center
+            w.wcs.cunit = [units]*self.dimensionality
+            w.wcs.ctype = proj_type
+        
+            if self.dimensionality == 2:
+                w.wcs.cdelt = [dx,dy]
+            elif self.dimensionality == 3:
+                w.wcs.cdelt = [dx,dy,dz]
+
+            self._set_wcs(w)
+
         else:
-            dx = scale[0]
-            dy = scale[1]
-            if self.dimensionality == 3: dz = scale[2]
-            
-        w.wcs.crval = center
-        w.wcs.cunit = [units]*self.dimensionality
-        w.wcs.ctype = proj_type
-        
-        if self.dimensionality == 2:
-            w.wcs.cdelt = [dx,dy]
-        elif self.dimensionality == 3:
-            w.wcs.cdelt = [dx,dy,dz]
 
-        self._set_wcs(w)
-            
+            self._set_wcs(wcs)
+
     def _set_wcs(self, wcs):
         """
         Set the WCS coordinate information for all images
@@ -212,11 +219,13 @@
         elif self.dimensionality == 3:
             return self.nx, self.ny, self.nz
 
-    def to_glue(self, label="yt"):
+    def to_glue(self, label="yt", data_collection=None):
         """
         Takes the data in the FITSImageBuffer and exports it to
         Glue (http://www.glueviz.org) for interactive
-        analysis. Optionally add a *label*. 
+        analysis. Optionally add a *label*. If you are already within
+        the Glue environment, you can pass a *data_collection* object,
+        otherwise Glue will be started.
         """
         from glue.core import DataCollection, Data
         from glue.core.coordinates import coordinates_from_header
@@ -228,10 +237,22 @@
         image.coords = coordinates_from_header(self.wcs.to_header())
         for k,v in field_dict.items():
             image.add_component(v, k)
-        dc = DataCollection([image])
+        if data_collection is None:
+            dc = DataCollection([image])
+            app = GlueApplication(dc)
+            app.start()
+        else:
+            data_collection.append(image)
 
-        app = GlueApplication(dc)
-        app.start()
+    def to_aplpy(self, **kwargs):
+        """
+        Use APLpy (http://aplpy.github.io) for plotting. Returns an `aplpy.FITSFigure`
+        instance. All keyword arguments are passed to the
+        `aplpy.FITSFigure` constructor.
+        """
+        import aplpy
+        return aplpy.FITSFigure(self, **kwargs)
+
 
         
 

diff -r 5c40f4d50013ec830e3c37296196b2754fee685c -r 732be20f7ed12269689d3e5040c7e9e2acbce121 yt/visualization/color_maps.py
--- a/yt/visualization/color_maps.py
+++ b/yt/visualization/color_maps.py
@@ -98,6 +98,32 @@
 
 add_cmap('black_green', cdict)
 
+# This one is a variant of a colormap commonly
+# used for X-ray observations by Maxim Markevitch
+
+cdict = {'red': ((0.0, 0.0, 0.0),
+                 (0.3, 0.0, 0.0),
+                 (0.352, 0.245, 0.245),
+                 (0.42, 0.5, 0.5),
+                 (0.51, 0.706, 0.706),
+                 (0.613, 0.882, 0.882),
+                 (0.742, 1.0, 1.0),
+                 (1.0, 1.0, 1.0)),
+         'green': ((0.0, 0.0, 0.0),
+                   (0.585, 0.0, 0.0),
+                   (0.613, 0.196, 0.196),
+                   (0.693, 0.48, 0.48),
+                   (0.785, 0.696, 0.696),
+                   (0.885, 0.882, 0.882),
+                   (1.0, 1.0, 1.0)),
+         'blue': ((0.0, 0.0, 0.0),
+                  (0.136, 0.0, 0.0),
+                  (0.136, 0.373, 0.373),
+                  (0.391, 1.0, 1.0),
+                  (1.0, 1.0, 1.0))}
+
+add_cmap("purple_mm", cdict)
+
 # This one comes from
 # http://permalink.gmane.org/gmane.comp.python.matplotlib.devel/10518
 # and is an implementation of http://arxiv.org/abs/1108.5083


https://bitbucket.org/yt_analysis/yt/commits/2e4fa7396887/
Changeset:   2e4fa7396887
Branch:      yt-3.0
User:        MatthewTurk
Date:        2014-03-17 15:55:39
Summary:     Changing "De" to "El" for electrons.
Affected #:  2 files

diff -r 732be20f7ed12269689d3e5040c7e9e2acbce121 -r 2e4fa7396887aa889de8e6a6de1ebc434503a815 yt/frontends/enzo/fields.py
--- a/yt/frontends/enzo/fields.py
+++ b/yt/frontends/enzo/fields.py
@@ -59,7 +59,7 @@
     'DI'      : 'D',
     'DII'     : 'D_p1',
     'HD'      : 'HD',
-    'Electron': 'De'
+    'Electron': 'El'
 }
 
 class EnzoFieldInfo(FieldInfoContainer):
@@ -145,7 +145,7 @@
                          if sp in known_species_names]
         def _electron_density(field, data):
             return data["Electron_Density"] * (me/mp)
-        self.add_field(("gas", "De_density"),
+        self.add_field(("gas", "El_density"),
                        function = _electron_density,
                        units = "g/cm**3")
         for sp in species_names:

diff -r 732be20f7ed12269689d3e5040c7e9e2acbce121 -r 2e4fa7396887aa889de8e6a6de1ebc434503a815 yt/utilities/periodic_table.py
--- a/yt/utilities/periodic_table.py
+++ b/yt/utilities/periodic_table.py
@@ -137,7 +137,7 @@
     (118, 294.0000000000, "Ununoctium", "Uuo"),
     # Now some special cases that are *not* elements
     (-1, 2.014102, "Deuterium", "D"),
-    (-1, 0.00054858, "Electron", "De"),
+    (-1, 0.00054858, "Electron", "El"),
 )
 
 class Element:


https://bitbucket.org/yt_analysis/yt/commits/e67d4515547e/
Changeset:   e67d4515547e
Branch:      yt-3.0
User:        MatthewTurk
Date:        2014-03-17 20:14:02
Summary:     Removing unused known_species_masses dict.
Affected #:  1 file

diff -r 2e4fa7396887aa889de8e6a6de1ebc434503a815 -r e67d4515547e2dddd0188571a10617d3be6c51f2 yt/frontends/enzo/fields.py
--- a/yt/frontends/enzo/fields.py
+++ b/yt/frontends/enzo/fields.py
@@ -31,22 +31,6 @@
 rho_units = "code_mass / code_length**3"
 vel_units = "code_velocity"
 
-known_species_masses = dict(
-  (sp, mh * v) for sp, v in [
-                ("HI", 1.0),
-                ("HII", 1.0),
-                ("De", 1.0),
-                ("HeI", 4.0),
-                ("HeII", 4.0),
-                ("HeIII", 4.0),
-                ("H2I", 2.0),
-                ("H2II", 2.0),
-                ("HM", 1.0),
-                ("DI", 2.0),
-                ("DII", 2.0),
-                ("HDI", 3.0),
-    ])
-
 known_species_names = {
     'HI'      : 'H',
     'HII'     : 'H_p1',


https://bitbucket.org/yt_analysis/yt/commits/a10dc67c0f35/
Changeset:   a10dc67c0f35
Branch:      yt-3.0
User:        brittonsmith
Date:        2014-03-20 12:45:01
Summary:     Merged in MatthewTurk/yt/yt-3.0 (pull request #719)

Adding some species handler helpers
Affected #:  6 files

diff -r 794514da3f86f4792efff526fdeeb9556e53b0b5 -r a10dc67c0f352dd88cc4dd8954e1dc1da210fec1 yt/fields/species_fields.py
--- a/yt/fields/species_fields.py
+++ b/yt/fields/species_fields.py
@@ -15,22 +15,56 @@
 
 import numpy as np
 
-from yt.fields.field_info_container import \
-    FieldInfoContainer, \
-    NullFunc, \
-    TranslationFunc, \
-    FieldInfo, \
-    ValidateParameter, \
-    ValidateDataField, \
-    ValidateProperty, \
-    ValidateSpatial, \
-    ValidateGridType
-from yt.fields.particle_fields import \
-    particle_deposition_functions, \
-    particle_vector_functions, \
-    standard_particle_fields
 from yt.utilities.physical_constants import \
     mh, \
-    mass_sun_cgs
+    mass_sun_cgs, \
+    amu_cgs
 from yt.funcs import *
+from yt.utilities.chemical_formulas import \
+    ChemicalFormula
 
+# See YTEP-0003 for details, but we want to ensure these fields are all
+# populated:
+#
+#   * _mass
+#   * _density
+#   * _fraction
+#   * _number_density
+#
+
+def _create_fraction_func(ftype, species):
+    def _frac(field, data):
+        return data[ftype, "%s_density" % species] \
+             / data[ftype, "density"]
+    return _frac
+
+def _create_mass_func(ftype, species):
+    def _mass(field, data):
+        return data[ftype, "%s_density" % species] \
+             * data["index", "cell_volume"]
+    return _mass
+
+def _create_number_density_func(ftype, species):
+    formula = ChemicalFormula(species)
+    weight = formula.weight # This is in AMU
+    weight *= amu_cgs
+    def _number_density(field, data):
+        return data[ftype, "%s_density" % species] \
+             / amu_cgs
+    return _number_density
+
+def add_species_field_by_density(registry, ftype, species):
+    """
+    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 = "")
+    registry.add_field((ftype, "%s_mass" % species),
+                        function = _create_mass_func(ftype, species),
+                        units = "g")
+    registry.add_field((ftype, "%s_number_density" % species),
+                        function = _create_number_density_func(ftype, species),
+                        units = "cm**-3")

diff -r 794514da3f86f4792efff526fdeeb9556e53b0b5 -r a10dc67c0f352dd88cc4dd8954e1dc1da210fec1 yt/frontends/enzo/fields.py
--- a/yt/frontends/enzo/fields.py
+++ b/yt/frontends/enzo/fields.py
@@ -20,9 +20,10 @@
     FieldInfoContainer
 from yt.units.yt_array import \
     YTArray
-
+from yt.fields.species_fields import \
+    add_species_field_by_density
 from yt.utilities.physical_constants import \
-    mh, \
+    mh, me, mp, \
     mass_sun_cgs
 
 b_units = "code_magnetic"
@@ -30,22 +31,20 @@
 rho_units = "code_mass / code_length**3"
 vel_units = "code_velocity"
 
-known_species_masses = dict(
-  (sp, mh * v) for sp, v in [
-                ("HI", 1.0),
-                ("HII", 1.0),
-                ("Electron", 1.0),
-                ("HeI", 4.0),
-                ("HeII", 4.0),
-                ("HeIII", 4.0),
-                ("H2I", 2.0),
-                ("H2II", 2.0),
-                ("HM", 1.0),
-                ("DI", 2.0),
-                ("DII", 2.0),
-                ("HDI", 3.0),
-    ])
-
+known_species_names = {
+    'HI'      : 'H',
+    'HII'     : 'H_p1',
+    'HeI'     : 'He',
+    'HeII'    : 'He_p1',
+    'HeIII'   : 'He_p2',
+    'H2I'     : 'H2',
+    'H2II'    : 'H2_p1',
+    'HM'      : 'H_m1',
+    'DI'      : 'D',
+    'DII'     : 'D_p1',
+    'HD'      : 'HD',
+    'Electron': 'El'
+}
 
 class EnzoFieldInfo(FieldInfoContainer):
     known_other_fields = (
@@ -72,6 +71,8 @@
         ("Density", (rho_units, ["density"], None)),
         ("Metal_Density", (rho_units, ["metal_density"], None)),
         ("SN_Colour", (rho_units, [], None)),
+        # Note: we do not alias Electron_Density to anything
+        ("Electron_Density", (rho_units, [], None)),
     )
 
     known_particle_fields = (
@@ -117,46 +118,24 @@
         self.add_output_field(("enzo", "%s_Density" % species),
                            take_log=True,
                            units="code_mass/code_length**3")
-        self.alias(("gas", "%s_density" % species),
+        yt_name = known_species_names[species]
+        self.alias(("gas", "%s_density" % yt_name),
                    ("enzo", "%s_Density" % species))
-        def _species_mass(field, data):
-            return data["gas", "%s_density" % species] \
-                 * data["cell_volume"]
-        self.add_field(("gas", "%s_mass" % species),
-                           function=_species_mass,
-                           units = "g")
-        def _species_fraction(field, data):
-            return data["gas", "%s_density" % species] \
-                 / data["gas","density"]
-        self.add_field(("gas", "%s_fraction" % species),
-                           function=_species_fraction,
-                           units = "")
-        def _species_number_density(field, data):
-            return data["gas", "%s_density" % species] \
-                / known_species_masses[species]
-        self.add_field(("gas", "%s_number_density" % species),
-                           function=_species_number_density,
-                           units = "1/cm**3")
+        add_species_field_by_density(self, "gas", yt_name)
 
     def setup_species_fields(self):
         species_names = [fn.rsplit("_Density")[0] for ft, fn in 
                          self.field_list if fn.endswith("_Density")]
         species_names = [sp for sp in species_names
-                         if sp in known_species_masses]
+                         if sp in known_species_names]
+        def _electron_density(field, data):
+            return data["Electron_Density"] * (me/mp)
+        self.add_field(("gas", "El_density"),
+                       function = _electron_density,
+                       units = "g/cm**3")
         for sp in species_names:
             self.add_species_field(sp)
-        def _number_density(_sp_list, masses):
-            def _num_dens_func(field, data):
-                num = data.pf.arr(np.zeros_like(data["density"], np.float64),
-                                  "1/cm**3")
-                for sp in _sp_list:
-                    num += data["%s_density" % sp] / masses[sp]
-                return num
-            return _num_dens_func
-        func = _number_density(species_names, known_species_masses)
-        self.add_field(("gas", "number_density"),
-                           function = func,
-                           units = "1 / cm**3")
+
 
     def setup_fluid_fields(self):
         # Now we conditionally load a few other things.

diff -r 794514da3f86f4792efff526fdeeb9556e53b0b5 -r a10dc67c0f352dd88cc4dd8954e1dc1da210fec1 yt/utilities/chemical_formulas.py
--- /dev/null
+++ b/yt/utilities/chemical_formulas.py
@@ -0,0 +1,44 @@
+"""
+Very basic chemical formula parser.
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# 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 numpy as np
+import string
+import re
+from .periodic_table import periodic_table
+
+class ChemicalFormula:
+    def __init__(self, formula_string):
+        # See YTEP-0003 for information on the format.
+        self.formula_string = formula_string
+        self.elements = []
+        if "_" in self.formula_string:
+            molecule, ionization = self.formula_string.split("_")
+            if ionization[0] == "p":
+                charge = int(ionization[1:])
+            elif ionization[0] == "m":
+                charge = -int(ionization[1:])
+            else:
+                raise NotImplementedError
+        else:
+            molecule = self.formula_string
+            charge = 0
+        self.charge = charge
+        for element, count in re.findall(r'([A-Z][a-z]*)(\d*)', molecule):
+            if count == '': count = 1
+            self.elements.append((periodic_table[element], int(count)))
+        self.weight = sum(n * e.weight for e, n in self.elements)
+
+    def __repr__(self):
+        return self.formula_string

diff -r 794514da3f86f4792efff526fdeeb9556e53b0b5 -r a10dc67c0f352dd88cc4dd8954e1dc1da210fec1 yt/utilities/periodic_table.py
--- /dev/null
+++ b/yt/utilities/periodic_table.py
@@ -0,0 +1,178 @@
+"""
+A simple periodic table.
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# 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 numpy as np
+import numbers
+import types
+
+_elements = (
+    (1, 1.0079400000, "Hydrogen", "H"),
+    (2, 4.0026020000, "Helium", "He"),
+    (3, 6.9410000000, "Lithium", "Li"),
+    (4, 9.0121820000, "Beryllium", "Be"),
+    (5, 10.8110000000, "Boron", "B"),
+    (6, 12.0107000000, "Carbon", "C"),
+    (7, 14.0067000000, "Nitrogen", "N"),
+    (8, 15.9994000000, "Oxygen", "O"),
+    (9, 18.9994000000, "Fluorine", "F"),
+    (10, 20.1797000000, "Neon", "Ne"),
+    (11, 22.9897692800, "Sodium", "Na"),
+    (12, 24.3050000000, "Magnesium", "Mg"),
+    (13, 26.9815386000, "Aluminium", "Al"),
+    (14, 28.0855000000, "Silicon", "Si"),
+    (15, 30.9737620000, "Phosphorus", "P"),
+    (16, 32.0650000000, "Sulphur", "S"),
+    (17, 35.4530000000, "Chlorine", "Cl"),
+    (18, 39.9480000000, "Argon", "Ar"),
+    (19, 39.0983000000, "Potassium", "K"),
+    (20, 40.0780000000, "Calcium", "Ca"),
+    (21, 44.9559120000, "Scandium", "Sc"),
+    (22, 47.8670000000, "Titanium", "Ti"),
+    (23, 50.9415000000, "Vanadium", "V"),
+    (24, 51.9961000000, "Chromium", "Cr"),
+    (25, 54.9380450000, "Manganese", "Mn"),
+    (26, 55.8450000000, "Iron", "Fe"),
+    (27, 58.9331950000, "Cobalt", "Co"),
+    (28, 58.6934000000, "Nickel", "Ni"),
+    (29, 63.5460000000, "Copper", "Cu"),
+    (30, 65.3800000000, "Zinc", "Zn"),
+    (31, 69.7230000000, "Gallium", "Ga"),
+    (32, 72.6400000000, "Germanium", "Ge"),
+    (33, 74.9216000000, "Arsenic", "As"),
+    (34, 78.9600000000, "Selenium", "Se"),
+    (35, 79.9040000000, "Bromine", "Br"),
+    (36, 83.7980000000, "Krypton", "Kr"),
+    (37, 85.4678000000, "Rubidium", "Rb"),
+    (38, 87.6200000000, "Strontium", "Sr"),
+    (39, 88.9058500000, "Yttrium", "Y"),
+    (40, 91.2240000000, "Zirkonium", "Zr"),
+    (41, 92.9063800000, "Niobium", "Nb"),
+    (42, 95.9600000000, "Molybdaenum", "Mo"),
+    (43, 98.0000000000, "Technetium", "Tc"),
+    (44, 101.0700000000, "Ruthenium", "Ru"),
+    (45, 102.9055000000, "Rhodium", "Rh"),
+    (46, 106.4200000000, "Palladium", "Pd"),
+    (47, 107.8682000000, "Silver", "Ag"),
+    (48, 112.4110000000, "Cadmium", "Cd"),
+    (49, 114.8180000000, "Indium", "In"),
+    (50, 118.7100000000, "Tin", "Sn"),
+    (51, 121.7600000000, "Antimony", "Sb"),
+    (52, 127.6000000000, "Tellurium", "Te"),
+    (53, 126.9044700000, "Iodine", "I"),
+    (54, 131.2930000000, "Xenon", "Xe"),
+    (55, 132.9054519000, "Cesium", "Cs"),
+    (56, 137.3270000000, "Barium", "Ba"),
+    (57, 138.9054700000, "Lanthanum", "La"),
+    (58, 140.1160000000, "Cerium", "Ce"),
+    (59, 140.9076500000, "Praseodymium", "Pr"),
+    (60, 144.2420000000, "Neodymium", "Nd"),
+    (61, 145.0000000000, "Promethium", "Pm"),
+    (62, 150.3600000000, "Samarium", "Sm"),
+    (63, 151.9640000000, "Europium", "Eu"),
+    (64, 157.2500000000, "Gadolinium", "Gd"),
+    (65, 158.9253500000, "Terbium", "Tb"),
+    (66, 162.5001000000, "Dysprosium", "Dy"),
+    (67, 164.9303200000, "Holmium", "Ho"),
+    (68, 167.2590000000, "Erbium", "Er"),
+    (69, 168.9342100000, "Thulium", "Tm"),
+    (70, 173.0540000000, "Ytterbium", "Yb"),
+    (71, 174.9668000000, "Lutetium", "Lu"),
+    (72, 178.4900000000, "Hafnium", "Hf"),
+    (73, 180.9478800000, "Tantalum", "Ta"),
+    (74, 183.8400000000, "Tungsten", "W"),
+    (75, 186.2070000000, "Rhenium", "Re"),
+    (76, 190.2300000000, "Osmium", "Os"),
+    (77, 192.2170000000, "Iridium", "Ir"),
+    (78, 192.0840000000, "Platinum", "Pt"),
+    (79, 196.9665690000, "Gold", "Au"),
+    (80, 200.5900000000, "Hydrargyrum", "Hg"),
+    (81, 204.3833000000, "Thallium", "Tl"),
+    (82, 207.2000000000, "Lead", "Pb"),
+    (83, 208.9804010000, "Bismuth", "Bi"),
+    (84, 210.0000000000, "Polonium", "Po"),
+    (85, 210.0000000000, "Astatine", "At"),
+    (86, 220.0000000000, "Radon", "Rn"),
+    (87, 223.0000000000, "Francium", "Fr"),
+    (88, 226.0000000000, "Radium", "Ra"),
+    (89, 227.0000000000, "Actinium", "Ac"),
+    (90, 232.0380600000, "Thorium", "Th"),
+    (91, 231.0358800000, "Protactinium", "Pa"),
+    (92, 238.0289100000, "Uranium", "U"),
+    (93, 237.0000000000, "Neptunium", "Np"),
+    (94, 244.0000000000, "Plutonium", "Pu"),
+    (95, 243.0000000000, "Americium", "Am"),
+    (96, 247.0000000000, "Curium", "Cm"),
+    (97, 247.0000000000, "Berkelium", "Bk"),
+    (98, 251.0000000000, "Californium", "Cf"),
+    (99, 252.0000000000, "Einsteinium", "Es"),
+    (100, 257.0000000000, "Fermium", "Fm"),
+    (101, 258.0000000000, "Mendelevium", "Md"),
+    (102, 259.0000000000, "Nobelium", "No"),
+    (103, 262.0000000000, "Lawrencium", "Lr"),
+    (104, 261.0000000000, "Rutherfordium", "Rf"),
+    (105, 262.0000000000, "Dubnium", "Db"),
+    (106, 266.0000000000, "Seaborgium", "Sg"),
+    (107, 264.0000000000, "Bohrium", "Bh"),
+    (108, 277.0000000000, "Hassium", "Hs"),
+    (109, 268.0000000000, "Meitnerium", "Mt"),
+    (110, 271.0000000000, "Ununnilium", "Ds"),
+    (111, 272.0000000000, "Unununium", "Rg"),
+    (112, 285.0000000000, "Ununbium", "Uub"),
+    (113, 284.0000000000, "Ununtrium", "Uut"),
+    (114, 289.0000000000, "Ununquadium", "Uuq"),
+    (115, 288.0000000000, "Ununpentium", "Uup"),
+    (116, 292.0000000000, "Ununhexium", "Uuh"),
+    (118, 294.0000000000, "Ununoctium", "Uuo"),
+    # Now some special cases that are *not* elements
+    (-1, 2.014102, "Deuterium", "D"),
+    (-1, 0.00054858, "Electron", "El"),
+)
+
+class Element:
+    def __init__(self, num, weight, name, symbol):
+        self.num = num
+        self.weight = weight
+        self.name = name
+        self.symbol = symbol
+
+    def __repr__(self):
+        return "Element: %s (%s)" % (self.symbol, self.name)
+
+class PeriodicTable:
+    def __init__(self):
+        self.elements_by_number = {}
+        self.elements_by_name = {}
+        self.elements_by_symbol = {}
+        for num, weight, name, symbol in _elements:
+            e = Element(num, weight, name, symbol)
+            self.elements_by_number[num] = e
+            self.elements_by_name[name] = e
+            self.elements_by_symbol[symbol] = e
+
+    def __getitem__(self, key):
+        if isinstance(key, (np.number, numbers.Number)):
+            d = self.elements_by_number
+        elif isinstance(key, types.StringTypes):
+            if len(key) <= 2:
+                d = self.elements_by_symbol
+            elif len(key) == 3 and key[0] == "U":
+                d = self.elements_by_symbol
+            else:
+                d = self.elements_by_name
+        else:
+            raise KeyError(key)
+        return d[key]
+
+periodic_table = PeriodicTable()

diff -r 794514da3f86f4792efff526fdeeb9556e53b0b5 -r a10dc67c0f352dd88cc4dd8954e1dc1da210fec1 yt/utilities/tests/test_chemical_formulas.py
--- /dev/null
+++ b/yt/utilities/tests/test_chemical_formulas.py
@@ -0,0 +1,23 @@
+from yt.testing import *
+from yt.utilities.chemical_formulas import ChemicalFormula
+from yt.utilities.periodic_table import periodic_table
+
+_molecules = (
+
+    ("H2O_p1", (("H", 2), ("O", 1)), 1),
+    ("H2O_m1", (("H", 2), ("O", 1)), -1),
+    ("H2O", (("H", 2), ("O", 1)), 0),
+    ("H2SO4", (("H", 2), ("S", 1), ("O", 4)), 0),
+    # Now a harder one
+    ("UuoMtUuq3", (("Uuo", 1), ("Mt", 1), ("Uuq", 3)), 0)
+)
+
+def test_formulas():
+    for formula, components, charge in _molecules:
+        f = ChemicalFormula(formula)
+        w = sum( n * periodic_table[e].weight for e, n in components)
+        yield assert_equal, f.charge, charge
+        yield assert_equal, f.weight, w
+        for (n, c1), (e, c2) in zip(components, f.elements):
+            yield assert_equal, n, e.symbol
+            yield assert_equal, c1, c2

diff -r 794514da3f86f4792efff526fdeeb9556e53b0b5 -r a10dc67c0f352dd88cc4dd8954e1dc1da210fec1 yt/utilities/tests/test_periodic_table.py
--- /dev/null
+++ b/yt/utilities/tests/test_periodic_table.py
@@ -0,0 +1,14 @@
+from yt.testing import *
+from yt.utilities.periodic_table import _elements, periodic_table
+
+def test_element_accuracy():
+    for num, w, name, sym in _elements:
+        e0 = periodic_table[num]
+        e1 = periodic_table[name]
+        e2 = periodic_table[sym]
+        yield assert_equal, id(e0), id(e1)
+        yield assert_equal, id(e0), id(e2)
+        yield assert_equal, e0.num, num
+        yield assert_equal, e0.weight, w
+        yield assert_equal, e0.name, name
+        yield assert_equal, e0.symbol, sym

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