diff --git a/storm_analysis/decon_storm/mlem_sparse.c b/storm_analysis/decon_storm/mlem_sparse.c
index 46398bb2..d580b809 100644
--- a/storm_analysis/decon_storm/mlem_sparse.c
+++ b/storm_analysis/decon_storm/mlem_sparse.c
@@ -509,7 +509,7 @@ void backward(decon *dec)
  * Calculate ej and xj w/ compression parameter 
  * following Lingenfelter L1 norm.
  *
- * cmpr - compression parameter
+ * compression - compression parameter
  */
 void backwardCompressed(decon *dec, double compression)
 {
@@ -539,7 +539,7 @@ void backwardCompressed(decon *dec, double compression)
  * following Lingenfelter L1 norm, but do not
  * update the background.
  *
- * cmpr - compression parameter
+ * compression - compression parameter
  */
 void backwardCompressedFixedBg(decon *dec, double compression)
 {
diff --git a/storm_analysis/multi_plane/psf_localizations.py b/storm_analysis/multi_plane/psf_localizations.py
index bbc9fe0b..bb3aeb4c 100644
--- a/storm_analysis/multi_plane/psf_localizations.py
+++ b/storm_analysis/multi_plane/psf_localizations.py
@@ -1,8 +1,8 @@
 #!/usr/bin/env python
 """
-Giving a mapping file (from multi_plane.mapper), and a 
-molecule list, generate molecule lists to use for the
-PSF extraction step.
+Given a mapping file (from multi_plane.mapper) and a localizations
+HDF5 file, generate the per channel localization files to use for
+the PSF extraction step.
 
 Hazen 05/17
 """
@@ -127,7 +127,7 @@ def psfLocalizations(h5_filename, mapping_filename, frame = 0, aoi_size = 8, min
     parser.add_argument('--map', dest='mapping', type=str, required=True,
                         help = "The name of the mapping file. This is the output of multi_plane.mapper.")
     parser.add_argument('--frame', dest='frame', type=int, required=False, default=0,
-                        help = "The frame in .bin file to get the localizations from. The default is 0.")
+                        help = "The frame of the localizations file to get the localizations from. The default is 0.")
     parser.add_argument('--aoi_size', dest='aoi_size', type=int, required=False, default=8,
                         help = "The size of the area of interest around the bead in pixels. The default is 8.")
     parser.add_argument('--min_height', dest='min_height', type=float, required=False, default = 0.0,
diff --git a/storm_analysis/sa_library/dao_fit.h b/storm_analysis/sa_library/dao_fit.h
index 04dd3340..52181c9f 100644
--- a/storm_analysis/sa_library/dao_fit.h
+++ b/storm_analysis/sa_library/dao_fit.h
@@ -32,8 +32,17 @@ typedef struct daoFit
   double width_min;             /* Minimum allowed width for Gaussian fitting. */
   double width_max;             /* Maximum allowed width for Gaussian fitting. */
   
-  double wx_z_params[5];        /* x width versus z parameters. */
-  double wy_z_params[5];        /* y width versus z parameters. */
+  /*
+   * x and y width versus z parameters, wo, c, d, A and B.
+   *
+   * Note that the calibration can have seven terms, and that
+   * sa_utilities/fitz.c uses all seven. The Z fitting model uses only
+   * these five, so daoInitializeZ() copies five and any C and D terms in
+   * the calibration are dropped without warning. See the wx/wy vs z
+   * parameter documentation in sa_library/parameters.py.
+   */
+  double wx_z_params[5];
+  double wy_z_params[5];
 } daoFit;
 
 
diff --git a/storm_analysis/sa_library/ia_utilities_c.py b/storm_analysis/sa_library/ia_utilities_c.py
index cc3c4228..307df2d3 100644
--- a/storm_analysis/sa_library/ia_utilities_c.py
+++ b/storm_analysis/sa_library/ia_utilities_c.py
@@ -197,8 +197,8 @@ def findMaxima(self, images, want_height = False):
 
         # Verify images are the correct size and C contiguous.
         for image in images:
-            assert (image.shape[0] == self.taken[0].shape[0]), "Unexpected image x size!"
-            assert (image.shape[1] == self.taken[0].shape[1]), "Unexpected image y size!"
+            assert (image.shape[0] == self.taken[0].shape[0]), "Unexpected image y size!"
+            assert (image.shape[1] == self.taken[0].shape[1]), "Unexpected image x size!"
             assert (image.flags['C_CONTIGUOUS']), "Image is not C contiguous!"
 
         # Create pointer array to the images.
diff --git a/storm_analysis/sa_library/matched_filter.c b/storm_analysis/sa_library/matched_filter.c
index f00c335b..d4c827fc 100644
--- a/storm_analysis/sa_library/matched_filter.c
+++ b/storm_analysis/sa_library/matched_filter.c
@@ -61,7 +61,7 @@ filter *initialize(double *, double, int, int, int);
 /*
  * cleanup()
  *
- * filter - A pointer to a filter structure.
+ * flt - A pointer to a filter structure.
  */
 void cleanup(filter *flt)
 {
diff --git a/storm_analysis/sa_library/multi_fit.c b/storm_analysis/sa_library/multi_fit.c
index ca81f46c..61faaa19 100644
--- a/storm_analysis/sa_library/multi_fit.c
+++ b/storm_analysis/sa_library/multi_fit.c
@@ -62,7 +62,6 @@ void mFitAddPeak(fitData *fit_data)
  * Calculate the Anscombe transform.
  *
  * x - Signal / function value.
- * var - (Gaussian) variance.
  */
 double mFitAnscombe(double x)
 {
@@ -387,7 +386,9 @@ int mFitCalcErrDWLS(fitData *fit_data)
 /*
  * mFitCalcErrFWLS()
  *
- * The data weighted least squares version of the error function.
+ * The fit weighted least squares version of the error function. This
+ * weights by the fit value where mFitCalcErrDWLS() weights by the
+ * measurement, so unlike DWLS it fails if the fit goes negative.
  *
  * fit_data - pointer to a fitData structure.
  *
@@ -611,7 +612,7 @@ int mFitDeltaConvergence(fitData *fit_data, int index)
     return 0;
   }
 
-  /* 0.01 delta */    
+  /* 0.0001 delta */    
   if(fabs(old_params[ZCENTER] - cur_params[ZCENTER]) > (100.0 * fit_data->tolerance)){
     return 0;
   }
diff --git a/storm_analysis/sa_library/parameters.py b/storm_analysis/sa_library/parameters.py
index 1007f423..4481693e 100644
--- a/storm_analysis/sa_library/parameters.py
+++ b/storm_analysis/sa_library/parameters.py
@@ -492,6 +492,16 @@ def __init__(self, **kwds):
 ww_doc_string = """wx/wy vs z parameters. Units are nanometers or dimensionless.
 
 See Huang, Science 2008 for a more detailed explanation.
+
+Note that not everything uses all seven terms. sa_utilities/fitz.c, which
+calculates z from wx and wy after fitting for the "3d" model, uses all of
+them. The "Z" fitting model and sa_utilities.fitz_c.calcSxSy() use only
+the first five, wo, c, d, A and B, and silently ignore C and D.
+
+daostorm_3d.z_calibration defaults to --fit_order 2, which produces
+exactly those five, so this only matters if you fit a calibration with
+fit_order 3 or 4. Such a calibration is followed in full by the "3d"
+model's z calculation and truncated everywhere else.
 """
 
 
diff --git a/storm_analysis/sa_utilities/fitz_c.py b/storm_analysis/sa_utilities/fitz_c.py
index 879feb4c..8dc6b00f 100644
--- a/storm_analysis/sa_utilities/fitz_c.py
+++ b/storm_analysis/sa_utilities/fitz_c.py
@@ -44,6 +44,10 @@
 def calcSxSy(wx_params, wy_params, z):
     """
     Return sigma x and sigma y given the z calibration parameters.
+
+    This uses the first five terms only, wo, c, d, A and B. A calibration
+    fitted with fit_order 3 or 4 also has C and D terms, and fitz.c does
+    use them, so the two do not agree on such a calibration.
     """
     zx = (z - wx_params[1])/wx_params[2]
     sx = 0.5 * wx_params[0] * numpy.sqrt(1.0 + zx*zx + wx_params[3]*zx*zx*zx + wx_params[4]*zx*zx*zx*zx)
diff --git a/storm_analysis/simulator/drift.py b/storm_analysis/simulator/drift.py
index e4ca8268..a2c3b73b 100644
--- a/storm_analysis/simulator/drift.py
+++ b/storm_analysis/simulator/drift.py
@@ -22,8 +22,9 @@ def __init__(self, sim_fp, x_size, y_size, h5_data):
 
 class DriftFromFile(Drift):
     """
-    Add drift from a file. X and Y are in units of pixels, Z is in
-    nanometers.
+    Add drift from a file. The file has one row per frame and three
+    columns, x, y and z. X and Y are in units of pixels and Z is in
+    microns, the same units as the localizations it is added to.
     """
     def __init__(self, sim_fp, x_size, y_size, h5_data, drift_file):
         super(DriftFromFile, self).__init__(sim_fp, x_size, y_size, h5_data)
diff --git a/storm_analysis/spliner/README.txt b/storm_analysis/spliner/README.txt
index d9e6ee2f..c240df5f 100644
--- a/storm_analysis/spliner/README.txt
+++ b/storm_analysis/spliner/README.txt
@@ -12,8 +12,8 @@ cubic_spline_c.py - The Python wrapper of the C spline library. Given
    the splines coefficients this can be used to calculate the spline
    and it's derivatives.
 
-find_peaks_fista.py - Peak finding using a compressed sensing approach
-   and the FISTA solver.
+find_peaks_decon.py - Peak finding using a compressed sensing
+   deconvolution approach, with either the FISTA or the ADMM solver.
 
 find_peaks_std.py - Peak finding using a similar approach to 3D-DAOSTORM,
    start with the brightest local maxima and the proceed to dimmer
@@ -25,12 +25,13 @@ hdf5_to_beads.py - Convert an HDF5 file to a beads file for PSF
 measure_psf.py - Measure the PSF of the microscope given a z stack movie
    and the location of the emitters in each frame.
    
-measure_psf_beads.py - This is similar to measure_psf except that it uses
-   a text for input rather than a molecule list file (.bin file).
+measure_psf_beads.py - This is similar to measure_psf except that it
+   takes the bead positions from a text file rather than from an HDF5
+   localization file.
 
 measure_psf_utils.py - Utility functions used for PSF measurement.
 
-offset_to_Z.py - Converts a .off file (from storm-control) to a z_offset
+offset_to_z.py - Converts a .off file (from storm-control) to a z_offset
    file that can be used for measuring the PSF.
 
 print_psf.py - Print out some information about a PSF file.
@@ -45,6 +46,8 @@ spline3D.py - 3D cubic spline in Python.
 
 spline_analysis.py - Run spline analysis on a single STORM movie.
 
+spline_info.py - Print out the details of a spline.
+
 spline_to_psf.py - Generate a PSF at a particular z value given a spline.
 
 
diff --git a/storm_analysis/spliner/cubic_spline.c b/storm_analysis/spliner/cubic_spline.c
index 60962d52..443284d8 100644
--- a/storm_analysis/spliner/cubic_spline.c
+++ b/storm_analysis/spliner/cubic_spline.c
@@ -3,13 +3,21 @@
  * 
  * Notes:
  *
- *  1. Nominally xsize, ysize and zsize could be different, but this whether
- *     or not this works correctly has not been tested.
+ *  1. Nominally xsize, ysize and zsize could be different, but whether or
+ *     not this works correctly has not been tested.
  *
  *  2. Unlike the Python version, this will not return the correct value if
  *     x is exactly at the maximum value that the spline covers.
  *
- *  3. This library is thread safe.. Pretty sure..
+ *  3. This library is thread safe provided that each thread has its own
+ *     splineData. Sharing one between threads is not safe: computeDelta2D()
+ *     and computeDelta3D() write their results into the delta_f, delta_dxf
+ *     and delta_dyf / delta_dzf arrays that live in the struct, and the
+ *     f/dxf/dyf/dzf At2D() and At3D() accessors read them back, so two
+ *     threads evaluating the same splineData will overwrite each other's
+ *     deltas between the write and the read. initSpline2D() and
+ *     initSpline3D() copy the coefficients into the struct, so separate
+ *     splineData share nothing and need no locking.
  *
  * Hazen 11/16
  *
@@ -200,7 +208,7 @@ double dxfAt2D(splineData *spline_data, int yc, int xc)
  * yc - The y coordinate (integer).
  * xc - The x coordinate (integer).
  *
- * Return - The derivative in y.
+ * Return - The derivative in x.
  */
 double dxfAt3D(splineData *spline_data, int zc, int yc, int xc)
 {
@@ -314,7 +322,7 @@ double dyfAt2D(splineData *spline_data, int yc, int xc)
 /*
  * dyfAt3D()
  *
- * Compute the derivative of the in spline y at x,y,z coordinate (integer). 
+ * Compute the derivative of the spline in y at x,y,z coordinate (integer). 
  * In order for this to work correctly computeDelta3D should have already 
  * been called.
  *
@@ -410,7 +418,7 @@ double dyfSpline3D(splineData *spline_data, double z, double y, double x)
 /*
  * dzfAt3D()
  *
- * Compute the derivative of the in spline z at x,y,z coordinate (integer). 
+ * Compute the derivative of the spline in z at x,y,z coordinate (integer). 
  * In order for this to work correctly computeDelta3D should have already 
  * been called.
  *
@@ -419,7 +427,7 @@ double dyfSpline3D(splineData *spline_data, double z, double y, double x)
  * yc - The y coordinate (integer).
  * xc - The x coordinate (integer).
  *
- * Return - The derivative in y.
+ * Return - The derivative in z.
  */
 double dzfAt3D(splineData *spline_data, int zc, int yc, int xc)
 {
diff --git a/storm_analysis/spliner/measure_psf.py b/storm_analysis/spliner/measure_psf.py
index d5b7a2d6..cd6fa29e 100644
--- a/storm_analysis/spliner/measure_psf.py
+++ b/storm_analysis/spliner/measure_psf.py
@@ -54,7 +54,7 @@ def measurePSF(movie_name, zfile_name, movie_h5_name, psf_name, want2d = False,
     # Create z scaling object.
     z_sclr = measurePSFUtils.ZScaler(z_range, z_step)
     
-    # Load dax file, z offset file and molecule list file.
+    # Load the movie, the z offset file and the localizations.
     dax_data = datareader.inferReader(movie_name)
     z_off = None
     if os.path.exists(zfile_name):
@@ -207,7 +207,7 @@ def measurePSF(movie_name, zfile_name, movie_h5_name, psf_name, want2d = False,
 
     import argparse
     
-    parser = argparse.ArgumentParser(description = 'Measure PSF given a movie, a list.bin file and (optionally) a z_offset file')
+    parser = argparse.ArgumentParser(description = 'Measure PSF given a movie, a localizations HDF5 file and (optionally) a z_offset file')
 
     parser.add_argument('--movie', dest='movie', type=str, required=True,
                         help = "The name of the movie to analyze, can be .dax, .tiff or .spe format.")
diff --git a/storm_analysis/spliner/measure_psf_utils.py b/storm_analysis/spliner/measure_psf_utils.py
index 91f87c0a..95e02fe5 100644
--- a/storm_analysis/spliner/measure_psf_utils.py
+++ b/storm_analysis/spliner/measure_psf_utils.py
@@ -25,7 +25,7 @@ def __init__(self, z_range, z_step):
         assert(z_step > 0.0), "The z step must be positive."
         assert(z_range >= z_step), "The z range must be greater than or equal to the step size."
     
-        # Assert that the z_step size is a multiple of the z_range.
+        # Assert that the z_range is a multiple of the z_step size.
         assert ((int(z_range*1.0e+3) % int(z_step*1.0e+3)) == 0), "The z range must be a multiple of the z step."
 
         self.z_mid = int(round(z_range/z_step))