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

commits-noreply at bitbucket.org commits-noreply at bitbucket.org
Mon Jul 14 11:02:46 PDT 2014


16 new commits in yt:

https://bitbucket.org/yt_analysis/yt/commits/1c42d116ca7b/
Changeset:   1c42d116ca7b
Branch:      yt-3.0
User:        drudd
Date:        2014-07-11 23:30:41
Summary:     Removed find_max from index object
Affected #:  4 files

diff -r c6b09204f3e728a7d3611ffc16a962ad037ddbee -r 1c42d116ca7be5d53936b759e31341925169a827 yt/frontends/artio/data_structures.py
--- a/yt/frontends/artio/data_structures.py
+++ b/yt/frontends/artio/data_structures.py
@@ -182,29 +182,6 @@
     def convert(self, unit):
         return self.parameter_file.conversion_factors[unit]
 
-    def find_max(self, field, finest_levels=3):
-        """
-        Returns (value, center) of location of maximum for a given field.
-        """
-        if (field, finest_levels) in self._max_locations:
-            return self._max_locations[(field, finest_levels)]
-        mv, pos = self.find_max_cell_location(field, finest_levels)
-        self._max_locations[(field, finest_levels)] = (mv, pos)
-        return mv, pos
-
-    def find_max_cell_location(self, field, finest_levels=3):
-        source = self.all_data()
-        if finest_levels is not False:
-            source.min_level = self.max_level - finest_levels
-        mylog.debug("Searching for maximum value of %s", field)
-        max_val, maxi, mx, my, mz = \
-            source.quantities["MaxLocation"](field)
-        mylog.info("Max Value is %0.5e at %0.16f %0.16f %0.16f",
-                   max_val, mx, my, mz)
-        self.pf.parameters["Max%sValue" % (field)] = max_val
-        self.pf.parameters["Max%sPos" % (field)] = "%s" % ((mx, my, mz),)
-        return max_val, np.array((mx, my, mz), dtype='float64')
-
     def _detect_output_fields(self):
         self.fluid_field_list = self._detect_fluid_fields()
         self.particle_field_list = self._detect_particle_fields()

diff -r c6b09204f3e728a7d3611ffc16a962ad037ddbee -r 1c42d116ca7be5d53936b759e31341925169a827 yt/frontends/boxlib/data_structures.py
--- a/yt/frontends/boxlib/data_structures.py
+++ b/yt/frontends/boxlib/data_structures.py
@@ -343,7 +343,6 @@
         self._parallel_locking = False
         self._data_file = None
         self._data_mode = None
-        self._max_locations = {}
 
     def _detect_output_fields(self):
         # This is all done in _parse_header_file

diff -r c6b09204f3e728a7d3611ffc16a962ad037ddbee -r 1c42d116ca7be5d53936b759e31341925169a827 yt/geometry/geometry_handler.py
--- a/yt/geometry/geometry_handler.py
+++ b/yt/geometry/geometry_handler.py
@@ -72,7 +72,6 @@
         self._parallel_locking = False
         self._data_file = None
         self._data_mode = None
-        self._max_locations = {}
         self.num_grids = None
 
     def _initialize_data_storage(self):

diff -r c6b09204f3e728a7d3611ffc16a962ad037ddbee -r 1c42d116ca7be5d53936b759e31341925169a827 yt/geometry/grid_geometry_handler.py
--- a/yt/geometry/grid_geometry_handler.py
+++ b/yt/geometry/grid_geometry_handler.py
@@ -201,31 +201,6 @@
         for item in ("Mpc", "pc", "AU", "cm"):
             print "\tWidth: %0.3e %s" % (dx.in_units(item), item)
 
-    def find_max(self, field, finest_levels = 3):
-        """
-        Returns (value, center) of location of maximum for a given field.
-        """
-        if (field, finest_levels) in self._max_locations:
-            return self._max_locations[(field, finest_levels)]
-        mv, pos = self.find_max_cell_location(field, finest_levels)
-        self._max_locations[(field, finest_levels)] = (mv, pos)
-        return mv, pos
-
-    def find_max_cell_location(self, field, finest_levels = 3):
-        if finest_levels is not False:
-            gi = (self.grid_levels >= self.max_level - finest_levels).ravel()
-            source = self.data_collection([0.0]*3, self.grids[gi])
-        else:
-            source = self.all_data()
-        mylog.debug("Searching for maximum value of %s", field)
-        max_val, maxi, mx, my, mz = \
-            source.quantities["MaxLocation"](field)
-        mylog.info("Max Value is %0.5e at %0.16f %0.16f %0.16f", 
-              max_val, mx, my, mz)
-        self.parameters["Max%sValue" % (field)] = max_val
-        self.parameters["Max%sPos" % (field)] = "%s" % ((mx,my,mz),)
-        return max_val, np.array((mx,my,mz), dtype='float64')
-
     def find_points(self, x, y, z) :
         """
         Returns the (objects, indices) of leaf grids containing a number of (x,y,z) points
@@ -241,40 +216,6 @@
         ind = pts.find_points_in_tree()
         return self.grids[ind], ind
 
-    def find_field_value_at_point(self, fields, coord):
-        r"""Find the value of fields at a coordinate.
-
-        Returns the values [field1, field2,...] of the fields at the given
-        (x, y, z) points. Returns a list of field values in the same order as
-        the input *fields*.
-
-        Parameters
-        ----------
-        fields : string or list of strings
-            The field(s) that will be returned.
-
-        coord : list or array of coordinates
-            The location for which field values will be returned.
-
-        Examples
-        --------
-        >>> pf.h.find_field_value_at_point(['Density', 'Temperature'],
-            [0.4, 0.3, 0.8])
-        [2.1489e-24, 1.23843e4]
-        """
-        this = self.find_points(*coord)[0][-1]
-        cellwidth = (this.RightEdge - this.LeftEdge) / this.ActiveDimensions
-        mark = np.zeros(3).astype('int')
-        # Find the index for the cell containing this point.
-        for dim in xrange(len(coord)):
-            mark[dim] = int((coord[dim] - this.LeftEdge[dim]) / cellwidth[dim])
-        out = []
-        fields = ensure_list(fields)
-        # Pull out the values and add it to the out list.
-        for field in fields:
-            out.append(this[field][mark[0], mark[1], mark[2]])
-        return out
-
     def get_grid_tree(self) :
 
         left_edge = self.pf.arr(np.zeros((self.num_grids, 3)),


https://bitbucket.org/yt_analysis/yt/commits/cd70d7a92798/
Changeset:   cd70d7a92798
Branch:      yt-3.0
User:        drudd
Date:        2014-07-11 23:31:41
Summary:     Removed find_max from oct_geometry_handler
Affected #:  1 file

diff -r 1c42d116ca7be5d53936b759e31341925169a827 -r cd70d7a927981945983947158280042be8b0c5ac yt/geometry/oct_geometry_handler.py
--- a/yt/geometry/oct_geometry_handler.py
+++ b/yt/geometry/oct_geometry_handler.py
@@ -49,27 +49,3 @@
 
     def convert(self, unit):
         return self.parameter_file.conversion_factors[unit]
-
-    def find_max(self, field, finest_levels = 3):
-        """
-        Returns (value, center) of location of maximum for a given field.
-        """
-        if (field, finest_levels) in self._max_locations:
-            return self._max_locations[(field, finest_levels)]
-        mv, pos = self.find_max_cell_location(field, finest_levels)
-        self._max_locations[(field, finest_levels)] = (mv, pos)
-        return mv, pos
-
-    def find_max_cell_location(self, field, finest_levels = 3):
-        source = self.all_data()
-        if finest_levels is not False:
-            source.min_level = self.max_level - finest_levels
-        mylog.debug("Searching for maximum value of %s", field)
-        max_val, maxi, mx, my, mz = \
-            source.quantities["MaxLocation"](field)
-        mylog.info("Max Value is %0.5e at %0.16f %0.16f %0.16f", 
-              max_val, mx, my, mz)
-        self.pf.parameters["Max%sValue" % (field,)] = max_val
-        self.pf.parameters["Max%sPos" % (field,)] = "%s" % ((mx,my,mz),)
-        return max_val, np.array((mx,my,mz), dtype='float64')
-


https://bitbucket.org/yt_analysis/yt/commits/8c4978612bab/
Changeset:   8c4978612bab
Branch:      yt-3.0
User:        drudd
Date:        2014-07-11 23:37:44
Summary:     Added YTPoint, and supporting PointSelector, YTSelectionContainer0D
Affected #:  5 files

diff -r cd70d7a927981945983947158280042be8b0c5ac -r 8c4978612babf319b248b5aec62d75e0ae520513 yt/data_objects/data_containers.py
--- a/yt/data_objects/data_containers.py
+++ b/yt/data_objects/data_containers.py
@@ -725,6 +725,12 @@
             self.index._identify_base_chunk(self)
         return self._current_chunk.fwidth
 
+class YTSelectionContainer0D(YTSelectionContainer):
+    _spatial = False
+    def __init__(self, pf, field_parameters):
+        super(YTSelectionContainer0D, self).__init__(
+            pf, field_parameters)
+
 class YTSelectionContainer1D(YTSelectionContainer):
     _spatial = False
     def __init__(self, pf, field_parameters):

diff -r cd70d7a927981945983947158280042be8b0c5ac -r 8c4978612babf319b248b5aec62d75e0ae520513 yt/data_objects/selection_data_containers.py
--- a/yt/data_objects/selection_data_containers.py
+++ b/yt/data_objects/selection_data_containers.py
@@ -21,7 +21,8 @@
 from yt.utilities.lib.alt_ray_tracers import cylindrical_ray_trace
 from yt.utilities.orientation import Orientation
 from .data_containers import \
-    YTSelectionContainer1D, YTSelectionContainer2D, YTSelectionContainer3D
+    YTSelectionContainer0D, YTSelectionContainer1D, \
+    YTSelectionContainer2D, YTSelectionContainer3D
 from yt.data_objects.derived_quantities import \
     DerivedQuantityCollection
 from yt.utilities.exceptions import \
@@ -34,6 +35,31 @@
 from yt.utilities.math_utils import get_rotation_matrix
 from yt.units.yt_array import YTQuantity
 
+
+class YTPointBase(YTSelectionContainer0D):
+    """
+    A 0-dimensional object defined by a single point
+
+    Parameters
+    ----------
+    p: array_like
+        A points defined within the domain.  If the domain is
+        periodic its position will be corrected to lie inside
+        the range [DLE,DRE) to ensure one and only one cell may
+        match that point
+
+    Examples
+    --------
+    >>> pf = load("DD0010/moving7_0010")
+    >>> c = [0.5,0.5,0.5]
+    >>> point = pf.point(c)
+    """
+    _type_name = "point"
+    _con_args = ('p',)
+    def __init__(self, p, pf = None, field_parameters = None):
+        super(YTPointBase, self).__init__(pf, field_parameters)
+        self.p = p
+
 class YTOrthoRayBase(YTSelectionContainer1D):
     """
     This is an orthogonal ray cast through the entire domain, at a specific
@@ -530,7 +556,7 @@
 
 class YTSphereBase(YTSelectionContainer3D):
     """
-    A sphere f points defined by a *center* and a *radius*.
+    A sphere of points defined by a *center* and a *radius*.
 
     Parameters
     ----------

diff -r cd70d7a927981945983947158280042be8b0c5ac -r 8c4978612babf319b248b5aec62d75e0ae520513 yt/data_objects/tests/test_points.py
--- /dev/null
+++ b/yt/data_objects/tests/test_points.py
@@ -0,0 +1,10 @@
+from yt.testing import *
+import numpy as np
+
+def setup():
+    from yt.config import ytcfg
+    ytcfg["yt","__withintesting"] = "True"
+
+def test_domain_point():
+    pf = fake_random_pf(16, fields = ("density"))
+    p = pf.point(pf.domain_center)

diff -r cd70d7a927981945983947158280042be8b0c5ac -r 8c4978612babf319b248b5aec62d75e0ae520513 yt/geometry/selection_routines.pyx
--- a/yt/geometry/selection_routines.pyx
+++ b/yt/geometry/selection_routines.pyx
@@ -534,6 +534,63 @@
                 self.domain_width[1],
                 self.domain_width[2])
 
+
+cdef class PointSelector(SelectorObject):
+    cdef np.float64_t p[3]
+
+    def __init__(self, dobj):
+        for i in range(3):
+            self.p[i] = _ensure_code(dobj.p[i])
+
+            # ensure the point lies in the domain
+            if self.periodicity[i]:
+                self.p[i] = np.fmod(self.p[i], self.domain_width[i])
+                if self.p[i] < 0.0:
+                    self.p[i] += self.domain_width[i]
+
+    @cython.boundscheck(False)
+    @cython.wraparound(False)
+    @cython.cdivision(True)
+    cdef int select_cell(self, np.float64_t pos[3], np.float64_t dds[3]) nogil:
+        if (pos[0] - 0.5*dds[0] <= self.p[0] < pos[0]+0.5*dds[0] and
+            pos[1] - 0.5*dds[1] <= self.p[1] < pos[1]+0.5*dds[1] and
+            pos[2] - 0.5*dds[2] <= self.p[2] < pos[2]+0.5*dds[2]):
+            return 1
+        else:
+            return 0
+
+    @cython.boundscheck(False)
+    @cython.wraparound(False)
+    @cython.cdivision(True)
+    cdef int select_sphere(self, np.float64_t pos[3], np.float64_t radius) nogil:
+        cdef int i
+        cdef np.float64_t dist, dist2 = 0
+        for i in range(3):
+            dist = self.difference(pos[i], self.p[i], i)
+            dist2 += dist*dist
+        if dist2 <= radius*radius: return 1
+        return 0
+
+    @cython.boundscheck(False)
+    @cython.wraparound(False)
+    @cython.cdivision(True)
+    cdef int select_bbox(self, np.float64_t left_edge[3],
+                               np.float64_t right_edge[3]) nogil:
+        cdef int i
+        # point definitely can only be in one cell
+        if (left_edge[0] <= self.p[0] < right_edge[0] and
+            left_edge[1] <= self.p[1] < right_edge[1] and
+            left_edge[2] <= self.p[2] < right_edge[2]):
+            return 1
+        else:
+            return 0
+
+    def _hash_vals(self):
+        return (self.p[0], self.p[1], self.p[2])
+
+point_selector = PointSelector
+
+
 cdef class SphereSelector(SelectorObject):
     cdef np.float64_t radius
     cdef np.float64_t radius2

diff -r cd70d7a927981945983947158280042be8b0c5ac -r 8c4978612babf319b248b5aec62d75e0ae520513 yt/utilities/tests/test_selectors.py
--- a/yt/utilities/tests/test_selectors.py
+++ b/yt/utilities/tests/test_selectors.py
@@ -1,7 +1,7 @@
 import numpy as np
 from yt.testing import \
     fake_random_pf, assert_equal, assert_array_less, \
-    YTArray
+    YTArray, fake_amr_pf
 from yt.utilities.math_utils import periodic_dist
 
 
@@ -9,6 +9,26 @@
     from yt.config import ytcfg
     ytcfg["yt", "__withintesting"] = "True"
 
+def test_point_selector():
+    # generate fake amr data
+    pf = fake_random_pf(64, nprocs=51)
+    assert(all(pf.periodicity))
+
+    dd = pf.h.all_data()
+    positions = np.array([dd[ax] for ax in 'xyz'])
+    delta = 0.5*np.array([dd['d'+ax] for ax in 'xyz'])
+    # ensure cell centers and corners always return one and 
+    # only one point object
+    for p in positions:
+        data = pf.point(p)
+        assert_equal(data["ones"].shape[0], 1) 
+    for p in positions - delta:
+        data = pf.point(p)
+        assert_equal(data["ones"].shape[0], 1)
+    for p in positions + delta:
+        data = pf.point(p)
+        assert_equal(data["ones"].shape[0], 1)
+ 
 def test_sphere_selector():
     # generate fake data with a number of non-cubical grids
     pf = fake_random_pf(64, nprocs=51)


https://bitbucket.org/yt_analysis/yt/commits/198f24efcdf5/
Changeset:   198f24efcdf5
Branch:      yt-3.0
User:        drudd
Date:        2014-07-12 00:13:24
Summary:     Removed test_find_points testcase (assuming find_points is okay to remove from index
Affected #:  1 file

diff -r 8c4978612babf319b248b5aec62d75e0ae520513 -r 198f24efcdf51c095bcc790525f59d4b6f4495f2 yt/utilities/lib/tests/test_grid_tree.py
--- a/yt/utilities/lib/tests/test_grid_tree.py
+++ b/yt/utilities/lib/tests/test_grid_tree.py
@@ -67,53 +67,3 @@
             grid_children = np.array([child.id - child._id_offset
                                       for child in grid.Children])
             yield assert_equal, grid_children, children[i]
-
-
-def test_find_points():
-    """Main test suite for MatchPoints"""
-    num_points = 100
-    randx = np.random.uniform(low=test_pf.domain_left_edge[0],
-                              high=test_pf.domain_right_edge[0],
-                              size=num_points)
-    randy = np.random.uniform(low=test_pf.domain_left_edge[1],
-                              high=test_pf.domain_right_edge[1],
-                              size=num_points)
-    randz = np.random.uniform(low=test_pf.domain_left_edge[2],
-                              high=test_pf.domain_right_edge[2],
-                              size=num_points)
-
-    point_grids, point_grid_inds = test_pf.index.find_points(randx, randy, randz)
-
-    grid_inds = np.zeros((num_points), dtype='int64')
-
-    for ind, ixx, iyy, izz in zip(range(num_points), randx, randy, randz):
-
-        pos = np.array([ixx, iyy, izz])
-        pt_level = -1
-
-        for grid in test_pf.index.grids:
-
-            if np.all(pos >= grid.LeftEdge) and \
-               np.all(pos <= grid.RightEdge) and \
-               grid.Level > pt_level:
-                pt_level = grid.Level
-                grid_inds[ind] = grid.id - grid._id_offset
-
-    yield assert_equal, point_grid_inds, grid_inds
-
-    # Test wheter find_points works for lists
-    point_grids, point_grid_inds = test_pf.index.find_points(randx.tolist(),
-                                                         randy.tolist(),
-                                                         randz.tolist())
-    yield assert_equal, point_grid_inds, grid_inds
-
-    # Test if find_points works for scalar
-    ind = random.randint(0, num_points - 1)
-    point_grids, point_grid_inds = test_pf.index.find_points(randx[ind],
-                                                         randy[ind],
-                                                         randz[ind])
-    yield assert_equal, point_grid_inds, grid_inds[ind]
-
-    # Test if find_points fails properly for non equal indices' array sizes
-    yield assert_raises, AssertionError, test_pf.index.find_points, \
-        [0], 1.0, [2, 3]


https://bitbucket.org/yt_analysis/yt/commits/84525632d209/
Changeset:   84525632d209
Branch:      yt-3.0
User:        drudd
Date:        2014-07-12 00:27:34
Summary:     Added crude implementation of find_values_at_points
Affected #:  2 files

diff -r 198f24efcdf51c095bcc790525f59d4b6f4495f2 -r 84525632d209088422ef3d90d7bdddc3bba03f30 setup.py
--- a/setup.py
+++ b/setup.py
@@ -193,7 +193,7 @@
 
 class my_build_src(build_src.build_src):
     def run(self):
-        self.run_command("build_forthon")
+        #self.run_command("build_forthon")
         build_src.build_src.run(self)
 
 
@@ -277,8 +277,8 @@
         configuration=configuration,
         zip_safe=False,
         data_files=REASON_FILES,
-        cmdclass={'build_py': my_build_py, 'build_forthon': BuildForthon,
-                  'build_src': my_build_src, 'install_data': my_install_data},
+        cmdclass={'build_py': my_build_py, #'build_forthon': BuildForthon,
+                  'build_src': my_build_src},# 'install_data': my_install_data},
     )
     return
 

diff -r 198f24efcdf51c095bcc790525f59d4b6f4495f2 -r 84525632d209088422ef3d90d7bdddc3bba03f30 yt/data_objects/static_output.py
--- a/yt/data_objects/static_output.py
+++ b/yt/data_objects/static_output.py
@@ -541,6 +541,31 @@
               min_val, mx, my, mz)
         return min_val, self.arr([mx, my, mz], 'code_length', dtype="float64")
 
+    def find_field_values_at_points(self, fields, coords):
+        """
+        Returns the values [field1, field2,...] of the fields at the given
+        [(x1, y1, z2), (x2, y2, z2),...] points.  Returns a list of field 
+        values in the same order as the input *fields*.
+
+        This is probably quite slow right now as it creates a new data object
+        for each point.  We'd probably need a YTPointCollection container to
+        optimize this.
+        """
+        fields = ensure_list(fields)
+        x, y, z = coords
+        x = ensure_list(x)
+        y = ensure_list(y)
+        z = ensure_list(z)
+        if not len(x) == len(y) == len(z):
+            raise AssertionError("Arrays of indices must be of the same size")
+       
+        out = [np.zeros(len(x), dtype=np.float64) for f in fields] 
+        for i in range(len(x)):
+            p = self.point([x[i],y[i],z[i]])
+            for j,field in enumerate(fields):
+                out[j][i] = p[field][0]
+        return out
+
     # Now all the object related stuff
     def all_data(self, find_max=False):
         if find_max: c = self.find_max("density")[1]


https://bitbucket.org/yt_analysis/yt/commits/617c99f73500/
Changeset:   617c99f73500
Branch:      yt-3.0
User:        drudd
Date:        2014-07-12 17:06:04
Summary:     Fixes to ARTIO when fileset lacks particles
Affected #:  1 file

diff -r c6b09204f3e728a7d3611ffc16a962ad037ddbee -r 617c99f73500d5c221d8acaf58f1cdfaee43e0c0 yt/frontends/artio/_artio_caller.pyx
--- a/yt/frontends/artio/_artio_caller.pyx
+++ b/yt/frontends/artio/_artio_caller.pyx
@@ -300,6 +300,10 @@
         cdef int num_fields = len(fields)
         cdef np.float64_t pos[3]
 
+        # since RuntimeErrors are not fatal, ensure no artio_particles* functions
+        # called if fileset lacks particles
+        if not self.has_particles: return
+
         data = {}
         accessed_species = np.zeros( self.num_species, dtype="int")
         selected_mass = [ None for i in range(self.num_species)]
@@ -587,15 +591,25 @@
         self.handle = artio_handle.handle
         self.oct_count = None
         self.root_mesh_data = NULL
-        self.pcount = <np.int64_t **> malloc(sizeof(np.int64_t*)
-            * artio_handle.num_species)
-        self.nvars[0] = artio_handle.num_species
-        self.nvars[1] = artio_handle.num_grid_variables
-        for i in range(artio_handle.num_species):
-            self.pcount[i] = <np.int64_t*> malloc(sizeof(np.int64_t)
-                * (self.sfc_end - self.sfc_start + 1))
-            for sfc in range(self.sfc_end - self.sfc_start + 1):
-                self.pcount[i][sfc] = 0
+        self.pcount = NULL
+
+        if artio_handle.has_particles:
+            self.pcount = <np.int64_t **> malloc(sizeof(np.int64_t*)
+                * artio_handle.num_species)
+            self.nvars[0] = artio_handle.num_species
+            for i in range(artio_handle.num_species):
+                self.pcount[i] = <np.int64_t*> malloc(sizeof(np.int64_t)
+                    * (self.sfc_end - self.sfc_start + 1))
+                for sfc in range(self.sfc_end - self.sfc_start + 1):
+                    self.pcount[i][sfc] = 0
+        else:
+            self.nvars[0] = 0
+
+        if artio_handle.has_grid:
+            self.nvars[1] = artio_handle.num_grid_variables
+        else:
+            self.nvars[1] = 0
+
         for i in range(3):
             self.dims[i] = domain_dimensions[i]
             self.DLE[i] = domain_left_edge[i]
@@ -604,12 +618,15 @@
 
     def __dealloc__(self):
         cdef int i
-        for i in range(self.nvars[0]):
-            free(self.pcount[i])
-        for i in range(self.nvars[1]):
-            free(self.root_mesh_data[i])
-        free(self.pcount)
-        free(self.root_mesh_data)
+        if self.artio_handle.has_particles:
+            for i in range(self.nvars[0]):
+                free(self.pcount[i])
+            free(self.pcount)
+        if self.artio_handle.has_grid:
+            if self.root_mesh_data != NULL:
+                for i in range(self.nvars[1]):
+                    free(self.root_mesh_data[i])
+                free(self.root_mesh_data)
 
     @cython.boundscheck(False)
     @cython.wraparound(False)
@@ -623,8 +640,7 @@
         cdef int *num_octs_per_level = <int *>malloc(
             (max_level + 1)*sizeof(int))
         cdef int num_species = self.artio_handle.num_species
-        cdef int *num_particles_per_species =  <int *>malloc(
-            sizeof(int)*num_species) 
+        cdef int *num_particles_per_species 
         cdef ARTIOOctreeContainer octree
         ngv = self.nvars[1]
         cdef float *grid_variables = <float *>malloc(
@@ -648,10 +664,10 @@
             status = artio_grid_read_root_cell_begin( self.handle,
                 sfc, dpos, grid_variables, &num_oct_levels,
                 num_octs_per_level)
+            check_artio_status(status)
             for i in range(ngv):
                 self.root_mesh_data[i][sfc - self.sfc_start] = \
                     grid_variables[i]
-            check_artio_status(status)
             if num_oct_levels > 0:
                 oc = 0
                 for level in range(num_oct_levels):
@@ -660,33 +676,39 @@
                 oct_count[sfc - self.sfc_start] = oc
                 octree.initialize_local_mesh(oc, num_oct_levels,
                     num_octs_per_level, sfc)
-            status = artio_grid_read_root_cell_end( self.handle )
+            status = artio_grid_read_root_cell_end(self.handle)
             check_artio_status(status)
-        status = artio_grid_clear_sfc_cache( self.handle)
+        status = artio_grid_clear_sfc_cache(self.handle)
         check_artio_status(status)
-        # Now particles
-        status = artio_particle_cache_sfc_range(self.handle, self.sfc_start,
+
+        if self.artio_handle.has_particles:
+            num_particles_per_species =  <int *>malloc(
+                    sizeof(int)*num_species)
+
+            # Now particles
+            status = artio_particle_cache_sfc_range(self.handle, self.sfc_start,
                                             self.sfc_end)
-        check_artio_status(status) 
-        for sfc in range(self.sfc_start, self.sfc_end + 1):
-            # Now particles
-            status = artio_particle_read_root_cell_begin( self.handle,
-                    sfc, num_particles_per_species )
+            check_artio_status(status) 
+            for sfc in range(self.sfc_start, self.sfc_end + 1):
+                # Now particles
+                status = artio_particle_read_root_cell_begin(self.handle,
+                        sfc, num_particles_per_species)
+                check_artio_status(status)
+
+                for i in range(num_species):
+                    self.pcount[i][sfc - self.sfc_start] = \
+                        num_particles_per_species[i]
+
+                status = artio_particle_read_root_cell_end(self.handle)
+                check_artio_status(status)
+
+            status = artio_particle_clear_sfc_cache(self.handle)
             check_artio_status(status)
 
-            for i in range(num_species):
-                self.pcount[i][sfc - self.sfc_start] = \
-                    num_particles_per_species[i]
-
-            status = artio_particle_read_root_cell_end( self.handle)
-            check_artio_status(status)
-
-        status = artio_particle_clear_sfc_cache( self.handle)
-        check_artio_status(status)
+            free(num_particles_per_species)
 
         free(grid_variables)
         free(num_octs_per_level)
-        free(num_particles_per_species)
         self.oct_count = oct_count
         self.doct_count = <np.int64_t *> oct_count.data
         self.root_mesh_handler = ARTIORootMeshContainer(self)
@@ -861,7 +883,7 @@
         cdef np.int64_t sfc_start, sfc_end
         sfc_start = self.domains[0].con_id
         sfc_end = self.domains[self.num_domains - 1].con_id
-        status = artio_grid_cache_sfc_range(handle, sfc_start, sfc_end )
+        status = artio_grid_cache_sfc_range(handle, sfc_start, sfc_end)
         check_artio_status(status) 
         cdef np.int64_t offset = 0 
         for si in range(self.num_domains):
@@ -884,7 +906,7 @@
                 status = artio_grid_read_level_end(handle)
                 check_artio_status(status)
                 lp += num_octs_per_level[level]
-            status = artio_grid_read_root_cell_end( handle )
+            status = artio_grid_read_root_cell_end(handle)
             check_artio_status(status)
             # Now we have all our sources.
             for j in range(nf):
@@ -943,6 +965,10 @@
     cdef np.ndarray[np.int64_t, ndim=1] npi64arr
     cdef np.ndarray[np.float64_t, ndim=1] npf64arr
 
+    if not artio_handle.has_particles:
+        raise RuntimeError("Attempted to read non-existent particles in ARTIO")
+        return
+
     # Now we set up our field pointers
     params = artio_handle.parameters
 


https://bitbucket.org/yt_analysis/yt/commits/34aa96803266/
Changeset:   34aa96803266
Branch:      yt-3.0
User:        drudd
Date:        2014-07-12 17:28:44
Summary:     Merged with upstream
Affected #:  35 files

diff -r 617c99f73500d5c221d8acaf58f1cdfaee43e0c0 -r 34aa96803266b7d40d590ab6174f713d2fcf7e6b .hgignore
--- a/.hgignore
+++ b/.hgignore
@@ -7,6 +7,7 @@
 rockstar.cfg
 yt_updater.log
 yt/frontends/artio/_artio_caller.c
+yt/analysis_modules/halo_finding/rockstar/rockstar_groupies.c
 yt/analysis_modules/halo_finding/rockstar/rockstar_interface.c
 yt/frontends/ramses/_ramses_reader.cpp
 yt/frontends/sph/smoothing_kernel.c

diff -r 617c99f73500d5c221d8acaf58f1cdfaee43e0c0 -r 34aa96803266b7d40d590ab6174f713d2fcf7e6b CREDITS
--- a/CREDITS
+++ b/CREDITS
@@ -2,15 +2,21 @@
 
 Contributors:   
                 Tom Abel (tabel at stanford.edu)
-                David Collins (dcollins at physics.ucsd.edu)
+                Gabriel Altay (gabriel.altay at gmail.com)
+                Kenza Arraki (karraki at gmail.com)
+                Alex Bogert (fbogert at ucsc.edu)
+                David Collins (dcollins4096 at gmail.com)
                 Brian Crosby (crosby.bd at gmail.com)
                 Andrew Cunningham (ajcunn at gmail.com)
+                Miguel de Val-Borro (miguel.deval at gmail.com)
                 Hilary Egan (hilaryye at gmail.com)
                 John Forces (jforbes at ucolick.org)
+                Sam Geen (samgeen at gmail.com)
                 Nathan Goldbaum (goldbaum at ucolick.org)
                 Markus Haider (markus.haider at uibk.ac.at)
                 Cameron Hummels (chummels at gmail.com)
                 Christian Karch (chiffre at posteo.de)
+                Ben W. Keller (kellerbw at mcmaster.ca)
                 Ji-hoon Kim (me at jihoonkim.org)
                 Steffen Klemer (sklemer at phys.uni-goettingen.de)
                 Kacper Kowalik (xarthisius.kk at gmail.com)
@@ -21,18 +27,23 @@
                 Chris Malone (chris.m.malone at gmail.com)
                 Josh Maloney (joshua.moloney at colorado.edu)
                 Chris Moody (cemoody at ucsc.edu)
+                Stuart Mumford (stuart at mumford.me.uk)
                 Andrew Myers (atmyers at astro.berkeley.edu)
                 Jill Naiman (jnaiman at ucolick.org)
+                Desika Narayanan (dnarayan at haverford.edu)
                 Kaylea Nelson (kaylea.nelson at yale.edu)
                 Jeff Oishi (jsoishi at gmail.com)
+                Brian O'Shea (bwoshea at gmail.com)
                 Jean-Claude Passy (jcpassy at uvic.ca)
+                John Regan (john.regan at helsinki.fi)
                 Mark Richardson (Mark.L.Richardson at asu.edu)
                 Thomas Robitaille (thomas.robitaille at gmail.com)
                 Anna Rosen (rosen at ucolick.org)
                 Douglas Rudd (drudd at uchicago.edu)
                 Anthony Scopatz (scopatz at gmail.com)
                 Noel Scudder (noel.scudder at stonybrook.edu)
-                Devin Silvia (devin.silvia at colorado.edu)
+                Pat Shriwise (shriwise at wisc.edu)
+                Devin Silvia (devin.silvia at gmail.com)
                 Sam Skillman (samskillman at gmail.com)
                 Stephen Skory (s at skory.us)
                 Britton Smith (brittonsmith at gmail.com)
@@ -42,8 +53,10 @@
                 Stephanie Tonnesen (stonnes at gmail.com)
                 Matthew Turk (matthewturk at gmail.com)
                 Rich Wagner (rwagner at physics.ucsd.edu)
+                Michael S. Warren (mswarren at gmail.com)
                 Andrew Wetzel (andrew.wetzel at yale.edu)
                 John Wise (jwise at physics.gatech.edu)
+                Michael Zingale (michael.zingale at stonybrook.edu)
                 John ZuHone (jzuhone at gmail.com)
 
 Several items included in the yt/extern directory were written by other

diff -r 617c99f73500d5c221d8acaf58f1cdfaee43e0c0 -r 34aa96803266b7d40d590ab6174f713d2fcf7e6b MANIFEST.in
--- a/MANIFEST.in
+++ b/MANIFEST.in
@@ -12,4 +12,3 @@
 prune tests
 graft yt/gui/reason/html/resources
 exclude clean.sh .hgchurn
-recursive-include yt/utilities/kdtree *.f90 *.v Makefile LICENSE

diff -r 617c99f73500d5c221d8acaf58f1cdfaee43e0c0 -r 34aa96803266b7d40d590ab6174f713d2fcf7e6b doc/source/analyzing/analysis_modules/running_halofinder.rst
--- a/doc/source/analyzing/analysis_modules/running_halofinder.rst
+++ b/doc/source/analyzing/analysis_modules/running_halofinder.rst
@@ -300,40 +300,11 @@
 Therefore Parallel HOP is not a direct substitution for
 normal HOP, but is very similar.
 
-.. _fkd_setup:
-
-Fortran kD Tree Setup
-^^^^^^^^^^^^^^^^^^^^^
-
-Parallel HOP will not build automatically with yt. Please follow the instructions
-below in order to setup Parallel HOP.
-
-  #. Download `Forthon <http://hifweb.lbl.gov/Forthon/>`_. Extract the files
-     (e.g. tar -zxvf Forthon.tgz) and cd into the new Forthon directory. 
-     Making sure you're using the same version of python you use with yt, invoke
-     ``python setup.py install``.
-  #. Change directory to your yt source. Starting from the top level, cd into
-     ``yt/utilities/kdtree``.
-  #. Inside that directory, you should see these files:
-  
-     .. code-block:: bash
-     
-        % ls
-        Makefile        fKD.f90         fKD_source.f90
-        __init__.py     fKD.v           test.py
-  
-  #. Type ``make``. If that is successful, there should be a file in the
-     directory named ``fKDpy.so``. If there are problems, please contact the
-     `yt-users email list <http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org>`_.
-  #. Go to the top level of the yt source directory, which from the ``kdtree``
-     directory is three levels up ``cd ../../..``, and invoke
-     ``python setup.py install``.
-  #. Parallel HOP should now work.
-     
-
 Running Parallel HOP
 ^^^^^^^^^^^^^^^^^^^^
 
+Note: This is probably broken now that the Fortran kdtree has been removed.
+
 In the simplest form, Parallel HOP is run very similarly to the other halo finders.
 In the example below, Parallel HOP will be run on a dataset with all the default
 values. Parallel HOP can be run in serial, but as mentioned above, it is

diff -r 617c99f73500d5c221d8acaf58f1cdfaee43e0c0 -r 34aa96803266b7d40d590ab6174f713d2fcf7e6b doc/source/analyzing/units/2)_Data_Selection_and_fields.ipynb
--- a/doc/source/analyzing/units/2)_Data_Selection_and_fields.ipynb
+++ b/doc/source/analyzing/units/2)_Data_Selection_and_fields.ipynb
@@ -1,7 +1,7 @@
 {
  "metadata": {
   "name": "",
-  "signature": "sha256:b7541e0167001c6dd74306c8490385ace7bdb0533a829286f0505c0b24c67f16"
+  "signature": "sha256:882b31591c60bfe6ad4cb0f8842953d2e94fb8a12ce742be831a65642eea72c9"
  },
  "nbformat": 3,
  "nbformat_minor": 0,
@@ -325,8 +325,7 @@
      "input": [
       "from astropy import units as u\n",
       "x = 42.0 * u.meter\n",
-      "y = YTQuantity(x)\n",
-      "y2 = YTQuantity.from_astropy(x) # Another way to create the quantity"
+      "y = YTQuantity.from_astropy(x) "
      ],
      "language": "python",
      "metadata": {},
@@ -337,8 +336,7 @@
      "collapsed": false,
      "input": [
       "print x, type(x)\n",
-      "print y, type(y)\n",
-      "print y2, type(y2)"
+      "print y, type(y)"
      ],
      "language": "python",
      "metadata": {},
@@ -349,8 +347,7 @@
      "collapsed": false,
      "input": [
       "a = np.random.random(size=10) * u.km/u.s\n",
-      "b = YTArray(a)\n",
-      "b2 = YTArray.from_astropy(a) # Another way to create the quantity"
+      "b = YTArray.from_astropy(a)"
      ],
      "language": "python",
      "metadata": {},
@@ -361,8 +358,7 @@
      "collapsed": false,
      "input": [
       "print a, type(a)\n",
-      "print b, type(b)\n",
-      "print b2, type(b2)"
+      "print b, type(b)"
      ],
      "language": "python",
      "metadata": {},
@@ -438,7 +434,7 @@
      "collapsed": false,
      "input": [
       "k1 = kboltz.to_astropy()\n",
-      "k2 = YTQuantity(kb)\n",
+      "k2 = YTQuantity.from_astropy(kb)\n",
       "print k1 == k2"
      ],
      "language": "python",
@@ -449,7 +445,7 @@
      "cell_type": "code",
      "collapsed": false,
      "input": [
-      "c = YTArray(a)\n",
+      "c = YTArray.from_astropy(a)\n",
       "d = c.to_astropy()\n",
       "print a == d"
      ],

diff -r 617c99f73500d5c221d8acaf58f1cdfaee43e0c0 -r 34aa96803266b7d40d590ab6174f713d2fcf7e6b doc/source/cookbook/fits_radio_cubes.ipynb
--- a/doc/source/cookbook/fits_radio_cubes.ipynb
+++ b/doc/source/cookbook/fits_radio_cubes.ipynb
@@ -98,7 +98,7 @@
      "cell_type": "code",
      "collapsed": false,
      "input": [
-      "wcs_slc.save()"
+      "slc.save()"
      ],
      "language": "python",
      "metadata": {},
@@ -462,4 +462,4 @@
    "metadata": {}
   }
  ]
-}
\ No newline at end of file
+}

diff -r 617c99f73500d5c221d8acaf58f1cdfaee43e0c0 -r 34aa96803266b7d40d590ab6174f713d2fcf7e6b doc/source/cookbook/simple_plots.rst
--- a/doc/source/cookbook/simple_plots.rst
+++ b/doc/source/cookbook/simple_plots.rst
@@ -118,6 +118,8 @@
 
 .. yt_cookbook:: show_hide_axes_colorbar.py
 
+.. _matplotlib-primitives:
+
 Accessing and Modifying Plots Directly
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 

diff -r 617c99f73500d5c221d8acaf58f1cdfaee43e0c0 -r 34aa96803266b7d40d590ab6174f713d2fcf7e6b doc/source/reference/api/api.rst
--- a/doc/source/reference/api/api.rst
+++ b/doc/source/reference/api/api.rst
@@ -15,6 +15,7 @@
    ~yt.visualization.plot_window.OffAxisSlicePlot
    ~yt.visualization.plot_window.ProjectionPlot
    ~yt.visualization.plot_window.OffAxisProjectionPlot
+   ~yt.visualization.plot_window.WindowPlotMPL
 
 ProfilePlot and PhasePlot
 ^^^^^^^^^^^^^^^^^^^^^^^^^
@@ -24,6 +25,7 @@
 
    ~yt.visualization.profile_plotter.ProfilePlot
    ~yt.visualization.profile_plotter.PhasePlot
+   ~yt.visualization.profile_plotter.PhasePlotMPL
 
 Fixed Resolution Pixelization
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

diff -r 617c99f73500d5c221d8acaf58f1cdfaee43e0c0 -r 34aa96803266b7d40d590ab6174f713d2fcf7e6b doc/source/visualizing/_cb_docstrings.inc
--- a/doc/source/visualizing/_cb_docstrings.inc
+++ b/doc/source/visualizing/_cb_docstrings.inc
@@ -1,370 +1,445 @@
+Arrow callback
+~~~~~~~~~~~~~~
+
 .. function:: annotate_arrow(self, pos, code_size, plot_args=None):
 
-   (This is a proxy for :class:`~yt.visualization.plot_modifications.ArrowCallback`.)
+   (This is a proxy for
+   :class:`~yt.visualization.plot_modifications.ArrowCallback`.)
 
-   This adds an arrow pointing at *pos* with size
-   *code_size* in code units.  *plot_args* is a dict fed to
+   This adds an arrow pointing at ``pos`` with size
+   ``code_size`` in code units.  ``plot_args`` is a dict fed to
    matplotlib with arrow properties.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   slc = SlicePlot(pf, 'z', 'density', width=(10,'kpc'), center='max')
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   slc = yt.SlicePlot(ds, 'z', 'density', width=(10,'kpc'), center='c')
    slc.annotate_arrow((0.5, 0.5, 0.5), (1, 'kpc'))
    slc.save()
 
--------------
+Clump Finder Callback
+~~~~~~~~~~~~~~~~~~~~~
 
 .. function:: annotate_clumps(self, clumps, plot_args=None):
 
-   (This is a proxy for :class:`~yt.visualization.plot_modifications.ClumpContourCallback`.)
+   (This is a proxy for
+   :class:`~yt.visualization.plot_modifications.ClumpContourCallback`.)
 
-   Take a list of *clumps* and plot them as a set of
+   Take a list of ``clumps`` and plot them as a set of
    contours.
 
 .. python-script::
 
-   from yt.mods import *
-   from yt.analysis_modules.level_sets.api import *
+   import yt
+   import numpy as np
+   from yt.analysis_modules.level_sets.api import \
+       Clump, find_clumps, get_lowest_clumps
 
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   data_source = pf.disk([0.5, 0.5, 0.5], [0., 0., 1.],
-                           (8., 'kpc'), (1., 'kpc'))
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   data_source = ds.disk([0.5, 0.5, 0.5], [0., 0., 1.],
+                         (8., 'kpc'), (1., 'kpc'))
 
    c_min = 10**np.floor(np.log10(data_source['density']).min()  )
    c_max = 10**np.floor(np.log10(data_source['density']).max()+1)
 
-   function = 'self.data[\'Density\'].size > 20'
+   function = 'self.data[\'density\'].size > 20'
    master_clump = Clump(data_source, None, 'density', function=function)
    find_clumps(master_clump, c_min, c_max, 2.0)
    leaf_clumps = get_lowest_clumps(master_clump)
 
-   prj = ProjectionPlot(pf, 2, 'density', center='c', width=(20,'kpc'))
+   prj = yt.ProjectionPlot(ds, 2, 'density', center='c', width=(20,'kpc'))
    prj.annotate_clumps(leaf_clumps)
    prj.save('clumps')
 
--------------
+Overplot Contours
+~~~~~~~~~~~~~~~~~
 
-.. function:: annotate_contour(self, field, ncont=5, factor=4, take_log=False, clim=None, plot_args=None):
+.. function:: annotate_contour(self, field, ncont=5, factor=4, take_log=False,
+                               clim=None, plot_args=None):
 
-   (This is a proxy for :class:`~yt.visualization.plot_modifications.ContourCallback`.)
+   (This is a proxy for
+   :class:`~yt.visualization.plot_modifications.ContourCallback`.)
 
-   Add contours in *field* to the plot.  *ncont* governs the
-   number of contours generated, *factor* governs the number
-   of points used in the interpolation, *take_log* governs
-   how it is contoured and *clim* gives the (upper, lower)
-   limits for contouring.
+   Add contours in ``field`` to the plot.  ``ncont`` governs the number of
+   contours generated, ``factor`` governs the number of points used in the
+   interpolation, ``take_log`` governs how it is contoured and ``clim`` gives
+   the (upper, lower) limits for contouring.
 
 .. python-script::
-   
-   from yt.mods import *
-   pf = load("Enzo_64/DD0043/data0043")
-   s = SlicePlot(pf, "x", ["density"], center="max")
+
+   import yt
+   ds = yt.load("Enzo_64/DD0043/data0043")
+   s = yt.SlicePlot(ds, "x", "density", center="max")
    s.annotate_contour("temperature")
    s.save()
 
--------------
+Overplot quivers
+~~~~~~~~~~~~~~~~
+
+Axis-Aligned data sources
+^^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. function:: annotate_quiver(self, field_x, field_y, factor, scale=None,
+                              scale_units=None, normalize=False):
+
+   (This is a proxy for
+   :class:`~yt.visualization.plot_modifications.QuiverCallback`.)
+
+   Adds a 'quiver' plot to any plot, using the ``field_x`` and ``field_y`` from
+   the associated data, skipping every ``factor`` datapoints ``scale`` is the
+   data units per arrow length unit using ``scale_units`` (see
+   matplotlib.axes.Axes.quiver for more info)
+
+.. python-script::
+
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   p = yt.ProjectionPlot(ds, 'z', 'density', center=[0.5, 0.5, 0.5],
+                         weight_field='density', width=(20, 'kpc'))
+   p.annotate_quiver('velocity_x', 'velocity_y', 16)
+   p.save()
+
+Off-axis Data Sources
+^^^^^^^^^^^^^^^^^^^^^
 
 .. function:: annotate_cquiver(self, field_x, field_y, factor):
 
-   (This is a proxy for :class:`~yt.visualization.plot_modifications.CuttingQuiverCallback`.)
+   (This is a proxy for
+   :class:`~yt.visualization.plot_modifications.CuttingQuiverCallback`.)
 
-   Get a quiver plot on top of a cutting plane, using
-   *field_x* and *field_y*, skipping every *factor*
-   datapoint in the discretization.
+   Get a quiver plot on top of a cutting plane, using ``field_x`` and
+   ``field_y``, skipping every ``factor`` datapoint in the discretization.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("Enzo_64/DD0043/data0043")
-   s = OffAxisSlicePlot(pf, [1,1,0], ["density"], center="c")
+   import yt
+   ds = yt.load("Enzo_64/DD0043/data0043")
+   s = yt.OffAxisSlicePlot(ds, [1,1,0], ["density"], center="c")
    s.annotate_cquiver('cutting_plane_velocity_x', 'cutting_plane_velocity_y', 10)
    s.zoom(1.5)
    s.save()
 
--------------
+Overplot grids
+~~~~~~~~~~~~~~
 
-.. function:: annotate_grids(self, alpha=1.0, min_pix=1, annotate=False, periodic=True):
+.. function:: annotate_grids(self, alpha=1.0, min_pix=1, annotate=False,
+                             periodic=True):
 
-   (This is a proxy for :class:`~yt.visualization.plot_modifications.GridBoundaryCallback`.)
+   (This is a proxy for
+   :class:`~yt.visualization.plot_modifications.GridBoundaryCallback`.)
 
-   Adds grid boundaries to a plot, optionally with
-   *alpha*-blending. Cuttoff for display is at *min_pix*
-   wide. *annotate* puts the grid id in the corner of the
-   grid.  (Not so great in projections...)
+   Adds grid boundaries to a plot, optionally with alpha-blending via the
+   ``alpha`` keyword. Cuttoff for display is at ``min_pix`` wide. ``annotate``
+   puts the grid id in the corner of the grid.  (Not so great in projections...)
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   slc = SlicePlot(pf, 'z', 'density', width=(10,'kpc'), center='max')
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   slc = yt.SlicePlot(ds, 'z', 'density', width=(10,'kpc'), center='max')
    slc.annotate_grids()
    slc.save()
 
--------------
+Overplot Halo Annotations
+~~~~~~~~~~~~~~~~~~~~~~~~~
 
-.. function:: annotate_halos(self, halo_catalog, col='white', alpha =1, width= None):
+.. function:: annotate_halos(self, halo_catalog, col='white', alpha =1,
+                             width=None):
 
-   (This is a proxy for :class:`~yt.visualization.plot_modifications.HaloCatalogCallback`.)
+   (This is a proxy for
+   :class:`~yt.visualization.plot_modifications.HaloCatalogCallback`.)
 
-   Accepts a :class:`yt.HaloCatalog` *HaloCatalog* and plots 
-   a circle at the location of each halo with the radius of
-   the circle corresponding to the virial radius of the halo.
-   If *width* is set to None (default) all halos are plotted.
-   Otherwise, only halos that fall within a slab with width
-   *width* centered on the center of the plot data. The 
-   color and transparency of the circles can be controlled with
-   *col* and *alpha* respectively.
+   Accepts a :class:`yt.HaloCatalog` and plots a circle at the location of each
+   halo with the radius of the circle corresponding to the virial radius of the
+   halo.  If ``width`` is set to None (default) all halos are plotted.
+   Otherwise, only halos that fall within a slab with width ``width`` centered
+   on the center of the plot data. The color and transparency of the circles can
+   be controlled with ``col`` and ``alpha`` respectively.
 
 .. python-script::
-   
-   from yt.mods import *
+
+   import yt
    from yt.analysis_modules.halo_analysis.halo_catalog import HaloCatalog
 
-   data_pf = load('Enzo_64/RD0006/RedshiftOutput0006')
-   halos_pf = load('rockstar_halos/halos_0.0.bin')
+   data_ds = yt.load('Enzo_64/RD0006/RedshiftOutput0006')
+   halos_ds = yt.load('rockstar_halos/halos_0.0.bin')
 
-   hc = HaloCatalog(halos_pf=halos_pf)
+   hc = HaloCatalog(halos_pf=halos_ds)
    hc.create()
 
-   prj = ProjectionPlot(data_pf, 'z', 'density')
+   prj = yt.ProjectionPlot(data_ds, 'z', 'density')
    prj.annotate_halos(hc)
    prj.save()
 
--------------
+Overplot a Straight Line
+~~~~~~~~~~~~~~~~~~~~~~~~
 
-.. function:: annotate_image_line(self, p1, p2, data_coords=False, plot_args=None):
+.. function:: annotate_image_line(self, p1, p2, data_coords=False,
+                                  plot_args=None):
 
-   (This is a proxy for :class:`~yt.visualization.plot_modifications.ImageLineCallback`.)
+   (This is a proxy for
+   :class:`~yt.visualization.plot_modifications.ImageLineCallback`.)
 
-   Plot from *p1* to *p2* (normalized image plane coordinates) with
-   *plot_args* fed into the plot.
+   Plot from ``p1`` to ``p2`` (normalized image plane coordinates) with
+  ``plot_args`` fed into the plot.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   p = ProjectionPlot(pf, 'z', 'density', center='m', width=(10, 'kpc'))
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   p = yt.ProjectionPlot(ds, 'z', 'density', center='m', width=(10, 'kpc'))
    p.annotate_image_line((0.3, 0.4), (0.8, 0.9), plot_args={'linewidth':5})
    p.save()
 
--------------
+Overplot a line plot
+~~~~~~~~~~~~~~~~~~~~
 
 .. function:: annotate_line(self, x, y, plot_args=None):
 
-   (This is a proxy for :class:`~yt.visualization.plot_modifications.LinePlotCallback`.)
+   (This is a proxy for
+   :class:`~yt.visualization.plot_modifications.LinePlotCallback`.)
 
-   Over plot *x* and *y* (in code units) with *plot_args* fed into the plot.
+   Over plot numpy arrays or lists ``x`` and ``y`` (in code units) with
+   ``plot_args`` fed into the plot.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   p = ProjectionPlot(pf, 'z', 'density', center='m', width=(10, 'kpc'))
-   p.annotate_line([-6, -4, -2, 0, 2, 4, 6], [3.6, 1.6, 0.4, 0, 0.4, 1.6, 3.6], plot_args={'linewidth':5})
+   import yt
+   import numpy as np
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   p = yt.ProjectionPlot(ds, 'z', 'density', center='m', width=(20, 'kpc'))
+   x = np.array([-6, -4, -2, 0, 2, 4, 6])
+   y = x**2/10
+   p.annotate_line(x, y, plot_args={'linewidth':5})
    p.save()
 
--------------
+Overplot Magnetic Field Quivers
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-.. function:: annotate_magnetic_field(self, factor=16, scale=None, scale_units=None, normalize=False):
+.. function:: annotate_magnetic_field(self, factor=16, scale=None,
+                                      scale_units=None, normalize=False):
 
-   (This is a proxy for :class:`~yt.visualization.plot_modifications.MagFieldCallback`.)
+   (This is a proxy for
+   :class:`~yt.visualization.plot_modifications.MagFieldCallback`.)
 
-   Adds a 'quiver' plot of magnetic field to the plot,
-   skipping all but every *factor* datapoint. *scale* is the
-   data units per arrow length unit using *scale_units* (see
-   matplotlib.axes.Axes.quiver for more info). if
-   *normalize* is True, the magnetic fields will be scaled
-   by their local (in-plane) length, allowing morphological
-   features to be more clearly seen for fields with
-   substantial variation in field strength.
+   Adds a 'quiver' plot of magnetic field to the plot, skipping all but every
+   ``factor`` datapoint. ``scale`` is the data units per arrow length unit using
+   ``scale_units`` (see matplotlib.axes.Axes.quiver for more info). if
+   ``normalize`` is ``True``, the magnetic fields will be scaled by their local
+   (in-plane) length, allowing morphological features to be more clearly seen
+   for fields with substantial variation in field strength.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("MHDSloshing/virgo_low_res.0054.vtk",
-             parameters={"TimeUnits":3.1557e13, "LengthUnits":3.0856e24,
-                         "DensityUnits":6.770424595218825e-27})
-   p = ProjectionPlot(pf, 'z', 'density', center='c', width=(300, 'kpc'))
+   import yt
+   ds = yt.load("MHDSloshing/virgo_low_res.0054.vtk",
+                parameters={"time_unit":(1, 'Myr'), "length_unit":(1, 'Mpc'),
+                            "mass_unit":(1e17, 'Msun')})
+   p = yt.ProjectionPlot(ds, 'z', 'density', center='c', width=(300, 'kpc'))
    p.annotate_magnetic_field()
    p.save()
 
--------------
+Annotate a Point With a Marker
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 .. function:: annotate_marker(self, pos, marker='x', plot_args=None):
 
-   (This is a proxy for :class:`~yt.visualization.plot_modifications.MarkerAnnotateCallback`.)
+   (This is a proxy for
+   :class:`~yt.visualization.plot_modifications.MarkerAnnotateCallback`.)
 
-   Adds text *marker* at *pos* in code coordinates.
-   *plot_args* is a dict that will be forwarded to the plot
+   Adds ``marker`` at ``pos`` in code coordinates.
+   ``plot_args`` is a dict that will be forwarded to the plot
    command.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   s = SlicePlot(pf, 'z', 'density', center='m', width=(10, 'kpc'))
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   s = yt.SlicePlot(ds, 'z', 'density', center='c', width=(10, 'kpc'))
    s.annotate_marker([0.5, 0.5, 0.5], plot_args={'s':10000})
-   s.save()   
+   s.save()
 
--------------
+Overplotting Particle Positions
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-.. function:: annotate_particles(self, width, p_size=1.0, col='k', marker='o', stride=1.0, ptype=None, stars_only=False, dm_only=False, minimum_mass=None):
+.. function:: annotate_particles(self, width, p_size=1.0, col='k', marker='o',
+                                 stride=1.0, ptype=None, stars_only=False,
+                                 dm_only=False, minimum_mass=None):
 
-   (This is a proxy for :class:`~yt.visualization.plot_modifications.ParticleCallback`.)
+   (This is a proxy for
+   :class:`~yt.visualization.plot_modifications.ParticleCallback`.)
 
-   Adds particle positions, based on a thick slab along
-   *axis* with a *width* along the line of sight.  *p_size*
-   controls the number of pixels per particle, and *col*
-   governs the color.  *ptype* will restrict plotted
-   particles to only those that are of a given type.
-   *minimum_mass* will require that the particles be of a
-   given mass, calculated via ParticleMassMsun, to be
-   plotted.
+   Adds particle positions, based on a thick slab along ``axis`` with a
+   ``width`` along the line of sight.  ``p_size`` controls the number of pixels
+   per particle, and ``col`` governs the color.  ``ptype`` will restrict plotted
+   particles to only those that are of a given type.  ``minimum_mass`` will
+   require that the particles be of a given mass minimum mass in solar units.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("Enzo_64/DD0043/data0043")
-   p = ProjectionPlot(pf, "x", "density", center='m', width=(10, 'Mpc'))
+   import yt
+   ds = yt.load("Enzo_64/DD0043/data0043")
+   p = yt.ProjectionPlot(ds, "x", "density", center='m', width=(10, 'Mpc'))
    p.annotate_particles((10, 'Mpc'))
    p.save()
 
--------------
+Annotate a point with text
+~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 .. function:: annotate_point(self, pos, text, text_args=None):
 
-   (This is a proxy for :class:`~yt.visualization.plot_modifications.PointAnnotateCallback`.)
+   (This is a proxy for
+   :class:`~yt.visualization.plot_modifications.PointAnnotateCallback`.)
 
-   This adds *text* at position *pos*, where *pos* is in
-   code-space. *text_args* is a dict fed to the text
-   placement code.
+   This adds ``text`` at position ``pos``, where ``pos`` is in
+   code-space. ``text_args`` is a dict fed to the text placement code.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   p = ProjectionPlot(pf, 'z', 'density', center='m', width=(10, 'kpc'))
-   p.annotate_point([0.5, 0.496, 0.5], "What's going on here?", text_args={'size':'xx-large', 'color':'w'})
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   p = yt.ProjectionPlot(ds, 'z', 'density', center='m', width=(10, 'kpc'))
+   p.annotate_point([0.5, 0.496, 0.5], "What's going on here?",
+                    text_args={'size':'xx-large', 'color':'w'})
    p.save()
 
--------------
+Overplot a circle on a plot
+~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-.. function:: annotate_quiver(self, field_x, field_y, factor, scale=None, scale_units=None, normalize=False):
+.. function:: annotate_sphere(self, center, radius, circle_args=None, text=None,
+                              text_args=None):
 
-   (This is a proxy for :class:`~yt.visualization.plot_modifications.QuiverCallback`.)
+   (This is a proxy for
+   :class:`~yt.visualization.plot_modifications.SphereCallback`.)
 
-   Adds a 'quiver' plot to any plot, using the *field_x* and
-   *field_y* from the associated data, skipping every
-   *factor* datapoints *scale* is the data units per arrow
-   length unit using *scale_units*  (see
-   matplotlib.axes.Axes.quiver for more info)
+   A sphere centered at ``center`` in code units with radius ``radius`` in code
+   units will be created, with optional ``circle_args``, ``text``, and
+   ``text_args``.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   p = ProjectionPlot(pf, 'z', 'density', center=[0.5, 0.5, 0.5], 
-                      weight_field='density', width=(20, 'kpc'))
-   p.annotate_quiver('velocity_x', 'velocity_y', 16)
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   p = yt.ProjectionPlot(ds, 'z', 'density', center='c', width=(20, 'kpc'))
+   p.annotate_sphere([0.5, 0.5, 0.5], (2, 'kpc'), {'fill':True})
    p.save()
 
--------------
+Overplot streamlines
+~~~~~~~~~~~~~~~~~~~~
 
-.. function:: annotate_sphere(self, center, radius, circle_args=None, text=None, text_args=None):
+.. function:: annotate_streamlines(self, field_x, field_y, factor=6.0, nx=16,
+                                   ny=16, xstart=(0, 1), ystart=(0, 1),
+                                   nsample=256, start_at_xedge=False,
+                                   start_at_yedge=False, plot_args=None):
 
-   (This is a proxy for :class:`~yt.visualization.plot_modifications.SphereCallback`.)
+   (This is a proxy for
+   :class:`~yt.visualization.plot_modifications.StreamlineCallback`.)
 
-   A sphere centered at *center* in code units with radius
-   *radius* in code units will be created, with optional
-   *circle_args*, *text*, and *text_args*.
+   Add streamlines to any plot, using the ``field_x`` and ``field_y`` from the
+   associated data, using ``nx`` and ``ny`` starting points that are bounded by
+   ``xstart`` and ``ystart``.  To begin streamlines from the left edge of the
+   plot, set ``start_at_xedge`` to ``True``; for the bottom edge, use
+   ``start_at_yedge``.  A line with the qmean vector magnitude will cover
+   1.0/``factor`` of the image.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   p = ProjectionPlot(pf, 'z', 'density', center='c', width=(20, 'kpc'))
-   p.annotate_sphere([0.5, 0.5, 0.5], (2, 'kpc'), {'fill':True})
-   p.save()
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   s = yt.SlicePlot(ds, 'z', 'density', center='c', width=(20, 'kpc'))
+   s.annotate_streamlines('velocity_x', 'velocity_y')
+   s.save()
 
--------------
+Add text
+~~~~~~~~
 
-.. function:: annotate_streamlines(self, field_x, field_y, factor=6.0, nx=16, ny=16, xstart=(0, 1), ystart=(0, 1), nsample=256, start_at_xedge=False, start_at_yedge=False, plot_args=None):
+.. function:: annotate_text(self, pos, text, data_coords=False, text_args=None):
 
-   (This is a proxy for :class:`~yt.visualization.plot_modifications.StreamlineCallback`.)
+   (This is a proxy for
+   :class:`~yt.visualization.plot_modifications.TextLabelCallback`.)
 
-   Add streamlines to any plot, using the *field_x* and
-   *field_y* from the associated data, using *nx* and *ny*
-   starting points that are bounded by *xstart* and
-   *ystart*.  To begin streamlines from the left edge of the
-   plot, set *start_at_xedge* to True; for the bottom edge,
-   use *start_at_yedge*.  A line with the qmean vector
-   magnitude will cover 1.0/*factor* of the image.
+   Accepts a position in (0..1, 0..1) of the image, some text and optionally
+   some text arguments. If data_coords is True, position will be in code units
+   instead of image coordinates.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   s = SlicePlot(pf, 'z', 'density', center='c', width=(20, 'kpc'))
-   s.annotate_streamlines('velocity_x', 'velocity_y')
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   s = yt.SlicePlot(ds, 'z', 'density', center='m', width=(10, 'kpc'))
+   s.annotate_text((0.5, 0.5), 'Sample text',
+                   text_args={'size':'xx-large', 'color':'w'})
    s.save()
 
--------------
+Add a title to the plot
+~~~~~~~~~~~~~~~~~~~~~~~
 
-.. function:: annotate_text(self, pos, text, data_coords=False, text_args=None):
+.. function:: annotate_title(self, title='Plot'):
 
-   (This is a proxy for :class:`~yt.visualization.plot_modifications.TextLabelCallback`.)
+   (This is a proxy for
+   :class:`~yt.visualization.plot_modifications.TitleCallback`.)
 
-   Accepts a position in (0..1, 0..1) of the image, some
-   text and optionally some text arguments. If data_coords
-   is True, position will be in code units instead of image
-   coordinates.
+   Accepts a ``title`` and adds it to the plot.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   s = SlicePlot(pf, 'z', 'density', center='m', width=(10, 'kpc'))
-   s.annotate_text((0.5, 0.5), 'Sample text', text_args={'size':'xx-large', 'color':'w'})
-   s.save()
-
--------------
-
-.. function:: annotate_title(self, title='Plot'):
-
-   (This is a proxy for :class:`~yt.visualization.plot_modifications.TitleCallback`.)
-
-   Accepts a *title* and adds it to the plot
-
-.. python-script::
-
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   p = ProjectionPlot(pf, 'z', 'density', center='c', width=(20, 'kpc'))
-   p.annotate_title('Density plot')
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   p = yt.ProjectionPlot(ds, 'z', 'density', center='c', width=(20, 'kpc'))
+   p.annotate_title('Density Plot')
    p.save()
 
--------------
+Overplot quivers for the velocity field
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-.. function:: annotate_velocity(self, factor=16, scale=None, scale_units=None, normalize=False):
+.. function:: annotate_velocity(self, factor=16, scale=None, scale_units=None,
+                                normalize=False):
 
-   (This is a proxy for :class:`~yt.visualization.plot_modifications.VelocityCallback`.)
+   (This is a proxy for
+   :class:`~yt.visualization.plot_modifications.VelocityCallback`.)
 
-   Adds a 'quiver' plot of velocity to the plot, skipping
-   all but every *factor* datapoint. *scale* is the data
-   units per arrow length unit using *scale_units* (see
-   matplotlib.axes.Axes.quiver for more info). if
-   *normalize* is True, the velocity fields will be scaled
-   by their local (in-plane) length, allowing morphological
-   features to be more clearly seen for fields with
-   substantial variation in field strength (normalize is not
+   Adds a 'quiver' plot of velocity to the plot, skipping all but every
+   ``factor`` datapoint. ``scale`` is the data units per arrow length unit using
+   ``scale_units`` (see matplotlib.axes.Axes.quiver for more info). if
+   ``normalize`` is ``True``, the velocity fields will be scaled by their local
+   (in-plane) length, allowing morphological features to be more clearly seen
+   for fields with substantial variation in field strength (normalize is not
    implemented and thus ignored for Cutting Planes).
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   p = SlicePlot(pf, 'z', 'density', center='m', width=(10, 'kpc'))
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   p = yt.SlicePlot(ds, 'z', 'density', center='m', width=(10, 'kpc'))
    p.annotate_velocity()
    p.save()
+
+Add a Timestamp Inset Box
+~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. function:: annotate_timestamp(x, y, units=None, format="{time:.3G} {units}", 
+                                 **kwargs, normalized=False, bbox_dict=None)
+
+   (This is a proxy for
+   :class:`~yt.visualization.plot_modifications.TimestampCallback`.)
+
+   Adds the current time to the plot at point given by *x* and *y*.  If *units*
+   is given ('s', 'ms', 'ns', etc), it will covert the time to this basis.  If
+   *units* is None, it will attempt to figure out the correct value by which to
+   scale.  The *format* keyword is a template string that will be evaluated and
+   displayed on the plot.  If *normalized* is true, *x* and *y* are interpreted
+   as normalized plot coordinates (0,0 is lower-left and 1,1 is upper-right)
+   otherwise *x* and *y* are assumed to be in plot coordinates. The *bbox_dict*
+   is an optional dict of arguments for the bbox that frames the timestamp, see
+   matplotlib's text annotation guide for more details. All other *kwargs* will
+   be passed to the text() method on the plot axes.  See matplotlib's text()
+   functions for more information.
+
+.. python-script::
+
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   p = yt.SlicePlot(ds, 'z', 'density', center='c', width=(20, 'kpc'))
+   p.annotate_timestamp()
+   p.save()

diff -r 617c99f73500d5c221d8acaf58f1cdfaee43e0c0 -r 34aa96803266b7d40d590ab6174f713d2fcf7e6b doc/source/visualizing/manual_plotting.rst
--- a/doc/source/visualizing/manual_plotting.rst
+++ b/doc/source/visualizing/manual_plotting.rst
@@ -35,11 +35,12 @@
 .. python-script::
    
    import pylab as P
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   import numpy as np
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
 
-   c = pf.h.find_max('density')[1]
-   proj = pf.proj('density', 0)
+   c = ds.find_max('density')[1]
+   proj = ds.proj('density', 0)
 
    width = (10, 'kpc') # we want a 1.5 mpc view
    res = [1000, 1000] # create an image with 1000x1000 pixels
@@ -64,22 +65,33 @@
 Line Plots
 ----------
 
-This is perhaps the simplest thing to do. ``yt`` provides a number of one dimensional objects, and these return a 1-D numpy array of their contents with direct dictionary access. As a simple example, take a :class:`~yt.data_objects.data_containers.AMROrthoRayBase` object, which can be created from a index by calling ``pf.ortho_ray(axis, center)``. 
+This is perhaps the simplest thing to do. ``yt`` provides a number of one
+dimensional objects, and these return a 1-D numpy array of their contents with
+direct dictionary access. As a simple example, take a
+:class:`~yt.data_objects.data_containers.AMROrthoRayBase` object, which can be
+created from a index by calling ``pf.ortho_ray(axis, center)``.
 
 .. python-script::
 
-   from yt.mods import *
+   import yt
+   import numpy as np
    import pylab as P
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   c = pf.h.find_max("density")[1]
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   c = ds.find_max("density")[1]
    ax = 0 # take a line cut along the x axis
-   ray = pf.ortho_ray(ax, (c[1], c[2])) # cutting through the y0,z0 such that we hit the max density
+
+   # cutting through the y0,z0 such that we hit the max density
+   ray = ds.ortho_ray(ax, (c[1], c[2]))
+
+   # Sort the ray values by 'x' so there are no discontinuities
+   # in the line plot
+   srt = np.argsort(ray['x'])
 
    P.subplot(211)
-   P.semilogy(np.array(ray['x']), np.array(ray['density']))
+   P.semilogy(np.array(ray['x'][srt]), np.array(ray['density'][srt]))
    P.ylabel('density')
    P.subplot(212)
-   P.semilogy(np.array(ray['x']), np.array(ray['temperature']))
+   P.semilogy(np.array(ray['x'][srt]), np.array(ray['temperature'][srt]))
    P.xlabel('x')
    P.ylabel('temperature')
 

diff -r 617c99f73500d5c221d8acaf58f1cdfaee43e0c0 -r 34aa96803266b7d40d590ab6174f713d2fcf7e6b doc/source/visualizing/plots.rst
--- a/doc/source/visualizing/plots.rst
+++ b/doc/source/visualizing/plots.rst
@@ -25,10 +25,10 @@
 If you need to take a quick look at a single simulation output, ``yt``
 provides the ``PlotWindow`` interface for generating annotated 2D
 visualizations of simulation data.  You can create a ``PlotWindow`` plot by
-supplying a parameter file, a list of fields to plot, and a plot center to
-create a :class:`~yt.visualization.plot_window.SlicePlot`,
+supplying a dataset, a list of fields to plot, and a plot center to
+create a :class:`~yt.visualization.plot_window.AxisAlignedSlicePlot`,
 :class:`~yt.visualization.plot_window.ProjectionPlot`, or
-:class:`~yt.visualization.plot_window.OffAxisSlicePlot`.
+:class:`~yt.visualization.plot_window.OffAxisProjectionPlot`.
 
 Plot objects use ``yt`` data objects to extract the maximum resolution
 data available to render a 2D image of a field. Whenever a
@@ -37,7 +37,8 @@
 is requested of it -- for instance, when the width or field is changed
 -- this high-resolution data is then pixelized and placed in a buffer
 of fixed size. This is accomplished behind the scenes using
-:class:`~yt.visualization.fixed_resolution.FixedResolutionBuffer`
+:class:`~yt.visualization.fixed_resolution.FixedResolutionBuffer`.
+
 ``PlotWindow`` expose the underlying matplotlib ``figure`` and
 ``axes`` objects, making it easy to customize your plots and 
 add new annotations.
@@ -45,28 +46,33 @@
 .. _slice-plots:
 
 Slice Plots
------------
+~~~~~~~~~~~
 
-The quickest way to plot a slice of a field through your data is to use
-:class:`~yt.visualization.plot_window.SlicePlot`.  Say we want to visualize a
-slice through the Density field along the z-axis centered on the center of the
-simulation box in a simulation dataset we've opened and stored in the parameter
-file object ``pf``.  This can be accomplished with the following command:
+The quickest way to plot a slice of a field through your data is via
+:class:`~yt.visualization.plot_window.SlicePlot`.  These plots are generally
+quicker than projections because they only need to read and process a slice
+through the dataset.
+
+The following script plots a slice through the density field along the z-axis
+centered on the center of the simulation box in a simulation dataset we've
+opened and stored in ``ds``:
 
 .. code-block:: python
 
-   >>> slc = SlicePlot(pf, 'z', 'density')
+   >>> slc = yt.SlicePlot(ds, 'z', 'density')
    >>> slc.save()
 
 These two commands will create a slice object and store it in a variable we've
-called ``slc``.  We then call the ``save()`` function that is associated with
-the slice object.  This automatically saves the plot in png image format with an
-automatically generated filename.  If you don't want the slice object to stick
-around, you can accomplish the same thing in one line:
+called ``slc``.  Since this plot is aligned with the simulation coordinate
+system, ``slc`` is an instance of
+:class:`~yt.visualization.plot_window.AxisAlignedSlicePlot`. We then call the
+``save()`` function, which automatically saves the plot in png image format with
+an automatically generated filename.  If you don't want the slice object to
+stick around, you can accomplish the same thing in one line:
 
 .. code-block:: python
    
-   >>> SlicePlot(pf, 'z', 'density').save()
+   >>> yt.SlicePlot(ds, 'z', 'density').save()
 
 It's nice to keep the slice object around if you want to modify the plot.  By
 default, the plot width will be set to the size of the simulation box.  To zoom
@@ -75,41 +81,58 @@
 
 .. code-block:: python
 
-   >>> slc = SlicePlot(pf, 'z', 'density')
+   >>> slc = yt.SlicePlot(ds, 'z', 'density')
    >>> slc.zoom(10)
    >>> slc.save('zoom')
 
 This will save a new plot to disk with a different filename - prepended with
-'zoom' instead of the name of the parameter file. If you want to set the width
+'zoom' instead of the name of the dataset. If you want to set the width
 manually, you can do that as well. For example, the following sequence of
 commands will create a slice, set the width of the plot to 10 kiloparsecs, and
 save it to disk.
 
 .. code-block:: python
 
-   >>> slc = SlicePlot(pf, 'z', 'density')
-   >>> slc.set_width((10,'kpc'))
+   >>> from yt.units import kpc
+   >>> slc = yt.SlicePlot(ds, 'z', 'density')
+   >>> slc.set_width(10*kpc)
    >>> slc.save('10kpc')
 
-The SlicePlot also optionally accepts the coordinate to center the plot on and
-the width of the plot:
+The plot width can be specified independently along the x and y direction by
+passing a tuple of widths.  An individual width can also be represented using a
+``(value, unit)`` tuple.  The following sequence of commands all equivalently
+set the width of the plot to 200 kiloparsecs in the ``x`` and ``y`` direction.
 
 .. code-block:: python
 
-   >>> SlicePlot(pf, 'z', 'density', center=[0.2, 0.3, 0.8], 
-   ...           width = (10,'kpc')).save()
+   >>> from yt.units import kpc
+   >>> slc.set_width(200*kpc)
+   >>> slc.set_width((200, 'kpc'))
+   >>> slc.set_width((200*kpc, 200*kpc))
 
-The center must be given in code units.  Optionally, you can supply 'c' or 'm'
-for the center.  These two choices will center the plot on the center of the
-simulation box and the coordinate of the maximum density cell, respectively.
+The ``SlicePlot`` also optionally accepts the coordinate to center the plot on
+and the width of the plot:
+
+.. code-block:: python
+
+   >>> yt.SlicePlot(ds, 'z', 'density', center=[0.2, 0.3, 0.8],
+   ...              width = (10,'kpc')).save()
+
+The plot center is relative to the simulation coordinate system.  If supplied
+without units, the center is assumed by in code units.  Optionally, you can
+supply 'c' or 'm' for the center.  These two choices will center the plot on the
+center of the simulation box and the coordinate of the maximum density cell,
+respectively.
 
 Here is an example that combines all of the options we just discussed.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   slc = SlicePlot(pf, 'z', 'density', center=[0.5, 0.5, 0.5], width=(20,'kpc'))
+   import yt
+   from yt.units import kpc
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   slc = yt.SlicePlot(ds, 'z', 'density', center=[0.5, 0.5, 0.5],
+                      width=(20,'kpc'))
    slc.save()
 
 The above example will display an annotated plot of a slice of the
@@ -118,14 +141,14 @@
 letter 'z', corresponding to the z-axis.  Finally, the image is saved to
 a png file.
 
-Conceptually, you can think of the SlicePlot as an adjustable window
+Conceptually, you can think of the plot object as an adjustable window
 into the data. For example:
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   slc = SlicePlot(pf, 'z', 'pressure', center='c')
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   slc = yt.SlicePlot(ds, 'z', 'pressure', center='c')
    slc.save()
    slc.zoom(30)
    slc.save('zoom')
@@ -138,32 +161,43 @@
 interesting region in the simulation and adjust the boundaries of the
 region to visualize on the fly.
 
-A slice object can also add annotations like a title, an overlying
-quiver plot, the location of grid boundaries, halo-finder annotations,
-and many other annotations, including user-customizable annotations.
-For example:
+See :class:`~yt.visualization.plot_window.AxisAlignedSlicePlot` for the 
+full class description.
+
+.. _off-axis-slices:
+
+Off Axis Slices
+~~~~~~~~~~~~~~~
+
+Off axis slice plots can be generated in much the same way as
+grid-aligned slices.  Off axis slices use
+:class:`~yt.data_objects.data_containers.AMRCuttingPlaneBase` to slice
+through simulation domains at an arbitrary oblique angle.  A
+:class:`~yt.visualization.plot_window.OffAxisSlicePlot` can be
+instantiated by specifying a dataset, the normal to the cutting
+plane, and the name of the fields to plot.  Just like an
+:class:`~yt.visualization.plot_window.AxisAlignedSlicePlot`, an
+``OffAxisSlicePlot`` can be created via the
+:class:`~yt.visualization.plot_window.SlicePlot` function. For example:
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   slc = SlicePlot(pf, 'z', 'density', width=(10,'kpc'))
-   slc.annotate_grids()
-   slc.save()
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   L = [1,1,0] # vector normal to cutting plane
+   north_vector = [-1,1,0]
+   cut = yt.SlicePlot(ds, L, 'density', width=(25, 'kpc'),
+                      north_vector=north_vector)
+   cut.save()
 
-will plot the density field in a 10 kiloparsec slice through the
-z-axis centered on the highest density point in the simulation domain.
-Before saving the plot, the script annotates it with the grid
-boundaries, which are drawn as thick black lines by default.
-
-Annotations are described in :ref:`callbacks`.  See
-:class:`~yt.visualization.plot_window.SlicePlot` for the full class
-description.
+In this case, a normal vector for the cutting plane is supplied in the second
+argument. Optionally, a `north_vector` can be specified to fix the orientation
+of the image plane.
 
 .. _projection-plots:
 
 Projection Plots
-----------------
+~~~~~~~~~~~~~~~~
 
 Using a fast adaptive projection, ``yt`` is able to quickly project
 simulation data along the coordinate axes.
@@ -174,15 +208,15 @@
 
 .. python-script::
  
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   prj = ProjectionPlot(pf, 2, 'density', width=(25, 'kpc'), 
-                        weight_field=None)
+   import yt
+   from yt.units import kpc
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   prj = yt.ProjectionPlot(ds, 2, 'temperature', width=25*kpc,
+                           weight_field='density')
    prj.save()
 
-will create a projection of Density field along the x axis, plot it,
-and then save it to a png image file.  The projection is only carried
-out to level 2 of the AMR index and no weighting is applied.
+will create a density-weighted projection of the temperature field along the x
+axis, plot it, and then save the plot to a png image file.
 
 Like :ref:`slice-plots`, annotations and modifications can be applied
 after creating the ``ProjectionPlot`` object.  Annotations are
@@ -190,43 +224,8 @@
 :class:`~yt.visualization.plot_window.ProjectionPlot` for the full
 class description.
 
-.. _off-axis-slices:
-
-Off Axis Slice Plots
---------------------
-
-Off axis slice plots can be generated in much the same way as
-grid-aligned slices.  Off axis slices use
-:class:`~yt.data_objects.data_containers.AMRCuttingPlaneBase` to slice
-through simulation domains at an arbitrary oblique angle.  A
-:class:`~yt.visualization.plot_window.OffAxisSlicePlot` can be
-instantiated by specifying a parameter file, the normal to the cutting
-plane, and the name of the fields to plot.  For example:
-
-.. python-script::
-
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   L = [1,1,0] # vector normal to cutting plane
-   north_vector = [-1,1,0]
-   cut = OffAxisSlicePlot(pf, L, 'density', width=(25, 'kpc'), 
-                          north_vector=north_vector)
-   cut.save()
-
-creates an off-axis slice in the plane perpendicular to ``L``,
-oriented such that ``north_vector`` is the up direction.  If ``L`` and
-``north_vector`` are not perpendicular, the component of
-``north_vector`` perpendicular to ``L`` is used. Like
-:ref:`slice-plots`, annotations and modifications can be applied after
-creating the ``OffAxisSlicePlot`` object.  Annotations are described
-in :ref:`callbacks`.  See
-:class:`~yt.visualization.plot_window.OffAxisSlicePlot` for the full
-class description.
-
-.. _off-axis-projections:
-
 Off Axis Projection Plots
--------------------------
+~~~~~~~~~~~~~~~~~~~~~~~~~
 
 Off axis projection plots .  Internally, off axis projections are
 created using :ref:`the-camera-interface` by applying the
@@ -246,15 +245,16 @@
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   import yt
+   import numpy as np
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
    L = [1,1,0] # vector normal to cutting plane
    north_vector = [-1,1,0]
    W = [0.02, 0.02, 0.02]
    c = [0.5, 0.5, 0.5]
    N = 512
-   image = off_axis_projection(pf, c, L, W, N, "density")
-   write_image(np.log10(image), "%s_offaxis_projection.png" % pf)
+   image = yt.off_axis_projection(ds, c, L, W, N, "density")
+   yt.write_image(np.log10(image), "%s_offaxis_projection.png" % ds)
 
 Here, ``W`` is the width of the projection in the x, y, *and* z
 directions.
@@ -268,32 +268,30 @@
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
    L = [1,1,0] # vector normal to cutting plane
    north_vector = [-1,1,0]
-   prj = OffAxisProjectionPlot(pf,L,'density',width=(25, 'kpc'), 
-                               north_vector=north_vector)
+   prj = yt.OffAxisProjectionPlot(ds,L,'density',width=(25, 'kpc'),
+                                  north_vector=north_vector)
    prj.save()
 
 OffAxisProjectionPlots can also be created with a number of
-keyword arguments, as described in the `api reference`__ for the
-class initializer.
-
-__ :class:`~yt.visualization.plot_window.OffAxisProjectionPlot`
+keyword arguments, as described in
+:class:`~yt.visualization.plot_window.OffAxisProjectionPlot`
 
 Plot Customization
 ------------------
 
 You can customize each of the four plot types above in identical ways.  We'll go
-over each of the customizations methos below.  For each of the examples below we
+over each of the customizations methods below.  For each of the examples below we
 will modify the following plot.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   slc = SlicePlot(pf, 'z', 'density', width=(10,'kpc'))
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   slc = yt.SlicePlot(ds, 'z', 'density', width=(10,'kpc'))
    slc.save()
 
 Panning and zooming
@@ -301,78 +299,84 @@
 
 There are three methods to dynamically pan around the data.  
 
-:class:`~yt.visualization.plot_window.SlicePlot.pan` accepts x and y deltas in code
-units.
+:class:`~yt.visualization.plot_window.AxisAlignedSlicePlot.pan` accepts x and y
+deltas.
 
 .. python-script::
 
-   from yt.mods import *
+   import yt
    from yt.units import kpc
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   slc = SlicePlot(pf, 'z', 'density', width=(10,'kpc'))
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   slc = yt.SlicePlot(ds, 'z', 'density', width=(10,'kpc'))
    slc.pan((2*kpc, 2*kpc))
    slc.save()
 
-:class:`~yt.visualization.plot_window.SlicePlot.pan_rel` accepts deltas in units relative
-to the field of view of the plot.  
+:class:`~yt.visualization.plot_window.AxisAlignedSlicePlot.pan_rel` accepts deltas 
+in units relative to the field of view of the plot.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   slc = SlicePlot(pf, 'z', 'density', width=(10,'kpc'))
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   slc = yt.SlicePlot(ds, 'z', 'density', width=(10,'kpc'))
    slc.pan_rel((0.1, -0.1))
    slc.save()
 
-:class:`~yt.visualization.plot_window.SlicePlot.zoom` accepts a factor to zoom in by.
+:class:`~yt.visualization.plot_window.AxisAlignedSlicePlot.zoom` accepts a factor to zoom in by.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   slc = SlicePlot(pf, 'z', 'density', width=(10,'kpc'))
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   slc = yt.SlicePlot(ds, 'z', 'density', width=(10,'kpc'))
    slc.zoom(2)
    slc.save()
 
 Set axes units
 ~~~~~~~~~~~~~~
 
-:class:`~yt.visualization.plot_window.SlicePlot.set_axes_unit` allows the customization of
+:class:`~yt.visualization.plot_window.AxisAlignedSlicePlot.set_axes_unit` allows the customization of
 the axes unit labels.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   slc = SlicePlot(pf, 'z', 'density', width=(10,'kpc'))
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   slc = yt.SlicePlot(ds, 'z', 'density', width=(10,'kpc'))
    slc.set_axes_unit('Mpc')
    slc.save()
 
+The same result could have been accomplished by explicitly setting the ``width``
+to ``(.01, 'Mpc')``.
+
 Set the plot center
 ~~~~~~~~~~~~~~~~~~~
 
-The :class:`~yt.visualization.plot_window.SlicePlot.set_center` function accepts a new
-center for the plot, in code units.  New centers must be two element tuples.
+The :class:`~yt.visualization.plot_window.AxisAlignedSlicePlot.set_center`
+function accepts a new center for the plot, in code units.  New centers must be
+two element tuples.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   slc = SlicePlot(pf, 'z', 'density', width=(10,'kpc'))
-   slc.set_center((0.5, 0.5))
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   slc = yt.SlicePlot(ds, 'z', 'density', width=(10,'kpc'))
+   slc.set_center((0.5, 0.503))
    slc.save()
 
 Fonts
 ~~~~~
 
-:class:`~yt.visualization.plot_window.SlicePlot.set_font` allows font costomization.
+:class:`~yt.visualization.plot_window.AxisAlignedSlicePlot.set_font` allows font
+costomization.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   slc = SlicePlot(pf, 'z', 'density', width=(10,'kpc'))
-   slc.set_font({'family': 'sans-serif', 'style': 'italic','weight': 'bold', 'size': 24})
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   slc = yt.SlicePlot(ds, 'z', 'density', width=(10,'kpc'))
+   slc.set_font({'family': 'sans-serif', 'style': 'italic',
+                 'weight': 'bold', 'size': 24})
    slc.save()
 
 Colormaps
@@ -383,117 +387,173 @@
 different fields tracked by the plot object.
 
 To change the colormap for the plot, call the
-:class:`~yt.visualization.plot_window.SlicePlot.set_cmap` function.  Use any of the
-colormaps listed in the :ref:`colormaps` section.
+:class:`~yt.visualization.plot_window.AxisAlignedSlicePlot.set_cmap` function.
+Use any of the colormaps listed in the :ref:`colormaps` section.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   slc = SlicePlot(pf, 'z', 'density', width=(10,'kpc'))
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   slc = yt.SlicePlot(ds, 'z', 'density', width=(10,'kpc'))
    slc.set_cmap('density', 'RdBu_r')
    slc.save()
 
-The :class:`~yt.visualization.plot_window.SlicePlot.set_log` function accepts a field name
-and a boolean.  If the boolean is :code:`True`, the colormap for the field will
-be log scaled.  If it is `False` the colormap will be linear.
+The :class:`~yt.visualization.plot_window.AxisAlignedSlicePlot.set_log` function
+accepts a field name and a boolean.  If the boolean is ``True``, the colormap
+for the field will be log scaled.  If it is ``False`` the colormap will be
+linear.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   slc = SlicePlot(pf, 'z', 'density', width=(10,'kpc'))
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   slc = yt.SlicePlot(ds, 'z', 'density', width=(10,'kpc'))
    slc.set_log('density', False)
    slc.save()
 
-Lastly, the :class:`~yt.visualization.plot_window.SlicePlot.set_zlim` function makes it
-possible to set a custom colormap range.
+Lastly, the :class:`~yt.visualization.plot_window.AxisAlignedSlicePlot.set_zlim`
+function makes it possible to set a custom colormap range.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   slc = SlicePlot(pf, 'z', 'density', width=(10,'kpc'))
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   slc = yt.SlicePlot(ds, 'z', 'density', width=(10,'kpc'))
    slc.set_zlim('density', 1e-30, 1e-25)
    slc.save()
 
+Annotations
+~~~~~~~~~~~
+
+A slice object can also add annotations like a title, an overlying
+quiver plot, the location of grid boundaries, halo-finder annotations,
+and many other annotations, including user-customizable annotations.
+For example:
+
+.. python-script::
+
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   slc = yt.SlicePlot(ds, 'z', 'density', width=(10,'kpc'))
+   slc.annotate_grids()
+   slc.save()
+
+will plot the density field in a 10 kiloparsec slice through the
+z-axis centered on the highest density point in the simulation domain.
+Before saving the plot, the script annotates it with the grid
+boundaries, which are drawn as thick black lines by default.
+
+Annotations are described in :ref:`callbacks`.
+
 Set the size of the plot
 ~~~~~~~~~~~~~~~~~~~~~~~~
 
 To set the size of the plot, use the
-:class:`~yt.visualization.plot_window.SlicePlot.set_window_size` function.  The argument
+:class:`~yt.visualization.plot_window.AxisAlignedSlicePlot.set_figure_size` function.  The argument
 is the size of the longest edge of the plot in inches.  View the full resolution
 image to see the difference more clearly.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   slc = SlicePlot(pf, 'z', 'density', width=(10,'kpc'))
-   slc.set_window_size(10)
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   slc = yt.SlicePlot(ds, 'z', 'density', width=(10,'kpc'))
+   slc.set_figure_size(10)
    slc.save()
 
 To change the resolution of the image, call the
-:class:`~yt.visualization.plot_window.SlicePlot.set_buff_size` function.
+:class:`~yt.visualization.plot_window.AxisAlignedSlicePlot.set_buff_size` function.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   slc = SlicePlot(pf, 'z', 'density', width=(10,'kpc'))
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   slc = yt.SlicePlot(ds, 'z', 'density', width=(10,'kpc'))
    slc.set_buff_size(1600)
    slc.save()
 
+Further customization via matplotlib
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Each ``PlotWindow`` object is really a container for plots - one plot for each
+field specified in the list of fields supplied when the plot object is
+created. The individual plots can be accessed via the ``plots`` dictionary
+attached to each ``PlotWindow`` object:
+
+.. code-block:: python
+
+   >>> slc = SlicePlot(ds, 2, ['density', 'temperature']
+   >>> dens_plot = slc.plots['density']
+
+In this example ``dens_plot`` is an instance of
+:class:`~yt.visualization.plot_window.WindowPlotMPL`, an object that wraps the
+matplotlib ``figure`` and ``axes`` objects.  We can access these matplotlib primitives via attributes of ``dens_plot``.
+
+.. code-block:: python
+
+   >>> figure = dens_plot.figure
+   >>> axes = dens_plot.axes
+   >>> colorbar_axes = dens_plot.cax
+
+These are the :matplotlib:class:`figure`, and :matplotlib:class:`axes` objects
+that control the actual drawing of the plot.  Arbitrary plot customizations
+are possible by manipulating these objects.  See :ref:`matplotlib-primitives` for
+an example.
+
 .. _how-to-make-1d-profiles:
 
 1D Profile Plots
 ----------------
 
-1D profiles are used to calculated the average or the sum of a given quantity
-with respect to a second quantity.  This means the "average density as a
-function of radius" or "the total mass within a given set of density bins."
+1D profiles are used to calculate the average or the sum of a given quantity
+with respect to a second quantity.  Two common examples are the "average density
+as a function of radius" or "the total mass within a given set of density bins."
 When created, they default to the average: in fact, they default to the average
 as weighted by the total cell mass.  However, this can be modified to take
 either the total value or the average with respect to a different quantity.
 
-Profiles operate on data objects; they will take the entire data contained in a
-sphere, a prism, an extracted region and so on, and they will calculate and use
-that as input to their calculation.  To make a 1D profile plot, create a
-(:class:`~yt.visualization.profile_plotter.ProfilePlot`) object, supplying the 
-data object, the field for binning, and a list of fields to be profiled.
+Profiles operate on :ref:`data objects <using-objects>`; they will take the
+entire data contained in a sphere, a prism, an extracted region and so on, and
+they will calculate and use that as input to their calculation.  To make a 1D
+profile plot, create a (:class:`~yt.visualization.profile_plotter.ProfilePlot`)
+object, supplying the data object, the field for binning, and a list of fields
+to be profiled.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   my_galaxy = pf.disk([0.5, 0.5, 0.5], [0.0, 0.0, 1.0], 0.01, 0.003)
-   plot = ProfilePlot(my_galaxy, "density", ["temperature"])
+   import yt
+   from yt.units import kpc
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   my_galaxy = ds.disk(ds.domain_center, [0.0, 0.0, 1.0], 10*kpc, 3*kpc)
+   plot = yt.ProfilePlot(my_galaxy, "density", ["temperature"])
    plot.save()
 
 This will create a :class:`yt.data_objects.data_containers.AMRCylinderBase`
-centered at [0.3, 0.5, 0.8], with a normal vector of [0.4, 0.5, 0.1], radius of
-0.01 and height of 0.001 and will then make a plot of the average (as a 
-function of the cell mass) temperature as a function of density.
+centered at [0.5, 0.5, 0.5], with a normal vector of [0.0, 0.0, 1.0], radius of
+10 kiloparsecs and height of 3 kiloparsecs and will then make a plot of the
+mass-weighted average temperature as a function of density for all of the gas
+contained in the cylinder.
 
-As another example, we create a sphere of radius 100 pc and plot total mass 
-in every equally-spaced temperature bin/
-
-We could also have allowed the plot collection to create a sphere for us, as
-well.  For instance:
+We could also have made a profile considering only the gas in a sphere.
+For instance:
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   my_sphere = pf.sphere([0.5, 0.5, 0.5], (100, "kpc"))
-   plot = ProfilePlot(my_sphere, "temperature", ["cell_mass"],
-                      weight_field=None)
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   my_sphere = ds.sphere([0.5, 0.5, 0.5], (100, "kpc"))
+   plot = yt.ProfilePlot(my_sphere, "temperature", ["cell_mass"],
+                         weight_field=None)
    plot.save()
 
-Note that because we have specified the weighting field to be none, it operates 
-as a local-bin accumulator.  We can also accumulate along the x-axis by setting 
-the **accumulation** keyword argument to True, which is useful for plots of 
-enclosed mass.
+Note that because we have specified the weighting field to be none, it operates
+the profile plot will display the accumulated cell mass as a function of
+temperature rather than the average.  We can also accumulate along the x-axis by
+setting the **accumulation** keyword argument to True, which is useful for plots
+of enclosed mass.
+
+Also note the use of a ``(value, unit)`` tuple. These can be used interchangably
+with units explicitly imported from ``yt.units``.
 
 You can also access the data generated by profiles directly, which can be
 useful for overplotting average quantities on top of phase plots, or for
@@ -506,32 +566,30 @@
 
    plot = ProfilePlot(my_sphere, "temperature", ["cell_mass"],
                       weight_field=None)
-   # print the x field
+   profile = plot.profiles[0]
+   # print the bin field, in this case temperature
    print plot.profiles[-1].x
-   # print the profiled temperature field
-   print plot.profiles[-1].field_data["temperature"]
+   # print the profiled cell_mass field
+   print plot.profiles[-1]["cell_mass"]
 
-Other options, such as the number of bins, are also configurable. See the 
-documentation for 
-The number of bins and other options and tweaks are 
-available for these methods.  See the documentation for 
-:class:`~yt.visualization.profile_plotter.ProfilePlot`
-for more information.
+Other options, such as the number of bins, are also configurable. See the
+documentation for :class:`~yt.visualization.profile_plotter.ProfilePlot` for
+more information.
 
 Overplotting Multiple 1D Profiles
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 It is often desirable to overplot multiple 1D profile to show evolution 
 with time.  This is supported with the ``from_profiles`` class method.  
-1D profiles are created with the :meth:`yt.data_objects.profiles.create_profile` 
+1D profiles are created with the :func:`yt.data_objects.profiles.create_profile` 
 method and then given to the ProfilePlot object.
 
 .. python-script::
 
-   from yt.mods import *
+   import yt
 
    # Create a time-series object.
-   es = simulation("enzo_tiny_cosmology/32Mpc_32.enzo", "Enzo")
+   es = yt.simulation("enzo_tiny_cosmology/32Mpc_32.enzo", "Enzo")
    es.get_time_series(redshifts=[5, 4, 3, 2, 1, 0])
 
 
@@ -540,23 +598,114 @@
    labels = []
 
    # Loop over each dataset in the time-series.
-   for pf in es:
+   for ds in es:
        # Create a data container to hold the whole dataset.
-       ad = pf.h.all_data()
+       ad = ds.all_data()
        # Create a 1d profile of density vs. temperature.
-       profiles.append(create_profile(ad, ["temperature"], 
-                                      fields=["cell_mass"],
-                                      weight_field=None,
-                                      accumulation=True))
+       profiles.append(yt.create_profile(ad, ["temperature"], 
+                                         fields=["cell_mass"],
+                                         weight_field=None,
+                                         accumulation=True))
        # Add labels
-       labels.append("z = %.2f" % pf.current_redshift)
+       labels.append("z = %.2f" % ds.current_redshift)
 
    # Create the profile plot from the list of profiles.
-   plot = ProfilePlot.from_profiles(profiles, labels=labels)
+   plot = yt.ProfilePlot.from_profiles(profiles, labels=labels)
 
    # Save the image.
    plot.save()
 
+Customizing axis limits
+~~~~~~~~~~~~~~~~~~~~~~~
+
+By default the x and y limits for ``ProfilePlot`` are determined using the
+:class:`~yt.data_objects.derived_quantities.Extrema` derived quantity.  If you
+want to create a plot with custom axis limits, you have two options.
+
+First, you can create a custom profile object using
+:func:`~yt.data_objects.profiles.create_profile`.  This function accepts a dictionary of ``(max, min)`` tuples keyed to field names.
+
+.. python-script::
+
+    import yt
+    import yt.units as u
+    ds = yt.load('IsolatedGalaxy/galaxy0030/galaxy0030')
+    sp = ds.sphere('m', 10*u.kpc)
+    profiles = yt.create_profile(sp, "temperature", "density",
+                                 weight_field=None, 
+                                 extrema={'temperature': (1e3, 1e7),
+                                          'density': (1e-26, 1e-22)})
+    plot = yt.ProfilePlot.from_profiles(profiles)
+    plot.save()
+
+You can also make use of the
+:meth:`yt.visualization.profile_plotter.ProfilePlot.set_xlim` and
+:meth:`yt.visualization.profile_plotter.ProfilePlot.set_ylim` functions to
+customize the axes limits of a plot that has already been created.  Note that
+calling ``set_xlim`` is much slower than calling ``set_ylim``.  This is because
+```set_xlim`` must recreate the profile object using the specified extrema.
+Creating a profile directly via ``create_profile`` might be significant faster.
+If you find that this operation is slow, consider creating your profile via
+``create_profile``.  Note that since there is only one bin field, ``set_xlim``
+does not accept a field name as the first argument.
+
+.. python-script::
+
+   import yt
+   import yt.units as u
+   ds = yt.load('IsolatedGalaxy/galaxy0030/galaxy0030')
+   sp = ds.sphere('m', 10*u.kpc)
+   plot = yt.ProfilePlot(sp, "temperature", "density", weight_field=None)
+   plot.set_xlim(1e3, 1e7)
+   plot.set_ylim("density", 1e-26, 1e-22)
+   plot.save()
+
+
+Customizing Units
+~~~~~~~~~~~~~~~~~
+
+Units for both the x and y axis can be controlled via the
+:meth:`~yt.visualization.profile_plotter.ProfilePlot.set_unit` method.
+Adjusting the plot units does not require recreating the histogram, so adjusting
+units will always be inexpensive, requiring only an in-place unit conversion.
+
+In the following example we create a a plot of the average density in solar
+masses per cubic parsec as a function of radius in kiloparsecs.
+
+.. python-script::
+
+    import yt
+    import yt.units as u
+    ds = yt.load('IsolatedGalaxy/galaxy0030/galaxy0030')
+    sp = ds.sphere('m', 10*u.kpc)
+    plot = yt.ProfilePlot(sp, "radius", "density", weight_field=None)
+    plot.set_unit("density", "msun/pc**3")
+    plot.set_unit("radius", "kpc")
+    plot.save()
+
+Linear and Logarithmic Scaling
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The axis scaling can be manipulated via the
+:meth:`~yt.visualization.profile_plotter.ProfilePlot.set_log` function.  This
+function accepts a field name and a boolean.  If the boolean is ``True``, the
+field is plotted in log scale.  If ``False``, the field is plotted in linear
+scale.
+
+In the following example we create a plot of the average x velocity as a
+function of radius.  Since the x component of the velocity vector can be
+negative, we set the scaling to be linear for this field.
+
+.. python-script::
+
+   import yt
+   import yt.units as u
+   ds = yt.load('IsolatedGalaxy/galaxy0030/galaxy0030')
+   sp = ds.sphere('m', 10*u.kpc)
+   plot = yt.ProfilePlot(sp, "radius", "x-velocity", weight_field=None)
+   plot.set_log("x-velocity", False)
+   plot.save()
+
 Altering Line Properties
 ~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -581,47 +730,72 @@
 2D Phase Plots
 --------------
 
-2D phase plots function in much the same was as 1D phase plots, but with a 
-:class:`~yt.visualization.profile_plotter.PhasePlot` object.  Much like 1D profiles, 
-2D profiles (phase plots) are best thought of as plotting a distribution of points, 
-either taking the average or the accumulation in a bin.  For example, to generate a 
-2D distribution of mass enclosed in density and temperature bins, you can do:
+2D phase plots function in much the same was as 1D phase plots, but with a
+:class:`~yt.visualization.profile_plotter.PhasePlot` object.  Much like 1D
+profiles, 2D profiles (phase plots) are best thought of as plotting a
+distribution of points, either taking the average or the accumulation in a bin.
+For example, to generate a 2D distribution of mass enclosed in density and
+temperature bins, you can do:
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   my_sphere = pf.sphere("c", (50, "kpc"))
-   plot = PhasePlot(my_sphere, "density", "temperature", ["cell_mass"],
-                    weight_field=None)
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   my_sphere = ds.sphere("c", (50, "kpc"))
+   plot = yt.PhasePlot(my_sphere, "density", "temperature", ["cell_mass"],
+                       weight_field=None)
    plot.save()
 
 If you would rather see the average value of a field as a function of two other
-fields, you can neglect supplying the *weight* parameter.  This would look
+fields, you can set the ``weight_field`` parameter to ``None``.  This would look
 something like:
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   my_sphere = pf.sphere("c", (50, "kpc"))
-   plot = PhasePlot(my_sphere, "density", "temperature", ["H_fraction"],
-                    weight_field="cell_mass")
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   my_sphere = ds.sphere("c", (50, "kpc"))
+   plot = yt.PhasePlot(my_sphere, "density", "temperature", ["H_fraction"],
+                       weight_field=None)
+   plot.save()
+
+Customizing Phase Plots
+~~~~~~~~~~~~~~~~~~~~~~~
+
+Similarly to 1D profile plots, ``PhasePlot`` can be customized via ``set_unit``,
+``set_xlim``, ``set_ylim``, and ``set_zlim``.  The following example illustrates
+how to manipulate these functions.
+
+.. python-script::
+
+   import yt
+   ds = yt.load("sizmbhloz-clref04SNth-rs9_a0.9011/sizmbhloz-clref04SNth-rs9_a0.9011.art")
+   center = ds.arr([64.0, 64.0, 64.0], 'code_length')
+   rvir = ds.quan(1e-1, "Mpccm/h")
+
+   sph = ds.sphere(center, rvir)
+   plot = yt.PhasePlot(sph, "density", "temperature", "cell_mass",
+                       weight_field=None)
+   plot.set_unit('density', 'Msun/pc**3')
+   plot.set_unit('cell_mass', 'Msun')
+   plot.set_xlim(1e-5,1e1)
+   plot.set_ylim(1,1e7)
    plot.save()
 
 Probability Distribution Functions and Accumulation
 ---------------------------------------------------
 
-Both 1D and 2D profiles which show the total of amount of some field, such as mass, 
-in a bin (done by setting the **weight_field** keyword to None) can be turned into 
-probability distribution functions (PDFs) by setting the **fractional** keyword to 
-True.  When set to True, the value in each bin is divided by the sum total from all 
-bins.  These can be turned into cumulative distribution functions (CDFs) by setting 
-the **accumulation** keyword to True.  This will make is so that the value in any bin 
-N is the cumulative sum of all bins from 0 to N.  The direction of the summation can be 
-rversed by setting **accumulation** to -True.  For PhasePlots, the accumulation can be 
-set independently for each axis by setting **accumulation** to a list of True/-True/False 
-values.
+Both 1D and 2D profiles which show the total of amount of some field, such as
+mass, in a bin (done by setting the ``weight_field`` keyword to ``None``) can be
+turned into probability distribution functions (PDFs) by setting the
+``fractional`` keyword to ``True``.  When set to ``True``, the value in each bin
+is divided by the sum total from all bins.  These can be turned into cumulative
+distribution functions (CDFs) by setting the ``accumulation`` keyword to
+``True``.  This will make is so that the value in any bin N is the cumulative
+sum of all bins from 0 to N.  The direction of the summation can be reversed by
+setting ``accumulation`` to ``-True``.  For ``PhasePlot``, the accumulation can
+be set independently for each axis by setting ``accumulation`` to a list of
+``True``/ ``-True`` /``False`` values.
 
 .. _interactive-plotting:
 
@@ -646,9 +820,9 @@
 
 .. notebook-cell::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   p = ProjectionPlot(pf, "x", "density", center='m', width=(10,'kpc'),
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   p = ProjectionPlot(ds, "x", "density", center='m', width=(10,'kpc'),
                       weight_field='density')
    p.show()
 
@@ -682,7 +856,7 @@
 .. code-block:: python
 
    >>> import yt.visualization.eps_writer as eps
-   >>> slc = SlicePlot(pf, 'z', 'density')
+   >>> slc = yt.SlicePlot(ds, 'z', 'density')
    >>> slc.set_width(25, 'kpc')
    >>> eps_fig = eps.single_plot(slc)
    >>> eps_fig.save_fig('zoom', format='eps')
@@ -706,9 +880,11 @@
 
 .. code-block:: python
 
+   >>> import yt
    >>> import yt.visualization.eps_writer as eps
-   >>> slc = SlicePlot(pf, 'z', ['density', 'temperature', 'Pressure',
-                       'VelocityMagnitude'])
+   >>>
+   >>> slc = yt.SlicePlot(ds, 'z', ['density', 'temperature', 'Pressure',
+                          'VelocityMagnitude'])
    >>> slc.set_width(25, 'kpc')
    >>> eps_fig = eps.multiplot_yt(2, 2, slc, bare_axes=True)
    >>> eps_fig.scale_line(0.2, '5 kpc')

This diff is so big that we needed to truncate the remainder.

https://bitbucket.org/yt_analysis/yt/commits/3b264449bb41/
Changeset:   3b264449bb41
Branch:      yt-3.0
User:        drudd
Date:        2014-07-12 17:36:43
Summary:     Merged
Affected #:  36 files

diff -r 84525632d209088422ef3d90d7bdddc3bba03f30 -r 3b264449bb411235ba00ab34a01df79261899f6b .hgignore
--- a/.hgignore
+++ b/.hgignore
@@ -7,6 +7,7 @@
 rockstar.cfg
 yt_updater.log
 yt/frontends/artio/_artio_caller.c
+yt/analysis_modules/halo_finding/rockstar/rockstar_groupies.c
 yt/analysis_modules/halo_finding/rockstar/rockstar_interface.c
 yt/frontends/ramses/_ramses_reader.cpp
 yt/frontends/sph/smoothing_kernel.c

diff -r 84525632d209088422ef3d90d7bdddc3bba03f30 -r 3b264449bb411235ba00ab34a01df79261899f6b CREDITS
--- a/CREDITS
+++ b/CREDITS
@@ -2,15 +2,21 @@
 
 Contributors:   
                 Tom Abel (tabel at stanford.edu)
-                David Collins (dcollins at physics.ucsd.edu)
+                Gabriel Altay (gabriel.altay at gmail.com)
+                Kenza Arraki (karraki at gmail.com)
+                Alex Bogert (fbogert at ucsc.edu)
+                David Collins (dcollins4096 at gmail.com)
                 Brian Crosby (crosby.bd at gmail.com)
                 Andrew Cunningham (ajcunn at gmail.com)
+                Miguel de Val-Borro (miguel.deval at gmail.com)
                 Hilary Egan (hilaryye at gmail.com)
                 John Forces (jforbes at ucolick.org)
+                Sam Geen (samgeen at gmail.com)
                 Nathan Goldbaum (goldbaum at ucolick.org)
                 Markus Haider (markus.haider at uibk.ac.at)
                 Cameron Hummels (chummels at gmail.com)
                 Christian Karch (chiffre at posteo.de)
+                Ben W. Keller (kellerbw at mcmaster.ca)
                 Ji-hoon Kim (me at jihoonkim.org)
                 Steffen Klemer (sklemer at phys.uni-goettingen.de)
                 Kacper Kowalik (xarthisius.kk at gmail.com)
@@ -21,18 +27,23 @@
                 Chris Malone (chris.m.malone at gmail.com)
                 Josh Maloney (joshua.moloney at colorado.edu)
                 Chris Moody (cemoody at ucsc.edu)
+                Stuart Mumford (stuart at mumford.me.uk)
                 Andrew Myers (atmyers at astro.berkeley.edu)
                 Jill Naiman (jnaiman at ucolick.org)
+                Desika Narayanan (dnarayan at haverford.edu)
                 Kaylea Nelson (kaylea.nelson at yale.edu)
                 Jeff Oishi (jsoishi at gmail.com)
+                Brian O'Shea (bwoshea at gmail.com)
                 Jean-Claude Passy (jcpassy at uvic.ca)
+                John Regan (john.regan at helsinki.fi)
                 Mark Richardson (Mark.L.Richardson at asu.edu)
                 Thomas Robitaille (thomas.robitaille at gmail.com)
                 Anna Rosen (rosen at ucolick.org)
                 Douglas Rudd (drudd at uchicago.edu)
                 Anthony Scopatz (scopatz at gmail.com)
                 Noel Scudder (noel.scudder at stonybrook.edu)
-                Devin Silvia (devin.silvia at colorado.edu)
+                Pat Shriwise (shriwise at wisc.edu)
+                Devin Silvia (devin.silvia at gmail.com)
                 Sam Skillman (samskillman at gmail.com)
                 Stephen Skory (s at skory.us)
                 Britton Smith (brittonsmith at gmail.com)
@@ -42,8 +53,10 @@
                 Stephanie Tonnesen (stonnes at gmail.com)
                 Matthew Turk (matthewturk at gmail.com)
                 Rich Wagner (rwagner at physics.ucsd.edu)
+                Michael S. Warren (mswarren at gmail.com)
                 Andrew Wetzel (andrew.wetzel at yale.edu)
                 John Wise (jwise at physics.gatech.edu)
+                Michael Zingale (michael.zingale at stonybrook.edu)
                 John ZuHone (jzuhone at gmail.com)
 
 Several items included in the yt/extern directory were written by other

diff -r 84525632d209088422ef3d90d7bdddc3bba03f30 -r 3b264449bb411235ba00ab34a01df79261899f6b MANIFEST.in
--- a/MANIFEST.in
+++ b/MANIFEST.in
@@ -12,4 +12,3 @@
 prune tests
 graft yt/gui/reason/html/resources
 exclude clean.sh .hgchurn
-recursive-include yt/utilities/kdtree *.f90 *.v Makefile LICENSE

diff -r 84525632d209088422ef3d90d7bdddc3bba03f30 -r 3b264449bb411235ba00ab34a01df79261899f6b doc/source/analyzing/analysis_modules/running_halofinder.rst
--- a/doc/source/analyzing/analysis_modules/running_halofinder.rst
+++ b/doc/source/analyzing/analysis_modules/running_halofinder.rst
@@ -300,40 +300,11 @@
 Therefore Parallel HOP is not a direct substitution for
 normal HOP, but is very similar.
 
-.. _fkd_setup:
-
-Fortran kD Tree Setup
-^^^^^^^^^^^^^^^^^^^^^
-
-Parallel HOP will not build automatically with yt. Please follow the instructions
-below in order to setup Parallel HOP.
-
-  #. Download `Forthon <http://hifweb.lbl.gov/Forthon/>`_. Extract the files
-     (e.g. tar -zxvf Forthon.tgz) and cd into the new Forthon directory. 
-     Making sure you're using the same version of python you use with yt, invoke
-     ``python setup.py install``.
-  #. Change directory to your yt source. Starting from the top level, cd into
-     ``yt/utilities/kdtree``.
-  #. Inside that directory, you should see these files:
-  
-     .. code-block:: bash
-     
-        % ls
-        Makefile        fKD.f90         fKD_source.f90
-        __init__.py     fKD.v           test.py
-  
-  #. Type ``make``. If that is successful, there should be a file in the
-     directory named ``fKDpy.so``. If there are problems, please contact the
-     `yt-users email list <http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org>`_.
-  #. Go to the top level of the yt source directory, which from the ``kdtree``
-     directory is three levels up ``cd ../../..``, and invoke
-     ``python setup.py install``.
-  #. Parallel HOP should now work.
-     
-
 Running Parallel HOP
 ^^^^^^^^^^^^^^^^^^^^
 
+Note: This is probably broken now that the Fortran kdtree has been removed.
+
 In the simplest form, Parallel HOP is run very similarly to the other halo finders.
 In the example below, Parallel HOP will be run on a dataset with all the default
 values. Parallel HOP can be run in serial, but as mentioned above, it is

diff -r 84525632d209088422ef3d90d7bdddc3bba03f30 -r 3b264449bb411235ba00ab34a01df79261899f6b doc/source/analyzing/units/2)_Data_Selection_and_fields.ipynb
--- a/doc/source/analyzing/units/2)_Data_Selection_and_fields.ipynb
+++ b/doc/source/analyzing/units/2)_Data_Selection_and_fields.ipynb
@@ -1,7 +1,7 @@
 {
  "metadata": {
   "name": "",
-  "signature": "sha256:b7541e0167001c6dd74306c8490385ace7bdb0533a829286f0505c0b24c67f16"
+  "signature": "sha256:882b31591c60bfe6ad4cb0f8842953d2e94fb8a12ce742be831a65642eea72c9"
  },
  "nbformat": 3,
  "nbformat_minor": 0,
@@ -325,8 +325,7 @@
      "input": [
       "from astropy import units as u\n",
       "x = 42.0 * u.meter\n",
-      "y = YTQuantity(x)\n",
-      "y2 = YTQuantity.from_astropy(x) # Another way to create the quantity"
+      "y = YTQuantity.from_astropy(x) "
      ],
      "language": "python",
      "metadata": {},
@@ -337,8 +336,7 @@
      "collapsed": false,
      "input": [
       "print x, type(x)\n",
-      "print y, type(y)\n",
-      "print y2, type(y2)"
+      "print y, type(y)"
      ],
      "language": "python",
      "metadata": {},
@@ -349,8 +347,7 @@
      "collapsed": false,
      "input": [
       "a = np.random.random(size=10) * u.km/u.s\n",
-      "b = YTArray(a)\n",
-      "b2 = YTArray.from_astropy(a) # Another way to create the quantity"
+      "b = YTArray.from_astropy(a)"
      ],
      "language": "python",
      "metadata": {},
@@ -361,8 +358,7 @@
      "collapsed": false,
      "input": [
       "print a, type(a)\n",
-      "print b, type(b)\n",
-      "print b2, type(b2)"
+      "print b, type(b)"
      ],
      "language": "python",
      "metadata": {},
@@ -438,7 +434,7 @@
      "collapsed": false,
      "input": [
       "k1 = kboltz.to_astropy()\n",
-      "k2 = YTQuantity(kb)\n",
+      "k2 = YTQuantity.from_astropy(kb)\n",
       "print k1 == k2"
      ],
      "language": "python",
@@ -449,7 +445,7 @@
      "cell_type": "code",
      "collapsed": false,
      "input": [
-      "c = YTArray(a)\n",
+      "c = YTArray.from_astropy(a)\n",
       "d = c.to_astropy()\n",
       "print a == d"
      ],

diff -r 84525632d209088422ef3d90d7bdddc3bba03f30 -r 3b264449bb411235ba00ab34a01df79261899f6b doc/source/cookbook/fits_radio_cubes.ipynb
--- a/doc/source/cookbook/fits_radio_cubes.ipynb
+++ b/doc/source/cookbook/fits_radio_cubes.ipynb
@@ -98,7 +98,7 @@
      "cell_type": "code",
      "collapsed": false,
      "input": [
-      "wcs_slc.save()"
+      "slc.save()"
      ],
      "language": "python",
      "metadata": {},
@@ -462,4 +462,4 @@
    "metadata": {}
   }
  ]
-}
\ No newline at end of file
+}

diff -r 84525632d209088422ef3d90d7bdddc3bba03f30 -r 3b264449bb411235ba00ab34a01df79261899f6b doc/source/cookbook/simple_plots.rst
--- a/doc/source/cookbook/simple_plots.rst
+++ b/doc/source/cookbook/simple_plots.rst
@@ -118,6 +118,8 @@
 
 .. yt_cookbook:: show_hide_axes_colorbar.py
 
+.. _matplotlib-primitives:
+
 Accessing and Modifying Plots Directly
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 

diff -r 84525632d209088422ef3d90d7bdddc3bba03f30 -r 3b264449bb411235ba00ab34a01df79261899f6b doc/source/reference/api/api.rst
--- a/doc/source/reference/api/api.rst
+++ b/doc/source/reference/api/api.rst
@@ -15,6 +15,7 @@
    ~yt.visualization.plot_window.OffAxisSlicePlot
    ~yt.visualization.plot_window.ProjectionPlot
    ~yt.visualization.plot_window.OffAxisProjectionPlot
+   ~yt.visualization.plot_window.WindowPlotMPL
 
 ProfilePlot and PhasePlot
 ^^^^^^^^^^^^^^^^^^^^^^^^^
@@ -24,6 +25,7 @@
 
    ~yt.visualization.profile_plotter.ProfilePlot
    ~yt.visualization.profile_plotter.PhasePlot
+   ~yt.visualization.profile_plotter.PhasePlotMPL
 
 Fixed Resolution Pixelization
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

diff -r 84525632d209088422ef3d90d7bdddc3bba03f30 -r 3b264449bb411235ba00ab34a01df79261899f6b doc/source/visualizing/_cb_docstrings.inc
--- a/doc/source/visualizing/_cb_docstrings.inc
+++ b/doc/source/visualizing/_cb_docstrings.inc
@@ -1,370 +1,445 @@
+Arrow callback
+~~~~~~~~~~~~~~
+
 .. function:: annotate_arrow(self, pos, code_size, plot_args=None):
 
-   (This is a proxy for :class:`~yt.visualization.plot_modifications.ArrowCallback`.)
+   (This is a proxy for
+   :class:`~yt.visualization.plot_modifications.ArrowCallback`.)
 
-   This adds an arrow pointing at *pos* with size
-   *code_size* in code units.  *plot_args* is a dict fed to
+   This adds an arrow pointing at ``pos`` with size
+   ``code_size`` in code units.  ``plot_args`` is a dict fed to
    matplotlib with arrow properties.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   slc = SlicePlot(pf, 'z', 'density', width=(10,'kpc'), center='max')
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   slc = yt.SlicePlot(ds, 'z', 'density', width=(10,'kpc'), center='c')
    slc.annotate_arrow((0.5, 0.5, 0.5), (1, 'kpc'))
    slc.save()
 
--------------
+Clump Finder Callback
+~~~~~~~~~~~~~~~~~~~~~
 
 .. function:: annotate_clumps(self, clumps, plot_args=None):
 
-   (This is a proxy for :class:`~yt.visualization.plot_modifications.ClumpContourCallback`.)
+   (This is a proxy for
+   :class:`~yt.visualization.plot_modifications.ClumpContourCallback`.)
 
-   Take a list of *clumps* and plot them as a set of
+   Take a list of ``clumps`` and plot them as a set of
    contours.
 
 .. python-script::
 
-   from yt.mods import *
-   from yt.analysis_modules.level_sets.api import *
+   import yt
+   import numpy as np
+   from yt.analysis_modules.level_sets.api import \
+       Clump, find_clumps, get_lowest_clumps
 
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   data_source = pf.disk([0.5, 0.5, 0.5], [0., 0., 1.],
-                           (8., 'kpc'), (1., 'kpc'))
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   data_source = ds.disk([0.5, 0.5, 0.5], [0., 0., 1.],
+                         (8., 'kpc'), (1., 'kpc'))
 
    c_min = 10**np.floor(np.log10(data_source['density']).min()  )
    c_max = 10**np.floor(np.log10(data_source['density']).max()+1)
 
-   function = 'self.data[\'Density\'].size > 20'
+   function = 'self.data[\'density\'].size > 20'
    master_clump = Clump(data_source, None, 'density', function=function)
    find_clumps(master_clump, c_min, c_max, 2.0)
    leaf_clumps = get_lowest_clumps(master_clump)
 
-   prj = ProjectionPlot(pf, 2, 'density', center='c', width=(20,'kpc'))
+   prj = yt.ProjectionPlot(ds, 2, 'density', center='c', width=(20,'kpc'))
    prj.annotate_clumps(leaf_clumps)
    prj.save('clumps')
 
--------------
+Overplot Contours
+~~~~~~~~~~~~~~~~~
 
-.. function:: annotate_contour(self, field, ncont=5, factor=4, take_log=False, clim=None, plot_args=None):
+.. function:: annotate_contour(self, field, ncont=5, factor=4, take_log=False,
+                               clim=None, plot_args=None):
 
-   (This is a proxy for :class:`~yt.visualization.plot_modifications.ContourCallback`.)
+   (This is a proxy for
+   :class:`~yt.visualization.plot_modifications.ContourCallback`.)
 
-   Add contours in *field* to the plot.  *ncont* governs the
-   number of contours generated, *factor* governs the number
-   of points used in the interpolation, *take_log* governs
-   how it is contoured and *clim* gives the (upper, lower)
-   limits for contouring.
+   Add contours in ``field`` to the plot.  ``ncont`` governs the number of
+   contours generated, ``factor`` governs the number of points used in the
+   interpolation, ``take_log`` governs how it is contoured and ``clim`` gives
+   the (upper, lower) limits for contouring.
 
 .. python-script::
-   
-   from yt.mods import *
-   pf = load("Enzo_64/DD0043/data0043")
-   s = SlicePlot(pf, "x", ["density"], center="max")
+
+   import yt
+   ds = yt.load("Enzo_64/DD0043/data0043")
+   s = yt.SlicePlot(ds, "x", "density", center="max")
    s.annotate_contour("temperature")
    s.save()
 
--------------
+Overplot quivers
+~~~~~~~~~~~~~~~~
+
+Axis-Aligned data sources
+^^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. function:: annotate_quiver(self, field_x, field_y, factor, scale=None,
+                              scale_units=None, normalize=False):
+
+   (This is a proxy for
+   :class:`~yt.visualization.plot_modifications.QuiverCallback`.)
+
+   Adds a 'quiver' plot to any plot, using the ``field_x`` and ``field_y`` from
+   the associated data, skipping every ``factor`` datapoints ``scale`` is the
+   data units per arrow length unit using ``scale_units`` (see
+   matplotlib.axes.Axes.quiver for more info)
+
+.. python-script::
+
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   p = yt.ProjectionPlot(ds, 'z', 'density', center=[0.5, 0.5, 0.5],
+                         weight_field='density', width=(20, 'kpc'))
+   p.annotate_quiver('velocity_x', 'velocity_y', 16)
+   p.save()
+
+Off-axis Data Sources
+^^^^^^^^^^^^^^^^^^^^^
 
 .. function:: annotate_cquiver(self, field_x, field_y, factor):
 
-   (This is a proxy for :class:`~yt.visualization.plot_modifications.CuttingQuiverCallback`.)
+   (This is a proxy for
+   :class:`~yt.visualization.plot_modifications.CuttingQuiverCallback`.)
 
-   Get a quiver plot on top of a cutting plane, using
-   *field_x* and *field_y*, skipping every *factor*
-   datapoint in the discretization.
+   Get a quiver plot on top of a cutting plane, using ``field_x`` and
+   ``field_y``, skipping every ``factor`` datapoint in the discretization.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("Enzo_64/DD0043/data0043")
-   s = OffAxisSlicePlot(pf, [1,1,0], ["density"], center="c")
+   import yt
+   ds = yt.load("Enzo_64/DD0043/data0043")
+   s = yt.OffAxisSlicePlot(ds, [1,1,0], ["density"], center="c")
    s.annotate_cquiver('cutting_plane_velocity_x', 'cutting_plane_velocity_y', 10)
    s.zoom(1.5)
    s.save()
 
--------------
+Overplot grids
+~~~~~~~~~~~~~~
 
-.. function:: annotate_grids(self, alpha=1.0, min_pix=1, annotate=False, periodic=True):
+.. function:: annotate_grids(self, alpha=1.0, min_pix=1, annotate=False,
+                             periodic=True):
 
-   (This is a proxy for :class:`~yt.visualization.plot_modifications.GridBoundaryCallback`.)
+   (This is a proxy for
+   :class:`~yt.visualization.plot_modifications.GridBoundaryCallback`.)
 
-   Adds grid boundaries to a plot, optionally with
-   *alpha*-blending. Cuttoff for display is at *min_pix*
-   wide. *annotate* puts the grid id in the corner of the
-   grid.  (Not so great in projections...)
+   Adds grid boundaries to a plot, optionally with alpha-blending via the
+   ``alpha`` keyword. Cuttoff for display is at ``min_pix`` wide. ``annotate``
+   puts the grid id in the corner of the grid.  (Not so great in projections...)
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   slc = SlicePlot(pf, 'z', 'density', width=(10,'kpc'), center='max')
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   slc = yt.SlicePlot(ds, 'z', 'density', width=(10,'kpc'), center='max')
    slc.annotate_grids()
    slc.save()
 
--------------
+Overplot Halo Annotations
+~~~~~~~~~~~~~~~~~~~~~~~~~
 
-.. function:: annotate_halos(self, halo_catalog, col='white', alpha =1, width= None):
+.. function:: annotate_halos(self, halo_catalog, col='white', alpha =1,
+                             width=None):
 
-   (This is a proxy for :class:`~yt.visualization.plot_modifications.HaloCatalogCallback`.)
+   (This is a proxy for
+   :class:`~yt.visualization.plot_modifications.HaloCatalogCallback`.)
 
-   Accepts a :class:`yt.HaloCatalog` *HaloCatalog* and plots 
-   a circle at the location of each halo with the radius of
-   the circle corresponding to the virial radius of the halo.
-   If *width* is set to None (default) all halos are plotted.
-   Otherwise, only halos that fall within a slab with width
-   *width* centered on the center of the plot data. The 
-   color and transparency of the circles can be controlled with
-   *col* and *alpha* respectively.
+   Accepts a :class:`yt.HaloCatalog` and plots a circle at the location of each
+   halo with the radius of the circle corresponding to the virial radius of the
+   halo.  If ``width`` is set to None (default) all halos are plotted.
+   Otherwise, only halos that fall within a slab with width ``width`` centered
+   on the center of the plot data. The color and transparency of the circles can
+   be controlled with ``col`` and ``alpha`` respectively.
 
 .. python-script::
-   
-   from yt.mods import *
+
+   import yt
    from yt.analysis_modules.halo_analysis.halo_catalog import HaloCatalog
 
-   data_pf = load('Enzo_64/RD0006/RedshiftOutput0006')
-   halos_pf = load('rockstar_halos/halos_0.0.bin')
+   data_ds = yt.load('Enzo_64/RD0006/RedshiftOutput0006')
+   halos_ds = yt.load('rockstar_halos/halos_0.0.bin')
 
-   hc = HaloCatalog(halos_pf=halos_pf)
+   hc = HaloCatalog(halos_pf=halos_ds)
    hc.create()
 
-   prj = ProjectionPlot(data_pf, 'z', 'density')
+   prj = yt.ProjectionPlot(data_ds, 'z', 'density')
    prj.annotate_halos(hc)
    prj.save()
 
--------------
+Overplot a Straight Line
+~~~~~~~~~~~~~~~~~~~~~~~~
 
-.. function:: annotate_image_line(self, p1, p2, data_coords=False, plot_args=None):
+.. function:: annotate_image_line(self, p1, p2, data_coords=False,
+                                  plot_args=None):
 
-   (This is a proxy for :class:`~yt.visualization.plot_modifications.ImageLineCallback`.)
+   (This is a proxy for
+   :class:`~yt.visualization.plot_modifications.ImageLineCallback`.)
 
-   Plot from *p1* to *p2* (normalized image plane coordinates) with
-   *plot_args* fed into the plot.
+   Plot from ``p1`` to ``p2`` (normalized image plane coordinates) with
+  ``plot_args`` fed into the plot.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   p = ProjectionPlot(pf, 'z', 'density', center='m', width=(10, 'kpc'))
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   p = yt.ProjectionPlot(ds, 'z', 'density', center='m', width=(10, 'kpc'))
    p.annotate_image_line((0.3, 0.4), (0.8, 0.9), plot_args={'linewidth':5})
    p.save()
 
--------------
+Overplot a line plot
+~~~~~~~~~~~~~~~~~~~~
 
 .. function:: annotate_line(self, x, y, plot_args=None):
 
-   (This is a proxy for :class:`~yt.visualization.plot_modifications.LinePlotCallback`.)
+   (This is a proxy for
+   :class:`~yt.visualization.plot_modifications.LinePlotCallback`.)
 
-   Over plot *x* and *y* (in code units) with *plot_args* fed into the plot.
+   Over plot numpy arrays or lists ``x`` and ``y`` (in code units) with
+   ``plot_args`` fed into the plot.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   p = ProjectionPlot(pf, 'z', 'density', center='m', width=(10, 'kpc'))
-   p.annotate_line([-6, -4, -2, 0, 2, 4, 6], [3.6, 1.6, 0.4, 0, 0.4, 1.6, 3.6], plot_args={'linewidth':5})
+   import yt
+   import numpy as np
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   p = yt.ProjectionPlot(ds, 'z', 'density', center='m', width=(20, 'kpc'))
+   x = np.array([-6, -4, -2, 0, 2, 4, 6])
+   y = x**2/10
+   p.annotate_line(x, y, plot_args={'linewidth':5})
    p.save()
 
--------------
+Overplot Magnetic Field Quivers
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-.. function:: annotate_magnetic_field(self, factor=16, scale=None, scale_units=None, normalize=False):
+.. function:: annotate_magnetic_field(self, factor=16, scale=None,
+                                      scale_units=None, normalize=False):
 
-   (This is a proxy for :class:`~yt.visualization.plot_modifications.MagFieldCallback`.)
+   (This is a proxy for
+   :class:`~yt.visualization.plot_modifications.MagFieldCallback`.)
 
-   Adds a 'quiver' plot of magnetic field to the plot,
-   skipping all but every *factor* datapoint. *scale* is the
-   data units per arrow length unit using *scale_units* (see
-   matplotlib.axes.Axes.quiver for more info). if
-   *normalize* is True, the magnetic fields will be scaled
-   by their local (in-plane) length, allowing morphological
-   features to be more clearly seen for fields with
-   substantial variation in field strength.
+   Adds a 'quiver' plot of magnetic field to the plot, skipping all but every
+   ``factor`` datapoint. ``scale`` is the data units per arrow length unit using
+   ``scale_units`` (see matplotlib.axes.Axes.quiver for more info). if
+   ``normalize`` is ``True``, the magnetic fields will be scaled by their local
+   (in-plane) length, allowing morphological features to be more clearly seen
+   for fields with substantial variation in field strength.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("MHDSloshing/virgo_low_res.0054.vtk",
-             parameters={"TimeUnits":3.1557e13, "LengthUnits":3.0856e24,
-                         "DensityUnits":6.770424595218825e-27})
-   p = ProjectionPlot(pf, 'z', 'density', center='c', width=(300, 'kpc'))
+   import yt
+   ds = yt.load("MHDSloshing/virgo_low_res.0054.vtk",
+                parameters={"time_unit":(1, 'Myr'), "length_unit":(1, 'Mpc'),
+                            "mass_unit":(1e17, 'Msun')})
+   p = yt.ProjectionPlot(ds, 'z', 'density', center='c', width=(300, 'kpc'))
    p.annotate_magnetic_field()
    p.save()
 
--------------
+Annotate a Point With a Marker
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 .. function:: annotate_marker(self, pos, marker='x', plot_args=None):
 
-   (This is a proxy for :class:`~yt.visualization.plot_modifications.MarkerAnnotateCallback`.)
+   (This is a proxy for
+   :class:`~yt.visualization.plot_modifications.MarkerAnnotateCallback`.)
 
-   Adds text *marker* at *pos* in code coordinates.
-   *plot_args* is a dict that will be forwarded to the plot
+   Adds ``marker`` at ``pos`` in code coordinates.
+   ``plot_args`` is a dict that will be forwarded to the plot
    command.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   s = SlicePlot(pf, 'z', 'density', center='m', width=(10, 'kpc'))
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   s = yt.SlicePlot(ds, 'z', 'density', center='c', width=(10, 'kpc'))
    s.annotate_marker([0.5, 0.5, 0.5], plot_args={'s':10000})
-   s.save()   
+   s.save()
 
--------------
+Overplotting Particle Positions
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-.. function:: annotate_particles(self, width, p_size=1.0, col='k', marker='o', stride=1.0, ptype=None, stars_only=False, dm_only=False, minimum_mass=None):
+.. function:: annotate_particles(self, width, p_size=1.0, col='k', marker='o',
+                                 stride=1.0, ptype=None, stars_only=False,
+                                 dm_only=False, minimum_mass=None):
 
-   (This is a proxy for :class:`~yt.visualization.plot_modifications.ParticleCallback`.)
+   (This is a proxy for
+   :class:`~yt.visualization.plot_modifications.ParticleCallback`.)
 
-   Adds particle positions, based on a thick slab along
-   *axis* with a *width* along the line of sight.  *p_size*
-   controls the number of pixels per particle, and *col*
-   governs the color.  *ptype* will restrict plotted
-   particles to only those that are of a given type.
-   *minimum_mass* will require that the particles be of a
-   given mass, calculated via ParticleMassMsun, to be
-   plotted.
+   Adds particle positions, based on a thick slab along ``axis`` with a
+   ``width`` along the line of sight.  ``p_size`` controls the number of pixels
+   per particle, and ``col`` governs the color.  ``ptype`` will restrict plotted
+   particles to only those that are of a given type.  ``minimum_mass`` will
+   require that the particles be of a given mass minimum mass in solar units.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("Enzo_64/DD0043/data0043")
-   p = ProjectionPlot(pf, "x", "density", center='m', width=(10, 'Mpc'))
+   import yt
+   ds = yt.load("Enzo_64/DD0043/data0043")
+   p = yt.ProjectionPlot(ds, "x", "density", center='m', width=(10, 'Mpc'))
    p.annotate_particles((10, 'Mpc'))
    p.save()
 
--------------
+Annotate a point with text
+~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 .. function:: annotate_point(self, pos, text, text_args=None):
 
-   (This is a proxy for :class:`~yt.visualization.plot_modifications.PointAnnotateCallback`.)
+   (This is a proxy for
+   :class:`~yt.visualization.plot_modifications.PointAnnotateCallback`.)
 
-   This adds *text* at position *pos*, where *pos* is in
-   code-space. *text_args* is a dict fed to the text
-   placement code.
+   This adds ``text`` at position ``pos``, where ``pos`` is in
+   code-space. ``text_args`` is a dict fed to the text placement code.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   p = ProjectionPlot(pf, 'z', 'density', center='m', width=(10, 'kpc'))
-   p.annotate_point([0.5, 0.496, 0.5], "What's going on here?", text_args={'size':'xx-large', 'color':'w'})
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   p = yt.ProjectionPlot(ds, 'z', 'density', center='m', width=(10, 'kpc'))
+   p.annotate_point([0.5, 0.496, 0.5], "What's going on here?",
+                    text_args={'size':'xx-large', 'color':'w'})
    p.save()
 
--------------
+Overplot a circle on a plot
+~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-.. function:: annotate_quiver(self, field_x, field_y, factor, scale=None, scale_units=None, normalize=False):
+.. function:: annotate_sphere(self, center, radius, circle_args=None, text=None,
+                              text_args=None):
 
-   (This is a proxy for :class:`~yt.visualization.plot_modifications.QuiverCallback`.)
+   (This is a proxy for
+   :class:`~yt.visualization.plot_modifications.SphereCallback`.)
 
-   Adds a 'quiver' plot to any plot, using the *field_x* and
-   *field_y* from the associated data, skipping every
-   *factor* datapoints *scale* is the data units per arrow
-   length unit using *scale_units*  (see
-   matplotlib.axes.Axes.quiver for more info)
+   A sphere centered at ``center`` in code units with radius ``radius`` in code
+   units will be created, with optional ``circle_args``, ``text``, and
+   ``text_args``.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   p = ProjectionPlot(pf, 'z', 'density', center=[0.5, 0.5, 0.5], 
-                      weight_field='density', width=(20, 'kpc'))
-   p.annotate_quiver('velocity_x', 'velocity_y', 16)
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   p = yt.ProjectionPlot(ds, 'z', 'density', center='c', width=(20, 'kpc'))
+   p.annotate_sphere([0.5, 0.5, 0.5], (2, 'kpc'), {'fill':True})
    p.save()
 
--------------
+Overplot streamlines
+~~~~~~~~~~~~~~~~~~~~
 
-.. function:: annotate_sphere(self, center, radius, circle_args=None, text=None, text_args=None):
+.. function:: annotate_streamlines(self, field_x, field_y, factor=6.0, nx=16,
+                                   ny=16, xstart=(0, 1), ystart=(0, 1),
+                                   nsample=256, start_at_xedge=False,
+                                   start_at_yedge=False, plot_args=None):
 
-   (This is a proxy for :class:`~yt.visualization.plot_modifications.SphereCallback`.)
+   (This is a proxy for
+   :class:`~yt.visualization.plot_modifications.StreamlineCallback`.)
 
-   A sphere centered at *center* in code units with radius
-   *radius* in code units will be created, with optional
-   *circle_args*, *text*, and *text_args*.
+   Add streamlines to any plot, using the ``field_x`` and ``field_y`` from the
+   associated data, using ``nx`` and ``ny`` starting points that are bounded by
+   ``xstart`` and ``ystart``.  To begin streamlines from the left edge of the
+   plot, set ``start_at_xedge`` to ``True``; for the bottom edge, use
+   ``start_at_yedge``.  A line with the qmean vector magnitude will cover
+   1.0/``factor`` of the image.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   p = ProjectionPlot(pf, 'z', 'density', center='c', width=(20, 'kpc'))
-   p.annotate_sphere([0.5, 0.5, 0.5], (2, 'kpc'), {'fill':True})
-   p.save()
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   s = yt.SlicePlot(ds, 'z', 'density', center='c', width=(20, 'kpc'))
+   s.annotate_streamlines('velocity_x', 'velocity_y')
+   s.save()
 
--------------
+Add text
+~~~~~~~~
 
-.. function:: annotate_streamlines(self, field_x, field_y, factor=6.0, nx=16, ny=16, xstart=(0, 1), ystart=(0, 1), nsample=256, start_at_xedge=False, start_at_yedge=False, plot_args=None):
+.. function:: annotate_text(self, pos, text, data_coords=False, text_args=None):
 
-   (This is a proxy for :class:`~yt.visualization.plot_modifications.StreamlineCallback`.)
+   (This is a proxy for
+   :class:`~yt.visualization.plot_modifications.TextLabelCallback`.)
 
-   Add streamlines to any plot, using the *field_x* and
-   *field_y* from the associated data, using *nx* and *ny*
-   starting points that are bounded by *xstart* and
-   *ystart*.  To begin streamlines from the left edge of the
-   plot, set *start_at_xedge* to True; for the bottom edge,
-   use *start_at_yedge*.  A line with the qmean vector
-   magnitude will cover 1.0/*factor* of the image.
+   Accepts a position in (0..1, 0..1) of the image, some text and optionally
+   some text arguments. If data_coords is True, position will be in code units
+   instead of image coordinates.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   s = SlicePlot(pf, 'z', 'density', center='c', width=(20, 'kpc'))
-   s.annotate_streamlines('velocity_x', 'velocity_y')
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   s = yt.SlicePlot(ds, 'z', 'density', center='m', width=(10, 'kpc'))
+   s.annotate_text((0.5, 0.5), 'Sample text',
+                   text_args={'size':'xx-large', 'color':'w'})
    s.save()
 
--------------
+Add a title to the plot
+~~~~~~~~~~~~~~~~~~~~~~~
 
-.. function:: annotate_text(self, pos, text, data_coords=False, text_args=None):
+.. function:: annotate_title(self, title='Plot'):
 
-   (This is a proxy for :class:`~yt.visualization.plot_modifications.TextLabelCallback`.)
+   (This is a proxy for
+   :class:`~yt.visualization.plot_modifications.TitleCallback`.)
 
-   Accepts a position in (0..1, 0..1) of the image, some
-   text and optionally some text arguments. If data_coords
-   is True, position will be in code units instead of image
-   coordinates.
+   Accepts a ``title`` and adds it to the plot.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   s = SlicePlot(pf, 'z', 'density', center='m', width=(10, 'kpc'))
-   s.annotate_text((0.5, 0.5), 'Sample text', text_args={'size':'xx-large', 'color':'w'})
-   s.save()
-
--------------
-
-.. function:: annotate_title(self, title='Plot'):
-
-   (This is a proxy for :class:`~yt.visualization.plot_modifications.TitleCallback`.)
-
-   Accepts a *title* and adds it to the plot
-
-.. python-script::
-
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   p = ProjectionPlot(pf, 'z', 'density', center='c', width=(20, 'kpc'))
-   p.annotate_title('Density plot')
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   p = yt.ProjectionPlot(ds, 'z', 'density', center='c', width=(20, 'kpc'))
+   p.annotate_title('Density Plot')
    p.save()
 
--------------
+Overplot quivers for the velocity field
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-.. function:: annotate_velocity(self, factor=16, scale=None, scale_units=None, normalize=False):
+.. function:: annotate_velocity(self, factor=16, scale=None, scale_units=None,
+                                normalize=False):
 
-   (This is a proxy for :class:`~yt.visualization.plot_modifications.VelocityCallback`.)
+   (This is a proxy for
+   :class:`~yt.visualization.plot_modifications.VelocityCallback`.)
 
-   Adds a 'quiver' plot of velocity to the plot, skipping
-   all but every *factor* datapoint. *scale* is the data
-   units per arrow length unit using *scale_units* (see
-   matplotlib.axes.Axes.quiver for more info). if
-   *normalize* is True, the velocity fields will be scaled
-   by their local (in-plane) length, allowing morphological
-   features to be more clearly seen for fields with
-   substantial variation in field strength (normalize is not
+   Adds a 'quiver' plot of velocity to the plot, skipping all but every
+   ``factor`` datapoint. ``scale`` is the data units per arrow length unit using
+   ``scale_units`` (see matplotlib.axes.Axes.quiver for more info). if
+   ``normalize`` is ``True``, the velocity fields will be scaled by their local
+   (in-plane) length, allowing morphological features to be more clearly seen
+   for fields with substantial variation in field strength (normalize is not
    implemented and thus ignored for Cutting Planes).
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   p = SlicePlot(pf, 'z', 'density', center='m', width=(10, 'kpc'))
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   p = yt.SlicePlot(ds, 'z', 'density', center='m', width=(10, 'kpc'))
    p.annotate_velocity()
    p.save()
+
+Add a Timestamp Inset Box
+~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. function:: annotate_timestamp(x, y, units=None, format="{time:.3G} {units}", 
+                                 **kwargs, normalized=False, bbox_dict=None)
+
+   (This is a proxy for
+   :class:`~yt.visualization.plot_modifications.TimestampCallback`.)
+
+   Adds the current time to the plot at point given by *x* and *y*.  If *units*
+   is given ('s', 'ms', 'ns', etc), it will covert the time to this basis.  If
+   *units* is None, it will attempt to figure out the correct value by which to
+   scale.  The *format* keyword is a template string that will be evaluated and
+   displayed on the plot.  If *normalized* is true, *x* and *y* are interpreted
+   as normalized plot coordinates (0,0 is lower-left and 1,1 is upper-right)
+   otherwise *x* and *y* are assumed to be in plot coordinates. The *bbox_dict*
+   is an optional dict of arguments for the bbox that frames the timestamp, see
+   matplotlib's text annotation guide for more details. All other *kwargs* will
+   be passed to the text() method on the plot axes.  See matplotlib's text()
+   functions for more information.
+
+.. python-script::
+
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   p = yt.SlicePlot(ds, 'z', 'density', center='c', width=(20, 'kpc'))
+   p.annotate_timestamp()
+   p.save()

diff -r 84525632d209088422ef3d90d7bdddc3bba03f30 -r 3b264449bb411235ba00ab34a01df79261899f6b doc/source/visualizing/manual_plotting.rst
--- a/doc/source/visualizing/manual_plotting.rst
+++ b/doc/source/visualizing/manual_plotting.rst
@@ -35,11 +35,12 @@
 .. python-script::
    
    import pylab as P
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   import numpy as np
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
 
-   c = pf.h.find_max('density')[1]
-   proj = pf.proj('density', 0)
+   c = ds.find_max('density')[1]
+   proj = ds.proj('density', 0)
 
    width = (10, 'kpc') # we want a 1.5 mpc view
    res = [1000, 1000] # create an image with 1000x1000 pixels
@@ -64,22 +65,33 @@
 Line Plots
 ----------
 
-This is perhaps the simplest thing to do. ``yt`` provides a number of one dimensional objects, and these return a 1-D numpy array of their contents with direct dictionary access. As a simple example, take a :class:`~yt.data_objects.data_containers.AMROrthoRayBase` object, which can be created from a index by calling ``pf.ortho_ray(axis, center)``. 
+This is perhaps the simplest thing to do. ``yt`` provides a number of one
+dimensional objects, and these return a 1-D numpy array of their contents with
+direct dictionary access. As a simple example, take a
+:class:`~yt.data_objects.data_containers.AMROrthoRayBase` object, which can be
+created from a index by calling ``pf.ortho_ray(axis, center)``.
 
 .. python-script::
 
-   from yt.mods import *
+   import yt
+   import numpy as np
    import pylab as P
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   c = pf.h.find_max("density")[1]
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   c = ds.find_max("density")[1]
    ax = 0 # take a line cut along the x axis
-   ray = pf.ortho_ray(ax, (c[1], c[2])) # cutting through the y0,z0 such that we hit the max density
+
+   # cutting through the y0,z0 such that we hit the max density
+   ray = ds.ortho_ray(ax, (c[1], c[2]))
+
+   # Sort the ray values by 'x' so there are no discontinuities
+   # in the line plot
+   srt = np.argsort(ray['x'])
 
    P.subplot(211)
-   P.semilogy(np.array(ray['x']), np.array(ray['density']))
+   P.semilogy(np.array(ray['x'][srt]), np.array(ray['density'][srt]))
    P.ylabel('density')
    P.subplot(212)
-   P.semilogy(np.array(ray['x']), np.array(ray['temperature']))
+   P.semilogy(np.array(ray['x'][srt]), np.array(ray['temperature'][srt]))
    P.xlabel('x')
    P.ylabel('temperature')
 

diff -r 84525632d209088422ef3d90d7bdddc3bba03f30 -r 3b264449bb411235ba00ab34a01df79261899f6b doc/source/visualizing/plots.rst
--- a/doc/source/visualizing/plots.rst
+++ b/doc/source/visualizing/plots.rst
@@ -25,10 +25,10 @@
 If you need to take a quick look at a single simulation output, ``yt``
 provides the ``PlotWindow`` interface for generating annotated 2D
 visualizations of simulation data.  You can create a ``PlotWindow`` plot by
-supplying a parameter file, a list of fields to plot, and a plot center to
-create a :class:`~yt.visualization.plot_window.SlicePlot`,
+supplying a dataset, a list of fields to plot, and a plot center to
+create a :class:`~yt.visualization.plot_window.AxisAlignedSlicePlot`,
 :class:`~yt.visualization.plot_window.ProjectionPlot`, or
-:class:`~yt.visualization.plot_window.OffAxisSlicePlot`.
+:class:`~yt.visualization.plot_window.OffAxisProjectionPlot`.
 
 Plot objects use ``yt`` data objects to extract the maximum resolution
 data available to render a 2D image of a field. Whenever a
@@ -37,7 +37,8 @@
 is requested of it -- for instance, when the width or field is changed
 -- this high-resolution data is then pixelized and placed in a buffer
 of fixed size. This is accomplished behind the scenes using
-:class:`~yt.visualization.fixed_resolution.FixedResolutionBuffer`
+:class:`~yt.visualization.fixed_resolution.FixedResolutionBuffer`.
+
 ``PlotWindow`` expose the underlying matplotlib ``figure`` and
 ``axes`` objects, making it easy to customize your plots and 
 add new annotations.
@@ -45,28 +46,33 @@
 .. _slice-plots:
 
 Slice Plots
------------
+~~~~~~~~~~~
 
-The quickest way to plot a slice of a field through your data is to use
-:class:`~yt.visualization.plot_window.SlicePlot`.  Say we want to visualize a
-slice through the Density field along the z-axis centered on the center of the
-simulation box in a simulation dataset we've opened and stored in the parameter
-file object ``pf``.  This can be accomplished with the following command:
+The quickest way to plot a slice of a field through your data is via
+:class:`~yt.visualization.plot_window.SlicePlot`.  These plots are generally
+quicker than projections because they only need to read and process a slice
+through the dataset.
+
+The following script plots a slice through the density field along the z-axis
+centered on the center of the simulation box in a simulation dataset we've
+opened and stored in ``ds``:
 
 .. code-block:: python
 
-   >>> slc = SlicePlot(pf, 'z', 'density')
+   >>> slc = yt.SlicePlot(ds, 'z', 'density')
    >>> slc.save()
 
 These two commands will create a slice object and store it in a variable we've
-called ``slc``.  We then call the ``save()`` function that is associated with
-the slice object.  This automatically saves the plot in png image format with an
-automatically generated filename.  If you don't want the slice object to stick
-around, you can accomplish the same thing in one line:
+called ``slc``.  Since this plot is aligned with the simulation coordinate
+system, ``slc`` is an instance of
+:class:`~yt.visualization.plot_window.AxisAlignedSlicePlot`. We then call the
+``save()`` function, which automatically saves the plot in png image format with
+an automatically generated filename.  If you don't want the slice object to
+stick around, you can accomplish the same thing in one line:
 
 .. code-block:: python
    
-   >>> SlicePlot(pf, 'z', 'density').save()
+   >>> yt.SlicePlot(ds, 'z', 'density').save()
 
 It's nice to keep the slice object around if you want to modify the plot.  By
 default, the plot width will be set to the size of the simulation box.  To zoom
@@ -75,41 +81,58 @@
 
 .. code-block:: python
 
-   >>> slc = SlicePlot(pf, 'z', 'density')
+   >>> slc = yt.SlicePlot(ds, 'z', 'density')
    >>> slc.zoom(10)
    >>> slc.save('zoom')
 
 This will save a new plot to disk with a different filename - prepended with
-'zoom' instead of the name of the parameter file. If you want to set the width
+'zoom' instead of the name of the dataset. If you want to set the width
 manually, you can do that as well. For example, the following sequence of
 commands will create a slice, set the width of the plot to 10 kiloparsecs, and
 save it to disk.
 
 .. code-block:: python
 
-   >>> slc = SlicePlot(pf, 'z', 'density')
-   >>> slc.set_width((10,'kpc'))
+   >>> from yt.units import kpc
+   >>> slc = yt.SlicePlot(ds, 'z', 'density')
+   >>> slc.set_width(10*kpc)
    >>> slc.save('10kpc')
 
-The SlicePlot also optionally accepts the coordinate to center the plot on and
-the width of the plot:
+The plot width can be specified independently along the x and y direction by
+passing a tuple of widths.  An individual width can also be represented using a
+``(value, unit)`` tuple.  The following sequence of commands all equivalently
+set the width of the plot to 200 kiloparsecs in the ``x`` and ``y`` direction.
 
 .. code-block:: python
 
-   >>> SlicePlot(pf, 'z', 'density', center=[0.2, 0.3, 0.8], 
-   ...           width = (10,'kpc')).save()
+   >>> from yt.units import kpc
+   >>> slc.set_width(200*kpc)
+   >>> slc.set_width((200, 'kpc'))
+   >>> slc.set_width((200*kpc, 200*kpc))
 
-The center must be given in code units.  Optionally, you can supply 'c' or 'm'
-for the center.  These two choices will center the plot on the center of the
-simulation box and the coordinate of the maximum density cell, respectively.
+The ``SlicePlot`` also optionally accepts the coordinate to center the plot on
+and the width of the plot:
+
+.. code-block:: python
+
+   >>> yt.SlicePlot(ds, 'z', 'density', center=[0.2, 0.3, 0.8],
+   ...              width = (10,'kpc')).save()
+
+The plot center is relative to the simulation coordinate system.  If supplied
+without units, the center is assumed by in code units.  Optionally, you can
+supply 'c' or 'm' for the center.  These two choices will center the plot on the
+center of the simulation box and the coordinate of the maximum density cell,
+respectively.
 
 Here is an example that combines all of the options we just discussed.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   slc = SlicePlot(pf, 'z', 'density', center=[0.5, 0.5, 0.5], width=(20,'kpc'))
+   import yt
+   from yt.units import kpc
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   slc = yt.SlicePlot(ds, 'z', 'density', center=[0.5, 0.5, 0.5],
+                      width=(20,'kpc'))
    slc.save()
 
 The above example will display an annotated plot of a slice of the
@@ -118,14 +141,14 @@
 letter 'z', corresponding to the z-axis.  Finally, the image is saved to
 a png file.
 
-Conceptually, you can think of the SlicePlot as an adjustable window
+Conceptually, you can think of the plot object as an adjustable window
 into the data. For example:
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   slc = SlicePlot(pf, 'z', 'pressure', center='c')
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   slc = yt.SlicePlot(ds, 'z', 'pressure', center='c')
    slc.save()
    slc.zoom(30)
    slc.save('zoom')
@@ -138,32 +161,43 @@
 interesting region in the simulation and adjust the boundaries of the
 region to visualize on the fly.
 
-A slice object can also add annotations like a title, an overlying
-quiver plot, the location of grid boundaries, halo-finder annotations,
-and many other annotations, including user-customizable annotations.
-For example:
+See :class:`~yt.visualization.plot_window.AxisAlignedSlicePlot` for the 
+full class description.
+
+.. _off-axis-slices:
+
+Off Axis Slices
+~~~~~~~~~~~~~~~
+
+Off axis slice plots can be generated in much the same way as
+grid-aligned slices.  Off axis slices use
+:class:`~yt.data_objects.data_containers.AMRCuttingPlaneBase` to slice
+through simulation domains at an arbitrary oblique angle.  A
+:class:`~yt.visualization.plot_window.OffAxisSlicePlot` can be
+instantiated by specifying a dataset, the normal to the cutting
+plane, and the name of the fields to plot.  Just like an
+:class:`~yt.visualization.plot_window.AxisAlignedSlicePlot`, an
+``OffAxisSlicePlot`` can be created via the
+:class:`~yt.visualization.plot_window.SlicePlot` function. For example:
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   slc = SlicePlot(pf, 'z', 'density', width=(10,'kpc'))
-   slc.annotate_grids()
-   slc.save()
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   L = [1,1,0] # vector normal to cutting plane
+   north_vector = [-1,1,0]
+   cut = yt.SlicePlot(ds, L, 'density', width=(25, 'kpc'),
+                      north_vector=north_vector)
+   cut.save()
 
-will plot the density field in a 10 kiloparsec slice through the
-z-axis centered on the highest density point in the simulation domain.
-Before saving the plot, the script annotates it with the grid
-boundaries, which are drawn as thick black lines by default.
-
-Annotations are described in :ref:`callbacks`.  See
-:class:`~yt.visualization.plot_window.SlicePlot` for the full class
-description.
+In this case, a normal vector for the cutting plane is supplied in the second
+argument. Optionally, a `north_vector` can be specified to fix the orientation
+of the image plane.
 
 .. _projection-plots:
 
 Projection Plots
-----------------
+~~~~~~~~~~~~~~~~
 
 Using a fast adaptive projection, ``yt`` is able to quickly project
 simulation data along the coordinate axes.
@@ -174,15 +208,15 @@
 
 .. python-script::
  
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   prj = ProjectionPlot(pf, 2, 'density', width=(25, 'kpc'), 
-                        weight_field=None)
+   import yt
+   from yt.units import kpc
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   prj = yt.ProjectionPlot(ds, 2, 'temperature', width=25*kpc,
+                           weight_field='density')
    prj.save()
 
-will create a projection of Density field along the x axis, plot it,
-and then save it to a png image file.  The projection is only carried
-out to level 2 of the AMR index and no weighting is applied.
+will create a density-weighted projection of the temperature field along the x
+axis, plot it, and then save the plot to a png image file.
 
 Like :ref:`slice-plots`, annotations and modifications can be applied
 after creating the ``ProjectionPlot`` object.  Annotations are
@@ -190,43 +224,8 @@
 :class:`~yt.visualization.plot_window.ProjectionPlot` for the full
 class description.
 
-.. _off-axis-slices:
-
-Off Axis Slice Plots
---------------------
-
-Off axis slice plots can be generated in much the same way as
-grid-aligned slices.  Off axis slices use
-:class:`~yt.data_objects.data_containers.AMRCuttingPlaneBase` to slice
-through simulation domains at an arbitrary oblique angle.  A
-:class:`~yt.visualization.plot_window.OffAxisSlicePlot` can be
-instantiated by specifying a parameter file, the normal to the cutting
-plane, and the name of the fields to plot.  For example:
-
-.. python-script::
-
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   L = [1,1,0] # vector normal to cutting plane
-   north_vector = [-1,1,0]
-   cut = OffAxisSlicePlot(pf, L, 'density', width=(25, 'kpc'), 
-                          north_vector=north_vector)
-   cut.save()
-
-creates an off-axis slice in the plane perpendicular to ``L``,
-oriented such that ``north_vector`` is the up direction.  If ``L`` and
-``north_vector`` are not perpendicular, the component of
-``north_vector`` perpendicular to ``L`` is used. Like
-:ref:`slice-plots`, annotations and modifications can be applied after
-creating the ``OffAxisSlicePlot`` object.  Annotations are described
-in :ref:`callbacks`.  See
-:class:`~yt.visualization.plot_window.OffAxisSlicePlot` for the full
-class description.
-
-.. _off-axis-projections:
-
 Off Axis Projection Plots
--------------------------
+~~~~~~~~~~~~~~~~~~~~~~~~~
 
 Off axis projection plots .  Internally, off axis projections are
 created using :ref:`the-camera-interface` by applying the
@@ -246,15 +245,16 @@
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   import yt
+   import numpy as np
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
    L = [1,1,0] # vector normal to cutting plane
    north_vector = [-1,1,0]
    W = [0.02, 0.02, 0.02]
    c = [0.5, 0.5, 0.5]
    N = 512
-   image = off_axis_projection(pf, c, L, W, N, "density")
-   write_image(np.log10(image), "%s_offaxis_projection.png" % pf)
+   image = yt.off_axis_projection(ds, c, L, W, N, "density")
+   yt.write_image(np.log10(image), "%s_offaxis_projection.png" % ds)
 
 Here, ``W`` is the width of the projection in the x, y, *and* z
 directions.
@@ -268,32 +268,30 @@
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
    L = [1,1,0] # vector normal to cutting plane
    north_vector = [-1,1,0]
-   prj = OffAxisProjectionPlot(pf,L,'density',width=(25, 'kpc'), 
-                               north_vector=north_vector)
+   prj = yt.OffAxisProjectionPlot(ds,L,'density',width=(25, 'kpc'),
+                                  north_vector=north_vector)
    prj.save()
 
 OffAxisProjectionPlots can also be created with a number of
-keyword arguments, as described in the `api reference`__ for the
-class initializer.
-
-__ :class:`~yt.visualization.plot_window.OffAxisProjectionPlot`
+keyword arguments, as described in
+:class:`~yt.visualization.plot_window.OffAxisProjectionPlot`
 
 Plot Customization
 ------------------
 
 You can customize each of the four plot types above in identical ways.  We'll go
-over each of the customizations methos below.  For each of the examples below we
+over each of the customizations methods below.  For each of the examples below we
 will modify the following plot.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   slc = SlicePlot(pf, 'z', 'density', width=(10,'kpc'))
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   slc = yt.SlicePlot(ds, 'z', 'density', width=(10,'kpc'))
    slc.save()
 
 Panning and zooming
@@ -301,78 +299,84 @@
 
 There are three methods to dynamically pan around the data.  
 
-:class:`~yt.visualization.plot_window.SlicePlot.pan` accepts x and y deltas in code
-units.
+:class:`~yt.visualization.plot_window.AxisAlignedSlicePlot.pan` accepts x and y
+deltas.
 
 .. python-script::
 
-   from yt.mods import *
+   import yt
    from yt.units import kpc
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   slc = SlicePlot(pf, 'z', 'density', width=(10,'kpc'))
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   slc = yt.SlicePlot(ds, 'z', 'density', width=(10,'kpc'))
    slc.pan((2*kpc, 2*kpc))
    slc.save()
 
-:class:`~yt.visualization.plot_window.SlicePlot.pan_rel` accepts deltas in units relative
-to the field of view of the plot.  
+:class:`~yt.visualization.plot_window.AxisAlignedSlicePlot.pan_rel` accepts deltas 
+in units relative to the field of view of the plot.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   slc = SlicePlot(pf, 'z', 'density', width=(10,'kpc'))
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   slc = yt.SlicePlot(ds, 'z', 'density', width=(10,'kpc'))
    slc.pan_rel((0.1, -0.1))
    slc.save()
 
-:class:`~yt.visualization.plot_window.SlicePlot.zoom` accepts a factor to zoom in by.
+:class:`~yt.visualization.plot_window.AxisAlignedSlicePlot.zoom` accepts a factor to zoom in by.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   slc = SlicePlot(pf, 'z', 'density', width=(10,'kpc'))
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   slc = yt.SlicePlot(ds, 'z', 'density', width=(10,'kpc'))
    slc.zoom(2)
    slc.save()
 
 Set axes units
 ~~~~~~~~~~~~~~
 
-:class:`~yt.visualization.plot_window.SlicePlot.set_axes_unit` allows the customization of
+:class:`~yt.visualization.plot_window.AxisAlignedSlicePlot.set_axes_unit` allows the customization of
 the axes unit labels.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   slc = SlicePlot(pf, 'z', 'density', width=(10,'kpc'))
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   slc = yt.SlicePlot(ds, 'z', 'density', width=(10,'kpc'))
    slc.set_axes_unit('Mpc')
    slc.save()
 
+The same result could have been accomplished by explicitly setting the ``width``
+to ``(.01, 'Mpc')``.
+
 Set the plot center
 ~~~~~~~~~~~~~~~~~~~
 
-The :class:`~yt.visualization.plot_window.SlicePlot.set_center` function accepts a new
-center for the plot, in code units.  New centers must be two element tuples.
+The :class:`~yt.visualization.plot_window.AxisAlignedSlicePlot.set_center`
+function accepts a new center for the plot, in code units.  New centers must be
+two element tuples.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   slc = SlicePlot(pf, 'z', 'density', width=(10,'kpc'))
-   slc.set_center((0.5, 0.5))
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   slc = yt.SlicePlot(ds, 'z', 'density', width=(10,'kpc'))
+   slc.set_center((0.5, 0.503))
    slc.save()
 
 Fonts
 ~~~~~
 
-:class:`~yt.visualization.plot_window.SlicePlot.set_font` allows font costomization.
+:class:`~yt.visualization.plot_window.AxisAlignedSlicePlot.set_font` allows font
+costomization.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   slc = SlicePlot(pf, 'z', 'density', width=(10,'kpc'))
-   slc.set_font({'family': 'sans-serif', 'style': 'italic','weight': 'bold', 'size': 24})
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   slc = yt.SlicePlot(ds, 'z', 'density', width=(10,'kpc'))
+   slc.set_font({'family': 'sans-serif', 'style': 'italic',
+                 'weight': 'bold', 'size': 24})
    slc.save()
 
 Colormaps
@@ -383,117 +387,173 @@
 different fields tracked by the plot object.
 
 To change the colormap for the plot, call the
-:class:`~yt.visualization.plot_window.SlicePlot.set_cmap` function.  Use any of the
-colormaps listed in the :ref:`colormaps` section.
+:class:`~yt.visualization.plot_window.AxisAlignedSlicePlot.set_cmap` function.
+Use any of the colormaps listed in the :ref:`colormaps` section.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   slc = SlicePlot(pf, 'z', 'density', width=(10,'kpc'))
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   slc = yt.SlicePlot(ds, 'z', 'density', width=(10,'kpc'))
    slc.set_cmap('density', 'RdBu_r')
    slc.save()
 
-The :class:`~yt.visualization.plot_window.SlicePlot.set_log` function accepts a field name
-and a boolean.  If the boolean is :code:`True`, the colormap for the field will
-be log scaled.  If it is `False` the colormap will be linear.
+The :class:`~yt.visualization.plot_window.AxisAlignedSlicePlot.set_log` function
+accepts a field name and a boolean.  If the boolean is ``True``, the colormap
+for the field will be log scaled.  If it is ``False`` the colormap will be
+linear.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   slc = SlicePlot(pf, 'z', 'density', width=(10,'kpc'))
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   slc = yt.SlicePlot(ds, 'z', 'density', width=(10,'kpc'))
    slc.set_log('density', False)
    slc.save()
 
-Lastly, the :class:`~yt.visualization.plot_window.SlicePlot.set_zlim` function makes it
-possible to set a custom colormap range.
+Lastly, the :class:`~yt.visualization.plot_window.AxisAlignedSlicePlot.set_zlim`
+function makes it possible to set a custom colormap range.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   slc = SlicePlot(pf, 'z', 'density', width=(10,'kpc'))
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   slc = yt.SlicePlot(ds, 'z', 'density', width=(10,'kpc'))
    slc.set_zlim('density', 1e-30, 1e-25)
    slc.save()
 
+Annotations
+~~~~~~~~~~~
+
+A slice object can also add annotations like a title, an overlying
+quiver plot, the location of grid boundaries, halo-finder annotations,
+and many other annotations, including user-customizable annotations.
+For example:
+
+.. python-script::
+
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   slc = yt.SlicePlot(ds, 'z', 'density', width=(10,'kpc'))
+   slc.annotate_grids()
+   slc.save()
+
+will plot the density field in a 10 kiloparsec slice through the
+z-axis centered on the highest density point in the simulation domain.
+Before saving the plot, the script annotates it with the grid
+boundaries, which are drawn as thick black lines by default.
+
+Annotations are described in :ref:`callbacks`.
+
 Set the size of the plot
 ~~~~~~~~~~~~~~~~~~~~~~~~
 
 To set the size of the plot, use the
-:class:`~yt.visualization.plot_window.SlicePlot.set_window_size` function.  The argument
+:class:`~yt.visualization.plot_window.AxisAlignedSlicePlot.set_figure_size` function.  The argument
 is the size of the longest edge of the plot in inches.  View the full resolution
 image to see the difference more clearly.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   slc = SlicePlot(pf, 'z', 'density', width=(10,'kpc'))
-   slc.set_window_size(10)
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   slc = yt.SlicePlot(ds, 'z', 'density', width=(10,'kpc'))
+   slc.set_figure_size(10)
    slc.save()
 
 To change the resolution of the image, call the
-:class:`~yt.visualization.plot_window.SlicePlot.set_buff_size` function.
+:class:`~yt.visualization.plot_window.AxisAlignedSlicePlot.set_buff_size` function.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   slc = SlicePlot(pf, 'z', 'density', width=(10,'kpc'))
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   slc = yt.SlicePlot(ds, 'z', 'density', width=(10,'kpc'))
    slc.set_buff_size(1600)
    slc.save()
 
+Further customization via matplotlib
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Each ``PlotWindow`` object is really a container for plots - one plot for each
+field specified in the list of fields supplied when the plot object is
+created. The individual plots can be accessed via the ``plots`` dictionary
+attached to each ``PlotWindow`` object:
+
+.. code-block:: python
+
+   >>> slc = SlicePlot(ds, 2, ['density', 'temperature']
+   >>> dens_plot = slc.plots['density']
+
+In this example ``dens_plot`` is an instance of
+:class:`~yt.visualization.plot_window.WindowPlotMPL`, an object that wraps the
+matplotlib ``figure`` and ``axes`` objects.  We can access these matplotlib primitives via attributes of ``dens_plot``.
+
+.. code-block:: python
+
+   >>> figure = dens_plot.figure
+   >>> axes = dens_plot.axes
+   >>> colorbar_axes = dens_plot.cax
+
+These are the :matplotlib:class:`figure`, and :matplotlib:class:`axes` objects
+that control the actual drawing of the plot.  Arbitrary plot customizations
+are possible by manipulating these objects.  See :ref:`matplotlib-primitives` for
+an example.
+
 .. _how-to-make-1d-profiles:
 
 1D Profile Plots
 ----------------
 
-1D profiles are used to calculated the average or the sum of a given quantity
-with respect to a second quantity.  This means the "average density as a
-function of radius" or "the total mass within a given set of density bins."
+1D profiles are used to calculate the average or the sum of a given quantity
+with respect to a second quantity.  Two common examples are the "average density
+as a function of radius" or "the total mass within a given set of density bins."
 When created, they default to the average: in fact, they default to the average
 as weighted by the total cell mass.  However, this can be modified to take
 either the total value or the average with respect to a different quantity.
 
-Profiles operate on data objects; they will take the entire data contained in a
-sphere, a prism, an extracted region and so on, and they will calculate and use
-that as input to their calculation.  To make a 1D profile plot, create a
-(:class:`~yt.visualization.profile_plotter.ProfilePlot`) object, supplying the 
-data object, the field for binning, and a list of fields to be profiled.
+Profiles operate on :ref:`data objects <using-objects>`; they will take the
+entire data contained in a sphere, a prism, an extracted region and so on, and
+they will calculate and use that as input to their calculation.  To make a 1D
+profile plot, create a (:class:`~yt.visualization.profile_plotter.ProfilePlot`)
+object, supplying the data object, the field for binning, and a list of fields
+to be profiled.
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   my_galaxy = pf.disk([0.5, 0.5, 0.5], [0.0, 0.0, 1.0], 0.01, 0.003)
-   plot = ProfilePlot(my_galaxy, "density", ["temperature"])
+   import yt
+   from yt.units import kpc
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   my_galaxy = ds.disk(ds.domain_center, [0.0, 0.0, 1.0], 10*kpc, 3*kpc)
+   plot = yt.ProfilePlot(my_galaxy, "density", ["temperature"])
    plot.save()
 
 This will create a :class:`yt.data_objects.data_containers.AMRCylinderBase`
-centered at [0.3, 0.5, 0.8], with a normal vector of [0.4, 0.5, 0.1], radius of
-0.01 and height of 0.001 and will then make a plot of the average (as a 
-function of the cell mass) temperature as a function of density.
+centered at [0.5, 0.5, 0.5], with a normal vector of [0.0, 0.0, 1.0], radius of
+10 kiloparsecs and height of 3 kiloparsecs and will then make a plot of the
+mass-weighted average temperature as a function of density for all of the gas
+contained in the cylinder.
 
-As another example, we create a sphere of radius 100 pc and plot total mass 
-in every equally-spaced temperature bin/
-
-We could also have allowed the plot collection to create a sphere for us, as
-well.  For instance:
+We could also have made a profile considering only the gas in a sphere.
+For instance:
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   my_sphere = pf.sphere([0.5, 0.5, 0.5], (100, "kpc"))
-   plot = ProfilePlot(my_sphere, "temperature", ["cell_mass"],
-                      weight_field=None)
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   my_sphere = ds.sphere([0.5, 0.5, 0.5], (100, "kpc"))
+   plot = yt.ProfilePlot(my_sphere, "temperature", ["cell_mass"],
+                         weight_field=None)
    plot.save()
 
-Note that because we have specified the weighting field to be none, it operates 
-as a local-bin accumulator.  We can also accumulate along the x-axis by setting 
-the **accumulation** keyword argument to True, which is useful for plots of 
-enclosed mass.
+Note that because we have specified the weighting field to be none, it operates
+the profile plot will display the accumulated cell mass as a function of
+temperature rather than the average.  We can also accumulate along the x-axis by
+setting the **accumulation** keyword argument to True, which is useful for plots
+of enclosed mass.
+
+Also note the use of a ``(value, unit)`` tuple. These can be used interchangably
+with units explicitly imported from ``yt.units``.
 
 You can also access the data generated by profiles directly, which can be
 useful for overplotting average quantities on top of phase plots, or for
@@ -506,32 +566,30 @@
 
    plot = ProfilePlot(my_sphere, "temperature", ["cell_mass"],
                       weight_field=None)
-   # print the x field
+   profile = plot.profiles[0]
+   # print the bin field, in this case temperature
    print plot.profiles[-1].x
-   # print the profiled temperature field
-   print plot.profiles[-1].field_data["temperature"]
+   # print the profiled cell_mass field
+   print plot.profiles[-1]["cell_mass"]
 
-Other options, such as the number of bins, are also configurable. See the 
-documentation for 
-The number of bins and other options and tweaks are 
-available for these methods.  See the documentation for 
-:class:`~yt.visualization.profile_plotter.ProfilePlot`
-for more information.
+Other options, such as the number of bins, are also configurable. See the
+documentation for :class:`~yt.visualization.profile_plotter.ProfilePlot` for
+more information.
 
 Overplotting Multiple 1D Profiles
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 It is often desirable to overplot multiple 1D profile to show evolution 
 with time.  This is supported with the ``from_profiles`` class method.  
-1D profiles are created with the :meth:`yt.data_objects.profiles.create_profile` 
+1D profiles are created with the :func:`yt.data_objects.profiles.create_profile` 
 method and then given to the ProfilePlot object.
 
 .. python-script::
 
-   from yt.mods import *
+   import yt
 
    # Create a time-series object.
-   es = simulation("enzo_tiny_cosmology/32Mpc_32.enzo", "Enzo")
+   es = yt.simulation("enzo_tiny_cosmology/32Mpc_32.enzo", "Enzo")
    es.get_time_series(redshifts=[5, 4, 3, 2, 1, 0])
 
 
@@ -540,23 +598,114 @@
    labels = []
 
    # Loop over each dataset in the time-series.
-   for pf in es:
+   for ds in es:
        # Create a data container to hold the whole dataset.
-       ad = pf.h.all_data()
+       ad = ds.all_data()
        # Create a 1d profile of density vs. temperature.
-       profiles.append(create_profile(ad, ["temperature"], 
-                                      fields=["cell_mass"],
-                                      weight_field=None,
-                                      accumulation=True))
+       profiles.append(yt.create_profile(ad, ["temperature"], 
+                                         fields=["cell_mass"],
+                                         weight_field=None,
+                                         accumulation=True))
        # Add labels
-       labels.append("z = %.2f" % pf.current_redshift)
+       labels.append("z = %.2f" % ds.current_redshift)
 
    # Create the profile plot from the list of profiles.
-   plot = ProfilePlot.from_profiles(profiles, labels=labels)
+   plot = yt.ProfilePlot.from_profiles(profiles, labels=labels)
 
    # Save the image.
    plot.save()
 
+Customizing axis limits
+~~~~~~~~~~~~~~~~~~~~~~~
+
+By default the x and y limits for ``ProfilePlot`` are determined using the
+:class:`~yt.data_objects.derived_quantities.Extrema` derived quantity.  If you
+want to create a plot with custom axis limits, you have two options.
+
+First, you can create a custom profile object using
+:func:`~yt.data_objects.profiles.create_profile`.  This function accepts a dictionary of ``(max, min)`` tuples keyed to field names.
+
+.. python-script::
+
+    import yt
+    import yt.units as u
+    ds = yt.load('IsolatedGalaxy/galaxy0030/galaxy0030')
+    sp = ds.sphere('m', 10*u.kpc)
+    profiles = yt.create_profile(sp, "temperature", "density",
+                                 weight_field=None, 
+                                 extrema={'temperature': (1e3, 1e7),
+                                          'density': (1e-26, 1e-22)})
+    plot = yt.ProfilePlot.from_profiles(profiles)
+    plot.save()
+
+You can also make use of the
+:meth:`yt.visualization.profile_plotter.ProfilePlot.set_xlim` and
+:meth:`yt.visualization.profile_plotter.ProfilePlot.set_ylim` functions to
+customize the axes limits of a plot that has already been created.  Note that
+calling ``set_xlim`` is much slower than calling ``set_ylim``.  This is because
+```set_xlim`` must recreate the profile object using the specified extrema.
+Creating a profile directly via ``create_profile`` might be significant faster.
+If you find that this operation is slow, consider creating your profile via
+``create_profile``.  Note that since there is only one bin field, ``set_xlim``
+does not accept a field name as the first argument.
+
+.. python-script::
+
+   import yt
+   import yt.units as u
+   ds = yt.load('IsolatedGalaxy/galaxy0030/galaxy0030')
+   sp = ds.sphere('m', 10*u.kpc)
+   plot = yt.ProfilePlot(sp, "temperature", "density", weight_field=None)
+   plot.set_xlim(1e3, 1e7)
+   plot.set_ylim("density", 1e-26, 1e-22)
+   plot.save()
+
+
+Customizing Units
+~~~~~~~~~~~~~~~~~
+
+Units for both the x and y axis can be controlled via the
+:meth:`~yt.visualization.profile_plotter.ProfilePlot.set_unit` method.
+Adjusting the plot units does not require recreating the histogram, so adjusting
+units will always be inexpensive, requiring only an in-place unit conversion.
+
+In the following example we create a a plot of the average density in solar
+masses per cubic parsec as a function of radius in kiloparsecs.
+
+.. python-script::
+
+    import yt
+    import yt.units as u
+    ds = yt.load('IsolatedGalaxy/galaxy0030/galaxy0030')
+    sp = ds.sphere('m', 10*u.kpc)
+    plot = yt.ProfilePlot(sp, "radius", "density", weight_field=None)
+    plot.set_unit("density", "msun/pc**3")
+    plot.set_unit("radius", "kpc")
+    plot.save()
+
+Linear and Logarithmic Scaling
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The axis scaling can be manipulated via the
+:meth:`~yt.visualization.profile_plotter.ProfilePlot.set_log` function.  This
+function accepts a field name and a boolean.  If the boolean is ``True``, the
+field is plotted in log scale.  If ``False``, the field is plotted in linear
+scale.
+
+In the following example we create a plot of the average x velocity as a
+function of radius.  Since the x component of the velocity vector can be
+negative, we set the scaling to be linear for this field.
+
+.. python-script::
+
+   import yt
+   import yt.units as u
+   ds = yt.load('IsolatedGalaxy/galaxy0030/galaxy0030')
+   sp = ds.sphere('m', 10*u.kpc)
+   plot = yt.ProfilePlot(sp, "radius", "x-velocity", weight_field=None)
+   plot.set_log("x-velocity", False)
+   plot.save()
+
 Altering Line Properties
 ~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -581,47 +730,72 @@
 2D Phase Plots
 --------------
 
-2D phase plots function in much the same was as 1D phase plots, but with a 
-:class:`~yt.visualization.profile_plotter.PhasePlot` object.  Much like 1D profiles, 
-2D profiles (phase plots) are best thought of as plotting a distribution of points, 
-either taking the average or the accumulation in a bin.  For example, to generate a 
-2D distribution of mass enclosed in density and temperature bins, you can do:
+2D phase plots function in much the same was as 1D phase plots, but with a
+:class:`~yt.visualization.profile_plotter.PhasePlot` object.  Much like 1D
+profiles, 2D profiles (phase plots) are best thought of as plotting a
+distribution of points, either taking the average or the accumulation in a bin.
+For example, to generate a 2D distribution of mass enclosed in density and
+temperature bins, you can do:
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   my_sphere = pf.sphere("c", (50, "kpc"))
-   plot = PhasePlot(my_sphere, "density", "temperature", ["cell_mass"],
-                    weight_field=None)
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   my_sphere = ds.sphere("c", (50, "kpc"))
+   plot = yt.PhasePlot(my_sphere, "density", "temperature", ["cell_mass"],
+                       weight_field=None)
    plot.save()
 
 If you would rather see the average value of a field as a function of two other
-fields, you can neglect supplying the *weight* parameter.  This would look
+fields, you can set the ``weight_field`` parameter to ``None``.  This would look
 something like:
 
 .. python-script::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   my_sphere = pf.sphere("c", (50, "kpc"))
-   plot = PhasePlot(my_sphere, "density", "temperature", ["H_fraction"],
-                    weight_field="cell_mass")
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   my_sphere = ds.sphere("c", (50, "kpc"))
+   plot = yt.PhasePlot(my_sphere, "density", "temperature", ["H_fraction"],
+                       weight_field=None)
+   plot.save()
+
+Customizing Phase Plots
+~~~~~~~~~~~~~~~~~~~~~~~
+
+Similarly to 1D profile plots, ``PhasePlot`` can be customized via ``set_unit``,
+``set_xlim``, ``set_ylim``, and ``set_zlim``.  The following example illustrates
+how to manipulate these functions.
+
+.. python-script::
+
+   import yt
+   ds = yt.load("sizmbhloz-clref04SNth-rs9_a0.9011/sizmbhloz-clref04SNth-rs9_a0.9011.art")
+   center = ds.arr([64.0, 64.0, 64.0], 'code_length')
+   rvir = ds.quan(1e-1, "Mpccm/h")
+
+   sph = ds.sphere(center, rvir)
+   plot = yt.PhasePlot(sph, "density", "temperature", "cell_mass",
+                       weight_field=None)
+   plot.set_unit('density', 'Msun/pc**3')
+   plot.set_unit('cell_mass', 'Msun')
+   plot.set_xlim(1e-5,1e1)
+   plot.set_ylim(1,1e7)
    plot.save()
 
 Probability Distribution Functions and Accumulation
 ---------------------------------------------------
 
-Both 1D and 2D profiles which show the total of amount of some field, such as mass, 
-in a bin (done by setting the **weight_field** keyword to None) can be turned into 
-probability distribution functions (PDFs) by setting the **fractional** keyword to 
-True.  When set to True, the value in each bin is divided by the sum total from all 
-bins.  These can be turned into cumulative distribution functions (CDFs) by setting 
-the **accumulation** keyword to True.  This will make is so that the value in any bin 
-N is the cumulative sum of all bins from 0 to N.  The direction of the summation can be 
-rversed by setting **accumulation** to -True.  For PhasePlots, the accumulation can be 
-set independently for each axis by setting **accumulation** to a list of True/-True/False 
-values.
+Both 1D and 2D profiles which show the total of amount of some field, such as
+mass, in a bin (done by setting the ``weight_field`` keyword to ``None``) can be
+turned into probability distribution functions (PDFs) by setting the
+``fractional`` keyword to ``True``.  When set to ``True``, the value in each bin
+is divided by the sum total from all bins.  These can be turned into cumulative
+distribution functions (CDFs) by setting the ``accumulation`` keyword to
+``True``.  This will make is so that the value in any bin N is the cumulative
+sum of all bins from 0 to N.  The direction of the summation can be reversed by
+setting ``accumulation`` to ``-True``.  For ``PhasePlot``, the accumulation can
+be set independently for each axis by setting ``accumulation`` to a list of
+``True``/ ``-True`` /``False`` values.
 
 .. _interactive-plotting:
 
@@ -646,9 +820,9 @@
 
 .. notebook-cell::
 
-   from yt.mods import *
-   pf = load("IsolatedGalaxy/galaxy0030/galaxy0030")
-   p = ProjectionPlot(pf, "x", "density", center='m', width=(10,'kpc'),
+   import yt
+   ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+   p = ProjectionPlot(ds, "x", "density", center='m', width=(10,'kpc'),
                       weight_field='density')
    p.show()
 
@@ -682,7 +856,7 @@
 .. code-block:: python
 
    >>> import yt.visualization.eps_writer as eps
-   >>> slc = SlicePlot(pf, 'z', 'density')
+   >>> slc = yt.SlicePlot(ds, 'z', 'density')
    >>> slc.set_width(25, 'kpc')
    >>> eps_fig = eps.single_plot(slc)
    >>> eps_fig.save_fig('zoom', format='eps')
@@ -706,9 +880,11 @@
 
 .. code-block:: python
 
+   >>> import yt
    >>> import yt.visualization.eps_writer as eps
-   >>> slc = SlicePlot(pf, 'z', ['density', 'temperature', 'Pressure',
-                       'VelocityMagnitude'])
+   >>>
+   >>> slc = yt.SlicePlot(ds, 'z', ['density', 'temperature', 'Pressure',
+                          'VelocityMagnitude'])
    >>> slc.set_width(25, 'kpc')
    >>> eps_fig = eps.multiplot_yt(2, 2, slc, bare_axes=True)
    >>> eps_fig.scale_line(0.2, '5 kpc')

This diff is so big that we needed to truncate the remainder.

https://bitbucket.org/yt_analysis/yt/commits/260316d16bbc/
Changeset:   260316d16bbc
Branch:      yt-3.0
User:        drudd
Date:        2014-07-12 17:47:28
Summary:     Add singular form of find_field_values_at_point which is a simple wrapper around dataset.point
Affected #:  1 file

diff -r 3b264449bb411235ba00ab34a01df79261899f6b -r 260316d16bbcdc5da0ab2b11e0c0766aa78c893e yt/data_objects/static_output.py
--- a/yt/data_objects/static_output.py
+++ b/yt/data_objects/static_output.py
@@ -541,6 +541,15 @@
               min_val, mx, my, mz)
         return min_val, self.arr([mx, my, mz], 'code_length', dtype="float64")
 
+    def find_field_values_at_point(self, fields, coord):
+        """
+        Returns the values [field1, field2,...] of the fields at the given
+        coordinates. Returns a list of field values in the same order as 
+        the input *fields*.
+        """
+        fields = ensure_list(fields)
+        out = [self.point(coords)[field] for field in fields]
+
     def find_field_values_at_points(self, fields, coords):
         """
         Returns the values [field1, field2,...] of the fields at the given


https://bitbucket.org/yt_analysis/yt/commits/eaff55614070/
Changeset:   eaff55614070
Branch:      yt-3.0
User:        drudd
Date:        2014-07-12 17:58:39
Summary:     Make Index.find_points a hidden (+optional) method for GridIndex
Affected #:  3 files

diff -r 260316d16bbcdc5da0ab2b11e0c0766aa78c893e -r eaff5561407033fb723f9e657084c9c649122582 yt/analysis_modules/particle_trajectories/particle_trajectories.py
--- a/yt/analysis_modules/particle_trajectories/particle_trajectories.py
+++ b/yt/analysis_modules/particle_trajectories/particle_trajectories.py
@@ -225,7 +225,8 @@
                     x = self["particle_position_x"][:,step].ndarray_view()
                     y = self["particle_position_y"][:,step].ndarray_view()
                     z = self["particle_position_z"][:,step].ndarray_view()
-                    particle_grids, particle_grid_inds = ds.index.find_points(x,y,z)
+                    # This will fail for non-grid index objects
+                    particle_grids, particle_grid_inds = ds.index._find_points(x,y,z)
                     for grid in particle_grids:
                         cube = grid.retrieve_ghost_zones(1, [fd])
                         CICSample_3(x,y,z,pfield,

diff -r 260316d16bbcdc5da0ab2b11e0c0766aa78c893e -r eaff5561407033fb723f9e657084c9c649122582 yt/frontends/stream/data_structures.py
--- a/yt/frontends/stream/data_structures.py
+++ b/yt/frontends/stream/data_structures.py
@@ -402,7 +402,7 @@
 def assign_particle_data(pf, pdata) :
 
     """
-    Assign particle data to the grids using find_points. This
+    Assign particle data to the grids using MatchPointsToGrids. This
     will overwrite any existing particle data, so be careful!
     """
     

diff -r 260316d16bbcdc5da0ab2b11e0c0766aa78c893e -r eaff5561407033fb723f9e657084c9c649122582 yt/geometry/grid_geometry_handler.py
--- a/yt/geometry/grid_geometry_handler.py
+++ b/yt/geometry/grid_geometry_handler.py
@@ -201,7 +201,8 @@
         for item in ("Mpc", "pc", "AU", "cm"):
             print "\tWidth: %0.3e %s" % (dx.in_units(item), item)
 
-    def find_points(self, x, y, z) :
+    def _find_
+    def _find_points(self, x, y, z) :
         """
         Returns the (objects, indices) of leaf grids containing a number of (x,y,z) points
         """
@@ -211,12 +212,12 @@
         if not len(x) == len(y) == len(z):
             raise AssertionError("Arrays of indices must be of the same size")
 
-        grid_tree = self.get_grid_tree()
+        grid_tree = self._get_grid_tree()
         pts = MatchPointsToGrids(grid_tree, len(x), x, y, z)
         ind = pts.find_points_in_tree()
         return self.grids[ind], ind
 
-    def get_grid_tree(self) :
+    def _get_grid_tree(self) :
 
         left_edge = self.pf.arr(np.zeros((self.num_grids, 3)),
                                'code_length')


https://bitbucket.org/yt_analysis/yt/commits/8f90a7d17d5b/
Changeset:   8f90a7d17d5b
Branch:      yt-3.0
User:        drudd
Date:        2014-07-12 17:59:07
Summary:     Typo
Affected #:  1 file

diff -r eaff5561407033fb723f9e657084c9c649122582 -r 8f90a7d17d5b51e2fdc7af6eeb34623762779952 yt/geometry/grid_geometry_handler.py
--- a/yt/geometry/grid_geometry_handler.py
+++ b/yt/geometry/grid_geometry_handler.py
@@ -201,7 +201,6 @@
         for item in ("Mpc", "pc", "AU", "cm"):
             print "\tWidth: %0.3e %s" % (dx.in_units(item), item)
 
-    def _find_
     def _find_points(self, x, y, z) :
         """
         Returns the (objects, indices) of leaf grids containing a number of (x,y,z) points


https://bitbucket.org/yt_analysis/yt/commits/2161c32b61ef/
Changeset:   2161c32b61ef
Branch:      yt-3.0
User:        drudd
Date:        2014-07-12 18:21:02
Summary:     Add shortcut for find_field_values_at_points defined on Index, change array order of coordinate parameter
Affected #:  1 file

diff -r 8f90a7d17d5b51e2fdc7af6eeb34623762779952 -r 2161c32b61efcb5ab6e78ca04ca8cdf1f98492f4 yt/data_objects/static_output.py
--- a/yt/data_objects/static_output.py
+++ b/yt/data_objects/static_output.py
@@ -556,23 +556,20 @@
         [(x1, y1, z2), (x2, y2, z2),...] points.  Returns a list of field 
         values in the same order as the input *fields*.
 
-        This is probably quite slow right now as it creates a new data object
-        for each point.  We'd probably need a YTPointCollection container to
-        optimize this.
+        This is quite slow right now as it creates a new data object for each
+        point.  If an optimized version exists on the Index object we'll use
+        that instead.
         """
+        if hasattr(self,"index") and \
+                hasattr(self.index,"_find_field_values_at_points"):
+            return self.index._find_field_values_at_points(fields,coords)
+
         fields = ensure_list(fields)
-        x, y, z = coords
-        x = ensure_list(x)
-        y = ensure_list(y)
-        z = ensure_list(z)
-        if not len(x) == len(y) == len(z):
-            raise AssertionError("Arrays of indices must be of the same size")
-       
-        out = [np.zeros(len(x), dtype=np.float64) for f in fields] 
-        for i in range(len(x)):
-            p = self.point([x[i],y[i],z[i]])
-            for j,field in enumerate(fields):
-                out[j][i] = p[field][0]
+        out = [np.zeros(len(coords), dtype=np.float64) for f in fields]
+        for i,coord in enumerate(coords):
+            data = self.point(coord)[fields]
+            for i in range(len(fields)):
+                out[j][i] = data[i]
         return out
 
     # Now all the object related stuff


https://bitbucket.org/yt_analysis/yt/commits/cb23bec6d341/
Changeset:   cb23bec6d341
Branch:      yt-3.0
User:        drudd
Date:        2014-07-12 18:26:13
Summary:     Changed from list return values
Affected #:  1 file

diff -r 2161c32b61efcb5ab6e78ca04ca8cdf1f98492f4 -r cb23bec6d341cea377082a98d917de4821db1452 yt/data_objects/static_output.py
--- a/yt/data_objects/static_output.py
+++ b/yt/data_objects/static_output.py
@@ -547,8 +547,7 @@
         coordinates. Returns a list of field values in the same order as 
         the input *fields*.
         """
-        fields = ensure_list(fields)
-        out = [self.point(coords)[field] for field in fields]
+        return self.point(coords)[fields]
 
     def find_field_values_at_points(self, fields, coords):
         """
@@ -565,11 +564,9 @@
             return self.index._find_field_values_at_points(fields,coords)
 
         fields = ensure_list(fields)
-        out = [np.zeros(len(coords), dtype=np.float64) for f in fields]
+        out = np.zeros((len(fields),len(coords)), dtype=np.float64)
         for i,coord in enumerate(coords):
-            data = self.point(coord)[fields]
-            for i in range(len(fields)):
-                out[j][i] = data[i]
+            out[:][i] = self.point(coord)[fields]
         return out
 
     # Now all the object related stuff


https://bitbucket.org/yt_analysis/yt/commits/7450fab2ea41/
Changeset:   7450fab2ea41
Branch:      yt-3.0
User:        drudd
Date:        2014-07-12 20:20:04
Summary:     Fixed failing tests due to moving of find_points to private method
Affected #:  2 files

diff -r cb23bec6d341cea377082a98d917de4821db1452 -r 7450fab2ea41bb71eaccaaa6921b559c92d481b3 yt/utilities/lib/tests/test_grid_tree.py
--- a/yt/utilities/lib/tests/test_grid_tree.py
+++ b/yt/utilities/lib/tests/test_grid_tree.py
@@ -50,7 +50,7 @@
 
 def test_grid_tree():
     """Main test suite for GridTree"""
-    grid_tree = test_pf.index.get_grid_tree()
+    grid_tree = test_pf.index._get_grid_tree()
     indices, levels, nchild, children = grid_tree.return_tree_info()
 
     grid_levels = [grid.Level for grid in test_pf.index.grids]

diff -r cb23bec6d341cea377082a98d917de4821db1452 -r 7450fab2ea41bb71eaccaaa6921b559c92d481b3 yt/utilities/particle_generator.py
--- a/yt/utilities/particle_generator.py
+++ b/yt/utilities/particle_generator.py
@@ -99,7 +99,7 @@
         Assigns grids to particles and sets up particle positions. *setup_fields* is
         a dict of fields other than the particle positions to set up. 
         """
-        particle_grids, particle_grid_inds = self.pf.index.find_points(x,y,z)
+        particle_grids, particle_grid_inds = self.pf.index._find_points(x,y,z)
         idxs = np.argsort(particle_grid_inds)
         self.particles[:,self.posx_index] = x[idxs]
         self.particles[:,self.posy_index] = y[idxs]


https://bitbucket.org/yt_analysis/yt/commits/253825ad96e4/
Changeset:   253825ad96e4
Branch:      yt-3.0
User:        drudd
Date:        2014-07-12 20:20:25
Summary:     Added back somewhat optimized version of find_field_values_at_points to grid_geometry_handler
Affected #:  1 file

diff -r 7450fab2ea41bb71eaccaaa6921b559c92d481b3 -r 253825ad96e484e07bf245bd481a178238fd656e yt/geometry/grid_geometry_handler.py
--- a/yt/geometry/grid_geometry_handler.py
+++ b/yt/geometry/grid_geometry_handler.py
@@ -201,6 +201,35 @@
         for item in ("Mpc", "pc", "AU", "cm"):
             print "\tWidth: %0.3e %s" % (dx.in_units(item), item)
 
+    def _find_field_values_at_points(self, fields, coords):
+        r"""Find the value of fields at a set of coordinates.
+
+        Returns the values [field1, field2,...] of the fields at the given
+        (x, y, z) points. Returns a numpy array of field values cross coords
+        """
+        coords = YTArray(ensure_numpy_array(coords),'code_length', registry=self.pf.unit_registry)
+        grids = self._find_points(coords[:,0], coords[:,1], coords[:,2])[0]
+        fields = ensure_list(fields)
+        mark = np.zeros(3, dtype=np.int)
+        out = []
+
+        # create point -> grid mapping
+        grid_index = {}
+        for coord_index, grid in enumerate(grids):
+            if not grid_index.has_key(grid):
+                grid_index[grid] = []
+            grid_index[grid].append(coord_index)
+
+        out = np.zeros((len(fields),len(coords)), dtype=np.float64)
+        for grid in grid_index:
+            cellwidth = (grid.RightEdge - grid.LeftEdge) / grid.ActiveDimensions
+            for field in fields:
+                for coord_index in grid_index[grid]:
+                    mark = ((coords[coord_index,:] - grid.LeftEdge) / cellwidth).astype('int')
+                    out[:,coord_index] = grid[field][mark[0],mark[1],mark[2]]
+        return out
+
+
     def _find_points(self, x, y, z) :
         """
         Returns the (objects, indices) of leaf grids containing a number of (x,y,z) points


https://bitbucket.org/yt_analysis/yt/commits/561d2e56371c/
Changeset:   561d2e56371c
Branch:      yt-3.0
User:        MatthewTurk
Date:        2014-07-14 20:02:36
Summary:     Merged in drudd/yt/yt-3.0 (pull request #1015)

Point search refactoring (take 2)
Affected #:  17 files

diff -r ef7c988b774a159c854d7a2fefc000e767007eda -r 561d2e56371ca0053ee093df7968f59ae978481a yt/analysis_modules/particle_trajectories/particle_trajectories.py
--- a/yt/analysis_modules/particle_trajectories/particle_trajectories.py
+++ b/yt/analysis_modules/particle_trajectories/particle_trajectories.py
@@ -225,7 +225,8 @@
                     x = self["particle_position_x"][:,step].ndarray_view()
                     y = self["particle_position_y"][:,step].ndarray_view()
                     z = self["particle_position_z"][:,step].ndarray_view()
-                    particle_grids, particle_grid_inds = ds.index.find_points(x,y,z)
+                    # This will fail for non-grid index objects
+                    particle_grids, particle_grid_inds = ds.index._find_points(x,y,z)
                     for grid in particle_grids:
                         cube = grid.retrieve_ghost_zones(1, [fd])
                         CICSample_3(x,y,z,pfield,

diff -r ef7c988b774a159c854d7a2fefc000e767007eda -r 561d2e56371ca0053ee093df7968f59ae978481a yt/data_objects/data_containers.py
--- a/yt/data_objects/data_containers.py
+++ b/yt/data_objects/data_containers.py
@@ -725,6 +725,12 @@
             self.index._identify_base_chunk(self)
         return self._current_chunk.fwidth
 
+class YTSelectionContainer0D(YTSelectionContainer):
+    _spatial = False
+    def __init__(self, pf, field_parameters):
+        super(YTSelectionContainer0D, self).__init__(
+            pf, field_parameters)
+
 class YTSelectionContainer1D(YTSelectionContainer):
     _spatial = False
     def __init__(self, pf, field_parameters):

diff -r ef7c988b774a159c854d7a2fefc000e767007eda -r 561d2e56371ca0053ee093df7968f59ae978481a yt/data_objects/selection_data_containers.py
--- a/yt/data_objects/selection_data_containers.py
+++ b/yt/data_objects/selection_data_containers.py
@@ -21,7 +21,8 @@
 from yt.utilities.lib.alt_ray_tracers import cylindrical_ray_trace
 from yt.utilities.orientation import Orientation
 from .data_containers import \
-    YTSelectionContainer1D, YTSelectionContainer2D, YTSelectionContainer3D
+    YTSelectionContainer0D, YTSelectionContainer1D, \
+    YTSelectionContainer2D, YTSelectionContainer3D
 from yt.data_objects.derived_quantities import \
     DerivedQuantityCollection
 from yt.utilities.exceptions import \
@@ -34,6 +35,31 @@
 from yt.utilities.math_utils import get_rotation_matrix
 from yt.units.yt_array import YTQuantity
 
+
+class YTPointBase(YTSelectionContainer0D):
+    """
+    A 0-dimensional object defined by a single point
+
+    Parameters
+    ----------
+    p: array_like
+        A points defined within the domain.  If the domain is
+        periodic its position will be corrected to lie inside
+        the range [DLE,DRE) to ensure one and only one cell may
+        match that point
+
+    Examples
+    --------
+    >>> pf = load("DD0010/moving7_0010")
+    >>> c = [0.5,0.5,0.5]
+    >>> point = pf.point(c)
+    """
+    _type_name = "point"
+    _con_args = ('p',)
+    def __init__(self, p, pf = None, field_parameters = None):
+        super(YTPointBase, self).__init__(pf, field_parameters)
+        self.p = p
+
 class YTOrthoRayBase(YTSelectionContainer1D):
     """
     This is an orthogonal ray cast through the entire domain, at a specific
@@ -529,7 +555,7 @@
 
 class YTSphereBase(YTSelectionContainer3D):
     """
-    A sphere f points defined by a *center* and a *radius*.
+    A sphere of points defined by a *center* and a *radius*.
 
     Parameters
     ----------

diff -r ef7c988b774a159c854d7a2fefc000e767007eda -r 561d2e56371ca0053ee093df7968f59ae978481a yt/data_objects/static_output.py
--- a/yt/data_objects/static_output.py
+++ b/yt/data_objects/static_output.py
@@ -541,6 +541,34 @@
               min_val, mx, my, mz)
         return min_val, self.arr([mx, my, mz], 'code_length', dtype="float64")
 
+    def find_field_values_at_point(self, fields, coord):
+        """
+        Returns the values [field1, field2,...] of the fields at the given
+        coordinates. Returns a list of field values in the same order as 
+        the input *fields*.
+        """
+        return self.point(coords)[fields]
+
+    def find_field_values_at_points(self, fields, coords):
+        """
+        Returns the values [field1, field2,...] of the fields at the given
+        [(x1, y1, z2), (x2, y2, z2),...] points.  Returns a list of field 
+        values in the same order as the input *fields*.
+
+        This is quite slow right now as it creates a new data object for each
+        point.  If an optimized version exists on the Index object we'll use
+        that instead.
+        """
+        if hasattr(self,"index") and \
+                hasattr(self.index,"_find_field_values_at_points"):
+            return self.index._find_field_values_at_points(fields,coords)
+
+        fields = ensure_list(fields)
+        out = np.zeros((len(fields),len(coords)), dtype=np.float64)
+        for i,coord in enumerate(coords):
+            out[:][i] = self.point(coord)[fields]
+        return out
+
     # Now all the object related stuff
     def all_data(self, find_max=False):
         if find_max: c = self.find_max("density")[1]

diff -r ef7c988b774a159c854d7a2fefc000e767007eda -r 561d2e56371ca0053ee093df7968f59ae978481a yt/data_objects/tests/test_points.py
--- /dev/null
+++ b/yt/data_objects/tests/test_points.py
@@ -0,0 +1,10 @@
+from yt.testing import *
+import numpy as np
+
+def setup():
+    from yt.config import ytcfg
+    ytcfg["yt","__withintesting"] = "True"
+
+def test_domain_point():
+    pf = fake_random_pf(16, fields = ("density"))
+    p = pf.point(pf.domain_center)

diff -r ef7c988b774a159c854d7a2fefc000e767007eda -r 561d2e56371ca0053ee093df7968f59ae978481a yt/frontends/artio/_artio_caller.pyx
--- a/yt/frontends/artio/_artio_caller.pyx
+++ b/yt/frontends/artio/_artio_caller.pyx
@@ -300,6 +300,10 @@
         cdef int num_fields = len(fields)
         cdef np.float64_t pos[3]
 
+        # since RuntimeErrors are not fatal, ensure no artio_particles* functions
+        # called if fileset lacks particles
+        if not self.has_particles: return
+
         data = {}
         accessed_species = np.zeros( self.num_species, dtype="int")
         selected_mass = [ None for i in range(self.num_species)]
@@ -587,15 +591,25 @@
         self.handle = artio_handle.handle
         self.oct_count = None
         self.root_mesh_data = NULL
-        self.pcount = <np.int64_t **> malloc(sizeof(np.int64_t*)
-            * artio_handle.num_species)
-        self.nvars[0] = artio_handle.num_species
-        self.nvars[1] = artio_handle.num_grid_variables
-        for i in range(artio_handle.num_species):
-            self.pcount[i] = <np.int64_t*> malloc(sizeof(np.int64_t)
-                * (self.sfc_end - self.sfc_start + 1))
-            for sfc in range(self.sfc_end - self.sfc_start + 1):
-                self.pcount[i][sfc] = 0
+        self.pcount = NULL
+
+        if artio_handle.has_particles:
+            self.pcount = <np.int64_t **> malloc(sizeof(np.int64_t*)
+                * artio_handle.num_species)
+            self.nvars[0] = artio_handle.num_species
+            for i in range(artio_handle.num_species):
+                self.pcount[i] = <np.int64_t*> malloc(sizeof(np.int64_t)
+                    * (self.sfc_end - self.sfc_start + 1))
+                for sfc in range(self.sfc_end - self.sfc_start + 1):
+                    self.pcount[i][sfc] = 0
+        else:
+            self.nvars[0] = 0
+
+        if artio_handle.has_grid:
+            self.nvars[1] = artio_handle.num_grid_variables
+        else:
+            self.nvars[1] = 0
+
         for i in range(3):
             self.dims[i] = domain_dimensions[i]
             self.DLE[i] = domain_left_edge[i]
@@ -604,12 +618,15 @@
 
     def __dealloc__(self):
         cdef int i
-        for i in range(self.nvars[0]):
-            free(self.pcount[i])
-        for i in range(self.nvars[1]):
-            free(self.root_mesh_data[i])
-        free(self.pcount)
-        free(self.root_mesh_data)
+        if self.artio_handle.has_particles:
+            for i in range(self.nvars[0]):
+                free(self.pcount[i])
+            free(self.pcount)
+        if self.artio_handle.has_grid:
+            if self.root_mesh_data != NULL:
+                for i in range(self.nvars[1]):
+                    free(self.root_mesh_data[i])
+                free(self.root_mesh_data)
 
     @cython.boundscheck(False)
     @cython.wraparound(False)
@@ -623,8 +640,7 @@
         cdef int *num_octs_per_level = <int *>malloc(
             (max_level + 1)*sizeof(int))
         cdef int num_species = self.artio_handle.num_species
-        cdef int *num_particles_per_species =  <int *>malloc(
-            sizeof(int)*num_species) 
+        cdef int *num_particles_per_species 
         cdef ARTIOOctreeContainer octree
         ngv = self.nvars[1]
         cdef float *grid_variables = <float *>malloc(
@@ -648,10 +664,10 @@
             status = artio_grid_read_root_cell_begin( self.handle,
                 sfc, dpos, grid_variables, &num_oct_levels,
                 num_octs_per_level)
+            check_artio_status(status)
             for i in range(ngv):
                 self.root_mesh_data[i][sfc - self.sfc_start] = \
                     grid_variables[i]
-            check_artio_status(status)
             if num_oct_levels > 0:
                 oc = 0
                 for level in range(num_oct_levels):
@@ -660,33 +676,39 @@
                 oct_count[sfc - self.sfc_start] = oc
                 octree.initialize_local_mesh(oc, num_oct_levels,
                     num_octs_per_level, sfc)
-            status = artio_grid_read_root_cell_end( self.handle )
+            status = artio_grid_read_root_cell_end(self.handle)
             check_artio_status(status)
-        status = artio_grid_clear_sfc_cache( self.handle)
+        status = artio_grid_clear_sfc_cache(self.handle)
         check_artio_status(status)
-        # Now particles
-        status = artio_particle_cache_sfc_range(self.handle, self.sfc_start,
+
+        if self.artio_handle.has_particles:
+            num_particles_per_species =  <int *>malloc(
+                    sizeof(int)*num_species)
+
+            # Now particles
+            status = artio_particle_cache_sfc_range(self.handle, self.sfc_start,
                                             self.sfc_end)
-        check_artio_status(status) 
-        for sfc in range(self.sfc_start, self.sfc_end + 1):
-            # Now particles
-            status = artio_particle_read_root_cell_begin( self.handle,
-                    sfc, num_particles_per_species )
+            check_artio_status(status) 
+            for sfc in range(self.sfc_start, self.sfc_end + 1):
+                # Now particles
+                status = artio_particle_read_root_cell_begin(self.handle,
+                        sfc, num_particles_per_species)
+                check_artio_status(status)
+
+                for i in range(num_species):
+                    self.pcount[i][sfc - self.sfc_start] = \
+                        num_particles_per_species[i]
+
+                status = artio_particle_read_root_cell_end(self.handle)
+                check_artio_status(status)
+
+            status = artio_particle_clear_sfc_cache(self.handle)
             check_artio_status(status)
 
-            for i in range(num_species):
-                self.pcount[i][sfc - self.sfc_start] = \
-                    num_particles_per_species[i]
-
-            status = artio_particle_read_root_cell_end( self.handle)
-            check_artio_status(status)
-
-        status = artio_particle_clear_sfc_cache( self.handle)
-        check_artio_status(status)
+            free(num_particles_per_species)
 
         free(grid_variables)
         free(num_octs_per_level)
-        free(num_particles_per_species)
         self.oct_count = oct_count
         self.doct_count = <np.int64_t *> oct_count.data
         self.root_mesh_handler = ARTIORootMeshContainer(self)
@@ -861,7 +883,7 @@
         cdef np.int64_t sfc_start, sfc_end
         sfc_start = self.domains[0].con_id
         sfc_end = self.domains[self.num_domains - 1].con_id
-        status = artio_grid_cache_sfc_range(handle, sfc_start, sfc_end )
+        status = artio_grid_cache_sfc_range(handle, sfc_start, sfc_end)
         check_artio_status(status) 
         cdef np.int64_t offset = 0 
         for si in range(self.num_domains):
@@ -884,7 +906,7 @@
                 status = artio_grid_read_level_end(handle)
                 check_artio_status(status)
                 lp += num_octs_per_level[level]
-            status = artio_grid_read_root_cell_end( handle )
+            status = artio_grid_read_root_cell_end(handle)
             check_artio_status(status)
             # Now we have all our sources.
             for j in range(nf):
@@ -943,6 +965,10 @@
     cdef np.ndarray[np.int64_t, ndim=1] npi64arr
     cdef np.ndarray[np.float64_t, ndim=1] npf64arr
 
+    if not artio_handle.has_particles:
+        raise RuntimeError("Attempted to read non-existent particles in ARTIO")
+        return
+
     # Now we set up our field pointers
     params = artio_handle.parameters
 

diff -r ef7c988b774a159c854d7a2fefc000e767007eda -r 561d2e56371ca0053ee093df7968f59ae978481a yt/frontends/artio/data_structures.py
--- a/yt/frontends/artio/data_structures.py
+++ b/yt/frontends/artio/data_structures.py
@@ -182,29 +182,6 @@
     def convert(self, unit):
         return self.parameter_file.conversion_factors[unit]
 
-    def find_max(self, field, finest_levels=3):
-        """
-        Returns (value, center) of location of maximum for a given field.
-        """
-        if (field, finest_levels) in self._max_locations:
-            return self._max_locations[(field, finest_levels)]
-        mv, pos = self.find_max_cell_location(field, finest_levels)
-        self._max_locations[(field, finest_levels)] = (mv, pos)
-        return mv, pos
-
-    def find_max_cell_location(self, field, finest_levels=3):
-        source = self.all_data()
-        if finest_levels is not False:
-            source.min_level = self.max_level - finest_levels
-        mylog.debug("Searching for maximum value of %s", field)
-        max_val, maxi, mx, my, mz = \
-            source.quantities["MaxLocation"](field)
-        mylog.info("Max Value is %0.5e at %0.16f %0.16f %0.16f",
-                   max_val, mx, my, mz)
-        self.pf.parameters["Max%sValue" % (field)] = max_val
-        self.pf.parameters["Max%sPos" % (field)] = "%s" % ((mx, my, mz),)
-        return max_val, np.array((mx, my, mz), dtype='float64')
-
     def _detect_output_fields(self):
         self.fluid_field_list = self._detect_fluid_fields()
         self.particle_field_list = self._detect_particle_fields()

diff -r ef7c988b774a159c854d7a2fefc000e767007eda -r 561d2e56371ca0053ee093df7968f59ae978481a yt/frontends/boxlib/data_structures.py
--- a/yt/frontends/boxlib/data_structures.py
+++ b/yt/frontends/boxlib/data_structures.py
@@ -343,7 +343,6 @@
         self._parallel_locking = False
         self._data_file = None
         self._data_mode = None
-        self._max_locations = {}
 
     def _detect_output_fields(self):
         # This is all done in _parse_header_file

diff -r ef7c988b774a159c854d7a2fefc000e767007eda -r 561d2e56371ca0053ee093df7968f59ae978481a yt/frontends/stream/data_structures.py
--- a/yt/frontends/stream/data_structures.py
+++ b/yt/frontends/stream/data_structures.py
@@ -402,7 +402,7 @@
 def assign_particle_data(pf, pdata) :
 
     """
-    Assign particle data to the grids using find_points. This
+    Assign particle data to the grids using MatchPointsToGrids. This
     will overwrite any existing particle data, so be careful!
     """
     

diff -r ef7c988b774a159c854d7a2fefc000e767007eda -r 561d2e56371ca0053ee093df7968f59ae978481a yt/geometry/geometry_handler.py
--- a/yt/geometry/geometry_handler.py
+++ b/yt/geometry/geometry_handler.py
@@ -72,7 +72,6 @@
         self._parallel_locking = False
         self._data_file = None
         self._data_mode = None
-        self._max_locations = {}
         self.num_grids = None
 
     def _initialize_data_storage(self):

diff -r ef7c988b774a159c854d7a2fefc000e767007eda -r 561d2e56371ca0053ee093df7968f59ae978481a yt/geometry/grid_geometry_handler.py
--- a/yt/geometry/grid_geometry_handler.py
+++ b/yt/geometry/grid_geometry_handler.py
@@ -201,32 +201,36 @@
         for item in ("Mpc", "pc", "AU", "cm"):
             print "\tWidth: %0.3e %s" % (dx.in_units(item), item)
 
-    def find_max(self, field, finest_levels = 3):
+    def _find_field_values_at_points(self, fields, coords):
+        r"""Find the value of fields at a set of coordinates.
+
+        Returns the values [field1, field2,...] of the fields at the given
+        (x, y, z) points. Returns a numpy array of field values cross coords
         """
-        Returns (value, center) of location of maximum for a given field.
-        """
-        if (field, finest_levels) in self._max_locations:
-            return self._max_locations[(field, finest_levels)]
-        mv, pos = self.find_max_cell_location(field, finest_levels)
-        self._max_locations[(field, finest_levels)] = (mv, pos)
-        return mv, pos
+        coords = YTArray(ensure_numpy_array(coords),'code_length', registry=self.pf.unit_registry)
+        grids = self._find_points(coords[:,0], coords[:,1], coords[:,2])[0]
+        fields = ensure_list(fields)
+        mark = np.zeros(3, dtype=np.int)
+        out = []
 
-    def find_max_cell_location(self, field, finest_levels = 3):
-        if finest_levels is not False:
-            gi = (self.grid_levels >= self.max_level - finest_levels).ravel()
-            source = self.data_collection([0.0]*3, self.grids[gi])
-        else:
-            source = self.all_data()
-        mylog.debug("Searching for maximum value of %s", field)
-        max_val, maxi, mx, my, mz = \
-            source.quantities["MaxLocation"](field)
-        mylog.info("Max Value is %0.5e at %0.16f %0.16f %0.16f", 
-              max_val, mx, my, mz)
-        self.parameters["Max%sValue" % (field)] = max_val
-        self.parameters["Max%sPos" % (field)] = "%s" % ((mx,my,mz),)
-        return max_val, np.array((mx,my,mz), dtype='float64')
+        # create point -> grid mapping
+        grid_index = {}
+        for coord_index, grid in enumerate(grids):
+            if not grid_index.has_key(grid):
+                grid_index[grid] = []
+            grid_index[grid].append(coord_index)
 
-    def find_points(self, x, y, z) :
+        out = np.zeros((len(fields),len(coords)), dtype=np.float64)
+        for grid in grid_index:
+            cellwidth = (grid.RightEdge - grid.LeftEdge) / grid.ActiveDimensions
+            for field in fields:
+                for coord_index in grid_index[grid]:
+                    mark = ((coords[coord_index,:] - grid.LeftEdge) / cellwidth).astype('int')
+                    out[:,coord_index] = grid[field][mark[0],mark[1],mark[2]]
+        return out
+
+
+    def _find_points(self, x, y, z) :
         """
         Returns the (objects, indices) of leaf grids containing a number of (x,y,z) points
         """
@@ -236,46 +240,12 @@
         if not len(x) == len(y) == len(z):
             raise AssertionError("Arrays of indices must be of the same size")
 
-        grid_tree = self.get_grid_tree()
+        grid_tree = self._get_grid_tree()
         pts = MatchPointsToGrids(grid_tree, len(x), x, y, z)
         ind = pts.find_points_in_tree()
         return self.grids[ind], ind
 
-    def find_field_value_at_point(self, fields, coord):
-        r"""Find the value of fields at a coordinate.
-
-        Returns the values [field1, field2,...] of the fields at the given
-        (x, y, z) points. Returns a list of field values in the same order as
-        the input *fields*.
-
-        Parameters
-        ----------
-        fields : string or list of strings
-            The field(s) that will be returned.
-
-        coord : list or array of coordinates
-            The location for which field values will be returned.
-
-        Examples
-        --------
-        >>> pf.h.find_field_value_at_point(['Density', 'Temperature'],
-            [0.4, 0.3, 0.8])
-        [2.1489e-24, 1.23843e4]
-        """
-        this = self.find_points(*coord)[0][-1]
-        cellwidth = (this.RightEdge - this.LeftEdge) / this.ActiveDimensions
-        mark = np.zeros(3).astype('int')
-        # Find the index for the cell containing this point.
-        for dim in xrange(len(coord)):
-            mark[dim] = int((coord[dim] - this.LeftEdge[dim]) / cellwidth[dim])
-        out = []
-        fields = ensure_list(fields)
-        # Pull out the values and add it to the out list.
-        for field in fields:
-            out.append(this[field][mark[0], mark[1], mark[2]])
-        return out
-
-    def get_grid_tree(self) :
+    def _get_grid_tree(self) :
 
         left_edge = self.pf.arr(np.zeros((self.num_grids, 3)),
                                'code_length')

diff -r ef7c988b774a159c854d7a2fefc000e767007eda -r 561d2e56371ca0053ee093df7968f59ae978481a yt/geometry/oct_geometry_handler.py
--- a/yt/geometry/oct_geometry_handler.py
+++ b/yt/geometry/oct_geometry_handler.py
@@ -49,27 +49,3 @@
 
     def convert(self, unit):
         return self.parameter_file.conversion_factors[unit]
-
-    def find_max(self, field, finest_levels = 3):
-        """
-        Returns (value, center) of location of maximum for a given field.
-        """
-        if (field, finest_levels) in self._max_locations:
-            return self._max_locations[(field, finest_levels)]
-        mv, pos = self.find_max_cell_location(field, finest_levels)
-        self._max_locations[(field, finest_levels)] = (mv, pos)
-        return mv, pos
-
-    def find_max_cell_location(self, field, finest_levels = 3):
-        source = self.all_data()
-        if finest_levels is not False:
-            source.min_level = self.max_level - finest_levels
-        mylog.debug("Searching for maximum value of %s", field)
-        max_val, maxi, mx, my, mz = \
-            source.quantities["MaxLocation"](field)
-        mylog.info("Max Value is %0.5e at %0.16f %0.16f %0.16f", 
-              max_val, mx, my, mz)
-        self.pf.parameters["Max%sValue" % (field,)] = max_val
-        self.pf.parameters["Max%sPos" % (field,)] = "%s" % ((mx,my,mz),)
-        return max_val, np.array((mx,my,mz), dtype='float64')
-

diff -r ef7c988b774a159c854d7a2fefc000e767007eda -r 561d2e56371ca0053ee093df7968f59ae978481a yt/geometry/selection_routines.pyx
--- a/yt/geometry/selection_routines.pyx
+++ b/yt/geometry/selection_routines.pyx
@@ -534,6 +534,63 @@
                 self.domain_width[1],
                 self.domain_width[2])
 
+
+cdef class PointSelector(SelectorObject):
+    cdef np.float64_t p[3]
+
+    def __init__(self, dobj):
+        for i in range(3):
+            self.p[i] = _ensure_code(dobj.p[i])
+
+            # ensure the point lies in the domain
+            if self.periodicity[i]:
+                self.p[i] = np.fmod(self.p[i], self.domain_width[i])
+                if self.p[i] < 0.0:
+                    self.p[i] += self.domain_width[i]
+
+    @cython.boundscheck(False)
+    @cython.wraparound(False)
+    @cython.cdivision(True)
+    cdef int select_cell(self, np.float64_t pos[3], np.float64_t dds[3]) nogil:
+        if (pos[0] - 0.5*dds[0] <= self.p[0] < pos[0]+0.5*dds[0] and
+            pos[1] - 0.5*dds[1] <= self.p[1] < pos[1]+0.5*dds[1] and
+            pos[2] - 0.5*dds[2] <= self.p[2] < pos[2]+0.5*dds[2]):
+            return 1
+        else:
+            return 0
+
+    @cython.boundscheck(False)
+    @cython.wraparound(False)
+    @cython.cdivision(True)
+    cdef int select_sphere(self, np.float64_t pos[3], np.float64_t radius) nogil:
+        cdef int i
+        cdef np.float64_t dist, dist2 = 0
+        for i in range(3):
+            dist = self.difference(pos[i], self.p[i], i)
+            dist2 += dist*dist
+        if dist2 <= radius*radius: return 1
+        return 0
+
+    @cython.boundscheck(False)
+    @cython.wraparound(False)
+    @cython.cdivision(True)
+    cdef int select_bbox(self, np.float64_t left_edge[3],
+                               np.float64_t right_edge[3]) nogil:
+        cdef int i
+        # point definitely can only be in one cell
+        if (left_edge[0] <= self.p[0] < right_edge[0] and
+            left_edge[1] <= self.p[1] < right_edge[1] and
+            left_edge[2] <= self.p[2] < right_edge[2]):
+            return 1
+        else:
+            return 0
+
+    def _hash_vals(self):
+        return (self.p[0], self.p[1], self.p[2])
+
+point_selector = PointSelector
+
+
 cdef class SphereSelector(SelectorObject):
     cdef np.float64_t radius
     cdef np.float64_t radius2

diff -r ef7c988b774a159c854d7a2fefc000e767007eda -r 561d2e56371ca0053ee093df7968f59ae978481a yt/utilities/lib/tests/test_grid_tree.py
--- a/yt/utilities/lib/tests/test_grid_tree.py
+++ b/yt/utilities/lib/tests/test_grid_tree.py
@@ -50,7 +50,7 @@
 
 def test_grid_tree():
     """Main test suite for GridTree"""
-    grid_tree = test_pf.index.get_grid_tree()
+    grid_tree = test_pf.index._get_grid_tree()
     indices, levels, nchild, children = grid_tree.return_tree_info()
 
     grid_levels = [grid.Level for grid in test_pf.index.grids]
@@ -67,53 +67,3 @@
             grid_children = np.array([child.id - child._id_offset
                                       for child in grid.Children])
             yield assert_equal, grid_children, children[i]
-
-
-def test_find_points():
-    """Main test suite for MatchPoints"""
-    num_points = 100
-    randx = np.random.uniform(low=test_pf.domain_left_edge[0],
-                              high=test_pf.domain_right_edge[0],
-                              size=num_points)
-    randy = np.random.uniform(low=test_pf.domain_left_edge[1],
-                              high=test_pf.domain_right_edge[1],
-                              size=num_points)
-    randz = np.random.uniform(low=test_pf.domain_left_edge[2],
-                              high=test_pf.domain_right_edge[2],
-                              size=num_points)
-
-    point_grids, point_grid_inds = test_pf.index.find_points(randx, randy, randz)
-
-    grid_inds = np.zeros((num_points), dtype='int64')
-
-    for ind, ixx, iyy, izz in zip(range(num_points), randx, randy, randz):
-
-        pos = np.array([ixx, iyy, izz])
-        pt_level = -1
-
-        for grid in test_pf.index.grids:
-
-            if np.all(pos >= grid.LeftEdge) and \
-               np.all(pos <= grid.RightEdge) and \
-               grid.Level > pt_level:
-                pt_level = grid.Level
-                grid_inds[ind] = grid.id - grid._id_offset
-
-    yield assert_equal, point_grid_inds, grid_inds
-
-    # Test wheter find_points works for lists
-    point_grids, point_grid_inds = test_pf.index.find_points(randx.tolist(),
-                                                         randy.tolist(),
-                                                         randz.tolist())
-    yield assert_equal, point_grid_inds, grid_inds
-
-    # Test if find_points works for scalar
-    ind = random.randint(0, num_points - 1)
-    point_grids, point_grid_inds = test_pf.index.find_points(randx[ind],
-                                                         randy[ind],
-                                                         randz[ind])
-    yield assert_equal, point_grid_inds, grid_inds[ind]
-
-    # Test if find_points fails properly for non equal indices' array sizes
-    yield assert_raises, AssertionError, test_pf.index.find_points, \
-        [0], 1.0, [2, 3]

diff -r ef7c988b774a159c854d7a2fefc000e767007eda -r 561d2e56371ca0053ee093df7968f59ae978481a yt/utilities/particle_generator.py
--- a/yt/utilities/particle_generator.py
+++ b/yt/utilities/particle_generator.py
@@ -99,7 +99,7 @@
         Assigns grids to particles and sets up particle positions. *setup_fields* is
         a dict of fields other than the particle positions to set up. 
         """
-        particle_grids, particle_grid_inds = self.pf.index.find_points(x,y,z)
+        particle_grids, particle_grid_inds = self.pf.index._find_points(x,y,z)
         idxs = np.argsort(particle_grid_inds)
         self.particles[:,self.posx_index] = x[idxs]
         self.particles[:,self.posy_index] = y[idxs]

diff -r ef7c988b774a159c854d7a2fefc000e767007eda -r 561d2e56371ca0053ee093df7968f59ae978481a yt/utilities/tests/test_selectors.py
--- a/yt/utilities/tests/test_selectors.py
+++ b/yt/utilities/tests/test_selectors.py
@@ -1,7 +1,7 @@
 import numpy as np
 from yt.testing import \
     fake_random_pf, assert_equal, assert_array_less, \
-    YTArray
+    YTArray, fake_amr_pf
 from yt.utilities.math_utils import periodic_dist
 
 
@@ -9,6 +9,26 @@
     from yt.config import ytcfg
     ytcfg["yt", "__withintesting"] = "True"
 
+def test_point_selector():
+    # generate fake amr data
+    pf = fake_random_pf(64, nprocs=51)
+    assert(all(pf.periodicity))
+
+    dd = pf.h.all_data()
+    positions = np.array([dd[ax] for ax in 'xyz'])
+    delta = 0.5*np.array([dd['d'+ax] for ax in 'xyz'])
+    # ensure cell centers and corners always return one and 
+    # only one point object
+    for p in positions:
+        data = pf.point(p)
+        assert_equal(data["ones"].shape[0], 1) 
+    for p in positions - delta:
+        data = pf.point(p)
+        assert_equal(data["ones"].shape[0], 1)
+    for p in positions + delta:
+        data = pf.point(p)
+        assert_equal(data["ones"].shape[0], 1)
+ 
 def test_sphere_selector():
     # generate fake data with a number of non-cubical grids
     pf = fake_random_pf(64, nprocs=51)

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