diff --git a/storm_analysis/sa_library/analysis_io.py b/storm_analysis/sa_library/analysis_io.py
index 2de99593..b8c36526 100644
--- a/storm_analysis/sa_library/analysis_io.py
+++ b/storm_analysis/sa_library/analysis_io.py
@@ -54,7 +54,7 @@ def loadCMOSCalibration(filename, verbose = False):
 
         rqe = numpy.ones_like(data[0])
         data = list(data) + [rqe]
-        return map(numpy.transpose, data)
+        return list(map(numpy.transpose, data))
 
     # Check for v1 format.
     elif (len(data) == 4):
diff --git a/storm_analysis/sa_library/dao_fit.c b/storm_analysis/sa_library/dao_fit.c
index 0be46b66..6c671d2b 100644
--- a/storm_analysis/sa_library/dao_fit.c
+++ b/storm_analysis/sa_library/dao_fit.c
@@ -983,7 +983,7 @@ void daoCalcPeakShape(fitData *fit_data)
   dao_peak = (daoPeak *)peak->peak_model;
 
   xo = -(peak->params[XCENTER] - (double)peak->xi) - fit_data->xoff;
-  yo = -(peak->params[YCENTER] - (double)peak->yi) - fit_data->xoff;
+  yo = -(peak->params[YCENTER] - (double)peak->yi) - fit_data->yoff;
 
   /* Calculate peak shape. */
   for(i=0;i<fit_data->fit_size_x;i++){
diff --git a/storm_analysis/sa_library/datareader.py b/storm_analysis/sa_library/datareader.py
index a2132e99..e55d83b7 100644
--- a/storm_analysis/sa_library/datareader.py
+++ b/storm_analysis/sa_library/datareader.py
@@ -167,8 +167,10 @@ def __init__(self, filename, verbose = False):
         self.inf_filename = dirname + os.path.splitext(os.path.basename(filename))[0] + ".inf"
 
         # defaults
+        self.bigendian = None
         self.image_height = None
         self.image_width = None
+        self.number_frames = None
 
         # extract the movie information from the associated inf file
         size_re = re.compile(r'frame dimensions = ([\d]+) x ([\d]+)')
@@ -217,17 +219,29 @@ def __init__(self, filename, verbose = False):
 
         # set defaults, probably correct, but warn the user 
         # that they couldn't be determined from the inf file.
-        if not self.image_height:
+        if self.image_height is None:
             print("Could not determine image size, assuming 256x256.")
             self.image_height = 256
             self.image_width = 256
 
+        # Not every .inf file records the endianness. Everything this package
+        # writes is little endian, so that is the better guess.
+        if self.bigendian is None:
+            print("Could not determine endianness of", self.inf_filename, "assuming little endian.")
+            self.bigendian = 0
+
         # open the dax file
-        if os.path.exists(filename):
-            self.fileptr = open(filename, "rb")
-        else:
-            if self.verbose:
-                print("dax data not found", filename)
+        if not os.path.exists(filename):
+            raise IOError("dax data not found '" + filename + "'")
+        self.fileptr = open(filename, "rb")
+
+        # A dax file is a bare block of 16 bit pixels with no header, so if the
+        # .inf file did not say how many frames there are the file size does.
+        if self.number_frames is None:
+            frame_size = 2 * self.image_height * self.image_width
+            self.number_frames = os.path.getsize(filename)//frame_size
+            print("Could not determine the number of frames in", self.inf_filename + ",",
+                  "the file size gives", self.number_frames)
 
     def loadAFrame(self, frame_number):
         """
diff --git a/storm_analysis/sa_library/fitting.py b/storm_analysis/sa_library/fitting.py
index ae1d6e82..c30cf27a 100644
--- a/storm_analysis/sa_library/fitting.py
+++ b/storm_analysis/sa_library/fitting.py
@@ -210,11 +210,11 @@ def peakMask(shape, parameters, margin):
         if parameters.hasAttr("x_start"):
             peak_mask[:,0:parameters.getAttr("x_start")+margin] = 0.0
         if parameters.hasAttr("x_stop"):
-            peak_mask[:,parameters.getAttr("x_stop")+margin:-1] = 0.0
+            peak_mask[:,parameters.getAttr("x_stop")+margin:] = 0.0
         if parameters.hasAttr("y_start"):
             peak_mask[0:parameters.getAttr("y_start")+margin,:] = 0.0
         if parameters.hasAttr("y_stop"):
-            peak_mask[parameters.getAttr("y_stop")+margin:-1,:] = 0.0
+            peak_mask[parameters.getAttr("y_stop")+margin:,:] = 0.0
         
     return peak_mask
 
diff --git a/storm_analysis/sa_library/multi_fit.c b/storm_analysis/sa_library/multi_fit.c
index 2aeb697c..54c94be5 100644
--- a/storm_analysis/sa_library/multi_fit.c
+++ b/storm_analysis/sa_library/multi_fit.c
@@ -777,12 +777,6 @@ void mFitGetPeakPropertyDouble(fitData *fit_data, double *values, char *what)
       values[i] = mFitPeakFgSum(fit_data, &fit_data->fit[i]);
     }
   }
-  else if (!strcmp(what, "fg_sum_sc")){
-    for(i=0;i<fit_data->nfit;i++){
-      values[i] = mFitPeakFgSum(fit_data, &fit_data->fit[i]);
-      //      values[i] = mFitPeakFgSumSensitivityCorrected(fit_data, &fit_data->fit[i]);
-    }
-  }
   else if (!strcmp(what, "height")){
     for(i=0;i<fit_data->nfit;i++){
       values[i] = fit_data->fit[i].params[HEIGHT];
diff --git a/storm_analysis/sa_library/multi_fit.h b/storm_analysis/sa_library/multi_fit.h
index 29afacf6..2c26fc91 100644
--- a/storm_analysis/sa_library/multi_fit.h
+++ b/storm_analysis/sa_library/multi_fit.h
@@ -178,7 +178,6 @@ void mFitNewImage(fitData *, double *);
 void mFitNewPeaks(fitData *, int);
 double mFitPeakBgSum(fitData *, peakData *);
 double mFitPeakFgSum(fitData *, peakData *);
-double mFitPeakFgSumSensitivityCorrected(fitData *, peakData *);
 double mFitPeakSum(fitData *, peakData *);
 void mFitRecenterPeaks(fitData *);
 void mFitRemoveErrorPeaks(fitData *);
diff --git a/storm_analysis/sa_library/sa_h5py.py b/storm_analysis/sa_library/sa_h5py.py
index 85a68cb4..4de1fefb 100644
--- a/storm_analysis/sa_library/sa_h5py.py
+++ b/storm_analysis/sa_library/sa_h5py.py
@@ -338,13 +338,28 @@ def getGroupName(self, frame_number):
         return "/fr_" + str(frame_number)
 
     def getLocalizations(self, drift_corrected = False, fields = None):
+        """
+        Return all of the localizations in the file as a single dictionary
+        of numpy arrays.
+
+        drift_corrected defaults to False here, so x, y and z come back
+        exactly as they are stored. Note that localizationsIterator()
+        defaults to True, so the two ways of reading a file do not give you
+        the same coordinates unless you say which you want.
+        """
         return self.getLocalizationsInFrameRange(0,
                                                  self.hdf5.attrs['movie_l'],
                                                  drift_corrected = drift_corrected,
                                                  fields = fields)
     
     def getLocalizationsInFrame(self, frame_number, drift_corrected = False, fields = None):
+        """
+        Return the localizations in a single frame.
 
+        This is where the drift correction is actually applied, by adding
+        the frame's dx, dy and dz attributes. drift_corrected defaults to
+        False, so the stored values are returned unchanged.
+        """
         locs = {}
         grp = self.getGroup(frame_number)
         
@@ -364,6 +379,9 @@ def getLocalizationsInFrame(self, frame_number, drift_corrected = False, fields
     def getLocalizationsInFrameRange(self, start, stop, drift_corrected = False, fields = None):
         """
         Return the localizations in the range start <= frame number < stop.
+
+        As with getLocalizations(), drift_corrected defaults to False and
+        the stored coordinates are returned unchanged.
         """
         assert(stop > start)
         locs = {}
@@ -528,6 +546,13 @@ def localizationsIterator(self, drift_corrected = True, fields = None, skip_empt
 
         If skip_empty is false you will get empty locs dictionaries for frames
         that have no localizations.
+
+        Note that drift_corrected defaults to True here, unlike
+        getLocalizations() and getLocalizationsInFrame(), which default to
+        False. Corrected coordinates are usually what you want when you are
+        just walking a finished file, but it does mean that switching
+        between loading and iterating silently changes whether dx, dy and dz
+        have been applied.
         """
         # This should be a zero length generator if there are no localizations.
         for i in range(self.getMovieLength()):
@@ -605,6 +630,10 @@ def tracksIterator(self, fields = None):
 
         The track center in x,y,z are the 'x', 'y' and 'z' fields. The other
         fields may need to be normalized by the track length.
+
+        There is no drift_corrected argument because tracks are built from
+        drift corrected localizations, so their coordinates are already
+        corrected and there is no way to get back the raw values.
         """
         # This should be a zero length generator if there are no tracks.
         if (not self.hasTracks()):
diff --git a/storm_analysis/simulator/astigmaticPSF.py b/storm_analysis/simulator/astigmaticPSF.py
index fca5343e..9fc7a77d 100644
--- a/storm_analysis/simulator/astigmaticPSF.py
+++ b/storm_analysis/simulator/astigmaticPSF.py
@@ -4,6 +4,9 @@
 of objects that emulates an astigmatic PSF in a format
 compatible with drawgaussians.
 
+x and y are in pixels and z is in microns, which is the convention the
+rest of the simulator uses and what h5_data['z'] holds.
+
 Hazen 01/16
 """
 
@@ -13,7 +16,9 @@
 
 psf_type = "astigmatic"
 
-# Parameters for astigmatic PSF.
+# Parameters for astigmatic PSF. calcSxSy() compares z against elements 1
+# and 2, so those are an offset of 150nm and a defocusing scale of 400nm,
+# both in microns like z.
 wx_params = numpy.array([2.0,  0.150, 0.40, 0.0, 0.0])
 wy_params = numpy.array([2.0, -0.150, 0.40, 0.0, 0.0])
 
@@ -21,7 +26,7 @@ def PSF(x, y, z, h):
     num_objects = x.size
     objects = numpy.zeros((num_objects, 5))
     for i in range(num_objects):
-        [sx, sy] = fitzC.calcSxSy(wx_params, wy_params, z[i] * 0.001)
+        [sx, sy] = fitzC.calcSxSy(wx_params, wy_params, z[i])
         objects[i,:] = [x[i], y[i], h[i], sx, sy]
 
     return objects
@@ -29,7 +34,7 @@ def PSF(x, y, z, h):
 def PSFIntegral(z, h):
     integral = numpy.zeros(z.size)
     for i in range(z.size):
-        [sx, sy] = fitzC.calcSxSy(wx_params, wy_params, z[i] * 0.001)
+        [sx, sy] = fitzC.calcSxSy(wx_params, wy_params, z[i])
         integral[i] = 2.0 * numpy.pi * h[i] * sx * sy
     return integral
 
diff --git a/storm_analysis/simulator/dhPSF.py b/storm_analysis/simulator/dhPSF.py
index 69822d3d..7c3daa9d 100644
--- a/storm_analysis/simulator/dhPSF.py
+++ b/storm_analysis/simulator/dhPSF.py
@@ -4,6 +4,9 @@
 of objects that emulates a double-helical PSF in a format
 compatible with drawgaussians.
 
+x and y are in pixels and z is in microns, which is the convention the
+rest of the simulator uses and what h5_data['z'] holds.
+
 Hazen 01/16
 """
 
@@ -11,9 +14,10 @@
 
 psf_type = "double-helix"
 
-# Double helix PSF range.
-z_min = -500.0
-z_max = 500.0
+# Double helix PSF range, in microns like z. The lobes rotate through 90
+# degrees between z_min and z_max.
+z_min = -0.5
+z_max = 0.5
 
 sigma = 1.0
 def PSF(x, y, z, h):
diff --git a/storm_analysis/spliner/cubic_spline_c.py b/storm_analysis/spliner/cubic_spline_c.py
index c6daad51..e87ef557 100644
--- a/storm_analysis/spliner/cubic_spline_c.py
+++ b/storm_analysis/spliner/cubic_spline_c.py
@@ -143,7 +143,7 @@ def getPSF(self, dy, dx):
         return psf
 
     def py_f(self, y, x):
-        return self.py_spline.f(x, y)
+        return self.py_spline.f(y, x)
 
 
 class CSpline3D(CSpline):
diff --git a/storm_analysis/test/test_analysis_io.py b/storm_analysis/test/test_analysis_io.py
index 87fb9cc5..7842a119 100644
--- a/storm_analysis/test/test_analysis_io.py
+++ b/storm_analysis/test/test_analysis_io.py
@@ -28,7 +28,17 @@ def test_cal_v0():
                      numpy.transpose(variance),
                      numpy.transpose(gain)], fp)
 
-    [o, v, g, r] = analysisIO.loadCMOSCalibration(cal_file)
+    cal = analysisIO.loadCMOSCalibration(cal_file)
+
+    # The doc string promises a list, and the v1 and v2 branches return one.
+    # v0 returned map(), which is not indexable and is exhausted by a single
+    # pass.
+    assert(isinstance(cal, list))
+    assert(cal[0].shape == cal_size)
+    assert(len(list(cal)) == 4)
+    assert(len(list(cal)) == 4)
+
+    [o, v, g, r] = cal
     assert(numpy.allclose(offset, o))
     assert(numpy.allclose(variance, v))
     assert(numpy.allclose(gain, g))
diff --git a/storm_analysis/test/test_dataio.py b/storm_analysis/test/test_dataio.py
index fa8f1d52..f829545d 100644
--- a/storm_analysis/test/test_dataio.py
+++ b/storm_analysis/test/test_dataio.py
@@ -211,10 +211,69 @@ def test_io_6():
     assert(ml == 1)
     assert(numpy.allclose(data, rd.loadAFrame(0)))
 
-    
+def test_io_7():
+    """
+    DaxReader on an incomplete .inf file, and on a missing .dax.
+
+    'number of frames' and the endianness had no fallback, so a .inf file
+    without them gave AttributeError from inside loadAFrame(). A missing
+    .dax constructed a reader that then failed with AttributeError on a
+    None file pointer.
+    """
+    import pytest
+
+    movie_h = 8
+    movie_w = 16
+    movie_l = 5
+
+    data = numpy.random.randint(0, 60000, (movie_h, movie_w)).astype(numpy.uint16)
+
+    movie_name = storm_analysis.getPathOutputTest("test_dataio_inf.dax")
+    inf_name = storm_analysis.getPathOutputTest("test_dataio_inf.inf")
+
+    wr = datawriter.inferWriter(movie_name)
+    for i in range(movie_l):
+        wr.addFrame(data)
+    wr.close()
+
+    with open(inf_name) as fp:
+        inf_lines = fp.readlines()
+
+    def writeInf(drop):
+        with open(inf_name, "w") as fp:
+            for line in inf_lines:
+                if (drop is None) or (not drop in line):
+                    fp.write(line)
+
+    ## No 'number of frames', the file size gives it instead.
+    writeInf("number of frames")
+    rd = datareader.inferReader(movie_name)
+    assert(rd.filmSize() == [movie_w, movie_h, movie_l])
+    assert(numpy.allclose(data, rd.loadAFrame(0)))
+    rd.close()
+
+    ## No endianness, little endian is assumed.
+    writeInf("data type")
+    rd = datareader.inferReader(movie_name)
+    assert(numpy.allclose(data, rd.loadAFrame(0)))
+    rd.close()
+
+    ## A complete .inf but no .dax, and the error names the file.
+    missing_name = storm_analysis.getPathOutputTest("test_dataio_missing.dax")
+    storm_analysis.removeFile(missing_name)
+    with open(storm_analysis.getPathOutputTest("test_dataio_missing.inf"), "w") as fp:
+        for line in inf_lines:
+            fp.write(line)
+
+    with pytest.raises(IOError) as err:
+        datareader.inferReader(missing_name)
+    assert("test_dataio_missing.dax" in str(err.value))
+
+
 if (__name__ == "__main__"):
     test_io_1()
     test_io_2()
     test_io_3()
     test_io_6()
+    test_io_7()
     
diff --git a/storm_analysis/test/test_fitting.py b/storm_analysis/test/test_fitting.py
index 93f8300e..da4b1a15 100644
--- a/storm_analysis/test/test_fitting.py
+++ b/storm_analysis/test/test_fitting.py
@@ -84,6 +84,33 @@ def test_fitting_peak_mask_square_aoi():
     assert (numpy.count_nonzero(mask[:10+margin,:]) == 0)
 
 
+def test_fitting_peak_mask_stop():
+    """
+    The x_stop and y_stop edges of a rectangular AOI.
+
+    Both slices ended at -1, which stops one element short, so the last
+    column and the last row kept the value 1.0 instead of being masked.
+    """
+    margin = 5
+    [x_stop, y_stop] = [70, 30]
+
+    p = params.ParametersDAO()
+    p.setAttr("x_stop", "int", x_stop)
+    p.setAttr("y_stop", "int", y_stop)
+
+    # More columns than rows, so the two axes cannot be confused.
+    shape = (60, 120)
+    mask = fitting.peakMask(shape, p, margin)
+
+    # Nothing at or beyond the stops survives, including the final
+    # column and the final row.
+    assert (numpy.count_nonzero(mask[:,x_stop+margin:]) == 0)
+    assert (numpy.count_nonzero(mask[y_stop+margin:,:]) == 0)
+
+    # And the region inside both stops is untouched.
+    assert (numpy.all(mask[:y_stop+margin,:x_stop+margin] == 1.0))
+
+
 def test_fitting_check_mode():
     """
     check_mode writes diagnostic images from the arbitrary PSF peak finder,
@@ -125,4 +152,5 @@ def test_fitting_check_mode():
     test_fitting_no_peak_finder()
     test_fitting_peak_mask_circular()
     test_fitting_peak_mask_square_aoi()
+    test_fitting_peak_mask_stop()
     test_fitting_check_mode()
diff --git a/storm_analysis/test/test_sa_h5py.py b/storm_analysis/test/test_sa_h5py.py
index e29e0b1b..ed039340 100644
--- a/storm_analysis/test/test_sa_h5py.py
+++ b/storm_analysis/test/test_sa_h5py.py
@@ -572,6 +572,54 @@ def test_sa_h5py_19():
             assert not elt in locs
     
         
+def test_sa_h5py_20():
+    """
+    Pin the drift_corrected defaults.
+
+    getLocalizations() and getLocalizationsInFrame() default to False and
+    return the stored coordinates. localizationsIterator() defaults to True
+    and returns corrected ones. The defaults suit each method's primary use,
+    but they disagree, so they are worth holding still: harmonizing them
+    would silently change what every existing caller gets back.
+    """
+    [dx, dy, dz] = [1.0, -2.0, 0.5]
+
+    peaks = {"x" : numpy.zeros(3),
+             "y" : numpy.zeros(3),
+             "z" : numpy.zeros(3)}
+
+    h5_name = storm_analysis.getPathOutputTest("test_sa_hdf5_drift.hdf5")
+    storm_analysis.removeFile(h5_name)
+
+    with saH5Py.SAH5Py(h5_name, is_existing = False, overwrite = True) as h5:
+        h5.setMovieInformation(100, 100, 1, "")
+        h5.addLocalizations(peaks, 0)
+        h5.setDriftCorrection(0, dx = dx, dy = dy, dz = dz)
+
+    with saH5Py.SAH5Py(h5_name) as h5:
+
+        # The getters hand back what is stored.
+        for locs in [h5.getLocalizations(),
+                     h5.getLocalizationsInFrame(0),
+                     h5.getLocalizationsInFrameRange(0, 1)]:
+            assert(numpy.allclose(locs["x"], numpy.zeros(3)))
+            assert(numpy.allclose(locs["y"], numpy.zeros(3)))
+            assert(numpy.allclose(locs["z"], numpy.zeros(3)))
+
+        # The iterator corrects.
+        for fnum, locs in h5.localizationsIterator():
+            assert(numpy.allclose(locs["x"], dx*numpy.ones(3)))
+            assert(numpy.allclose(locs["y"], dy*numpy.ones(3)))
+            assert(numpy.allclose(locs["z"], dz*numpy.ones(3)))
+
+        # And both are explicitly overridable.
+        locs = h5.getLocalizations(drift_corrected = True)
+        assert(numpy.allclose(locs["x"], dx*numpy.ones(3)))
+
+        for fnum, locs in h5.localizationsIterator(drift_corrected = False):
+            assert(numpy.allclose(locs["x"], numpy.zeros(3)))
+
+
 if (__name__ == "__main__"):
     test_sa_h5py_1()
     test_sa_h5py_2()
@@ -592,4 +640,5 @@ def test_sa_h5py_19():
     test_sa_h5py_17()
     test_sa_h5py_18()
     test_sa_h5py_19()
+    test_sa_h5py_20()
     
diff --git a/storm_analysis/test/test_simulator.py b/storm_analysis/test/test_simulator.py
index 8d474cd0..9f3e15ba 100644
--- a/storm_analysis/test/test_simulator.py
+++ b/storm_analysis/test/test_simulator.py
@@ -7,8 +7,10 @@
 
 import storm_analysis.spliner.spline_to_psf as splineToPSF
 
+import storm_analysis.simulator.astigmaticPSF as astigmaticPSF
 import storm_analysis.simulator.background as background
 import storm_analysis.simulator.camera as camera
+import storm_analysis.simulator.dhPSF as dhPSF
 import storm_analysis.simulator.photophysics as photophysics
 import storm_analysis.simulator.psf as psf
 import storm_analysis.simulator.simulate as simulate
@@ -149,6 +151,87 @@ def test_simulate_3():
     sim.simulate(dax_name, bin_name, 5)
 
     
+def test_astigmatic_psf():
+    """
+    astigmaticPSF takes z in microns, like everything else in the simulator.
+
+    It used to multiply z by 0.001 before calling calcSxSy(), which says the
+    incoming z is nanometers. Fed the simulator's own convention that put
+    every z within a nanometer of focus, so it returned a round,
+    z-independent PSF, sx = sy = 1.068 px at every z.
+    """
+    z = numpy.array([-0.4, -0.2, 0.0, 0.2, 0.4])
+    x = numpy.ones(z.size)
+    y = numpy.ones(z.size)
+    h = numpy.ones(z.size)*100.0
+
+    objects = astigmaticPSF.PSF(x, y, z, h)
+    sx = objects[:,3]
+    sy = objects[:,4]
+
+    # Round at focus.
+    assert(abs(sx[2] - sy[2]) < 1.0e-9)
+
+    # Elongated in x below focus, in y above it, and it is a real effect and
+    # not rounding. wx_params puts the x waist at +150nm and the y waist at
+    # -150nm, with a 400nm defocusing scale.
+    assert(sx[0] > 1.4*sy[0])
+    assert(sy[4] > 1.4*sx[4])
+
+    # Symmetric about focus.
+    assert(abs(sx[0] - sy[4]) < 1.0e-9)
+    assert(abs(sy[0] - sx[4]) < 1.0e-9)
+
+    # And the width actually varies with z, which is what a Z fitter needs.
+    assert(sx[0] > 1.5*sx[2])
+
+    # PSFIntegral uses the same widths.
+    integral = astigmaticPSF.PSFIntegral(z, h)
+    assert(numpy.allclose(integral, 2.0*numpy.pi*h*sx*sy))
+
+
+def test_dh_psf():
+    """
+    dhPSF takes z in microns too.
+
+    z_min and z_max were -500 and 500, which are nanometers, so feeding it
+    the simulator's own z froze the helix angle at 45 degrees across the
+    whole range. A double helix PSF that does not rotate carries no z
+    information at all, the same degeneracy astigmaticPSF had.
+    """
+    [xc, yc] = [10.0, 20.0]
+    z = numpy.array([-0.5, 0.0, 0.5])
+
+    objects = dhPSF.PSF(numpy.ones(z.size)*xc,
+                        numpy.ones(z.size)*yc,
+                        z,
+                        numpy.ones(z.size)*100.0)
+
+    # Two lobes per emitter, ordered [emitter0 +, emitter0 -, emitter1 +, ..].
+    assert(objects.shape[0] == 2*z.size)
+
+    lobe = objects[::2,:]
+    dx = lobe[:,0] - xc
+    dy = lobe[:,1] - yc
+
+    r = 2.0*dhPSF.sigma
+
+    # At z_min the lobes lie along x, at z_max along y, and at focus the
+    # angle is halfway between.
+    assert(numpy.allclose([dx[0], dy[0]], [r, 0.0], atol = 1.0e-9))
+    assert(numpy.allclose([dx[2], dy[2]], [0.0, r], atol = 1.0e-9))
+    assert(abs(dx[1] - dy[1]) < 1.0e-9)
+
+    # The angle has to actually sweep, which is the whole point.
+    angles = numpy.degrees(numpy.arctan2(dy, dx))
+    assert(numpy.allclose(angles, [0.0, 45.0, 90.0], atol = 1.0e-6))
+
+    # The second lobe of each pair is opposite the first.
+    other = objects[1::2,:]
+    assert(numpy.allclose(other[:,0] - xc, -dx))
+    assert(numpy.allclose(other[:,1] - yc, -dy))
+
+
 if (__name__ == "__main__"):
     test_psf_spline2D_1()
     test_psf_spline2D_2()
@@ -157,3 +240,5 @@ def test_simulate_3():
     test_simulate_1()
     test_simulate_2()
     test_simulate_3()
+    test_astigmatic_psf()
+    test_dh_psf()
diff --git a/storm_analysis/test/test_spline.py b/storm_analysis/test/test_spline.py
index 23a192a4..172bf7b5 100644
--- a/storm_analysis/test/test_spline.py
+++ b/storm_analysis/test/test_spline.py
@@ -171,6 +171,54 @@ def test_psf_3D_dz():
         assert (abs(py_spline.dzf(x, y, z) - c_spline.dzf(x, y, z)) < 1.0e-6)
 
 
+def test_psf_2D_py_f():
+    """
+    py_f() is the Python evaluation of the same point that f() evaluates in
+    C, so the two have to agree. CSpline2D.py_f() passed its arguments on in
+    the opposite order, which evaluates the transposed point. A PSF spline is
+    not symmetric in x and y, so the answers differ.
+    """
+    if (sys.version_info < (3, 0)):
+        return
+
+    spline_filename = storm_analysis.getData("test/data/test_spliner_psf_2d.spline")
+    with open(spline_filename, "rb") as fp:
+        spline_data = pickle.load(fp)
+
+    py_spline = spline2D.Spline2D(spline_data["spline"], spline_data["coeff"])
+    c_spline = cubicSplineC.CSpline2D(py_spline)
+
+    size = py_spline.getSize() - 1.0e-6
+
+    for i in range(reps):
+        x = random.uniform(1.0e-6, size)
+        y = random.uniform(1.0e-6, size)
+        assert (abs(c_spline.py_f(y, x) - c_spline.f(y, x)) < 1.0e-6)
+
+def test_psf_3D_py_f():
+    """
+    The same check on CSpline3D, which has always passed its arguments
+    through in the right order.
+    """
+    if (sys.version_info < (3, 0)):
+        return
+
+    spline_filename = storm_analysis.getData("test/data/test_spliner_psf.spline")
+    with open(spline_filename, "rb") as fp:
+        spline_data = pickle.load(fp)
+
+    py_spline = spline3D.Spline3D(spline_data["spline"], spline_data["coeff"])
+    c_spline = cubicSplineC.CSpline3D(py_spline)
+
+    size = py_spline.getSize() - 1.0e-6
+
+    for i in range(reps):
+        x = random.uniform(1.0e-6, size)
+        y = random.uniform(1.0e-6, size)
+        z = random.uniform(1.0e-6, size)
+        assert (abs(c_spline.py_f(z, y, x) - c_spline.f(z, y, x)) < 1.0e-6)
+
+
 if (__name__ == "__main__"):
     test_psf_2D_f()
     test_psf_2D_dx()
@@ -179,3 +227,5 @@ def test_psf_3D_dz():
     test_psf_3D_dx()
     test_psf_3D_dy()
     test_psf_3D_dz()
+    test_psf_2D_py_f()
+    test_psf_3D_py_f()