diff --git a/CMakeLists.txt b/CMakeLists.txt
index a14956616..5d7e5cb8c 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -49,7 +49,7 @@ endif()
 option(ENABLE_NBODY "Enable EXP n-body" ON)
 option(ENABLE_PYEXP "Enable the Python bindings" ON)
 option(ENABLE_PYEXP_ONLY "Python bindings and support libraries only" OFF)
-option(ENABLE_UTILS "Enable build of the EXP standalone utilities" OFF)
+option(ENABLE_UTILS "Enable build of the EXP standalone utilities" ON)
 option(ENABLE_PNG "Enable PNG graphics support" FALSE)
 option(ENABLE_CUDA "Enable CUDA" FALSE)
 option(ENABLE_SLURM "Enable SLURM checkpointing support" FALSE)
diff --git a/expui/BasisFactory.H b/expui/BasisFactory.H
index a1db14173..ca71e4ab6 100644
--- a/expui/BasisFactory.H
+++ b/expui/BasisFactory.H
@@ -184,13 +184,13 @@ namespace BasisClasses
     { return coefret; }
 
     //! Accumulate new coefficients from full phase space
-    virtual void accumulate(double mass,
-			    double x, double y, double z,
-			    double u, double v, double w) = 0;
+    virtual void accumulate(double x, double y, double z,
+			    double u, double v, double w,
+			    double mass, unsigned long int indx) = 0;
     
     //! Accumulate new coefficients from coordinates only
     virtual void accumulate(double x, double y, double z,
-			    double mass) = 0;
+			    double mass, unsigned long int indx) = 0;
     
     //! Get mass on grid
     double getMass(void) { return totalMass; }
diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H
index 7b6699d35..133f2a31e 100644
--- a/expui/BiorthBasis.H
+++ b/expui/BiorthBasis.H
@@ -3,11 +3,18 @@
 
 #include <functional>
 #include <tuple>
+#include <map>
+#include <set>
 
 #include <Eigen/Eigen>
 #include <unsupported/Eigen/CXX11/Tensor> // For 3d rectangular grids
 #include <yaml-cpp/yaml.h>
 
+#include <highfive/H5File.hpp>
+#include <highfive/H5DataSet.hpp>
+#include <highfive/H5DataSpace.hpp>
+#include <highfive/H5Attribute.hpp>
+
 #include "DiskDensityFunc.H"
 #include "ParticleReader.H"
 #include "Coefficients.H"
@@ -96,8 +103,52 @@ namespace BasisClasses
     RowMatrixXd varray;
     //@}
 
+    //! Coefficient variance computation enabled
+    bool pcavar = false;
+
+    //! Data type for covariance
+    bool floatType = false;
+
+    //@{
+    //! Sample counts and masses for covariance computation
+    Eigen::VectorXi sampleCounts;
+    Eigen::VectorXd sampleMasses;
+    //@}
+
+    //@{
+    //! HDF5 compression settings
+    unsigned H5compress  = 5;
+    unsigned H5chunk     = 1024*1024; // (1 MB)
+    bool     H5shuffle   = true;
+    bool     H5szip      = false;
+    //@}
+
+    //! Round time key to emulated fixed-point arithmetic
+    double roundTime(double time)
+    {
+      // Eight decimal places should be enough here...
+      const double multiplier = 1.0e+08;  // std::pow(10.0, 8);
+      return std::floor(time * multiplier + 0.5) / multiplier;
+    }
+
+    //! Write H5 covariance; returns updated group count
+    virtual unsigned writeCovarH5(HighFive::Group& group, unsigned count, double time);
+    
+    //! Write H5 parameters for coefficient covariance
+    virtual void writeCovarH5Params(HighFive::File& file) {}
+
+    //! Add coefficient covariance data to an HDF5 file
+    virtual void extendCoefCovariance(const std::string& filename, double time);
+
+    //! Variance file versioning (Version 1.0 has more granularity and
+    //! poor compression)
+    inline static const std::string CovarianceFileVersion = "1.1";
+
   public:
     
+    //! Basis identifier string
+    std::string BasisID = "BiorthBasis";
+
     //! Constructor from YAML node
     BiorthBasis(const YAML::Node& conf,
 		const std::string& name="BiorthBasis") : Basis(conf, name) {}
@@ -166,13 +217,14 @@ namespace BasisClasses
     }
     
     //! Accumulate new coefficients
-    virtual void accumulate(double x, double y, double z, double mass) = 0;
+    virtual void accumulate(double x, double y, double z, double mass,
+			    unsigned long indx = 0) = 0;
     
     //! Accumulate new coefficients
-    virtual void accumulate(double mass,
-			    double x, double y, double z,
-			    double u, double v, double w)
-    { accumulate(x, y, z, mass); }
+    virtual void accumulate(double x, double y, double z,
+			    double u, double v, double w,
+			    double mass, unsigned long indx = 0)
+    { accumulate(x, y, z, mass, indx); }
     
     //! Get mass on grid
     double getMass(void) { return totalMass; }
@@ -235,6 +287,51 @@ namespace BasisClasses
       return varray;
     }
 
+    //! Return a vector of tuples of basis functions and the
+    //! covariance matrix for subsamples of particles
+    using CoefCovarType = std::tuple<Eigen::VectorXcd, Eigen::MatrixXcd>;
+    virtual std::vector<std::vector<CoefCovarType>>
+    getCoefCovariance()
+    {
+      // Must be overriden; base implementation throws error
+      throw std::runtime_error("BiorthBasis::getCoefCovariance: Not implemented for this basis");
+    }
+
+    //! Get sample counts for the covariance computation
+    virtual std::tuple<Eigen::VectorXi, Eigen::VectorXd>
+    getCovarSamples()
+    {
+      return std::make_tuple(sampleCounts, sampleMasses);
+    }
+
+    //! Write coefficient covariance data to an HDF5 file
+    void writeCoefCovariance(const std::string& compname,
+			     const std::string& runtag,
+			     double time=0.0);
+
+    //! Choose between float and double storage for covariance
+    void setCovarFloatType(bool flt) { floatType = flt; }
+
+    //! HDF5 compression settings
+    void setCovarH5Compress(unsigned level, unsigned chunksize, bool shuffle, bool szip=false)
+    {
+      H5compress = level;
+      H5chunk    = chunksize;
+      H5shuffle  = shuffle;
+      H5szip     = szip;
+    }
+
+    //! Make covariance after accumulation
+    virtual void makeCoefCovariance(void) {}
+
+    //! Enable covariance computation with optional sample time
+    virtual void enableCoefCovariance(bool pcavar, int sampT_in)
+    {
+      // Must be overriden; base implementation throws error
+      throw std::runtime_error("BiorthBasis::enableCoefCovariance: "
+			       "Not implemented for this basis");
+    }
+
     //! Evaluate acceleration in Cartesian coordinates in centered
     //! coordinate system for a collections of points
     virtual RowMatrixXd& getAccel(Eigen::Ref<const RowMatrixXd> pos)
@@ -349,6 +446,17 @@ namespace BasisClasses
     virtual void computeAccel(double x, double y, double z, 
 			      Eigen::Ref<Eigen::Vector3d> vstore);
 
+    //@{
+    //! Covariance structures.  First index is T, second is harmonic
+    //! (l,m).
+    std::vector<std::vector<Eigen::VectorXcd>> meanV;
+    std::vector<std::vector<Eigen::MatrixXcd>> covrV;
+    void init_covariance();
+    void zero_covariance();
+    void writeCovarH5Params(HighFive::File& file);
+    int sampT = 100;
+    //@}
+
   public:
     
     //! Constructor from YAML node
@@ -377,7 +485,8 @@ namespace BasisClasses
     void make_coefs(void);
     
     //! Accumulate new coefficients
-    virtual void accumulate(double x, double y, double z, double mass);
+    virtual void accumulate(double x, double y, double z, double mass,
+			    unsigned long indx = 0);
     
     //! Return current maximum harmonic order in expansion
     int getLmax() { return lmax; }
@@ -404,6 +513,17 @@ namespace BasisClasses
       return ret;
     }
 
+    //! Return a vector of tuples of basis functions and the
+    //! covariance matrix for subsamples of particles
+    std::vector<std::vector<CoefCovarType>>
+    getCoefCovariance();
+
+    virtual void enableCoefCovariance(bool pcavar_in, int sampT_in=100)
+    {
+      pcavar = pcavar_in;
+      sampT  = sampT_in;
+      if (pcavar) init_covariance();
+    }
     //! Create and the coefficients from a function callback with the
     //! provided time value.  Set potential=true if function is a
     //! potential field.  Density is assumed if false (default).
@@ -632,7 +752,8 @@ namespace BasisClasses
     void make_coefs(void);
     
     //! Accumulate new coefficients
-    virtual void accumulate(double x, double y, double z, double mass);
+    virtual void accumulate(double x, double y, double z, double mass,
+			    unsigned long indx = 0);
     
     //! Return current maximum harmonic order in expansion
     int getMmax() { return mmax; }
@@ -789,7 +910,8 @@ namespace BasisClasses
     void make_coefs(void);
     
     //! Accumulate new coefficients
-    virtual void accumulate(double x, double y, double z, double mass);
+    virtual void accumulate(double x, double y, double z, double mass,
+			    unsigned long indx = 0);
     
     //! Return current maximum harmonic order in expansion
     int getMmax() { return mmax; }
@@ -839,8 +961,7 @@ namespace BasisClasses
     int lmaxfid, nmaxfid, mmax, mlim, nmax;
     int ncylodd, ncylnx, ncylny, ncylr, cmap, cmapR, cmapZ, vflag;
     int rnum, pnum, tnum;
-    double rmin, rmax, rcylmin, rcylmax;
-    double acyl, hcyl;
+    double rcylmin, rcylmax, acyl, hcyl;
     bool expcond, logarithmic, density, EVEN_M, sech2 = false;
     
     std::vector<Eigen::MatrixXd> potd, dpot, dpt2, dend;
@@ -849,8 +970,7 @@ namespace BasisClasses
     bool coefs_defined;
     Eigen::MatrixXd factorial;
     Eigen::MatrixXd expcoef;
-    double rscl;
-    int used;
+    int used, sampT = 1;
     
     std::string cachename;
     bool oldcache = false;
@@ -916,6 +1036,9 @@ namespace BasisClasses
     void computeAccel(double x, double y, double z,
 		      Eigen::Ref<Eigen::Vector3d> vstore);
 
+    //! Needed for covariance writing
+    void writeCovarH5Params(HighFive::File& file);
+
   public:
     
     //! Constructor from YAML node
@@ -934,9 +1057,6 @@ namespace BasisClasses
       return EmpCylSL::cacheInfo(cachefile);
     }
 
-    //! Prescaling factor
-    void set_scale(const double scl) { rscl = scl; }
-    
     //! Zero out coefficients to prepare for a new expansion
     void reset_coefs(void);
     
@@ -944,7 +1064,8 @@ namespace BasisClasses
     void make_coefs(void);
     
     //! Accumulate new coefficients
-    virtual void accumulate(double x, double y, double z, double mass);
+    virtual void accumulate(double x, double y, double z, double mass,
+			    unsigned long indx = 0);
     
     //! Return current maximum harmonic order in expansion
     int getMmax() { return mmax; }
@@ -975,6 +1096,31 @@ namespace BasisClasses
       return ret;
     }
 
+    //! Return a vector of tuples of basis functions and the
+    //! covariance matrix for subsamples of particles
+    std::vector<std::vector<CoefCovarType>>
+    getCoefCovariance()
+    {
+      return sl->getCoefCovariance();
+    }
+
+    std::tuple<Eigen::VectorXi, Eigen::VectorXd>
+    getCovarSamples()
+    {
+      // Copy to base class members
+      std::tie(sampleCounts, sampleMasses) = sl->getCovarSamples();
+      return {sampleCounts, sampleMasses};
+    }
+
+    virtual void enableCoefCovariance(bool pcavar_in, int sampT_in=100)
+    {
+      pcavar = pcavar_in;
+      sampT  = sampT_in;
+      if (pcavar)  {
+	sl->setSampT(std::max(1, sampT));
+	sl->set_covar(true);
+      }
+    }
     //! Create and the coefficients from a function callback with the
     //! provided time value.  Set potential=true if function is a
     //! potential field.  Density is assumed if false (default).
@@ -1040,7 +1186,7 @@ namespace BasisClasses
 
     coefType expcoef;
     
-    //! Notal mass on grid
+    //! Total mass on grid
     double totalMass;
 
     //! Number of particles
@@ -1103,7 +1249,8 @@ namespace BasisClasses
     void make_coefs(void);
     
     //! Accumulate new coefficients
-    virtual void accumulate(double x, double y, double z, double mass);
+    virtual void accumulate(double x, double y, double z, double mass,
+			    unsigned long indx = 0);
     
     //! Return current maximum harmonic order in expansion
     Eigen::Vector3i getNmax() { return {nmaxx, nmaxy, nmaxz}; }
@@ -1171,6 +1318,9 @@ namespace BasisClasses
     //! Maximum expansion order for construction
     int nmaxx, nmaxy, nmaxz;
 
+    //! Total number of modes
+    unsigned Itot;
+
     //! Number of integration knots for orthogonal check
     int knots;
 
@@ -1181,7 +1331,7 @@ namespace BasisClasses
 
     coefType expcoef;
     
-    //! Notal mass on grid
+    //! Total mass on grid
     double totalMass;
 
     //! Number of particles
@@ -1222,6 +1372,17 @@ namespace BasisClasses
     void computeAccel(double x, double y, double z, 
 		      Eigen::Ref<Eigen::Vector3d> vstore);
 
+    //@{
+    //! Covariance structures.  First index is T, second is the
+    //! flattened 3-d k vector
+    std::vector<Eigen::VectorXcd> meanV;
+    std::vector<Eigen::MatrixXcd> covrV;
+    void init_covariance();
+    void zero_covariance();
+    void writeCovarH5Params(HighFive::File& file);
+    bool covar = false;
+    int  sampT = 100;
+    //@}
 
   public:
     
@@ -1241,7 +1402,8 @@ namespace BasisClasses
     void make_coefs(void);
     
     //! Accumulate new coefficients
-    virtual void accumulate(double x, double y, double z, double mass);
+    virtual void accumulate(double x, double y, double z, double mass,
+			    unsigned long indx = 0);
     
     //! Return current maximum harmonic order in expansion
     Eigen::Vector3i getNmax() { return {nmaxx, nmaxy, nmaxz}; }
@@ -1270,8 +1432,101 @@ namespace BasisClasses
       return true;
     }
 
+    //! Return a vector of tuples of basis functions and the
+    //! covariance matrix for subsamples of particles
+    std::vector<std::vector<CoefCovarType>> getCoefCovariance();
+
+    //! Get sample count and mass arrays
+    std::tuple<Eigen::VectorXi, Eigen::VectorXd>
+    getCovarSamples()
+    {
+      return {sampleCounts, sampleMasses};
+    }
+
+    //! Enable coefficient covariance computation
+    virtual void enableCoefCovariance(bool pcavar_in, int sampT_in=100)
+    {
+      pcavar = pcavar_in;
+      sampT  = sampT_in;
+      if (pcavar) {
+	init_covariance();
+	zero_covariance();
+      }
+    }
+
+    //! Turn on/off covariance computation
+    void setCovar(bool flag=true)  { covar = flag; }
+
+    //! Return wave-number indices from flattened index
+    std::tuple<int, int, int> index3D(unsigned k);
+
+    //! Get flattened index from wave-number indicies
+    unsigned index1D(int kx, int ky, int kz);
   };
   
+  class CovarianceReader
+  {
+  public:
+    
+    using CoefCovarType = std::tuple<Eigen::VectorXcd, Eigen::MatrixXcd>;
+
+  protected:
+
+    std::vector<double> times;
+    std::vector<Eigen::VectorXi> sampleCounts;
+    std::vector<Eigen::VectorXd> sampleMasses;
+    std::vector<std::vector<std::vector<CoefCovarType>>> covarData;
+
+    const double fixedPointPrecision = 1.0e08;
+    std::map<int, unsigned> timeMap;
+    unsigned timeIndex(double time)
+    {
+      int itime = static_cast<int>(time * fixedPointPrecision + 0.5);
+      auto it = timeMap.find(itime);
+      if (it == timeMap.end())
+	throw std::runtime_error("CovarianceReader: time not found");
+      return it->second;
+    }
+
+    std::string basisID;
+
+  public:
+
+    //! Constructor from HDF5 file
+    CovarianceReader(const std::string& filename, int stride=1);
+
+    //! Get time list
+    std::vector<double> Times() { return times; }
+
+    //@{
+    //! Return a vector of tuples of basis functions and the
+    //! covariance matrix for subsamples of particles
+    std::vector<std::vector<CoefCovarType>>
+    getCoefCovariance(unsigned index) { return covarData[index]; }
+
+    std::vector<std::vector<CoefCovarType>>
+    getCoefCovariance(double time) { return covarData[timeIndex(time)]; }
+    //@}
+
+
+    //@{
+    //! Get sample counts for the covariance computation
+    std::tuple<Eigen::VectorXi, Eigen::VectorXd>
+    getCovarSamples(unsigned index)
+    {
+      return {sampleCounts[index], sampleMasses[index]};
+    }
+    
+    std::tuple<Eigen::VectorXi, Eigen::VectorXd>
+    getCovarSamples(double time)
+    {
+      unsigned index = timeIndex(time);
+      return {sampleCounts[index], sampleMasses[index]};
+    }
+    //@}
+    
+    std::string basisIDname() { return basisID; }
+  };
 
   //! Time-dependent potential-density model
   using BasisCoef = std::tuple<std::shared_ptr<Basis>, std::shared_ptr<CoefClasses::Coefs>>;
diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc
index 60f6faae8..9a1469078 100644
--- a/expui/BiorthBasis.cc
+++ b/expui/BiorthBasis.cc
@@ -262,6 +262,10 @@ namespace BasisClasses
       if (conf["EVEN_L"])    EVEN_L    = conf["EVEN_L"].as<bool>();
       if (conf["EVEN_M"])    EVEN_M    = conf["EVEN_M"].as<bool>();
       if (conf["M0_ONLY"])   M0_only   = conf["M0_ONLY"].as<bool>();
+      if (conf["pcavar"])    pcavar    = conf["pcavar"].as<bool>();
+      if (conf["subsamp"])   sampT     = conf["subsamp"].as<int>();
+
+      sampT = std::max(1, sampT); // Sanity
     } 
     catch (YAML::Exception & error) {
       if (myid==0) std::cout << "Error parsing parameter stanza for <"
@@ -291,15 +295,23 @@ namespace BasisClasses
     for (auto & v : dpt2) v.resize(lmax+1, nmax);
     for (auto & v : dend) v.resize(lmax+1, nmax);
 
+    // Wasteful but simple factorial table.  Could be done with a
+    // triangular indexing scheme...
+    //
     for (auto & v : legs  ) v.resize(lmax+1, lmax+1);
     for (auto & v : dlegs ) v.resize(lmax+1, lmax+1);
     for (auto & v : d2legs) v.resize(lmax+1, lmax+1);
 
+    // Allocate coefficient storage; stores real and imaginary parts as reals
+    //
     expcoef.resize((lmax+1)*(lmax+1), nmax);
     expcoef.setZero();
       
     work.resize(nmax);
       
+    // Wasteful but simple factorial table.  Could be done with a
+    // triangular indexing scheme...
+    //
     factorial.resize(lmax+1, lmax+1);
       
     for (int l=0; l<=lmax; l++) {
@@ -312,13 +324,57 @@ namespace BasisClasses
 
     used = 0;
 
-    // Set spherical coordindates
+    // Initialize covariance
+    //
+    if (pcavar) init_covariance();
+
+    // Set spherical coordinates
     //
     coordinates = Coord::Spherical;
+    BasisID = "Spherical";
   }
   
+  void Spherical::init_covariance()
+  {
+    if (pcavar) {
+      // Triangular for l,m>=0
+      int Ltot = (lmax+1)*(lmax+2)/2;
+	
+      meanV.resize(sampT);
+      for (auto& v : meanV) {
+	v.resize(Ltot);
+	for (auto& vec : v) vec.resize(nmax);
+      }
+
+      covrV.resize(sampT);
+      for (auto& v : covrV) {
+	v.resize(Ltot);
+	for (auto& mat : v) mat.resize(nmax, nmax);
+      }
+
+      sampleCounts.resize(sampT);
+      sampleMasses.resize(sampT);
+      
+      zero_covariance();
+    }
+  }
+
+  void Spherical::zero_covariance()
+  {
+    for (int T=0; T<sampT; T++) {
+      for (auto& v : meanV[T]) v.setZero();
+      for (auto& v : covrV[T]) v.setZero();
+    }
+
+    sampleCounts.setZero();
+    sampleMasses.setZero();
+  }
+
   void SphericalSL::initialize()
   {
+    // Identifier
+    //
+    BasisID = "SphereSL";
 
     // Assign some defaults
     //
@@ -380,6 +436,10 @@ namespace BasisClasses
   
   void Bessel::initialize()
   {
+    // Identifier
+    //
+    BasisID = "Bessel";
+
     try {
       if (conf["rnum"])
 	rnum = conf["rnum"].as<int>();
@@ -409,6 +469,7 @@ namespace BasisClasses
     if (expcoef.rows()>0 && expcoef.cols()>0) expcoef.setZero();
     totalMass = 0.0;
     used = 0;
+    if (pcavar) zero_covariance();
   }
   
   
@@ -422,11 +483,11 @@ namespace BasisClasses
     cf->time   = time;
     cf->normed = true;
 
-    // Angular storage dimension
+    // Angular storage dimension; triangular number for complex coefs
     int ldim = (lmax+1)*(lmax+2)/2;
 
     // Allocate the coefficient storage
-    cf->store.resize((lmax+1)*(lmax+2)/2*nmax);
+    cf->store.resize(ldim*nmax);
 
     // Make the coefficient map
     cf->coefs = std::make_shared<CoefClasses::SphStruct::coefType>
@@ -507,7 +568,8 @@ namespace BasisClasses
       coefctr = {0.0, 0.0, 0.0};
   }
 
-  void Spherical::accumulate(double x, double y, double z, double mass)
+  void Spherical::accumulate(double x, double y, double z, double mass,
+			     unsigned long indx)
   {
     double fac, fac1, fac2, fac4;
     double norm = -4.0*M_PI;
@@ -528,6 +590,7 @@ namespace BasisClasses
     
     if (r < rmin or r > rmax) return;
     
+    // Update counters
     used++;
     totalMass += mass;
     
@@ -535,17 +598,25 @@ namespace BasisClasses
     
     legendre_R(lmax, costh, legs[tid]);
     
+    // Sample index for pcavar
+    int T = 0;
+    if (pcavar) {
+      T = used % sampT;
+      sampleCounts(T) += 1;
+      sampleMasses(T) += mass;
+    }
+    
     // L loop
-    for (int l=0, loffset=0; l<=lmax; loffset+=(2*l+1), l++) {
+    for (int l=0, loffset=0, L=0; l<=lmax; loffset+=(2*l+1), l++) {
       
       Eigen::VectorXd workE;
-      int esize = (l+1)*nmax;
       
       // M loop
-      for (int m=0, moffset=0, moffE=0; m<=l; m++) {
+      for (int m=0, moffset=0; m<=l; m++, L++) {
 	
+	fac = factorial(l, m) * legs[tid](l, m);
+
 	if (m==0) {
-	  fac = factorial(l, m) * legs[tid](l, m);
 	  for (int n=0; n<nmax; n++) {
 	    fac4 = potd[tid](l, n)*fac;
 	    expcoef(loffset+moffset, n) += fac4 * norm * mass;
@@ -554,9 +625,9 @@ namespace BasisClasses
 	  moffset++;
 	}
 	else {
-	  fac  = factorial(l, m) * legs[tid](l, m);
 	  fac1 = fac*cos(phi*m);
 	  fac2 = fac*sin(phi*m);
+
 	  for (int n=0; n<nmax; n++) {
 	    fac4 = potd[tid](l, n);
 	    expcoef(loffset+moffset  , n) += fac1 * fac4 * norm * mass;
@@ -565,24 +636,62 @@ namespace BasisClasses
 	  
 	  moffset+=2;
 	}
+
+	if (pcavar) {
+	  Eigen::VectorXcd g = std::exp(std::complex<double>(0.0, m*phi)) * 
+	    potd[tid].row(l).transpose() * fac * norm;
+	  
+	  meanV[T][L].noalias() += g * mass;
+	  covrV[T][L].noalias() += g * g.adjoint() * mass;
+	}
+
       }
+      // END: m loop
     }
-    
+    // END: l loop
   }
   
   void Spherical::make_coefs()
   {
+    // MPI reduction of coefficients
     if (use_mpi) {
       
+      // Square of total number of angular coefficients in real form
+      int Ltot = (lmax+1)*(lmax+1);
+
       MPI_Allreduce(MPI_IN_PLACE, &used, 1, MPI_INT,
 		    MPI_SUM, MPI_COMM_WORLD);
       
-      for (int l=0; l<(lmax+1)*(lmax+1); l++) {
+      for (int l=0; l<Ltot; l++) {
 	work = expcoef.row(l);
 	MPI_Allreduce(MPI_IN_PLACE, work.data(), nmax, MPI_DOUBLE,
 		      MPI_SUM, MPI_COMM_WORLD);
 	expcoef.row(l) = work;
       }
+      
+      if (pcavar) {
+
+	// Triangular number of angular coefficients in complex form
+	int ltot = (lmax+1)*(lmax+2)/2;
+
+	MPI_Allreduce(MPI_IN_PLACE, sampleCounts.data(), sampleCounts.size(),
+		      MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+	
+	MPI_Allreduce(MPI_IN_PLACE, sampleMasses.data(), sampleMasses.size(),
+		      MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+
+	for (int T=0; T<sampT; T++) {
+
+	  for (int l=0; l<ltot; l++) {
+	    
+	    MPI_Allreduce(MPI_IN_PLACE, meanV[T][l].data(), meanV[T][l].size(),
+			  MPI_DOUBLE_COMPLEX, MPI_SUM, MPI_COMM_WORLD);
+
+	    MPI_Allreduce(MPI_IN_PLACE, covrV[T][l].data(), covrV[T][l].size(),
+			  MPI_DOUBLE_COMPLEX, MPI_SUM, MPI_COMM_WORLD);
+	  }
+	}
+      }
     }
   }
   
@@ -902,6 +1011,28 @@ namespace BasisClasses
     return ret;
   }
 
+  
+  /** Return a vector of tuples of basis functions and the covariance
+      matrix for subsamples of particles */
+  std::vector<std::vector<BiorthBasis::CoefCovarType>>
+  Spherical::getCoefCovariance()
+  {
+    std::vector<std::vector<BiorthBasis::CoefCovarType>> ret;
+   if (pcavar) {
+     ret.resize(sampT);
+     for (int T=0; T<sampT; T++) {
+       int ltot = (lmax+1)*(lmax+2)/2;
+       ret[T].resize(ltot);
+       for (int l=0; l<ltot; l++) {
+	 std::get<0>(ret[T][l]) = meanV[T][l];
+	 std::get<1>(ret[T][l]) = covrV[T][l];
+       }
+     }
+   }
+    
+    return ret;
+  }
+
 
 #define MINEPS 1.0e-10
   
@@ -1238,6 +1369,10 @@ namespace BasisClasses
 
   void Cylindrical::initialize()
   {
+    // Basis identifier
+    //
+    BasisID = "Cylindrical";
+
     // Assign some defaults
     //
     rcylmin     = 0.001;
@@ -1350,9 +1485,9 @@ namespace BasisClasses
 
       if (conf["aratio"    ])      aratio = conf["aratio"    ].as<double>();
       if (conf["hratio"    ])      hratio = conf["hratio"    ].as<double>();
-      if (conf["dweight"   ])      dweight = conf["dweight"   ].as<double>();
-      if (conf["Mfac"   ])         Mfac = conf["Mfac"   ].as<double>();
-      if (conf["HERNA"   ])        HERNA = conf["HERNA"   ].as<double>();
+      if (conf["dweight"   ])      dweight = conf["dweight"  ].as<double>();
+      if (conf["Mfac"      ])      Mfac   = conf["Mfac"      ].as<double>();
+      if (conf["HERNA"     ])      HERNA  = conf["HERNA"     ].as<double>();
       if (conf["rwidth"    ])      rwidth = conf["rwidth"    ].as<double>();
       if (conf["ashift"    ])      ashift = conf["ashift"    ].as<double>();
       if (conf["rfactor"   ])     rfactor = conf["rfactor"   ].as<double>();
@@ -1362,6 +1497,11 @@ namespace BasisClasses
       if (conf["dtype"     ])       dtype = conf["dtype"     ].as<std::string>();
       if (conf["vflag"     ])       vflag = conf["vflag"     ].as<int>();
       if (conf["pyname"    ])      pyname = conf["pyname"    ].as<std::string>();
+      if (conf["pcavar"]    )      pcavar = conf["pcavar"    ].as<bool>();
+      if (conf["subsamp"]   )      sampT  = conf["subsamp"   ].as<int>();
+
+      // Sanity
+      sampT = std::max(1, sampT);
 
       // Deprecation warning
       if (conf["density"   ]) {
@@ -1411,7 +1551,7 @@ namespace BasisClasses
     EmpCylSL::CMAPZ       = cmapZ;
     EmpCylSL::logarithmic = logarithmic;
     EmpCylSL::VFLAG       = vflag;
-    
+
     // Check for non-null cache file name.  This must be specified
     // to prevent recomputation and unexpected behavior.
     //
@@ -1433,7 +1573,11 @@ namespace BasisClasses
     //
     if (mlim>=0)  sl->set_mlim(mlim);
     if (EVEN_M)   sl->setEven(EVEN_M);
-      
+    if (pcavar)  {
+      sl->setSampT(sampT);
+      sl->set_covar(true);
+    }
+    
     // Cache override for old Eigen cache
     //
     if (oldcache) sl->AllowOldCache();
@@ -1662,11 +1806,12 @@ namespace BasisClasses
     }
   }
   
-  void Cylindrical::accumulate(double x, double y, double z, double mass)
+  void Cylindrical::accumulate(double x, double y, double z, double mass,
+			       unsigned long indx)
   {
     double R   = sqrt(x*x + y*y);
     double phi = atan2(y, x);
-    sl->accumulate(R, z, phi, mass, 0, 0);
+    sl->accumulate(R, z, phi, mass, indx, 0, 0, pcavar);
   }
   
   void Cylindrical::reset_coefs(void)
@@ -1729,7 +1874,7 @@ namespace BasisClasses
 
   void Cylindrical::make_coefs(void)
   {
-    sl->make_coefficients();
+    sl->make_coefficients(pcavar);
   }
   
   
@@ -1834,6 +1979,9 @@ namespace BasisClasses
 
   void FlatDisk::initialize()
   {
+    // Basis identifier
+    //
+    BasisID = "FlatDisk";
 
     // Assign some defaults
     //
@@ -2039,7 +2187,8 @@ namespace BasisClasses
       coefctr = {0.0, 0.0, 0.0};
   }
 
-  void FlatDisk::accumulate(double x, double y, double z, double mass)
+  void FlatDisk::accumulate(double x, double y, double z, double mass,
+			    unsigned long indx)
   {
     // Normalization factors
     //
@@ -2058,6 +2207,7 @@ namespace BasisClasses
 
     if (R < ortho->getRtable() and fabs(z) < ortho->getRtable()) {
     
+      // Update counters
       used++;
       totalMass += mass;
     
@@ -2413,6 +2563,9 @@ namespace BasisClasses
 
   void CBDisk::initialize()
   {
+    // Basis identifier
+    //
+    BasisID = "CBDisk";
 
     // Assign some defaults
     //
@@ -2840,7 +2993,8 @@ namespace BasisClasses
       coefctr = {0.0, 0.0, 0.0};
   }
 
-  void CBDisk::accumulate(double x, double y, double z, double mass)
+  void CBDisk::accumulate(double x, double y, double z, double mass,
+			  unsigned long int idx)
   {
     // Normalization factors
     //
@@ -2857,6 +3011,7 @@ namespace BasisClasses
     // Get thread id
     int tid = omp_get_thread_num();
 
+    // Update counters
     used++;
     totalMass += mass;
     
@@ -3148,6 +3303,10 @@ namespace BasisClasses
 
   void Slab::initialize()
   {
+    // Basis identifier
+    //
+    BasisID = "Slab";
+
     nminx = 0;
     nminy = 0;
 
@@ -3294,7 +3453,8 @@ namespace BasisClasses
     coefctr = {0.0, 0.0, 0.0};
   }
 
-  void Slab::accumulate(double x, double y, double z, double mass)
+  void Slab::accumulate(double x, double y, double z, double mass,
+			unsigned long int idx)
   {
     // Truncate to slab with sides in [0,1]
     if (x<0.0)
@@ -3307,6 +3467,7 @@ namespace BasisClasses
     else
       y -= std::floor( y);
     
+    // Update counters
     used++;
 
     // Storage for basis evaluation
@@ -3662,24 +3823,38 @@ namespace BasisClasses
     "knots",
     "verbose",
     "check",
-    "method"
+    "method",
+    "pcavar,"
+    "subsamp"
   };
 
   Cube::Cube(const YAML::Node& CONF) : BiorthBasis(CONF, "cube")
   {
     initialize();
+
+    // Initialize covariance
+    //
+    if (pcavar) init_covariance();
   }
 
   Cube::Cube(const std::string& confstr) : BiorthBasis(confstr, "cube")
   {
     initialize();
+
+    // Initialize covariance
+    //
+    if (pcavar) init_covariance();
   }
 
   void Cube::initialize()
   {
-    nminx = std::numeric_limits<int>::max();
-    nminy = std::numeric_limits<int>::max();
-    nminz = std::numeric_limits<int>::max();
+    // Basis identifier
+    //
+    BasisID = "Cube";
+
+    nminx = 0;
+    nminy = 0;
+    nminz = 0;
 
     nmaxx = 6;
     nmaxy = 6;
@@ -3705,17 +3880,20 @@ namespace BasisClasses
     // Assign values from YAML
     //
     try {
-      if (conf["nminx"])      nminx = conf["nminx"].as<int>();
-      if (conf["nminy"])      nminy = conf["nminy"].as<int>();
-      if (conf["nminz"])      nminz = conf["nminz"].as<int>();
+      if (conf["nminx"])      nminx  = conf["nminx"  ].as<int>();
+      if (conf["nminy"])      nminy  = conf["nminy"  ].as<int>();
+      if (conf["nminz"])      nminz  = conf["nminz"  ].as<int>();
       
-      if (conf["nmaxx"])      nmaxx = conf["nmaxx"].as<int>();
-      if (conf["nmaxy"])      nmaxy = conf["nmaxy"].as<int>();
-      if (conf["nmaxz"])      nmaxz = conf["nmaxz"].as<int>();
+      if (conf["nmaxx"])      nmaxx  = conf["nmaxx"  ].as<int>();
+      if (conf["nmaxy"])      nmaxy  = conf["nmaxy"  ].as<int>();
+      if (conf["nmaxz"])      nmaxz  = conf["nmaxz"  ].as<int>();
       
-      if (conf["knots"])      knots = conf["knots"].as<int>();
+      if (conf["knots"])      knots  = conf["knots"  ].as<int>();
 
-      if (conf["check"])      check = conf["check"].as<bool>();
+      if (conf["check"])      check  = conf["check"  ].as<bool>();
+
+      if (conf["pcavar"])     pcavar = conf["pcavar" ].as<bool>();
+      if (conf["subsamp"])    sampT  = conf["subsamp"].as<int>();
     } 
     catch (YAML::Exception & error) {
       if (myid==0) std::cout << "Error parsing parameter stanza for <"
@@ -3740,10 +3918,17 @@ namespace BasisClasses
     //
     int nthrds = omp_get_max_threads();
 
+    // Total number of wavenumbers
+    //
+    Itot = (2*nmaxx + 1) * (2*nmaxy + 1) * (2*nmaxz + 1);
+
     expcoef.resize(2*nmaxx+1, 2*nmaxy+1, 2*nmaxz+1);
     expcoef.setZero();
       
+    // Counters
+    //
     used = 0;
+    totalMass = 0.0;
 
     // Set cartesian coordindates
     //
@@ -3755,6 +3940,7 @@ namespace BasisClasses
     expcoef.setZero();
     totalMass = 0.0;
     used = 0;
+    if (pcavar) zero_covariance();
   }
   
   
@@ -3807,7 +3993,8 @@ namespace BasisClasses
     coefctr = {0.0, 0.0, 0.0};
   }
 
-  void Cube::accumulate(double x, double y, double z, double mass)
+  void Cube::accumulate(double x, double y, double z, double mass,
+			unsigned long int indx)
   {
     // Truncate to cube with sides in [0,1]
     if (x<0.0)
@@ -3825,6 +4012,9 @@ namespace BasisClasses
     else
       z -= std::floor( z);
     
+    // Update counters
+    used++;
+    totalMass += mass;
     
     // Recursion multipliers
     Eigen::Vector3cd step
@@ -3836,6 +4026,12 @@ namespace BasisClasses
        std::exp(-kfac*(y*nmaxy)),
        std::exp(-kfac*(z*nmaxz))};
     
+    Eigen::VectorXcd g;
+    if (pcavar) {
+      g.resize(Itot);
+      g.setZero();
+    }
+
     Eigen::Vector3cd curr(init);
     for (int ix=0; ix<=2*nmaxx; ix++, curr(0)*=step(0)) {
       curr(1) = init(1);
@@ -3850,13 +4046,30 @@ namespace BasisClasses
 	  int jj = iy-nmaxy;
 	  int kk = iz-nmaxz;
 
+	  // Limit to minimum wave number
+	  if (abs(ii)<nminx || abs(jj)<nminy || abs(kk)<nminz) continue;
+
 	  // Normalization
 	  double norm = 1.0/sqrt(M_PI*(ii*ii + jj*jj + kk*kk));;
 
 	  expcoef(ix, iy, iz) += - mass * curr(0)*curr(1)*curr(2) * norm;
+
+	  if (pcavar)
+	    g[index1D(ix, iy, iz)] = - curr(0)*curr(1)*curr(2) * norm;
 	}
       }
     }
+
+    if (pcavar) {
+      // Sample index for pcavar
+      int T = 0;
+      T = used % sampT;
+      sampleCounts(T) += 1;
+      sampleMasses(T) += mass;
+
+      meanV[T].noalias() += g * mass;
+      if (covar) covrV[T].noalias() += g * g.adjoint() * mass;
+    }
   }
   
   void Cube::make_coefs()
@@ -3866,11 +4079,56 @@ namespace BasisClasses
       MPI_Allreduce(MPI_IN_PLACE, &used, 1, MPI_INT,
 		    MPI_SUM, MPI_COMM_WORLD);
       
+      MPI_Allreduce(MPI_IN_PLACE, &totalMass, 1, MPI_DOUBLE,
+		    MPI_SUM, MPI_COMM_WORLD);
+      
       MPI_Allreduce(MPI_IN_PLACE, expcoef.data(), expcoef.size(), MPI_DOUBLE_COMPLEX,
 		    MPI_SUM, MPI_COMM_WORLD);
+
+
+      if (pcavar) {
+
+	MPI_Allreduce(MPI_IN_PLACE, sampleCounts.data(), sampleCounts.size(),
+		      MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+	
+	MPI_Allreduce(MPI_IN_PLACE, sampleMasses.data(), sampleMasses.size(),
+		      MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+
+	for (int T=0; T<sampT; T++) {
+
+	  MPI_Allreduce(MPI_IN_PLACE, meanV[T].data(), meanV[T].size(),
+			MPI_DOUBLE_COMPLEX, MPI_SUM, MPI_COMM_WORLD);
+	  
+	  if (covar)
+	    MPI_Allreduce(MPI_IN_PLACE, covrV[T].data(), covrV[T].size(),
+			  MPI_DOUBLE_COMPLEX, MPI_SUM, MPI_COMM_WORLD);
+	}
+	// END: sample loop
+      }
+      // END: pcavar
     }
+    // END: using mpi
   }
   
+  /** Return a vector of tuples of basis functions and the covariance
+      matrix for subsamples of particles for Cube type */
+  std::vector<std::vector<BiorthBasis::CoefCovarType>>
+  Cube::getCoefCovariance()
+  {
+    std::vector<std::vector<BiorthBasis::CoefCovarType>> ret;
+    if (pcavar) {
+      ret.resize(sampT);
+      for (int T=0; T<sampT; T++) {
+	ret[T].resize(1);
+	std::get<0>(ret[T][0]) = meanV[T];
+	if (covar) std::get<1>(ret[T][0]) = covrV[T];
+      }
+    }
+    
+    return ret;
+  }
+
+
   std::vector<double> Cube::crt_eval(double x, double y, double z)
   {
     // Get thread id
@@ -3977,6 +4235,91 @@ namespace BasisClasses
     return ret;
   }
   
+  void Cube::init_covariance()
+  {
+    if (pcavar) {
+
+      meanV.resize(sampT);
+      for (auto& v : meanV) {
+	v.resize(Itot);
+      }
+
+      if (covar) {
+	covrV.resize(sampT);
+	for (auto& v : covrV) {
+	  v.resize(Itot, Itot);
+	}
+      } else covrV.clear();
+
+      sampleCounts.resize(sampT);
+      sampleMasses.resize(sampT);
+      
+      zero_covariance();
+    }
+  }
+
+
+  void Cube::zero_covariance()
+  {
+    for (int T=0; T<sampT; T++) {
+      meanV[T].setZero();
+      if (covar) covrV[T].setZero();
+    }
+
+    sampleCounts.setZero();
+    sampleMasses.setZero();
+  }
+
+
+  unsigned Cube::index1D(int kx, int ky, int kz)
+  {
+    if (kx < 0 or kx > 2*nmaxx) {
+      std::ostringstream sout;
+      sout << "Cube::index1D: x index [" << kx << "] must be in [0, "
+	   << 2*nmaxx << "]";
+      throw std::runtime_error(sout.str());
+    }
+
+    if (ky < 0 or ky > 2*nmaxy) {
+      std::ostringstream sout;
+      sout << "Cube::index1D: y index [" << ky << "] must be in [0, "
+	   << 2*nmaxy << "]";
+      throw std::runtime_error(sout.str());
+    }
+
+    if (kz < 0 or kx > 2*nmaxz) {
+      std::ostringstream sout;
+      sout << "Cube::index1D: z index [" << kz << "] must be in [0, "
+	   << 2*nmaxz << "]";
+      throw std::runtime_error(sout.str());
+    }
+
+    return
+      kx*(2*nmaxy+1)*(2*nmaxz+1) +
+      ky*(2*nmaxz+1) +
+      kz;
+  }
+
+  std::tuple<int, int, int> Cube::index3D(unsigned indx)
+  {
+    // Sanity check
+    //
+    if (indx >= Itot) {
+      std::ostringstream sout;
+      sout << "Cube::index3D: index [" << indx << "] must be in 0 <= indx < " << Itot;
+      throw std::runtime_error(sout.str());
+    }
+
+    // Compute the 3d index
+    //
+    int ix = indx/((2*nmaxy+1)*(2*nmaxz+1));
+    int iy = (indx - ix*(2*nmaxy+1)*(2*nmaxz+1))/(2*nmaxz+1);
+    int iz = indx - ix*(2*nmaxy+1)*(2*nmaxz+1)/(2*nmaxz+1) - iy*(2*nmaxz+1);
+  
+    return {ix, iy, iz};
+  }
+
+
   // Generate coeffients from a particle reader
   CoefClasses::CoefStrPtr BiorthBasis::createFromReader
   (PR::PRptr reader, Eigen::Vector3d ctr, RowMatrix3d rot)
@@ -4037,7 +4380,7 @@ namespace BasisClasses
       }
 
       if (use) {
-	accumulate(pp(0), pp(1), pp(2), p->mass);
+	accumulate(pp(0), pp(1), pp(2), p->mass, p->indx);
       }
     }
     make_coefs();
@@ -4156,7 +4499,7 @@ namespace BasisClasses
 	  coefindx++;
 	  
 	  if (use) {
-	    accumulate(pp(0), pp(1), pp(2), m(n));
+	    accumulate(pp(0), pp(1), pp(2), m(n), n);
 	  }
 	}
       }
@@ -4194,7 +4537,7 @@ namespace BasisClasses
 	  coefindx++;
 	  
 	  if (use) {
-	    accumulate(pp(0), pp(1), pp(2), m(n));
+	    accumulate(pp(0), pp(1), pp(2), m(n), n);
 	  }
 	}
       }
@@ -4653,6 +4996,556 @@ namespace BasisClasses
     return {times, ret};
   }
 
+  void Spherical::writeCovarH5Params(HighFive::File& file)
+  {
+    file.createAttribute<int>("lmax", HighFive::DataSpace::From(lmax)).write(lmax);
+    file.createAttribute<int>("nmax", HighFive::DataSpace::From(nmax)).write(nmax);
+    file.createAttribute<double>("scale", HighFive::DataSpace::From(scale)).write(scale);
+    file.createAttribute<double>("rmin", HighFive::DataSpace::From(rmin)).write(rmin);
+    file.createAttribute<double>("rmax", HighFive::DataSpace::From(rmax)).write(rmax);
+  }
+  
+  void Cylindrical::writeCovarH5Params(HighFive::File& file)
+  {
+    file.createAttribute<int>("mmax", HighFive::DataSpace::From(mmax)).write(mmax);
+    file.createAttribute<int>("nmax", HighFive::DataSpace::From(nmax)).write(nmax);
+    file.createAttribute<double>("rcylmin", HighFive::DataSpace::From(rcylmin)).write(rcylmin);
+    file.createAttribute<double>("rcylmax", HighFive::DataSpace::From(rcylmax)).write(rcylmax);
+    file.createAttribute<double>("acyl", HighFive::DataSpace::From(acyl)).write(acyl);
+    file.createAttribute<double>("hcyl", HighFive::DataSpace::From(hcyl)).write(hcyl);
+  }
+  
+  void Cube::writeCovarH5Params(HighFive::File& file)
+  {
+    file.createAttribute<int>("nminx", HighFive::DataSpace::From(nminx)).write(nminx);
+    file.createAttribute<int>("nminy", HighFive::DataSpace::From(nminy)).write(nminy);
+    file.createAttribute<int>("nminz", HighFive::DataSpace::From(nminz)).write(nminz);
+    file.createAttribute<int>("nmaxx", HighFive::DataSpace::From(nmaxx)).write(nmaxx);
+    file.createAttribute<int>("nmaxy", HighFive::DataSpace::From(nmaxy)).write(nmaxy);
+    file.createAttribute<int>("nmaxz", HighFive::DataSpace::From(nmaxz)).write(nmaxz);
+  }
+  
+  unsigned BiorthBasis::writeCovarH5(HighFive::Group& snaps, unsigned count, double time)
+  {
+    std::ostringstream stim;
+    stim << std::setw(8) << std::setfill('0') << std::right << count++;
+      
+    // Make a new group for this time
+    //
+    HighFive::Group stanza = snaps.createGroup(stim.str());
+      
+    // Add a time attribute
+    //
+    time = roundTime(time);
+    stanza.createAttribute<double>("Time", HighFive::DataSpace::From(time)).write(time);
+      
+    // Enable compression
+    //
+    auto dcpl1 = HighFive::DataSetCreateProps{}; // sample stats
+    auto dcpl2 = HighFive::DataSetCreateProps{}; // coefficients
+    auto dcpl3 = HighFive::DataSetCreateProps{}; // covariance
+
+    // Properties for sample stats
+    if (H5compress) {
+      unsigned int csz = sampleCounts.size();
+      dcpl1.add(HighFive::Chunking({csz, 1}));
+      if (H5shuffle) dcpl1.add(HighFive::Shuffle());
+      dcpl1.add(HighFive::Deflate(H5compress));
+    }
+
+    // Add the sample statistics
+    //
+    HighFive::DataSet s1data = stanza.createDataSet("sampleCounts", sampleCounts, dcpl1);
+    HighFive::DataSet s2data = stanza.createDataSet("sampleMasses", sampleMasses, dcpl1);
+
+    // Covariance data
+    //
+    auto covar = getCoefCovariance();
+
+    // Number of samples
+    //
+    unsigned sampleSize   = covar.size();
+    unsigned ltot         = covar[0].size();
+    unsigned nmax         = std::get<0>(covar[0][0]).rows();
+    unsigned diagonalSize = nmax*(nmax + 1)/2;
+
+    // Add data dimensions
+    //
+    stanza.createAttribute<unsigned>
+      ("sampleSize", HighFive::DataSpace::From(sampleSize)).write(sampleSize);
+
+    stanza.createAttribute<unsigned>
+      ("angularSize", HighFive::DataSpace::From(ltot)).write(ltot);
+
+    stanza.createAttribute<unsigned>
+      ("rankSize", HighFive::DataSpace::From(nmax)).write(nmax);
+      
+    if (H5compress) {
+      // Szip parameters
+      const int options_mask = H5_SZIP_NN_OPTION_MASK;
+      const int pixels_per_block = 8;
+
+      // Properties for coefficients
+      //
+      unsigned int csz2 = nmax * ltot * sampleSize;
+      HighFive::Chunking data_dims2{std::min<unsigned>(csz2, H5chunk), 1};
+
+      dcpl2.add(data_dims2);
+      if (H5shuffle) dcpl2.add(HighFive::Shuffle());
+      if (H5szip) {
+	dcpl2.add(HighFive::Szip(options_mask, pixels_per_block));
+      } else {
+	dcpl2.add(HighFive::Deflate(H5compress));
+      }
+
+      // Properties for  covariance
+      //
+      unsigned int csz3 = ltot * diagonalSize * sampleSize;
+      HighFive::Chunking data_dims3{std::min<unsigned>(csz3, H5chunk), 1};
+
+      dcpl3.add(data_dims3);
+      if (H5shuffle) dcpl3.add(HighFive::Shuffle());
+      if (H5szip) {
+	dcpl3.add(HighFive::Szip(options_mask, pixels_per_block));
+      } else {
+	dcpl3.add(HighFive::Deflate(H5compress));
+      }
+    }
+
+    // Check for existence of a covariance matrix (only Cube can
+    // toggle this so far)
+    //
+    int varsz = std::get<1>(covar[0][0]).size();
+
+    // Pack the coefficient data
+    //
+    if (floatType) {
+      // Create a vector of doubles for the real and imaginary parts
+      Eigen::VectorXf real_part(nmax*ltot*sampleSize);
+      Eigen::VectorXf imag_part(nmax*ltot*sampleSize);
+
+      for (size_t T=0, c=0; T<sampleCounts.size(); T++) {
+	for (size_t l=0; l<ltot; l++) {
+	  for (size_t n=0; n<nmax; n++, c++) {
+	    real_part(c) = std::real(std::get<0>(covar[T][l])(n));
+	    imag_part(c) = std::imag(std::get<0>(covar[T][l])(n));
+	  }
+	}
+      }
+
+      // Create two separate, compressed datasets
+      stanza.createDataSet("coefficients_real", real_part, dcpl2);
+      stanza.createDataSet("coefficients_imag", imag_part, dcpl2);
+      
+      // Pack the covariance data in an upper triangular format
+      //
+      if (varsz) {
+
+	real_part.resize(ltot*diagonalSize*sampleSize);
+	imag_part.resize(ltot*diagonalSize*sampleSize);
+
+	for (size_t T=0, c=0; T<sampleCounts.size(); T++) {
+	  for (size_t l=0; l<ltot; l++) {
+	    for (size_t n1=0; n1<nmax; n1++) {
+	      for (size_t n2=n1; n2<nmax; n2++, c++) {
+		real_part(c) = std::real(std::get<1>(covar[T][l])(n1, n2));
+		imag_part(c) = std::imag(std::get<1>(covar[T][l])(n1, n2));
+	      }
+	    }
+	  }
+	}
+	
+	// Create two separate, compressed datasets
+	stanza.createDataSet("covariance_real", real_part, dcpl3);
+	stanza.createDataSet("covariance_imag", imag_part, dcpl3);
+      }
+      
+    } else {
+      Eigen::VectorXd real_part(ltot*nmax*sampleSize);
+      Eigen::VectorXd imag_part(ltot*nmax*sampleSize);
+
+      for (size_t T=0, c=0; T<sampleCounts.size(); T++) {
+	for (size_t l=0; l<ltot; l++) {
+	  for (size_t n=0; n<nmax; n++, c++) {
+	    real_part(c) = std::real(std::get<0>(covar[T][l])(n));
+	    imag_part(c) = std::imag(std::get<0>(covar[T][l])(n));
+	  }
+	}
+      }
+
+      // Create two separate, compressed datasets
+      //
+      stanza.createDataSet("coefficients_real", real_part, dcpl2);
+      stanza.createDataSet("coefficients_imag", imag_part, dcpl2);
+      
+      // Pack the covariance data in an upper triangular format
+      //
+      if (varsz) {
+	real_part.resize(ltot*diagonalSize*sampleSize);
+	imag_part.resize(ltot*diagonalSize*sampleSize);
+
+	for (size_t T=0, c=0; T<sampleCounts.size(); T++) {
+	  for (size_t l=0; l<ltot; l++) {
+	    for (size_t n1=0; n1<nmax; n1++) {
+	      for (size_t n2=n1; n2<nmax; n2++, c++) {
+		real_part(c) = std::real(std::get<1>(covar[T][l])(n1, n2));
+		imag_part(c) = std::imag(std::get<1>(covar[T][l])(n1, n2));
+	      }
+	    }
+	  }
+	}
+	
+	// Create two separate, compressed datasets
+	//
+	stanza.createDataSet("covariance_real", real_part, dcpl3);
+	stanza.createDataSet("covariance_imag", imag_part, dcpl3);
+      }
+    }
+    // END: sample loop
+
+    return count;
+  }
+  
+  void BiorthBasis::writeCoefCovariance(const std::string& compname, const std::string& runtag, double time)
+  {
+    // Check that variance computation is on
+    //
+    if (not pcavar) {
+      std::cout << "BiorthBasis::writeCoefCovariance: covariance computation is disabled.  "
+		<< "Set 'pcavar: true' to enable." << std::endl;
+      return;
+    }
+
+    // Only root process writes
+    //
+    if (myid) return;
+
+    // Check that there is something to write
+    //
+    int totalCount = 0;
+    std::tie(sampleCounts, sampleMasses) = getCovarSamples();
+    totalCount += sampleCounts.sum();
+
+    if (totalCount==0) {
+      std::cout << "BiorthBasis::writeCoefCovariance: no data" << std::endl;
+      return;
+    }
+
+    // Round time
+    //
+    time = roundTime(time);
+
+    // The H5 filename
+    //
+    std::string fname = "coefcovar." + compname + "." + runtag + ".h5";
+
+    // Check if file exists?
+    //
+    try {
+      // Open the HDF5 file in read-write mode, creating if it doesn't
+      // exist
+      HighFive::File file(fname,
+			  HighFive::File::ReadWrite |
+			  HighFive::File::Create);
+
+      // Check for version string
+      std::string path = "CovarianceFileVersion"; 
+
+      // Check for valid HDF file by attribute
+      if (file.hasAttribute(path)) {
+	extendCoefCovariance(fname, time);
+	return;
+      }
+
+      // Write the Version string
+      //
+      file.createAttribute<std::string>("CovarianceFileVersion", HighFive::DataSpace::From(CovarianceFileVersion)).write(CovarianceFileVersion);
+
+      // Write the basis identifier string
+      //
+      file.createAttribute<std::string>("BasisID", HighFive::DataSpace::From(BasisID)).write(BasisID);
+      
+      // Write the data type size
+      //
+      int sz = 8; if (floatType) sz = 4;
+      file.createAttribute<int>("FloatSize", HighFive::DataSpace::From(sz)).write(sz);
+
+      // Write the specific parameters
+      //
+      writeCovarH5Params(file);
+      
+      // Group count variable
+      //
+      unsigned count = 0;
+      HighFive::DataSet dataset = file.createDataSet("count", count);
+      
+      // Create a new group for coefficient snapshots
+      //
+      HighFive::Group group = file.createGroup("snapshots");
+      
+      // Write the coefficients
+      //
+      count = writeCovarH5(group, count, time);
+      
+      // Update the count
+      //
+      dataset.write(count);
+      
+    } catch (const HighFive::Exception& err) {
+      // Handle HighFive specific errors (e.g., file not found)
+      throw std::runtime_error
+	(std::string("BiorthBasis::writeCoefCovariance HighFive Error: ") + err.what());
+    } catch (const std::exception& err) {
+      // Handle other general exceptions
+      throw std::runtime_error
+	(std::string("BiorthBasis::writeCoefCovariance Error: ") + err.what());
+    }
+  }
+  
+  void BiorthBasis::extendCoefCovariance(const std::string& fname, double time)
+  {
+    try {
+      // Open an hdf5 file
+      //
+      HighFive::File file(fname, HighFive::File::ReadWrite);
+      
+      // Get the dataset
+      HighFive::DataSet dataset = file.getDataSet("count");
+      
+      unsigned count;
+      dataset.read(count);
+      
+      HighFive::Group group = file.getGroup("snapshots");
+      
+      // Write the coefficients
+      //
+      count = writeCovarH5(group, count, time);
+      
+      // Update the count
+      //
+      dataset.write(count);
+      
+    } catch (HighFive::Exception& err) {
+      throw std::runtime_error
+	(std::string("BiorthBasis::extendCoefCovariance: HighFive error: ") + err.what());
+    }
+  }
+
+  // Read covariance data
+  CovarianceReader::CovarianceReader(const std::string& filename, int stride)
+  {
+    try {
+      // Open an existing hdf5 file for reading
+      //
+      HighFive::File file(filename, HighFive::File::ReadOnly);
+      
+      // Write the Version string
+      //
+      std::string version;
+      file.getAttribute("CovarianceFileVersion").read(version);
+      // Check for alpha version
+      if (version == std::string("1.0")) {
+	throw std::runtime_error("CovarianceReader: this is an early alpha test version. Please remake your files");
+      }
+      // Test for current version
+      if (version != std::string("1.1")) {
+	throw std::runtime_error(std::string("CovarianceReader: unsupported file version, ") + version);
+      }
+
+      // Read the basis identifier string
+      //
+      file.getAttribute("BasisID").read(basisID);
+      
+      // Get the float size
+      int sz = 8;
+      file.getAttribute("FloatSize").read(sz);
+      if (sz != 4 and sz != 8) {
+	std::ostringstream sout;
+	sout << "CovarianceReader: unsupported float size, " << sz;
+	throw std::runtime_error(sout.str());
+      }
+
+      int lmax, nmax, ltot;
+
+      // Current implemented spherical types
+      const std::set<std::string> sphereType = {"Spherical", "SphereSL", "Bessel"};
+
+      // Currently implemented cylindrical types
+      const std::set<std::string> cylinderType = {"Cylindrical"};
+
+      if (sphereType.find(basisID) != sphereType.end()) {
+	file.getAttribute("lmax").read(lmax);
+	file.getAttribute("nmax").read(nmax);
+	ltot = (lmax+1)*(lmax+2)/2;
+      } else if (cylinderType.find(basisID) != cylinderType.end()) {
+	file.getAttribute("mmax").read(lmax);
+	file.getAttribute("nmax").read(nmax);
+	ltot = lmax + 1;
+      } else if (basisID == "Cube") {
+	int nmaxx, nmaxy, nmaxz;
+	file.getAttribute("nmaxx").read(nmaxx);
+	file.getAttribute("nmaxy").read(nmaxy);
+	file.getAttribute("nmaxz").read(nmaxz);
+	ltot = (2*nmaxx + 1) * (2*nmaxy + 1) * (2*nmaxz + 1);
+      } else {
+	throw std::runtime_error(std::string("CovarianceReader: unknown or unimplemented covariance for basis type, ") + basisID);
+      }
+
+      // Group count variable
+      //
+      unsigned count = 0;
+      file.getDataSet("count").read(count);
+
+      // Open the snapshot group
+      //
+      auto snaps = file.getGroup("snapshots");
+      
+      for (unsigned n=0; n<count; n+=stride) {
+
+	std::ostringstream sout;
+	sout << std::setw(8) << std::setfill('0') << std::right << n;
+      
+	auto stanza = snaps.getGroup(sout.str());
+      
+	double Time;
+	stanza.getAttribute("Time").read(Time);
+
+	int itime = static_cast<int>(Time * fixedPointPrecision + 0.5);
+	timeMap[itime] = times.size();
+	times.push_back(Time);
+
+	// Get sample properties
+	//
+	sampleCounts.push_back(Eigen::VectorXi());
+	stanza.getDataSet("sampleCounts").read(sampleCounts.back());
+	
+	sampleMasses.push_back(Eigen::VectorXd());
+	stanza.getDataSet("sampleMasses").read(sampleMasses.back());
+
+	// Get data attributes
+	//
+	int nT, lSize, rank;
+	stanza.getAttribute("sampleSize") .read(nT);
+	stanza.getAttribute("angularSize").read(lSize);
+	stanza.getAttribute("rankSize")   .read(rank);
+
+	// Allocate sample vector for current time
+	covarData.push_back(std::vector<std::vector<CoefCovarType>>(nT));
+
+	// Storage
+	Eigen::VectorXcd data0, data1;
+
+	// Get the flattened coefficient array
+	if (sz==4) {
+	  // Get the real and imaginary parts
+	  Eigen::VectorXf data_real =
+	    stanza.getDataSet("coefficients_real").read<Eigen::VectorXf>();
+
+	  Eigen::VectorXf data_imag =
+	    stanza.getDataSet("coefficients_imag").read<Eigen::VectorXf>();
+	  
+	  // Resize the complex array and assign
+	  data0.resize(data_real.size());
+	  data0.real() = data_real.cast<double>();
+	  data0.imag() = data_imag.cast<double>();
+
+	  // Check for existence of covariance
+	  //
+	  if (stanza.exist("covariance_real")) {
+
+	    data_real =
+	      stanza.getDataSet("covariance_real").read<Eigen::VectorXf>();
+
+	    data_imag =
+	      stanza.getDataSet("covariance_imag").read<Eigen::VectorXf>();
+	  
+	    // Resize the complex array and assign
+	    data1.resize(data_real.size());
+	    data1.real() = data_real.cast<double>();
+	    data1.imag() = data_imag.cast<double>();
+	  }
+	} else {
+	  // Get the real and imaginary parts
+	  Eigen::VectorXd data_real =
+	    stanza.getDataSet("coefficients_real").read<Eigen::VectorXd>();
+
+	  Eigen::VectorXd data_imag =
+	    stanza.getDataSet("coefficients_imag").read<Eigen::VectorXd>();
+	  
+	  // Resize the complex array and assign
+	  data0.resize(data_real.size());
+	  data0.real() = data_real;
+	  data0.imag() = data_imag;
+
+	  // Check for existence of covariance
+	  //
+	  if (stanza.exist("covariance_real")) {
+
+	    // Get the real and imaginary parts
+	    data_real =
+	      stanza.getDataSet("covariance_real").read<Eigen::VectorXd>();
+	    
+	    data_imag =
+	      stanza.getDataSet("covariance_imag").read<Eigen::VectorXd>();
+	    
+	    // Resize the complex array and assign
+	    data1.resize(data_real.size());
+	    data1.real() = data_real;
+	    data1.imag() = data_imag;
+	  }
+	}
+
+	// Positions in data stanzas
+	int sCof = 0, sCov = 0;
+
+	// Loop through all indices and repack
+	for (int T=0; T<nT; T++) {
+
+	  // Data element for this time
+	  std::vector<CoefCovarType> elem(lSize);
+	  for (auto & e : elem) {
+	    // Coefficients
+	    std::get<0>(e).resize(rank);
+	    // Covariance matrix
+	    if (data1.size()) std::get<1>(e).resize(rank, rank);
+	  }
+
+	  // Pack the coefficient data
+	  int c = 0;
+	  for (size_t l=0; l<lSize; l++) {
+	    for (size_t n=0; n<rank; n++) {
+	      std::get<0>(elem[l])(n) = data0(sCof + c++);
+	    }
+	  }
+	  sCof += c;
+
+	  // Pack the covariance data
+	  c = 0;
+	  for (size_t l=0; l<lSize; l++) {
+
+	    if (data1.size()) {
+	      for (size_t n1=0; n1<rank; n1++) {
+		for (size_t n2=n1; n2<rank; n2++) {
+		  std::get<1>(elem[l])(n1, n2) = data1(sCov + c++);
+		  if (n1 != n2)
+		    std::get<1>(elem[l])(n2, n1) = std::get<1>(elem[l])(n1, n2);
+		}
+	      }
+	    }
+	  }
+	  sCov += c;
+
+	  // Add the data
+	  covarData.back()[T] = std::move(elem);
+	}
+	// END: sample loop
+      }
+      // END: snapshot loop
+    } catch (HighFive::Exception& err) {
+      std::cerr << err.what() << std::endl;
+    }
+  }
+
   CoefClasses::CoefStrPtr Spherical::makeFromFunction
   (std::function<double(double, double, double, double time)> func,
    std::map<std::string, double>& params, double time, bool potential)
diff --git a/expui/CMakeLists.txt b/expui/CMakeLists.txt
index aeabee30d..fa993f75e 100644
--- a/expui/CMakeLists.txt
+++ b/expui/CMakeLists.txt
@@ -38,7 +38,7 @@ endif()
 set(expui_SOURCES BasisFactory.cc BiorthBasis.cc FieldBasis.cc
   CoefContainer.cc CoefStruct.cc FieldGenerator.cc expMSSA.cc
   Coefficients.cc KMeans.cc Centering.cc ParticleIterator.cc
-  Koopman.cc BiorthBess.cc SvdSignChoice.cc)
+  Koopman.cc BiorthBess.cc SvdSignChoice.cc KDdensity.cc)
 add_library(expui ${expui_SOURCES})
 set_target_properties(expui PROPERTIES OUTPUT_NAME expui)
 target_include_directories(expui PUBLIC ${common_INCLUDE})
diff --git a/expui/Coefficients.cc b/expui/Coefficients.cc
index 6f91a1129..f79ad70dd 100644
--- a/expui/Coefficients.cc
+++ b/expui/Coefficients.cc
@@ -2543,6 +2543,8 @@ namespace CoefClasses
       ret = std::make_shared<SphCoefs>();
     } else if (dynamic_cast<CylStruct*>(coef.get())) {
       ret = std::make_shared<CylCoefs>();
+    } else if (dynamic_cast<CubeStruct*>(coef.get())) {
+      ret = std::make_shared<CubeCoefs>();
     } else if (dynamic_cast<TblStruct*>(coef.get())) {
       ret = std::make_shared<TableData>();
     } else if (dynamic_cast<SphFldStruct*>(coef.get())) {
diff --git a/expui/FieldBasis.H b/expui/FieldBasis.H
index 605d7f0ee..bcb29ce98 100644
--- a/expui/FieldBasis.H
+++ b/expui/FieldBasis.H
@@ -156,12 +156,13 @@ namespace BasisClasses
     (Eigen::VectorXd& m, RowMatrixXd& p, bool roundrobin=true, bool posvelrows=false);
     
     //! Accumulate new coefficients
-    virtual void accumulate(double mass,
-			    double x, double y, double z,
-			    double u, double v, double w);
+    virtual void accumulate(double x, double y, double z,
+			    double u, double v, double w,
+			    double mass, unsigned long int indx = 0);
 			    
     [[deprecated("not relevant for this class")]]
-    virtual void accumulate(double x, double y, double z, double mass)
+    virtual void accumulate(double x, double y, double z, double mass,
+			    unsigned long indx = 0)
     {
     }
 
diff --git a/expui/FieldBasis.cc b/expui/FieldBasis.cc
index 47c0ec0f2..7f193fc03 100644
--- a/expui/FieldBasis.cc
+++ b/expui/FieldBasis.cc
@@ -316,9 +316,9 @@ namespace BasisClasses
     }
   }
 
-  void FieldBasis::accumulate(double mass,
-			      double x, double y, double z,
-			      double u, double v, double w)
+  void FieldBasis::accumulate(double x, double y, double z,
+			      double u, double v, double w,
+			      double mass, unsigned long int indx)
   {
     constexpr double fac2 = 0.5*M_2_SQRTPI/M_SQRT2; // 1/sqrt(2*pi)=0.3989422804
 
@@ -580,13 +580,13 @@ namespace BasisClasses
 	use = true;
       }
 
-      if (use) accumulate(p->mass,
-			  p->pos[0]-ctr[0],
+      if (use) accumulate(p->pos[0]-ctr[0],
 			  p->pos[1]-ctr[1],
 			  p->pos[2]-ctr[2],
 			  p->vel[0],
 			  p->vel[1],
-			  p->vel[2]);
+			  p->vel[2],
+			  p->mass, p->indx);
 
     }
     make_coefs();
@@ -672,13 +672,13 @@ namespace BasisClasses
 	  }
 	  coefindx++;
 	  
-	  if (use) accumulate(m(n),
-			      p(0, n)-coefctr[0],
+	  if (use) accumulate(p(0, n)-coefctr[0],
 			      p(1, n)-coefctr[1],
 			      p(2, n)-coefctr[2],
 			      p(3, n),
 			      p(4, n),
-			      p(5, n));
+			      p(5, n),
+			      m(n), n);
 	}
       }
       
@@ -707,13 +707,12 @@ namespace BasisClasses
 	  }
 	  coefindx++;
 	  
-	  if (use) accumulate(m(n),
-			      p(n, 0)-coefctr[0],
+	  if (use) accumulate(p(n, 0)-coefctr[0],
 			      p(n, 1)-coefctr[1],
 			      p(n, 2)-coefctr[2], 
 			      p(n, 3),
 			      p(n, 4),
-			      p(n, 5));
+			      p(n, 5), m(n), n);
 	}
       }
     }
diff --git a/expui/KDdensity.H b/expui/KDdensity.H
new file mode 100644
index 000000000..7afcd9c40
--- /dev/null
+++ b/expui/KDdensity.H
@@ -0,0 +1,58 @@
+#ifndef _KDdensity_H_
+#define _KDdensity_H_
+
+#include <functional>
+#include <cstdint>
+
+#include <Eigen/Eigen>
+
+#include "ParticleReader.H"
+#include "permutation.H"
+#include "localmpi.H"
+#include "KDtree.H"
+
+namespace Utility
+{
+  class KDdensity
+  {
+  protected:
+    using point3 = KDtree::point <double, 3>;
+    using tree3  = KDtree::kdtree<double, 3>;
+
+    //! Point list
+    std::vector<point3> points;
+
+    //! Density array
+    std::vector<double> KDdens;
+
+    //! Accumulated mass of phase space
+    double KDmass = 0.0;
+
+    //! Number of points per ball
+    int Ndens = 32;
+
+    //! The KD tree instance
+    std::shared_ptr<tree3> kdtree_;
+
+    //! The point index to particle ID map
+    std::map<unsigned long, unsigned> indx;
+    
+  public:
+    using RowMatrixXd = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
+
+    //! Create KD tree from a particle reader
+    KDdensity(PR::PRptr reader, int Ndens=32);
+
+    //! Create KD tree from a arrays of masses and positions
+    KDdensity(const Eigen::VectorXd& m, const RowMatrixXd& d, int Ndens=32);
+
+    //! Get density at a given point
+    double getDensityAtPoint(double x, double y, double z);
+
+    //! Get density for particle index
+    double getDensityByIndex(unsigned long index);
+
+  };
+}
+
+#endif
diff --git a/expui/KDdensity.cc b/expui/KDdensity.cc
new file mode 100644
index 000000000..caf02bf4a
--- /dev/null
+++ b/expui/KDdensity.cc
@@ -0,0 +1,138 @@
+#include "KDdensity.H"
+
+namespace Utility
+{
+  KDdensity::KDdensity(PR::PRptr reader, int Ndens) : Ndens(Ndens)
+  {
+    // Check for MPI and initialize if needed
+    //
+    int flag;
+    MPI_Initialized(&flag);
+    bool mpi_enabled = (flag != 0);
+
+    // Get number of particles
+    //
+    int nbod = reader->CurrentNumber();
+
+    // Every node needs to make the tree (a parallel share could be
+    // implemented in KDtree.H).  Build the ID map.
+    //
+    KDmass = 0.0;
+    for (auto part=reader->firstParticle(); part!=0; part=reader->nextParticle()) {
+      KDmass += part->mass;
+      points.push_back(point3({part->pos[0], part->pos[1], part->pos[2]}, part->mass));
+      indx[part->indx] = points.size()-1;
+    }
+    
+    // Build the k-d tree
+    //
+    kdtree_ = std::make_shared<tree3>(points.begin(), points.end());
+
+
+    // Share the density computation among the nodes
+    //
+    KDdens.resize(nbod, 0.0);
+    unsigned long badVol = 0;
+    for (int k=0; k<points.size(); k++) {
+      if (k % numprocs == myid) {
+	auto ret = kdtree_->nearestN(points[k], Ndens);
+	double volume = 4.0*M_PI/3.0*std::pow(std::get<2>(ret), 3.0);
+	if (volume>0.0 and KDmass>0.0)
+	  KDdens[k] = std::get<1>(ret)/volume;
+	else badVol++;
+      }
+    }
+
+    if (mpi_enabled)
+      MPI_Allreduce(MPI_IN_PLACE, KDdens.data(), nbod,
+		    MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+
+    if (myid==0) {
+      std::cout << "---- KDdensity: finished density estimate with " << badVol
+		<< " undetermined densities, mass=" << KDmass << std::endl;
+      std::cout << "---- KDdensity: a few densities are " << KDdens.front()
+		<< ", " << KDdens[nbod/2] << ", " << KDdens.back() << std::endl;
+      std::cout << "---- KDdensity: bad density count is " << badVol
+		<< " of " << nbod << " particles." << std::endl;
+    }
+  }
+
+  KDdensity::KDdensity(const Eigen::VectorXd& m,
+		       const KDdensity::RowMatrixXd& pos,
+		       int Ndens) : Ndens(Ndens)
+  {
+    // Check for MPI and initialize if needed
+    //
+    int flag;
+    MPI_Initialized(&flag);
+    bool mpi_enabled = (flag != 0);
+
+    // Get number of particles
+    //
+    int nbod = m.size();
+
+    // Sanity check
+    //
+    if (pos.rows() != nbod or pos.cols() != 3)
+      throw std::runtime_error("KDdensity: position array has wrong dimensions");
+
+    // Every node needs to make the tree (a parallel share could be
+    // implemented in KDtree.H).  Build the ID map.
+    //
+    KDmass = 0.0;
+    for (int i=0; i<nbod; i++) {
+      KDmass += m(i);
+      points.push_back(point3({pos(i, 0), pos(i, 1), pos(i, 2)}, m(i)));
+      indx[i] = i;
+    }
+    
+    // Build the k-d tree
+    //
+    kdtree_ = std::make_shared<tree3>(points.begin(), points.end());
+
+
+    // Share the density computation among the nodes
+    //
+    KDdens.resize(nbod, 0.0);
+    unsigned long badVol = 0;
+    for (int k=0; k<points.size(); k++) {
+      if (k % numprocs == myid) {
+	auto ret = kdtree_->nearestN(points[k], Ndens);
+	double volume = 4.0*M_PI/3.0*std::pow(std::get<2>(ret), 3.0);
+	if (volume>0.0 and KDmass>0.0)
+	  KDdens[k] = std::get<1>(ret)/volume;
+	else badVol++;
+      }
+    }
+
+    if (mpi_enabled)
+      MPI_Allreduce(MPI_IN_PLACE, KDdens.data(), nbod,
+		    MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+
+    if (myid==0) {
+      std::cout << "---- KDdensity: finished density estimate with " << badVol
+		<< " undetermined densities, mass=" << KDmass << std::endl;
+      std::cout << "---- KDdensity: a few densities are " << KDdens.front()
+		<< ", " << KDdens[nbod/2] << ", " << KDdens.back() << std::endl;
+      std::cout << "---- KDdensity: bad density count is " << badVol
+		<< " of " << nbod << " particles." << std::endl;
+    }
+  }
+
+  double KDdensity::getDensityAtPoint(double x, double y, double z)
+  {
+    point3 query_point({x, y, z}, 0.0);
+    auto ret = kdtree_->nearestN(query_point, Ndens);
+    double volume = 4.0*M_PI/3.0*std::pow(std::get<2>(ret), 3.0);
+    if (volume>0.0 and KDmass>0.0)
+      return std::get<1>(ret)/volume;
+    else
+      return 0.0;
+  }
+
+
+  double KDdensity::getDensityByIndex(unsigned long index)
+  {
+    return KDdens[indx[index]];
+  }
+}
diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc
index c1955422c..b99ddbb3e 100644
--- a/exputil/EmpCylSL.cc
+++ b/exputil/EmpCylSL.cc
@@ -83,7 +83,7 @@ EmpCylSL::EmpCylSL()
 {
   NORDER     = 0;
   eof_made   = false;
-  defSampT   = 1;
+  defSampT   = 0;
   sampT      = 1;
   tk_type    = None;
   EVEN_M     = false;
@@ -204,7 +204,7 @@ EmpCylSL::EmpCylSL(int nmax, int lmax, int mmax, int nord,
   eof_made   = false;
 
   sampT        = 1;
-  defSampT     = 1;
+  defSampT     = 0;
   tk_type      = None;
 
   cylmass      = 0.0;
@@ -332,7 +332,7 @@ EmpCylSL::EmpCylSL(int mlim, std::string cachename)
   eof_made   = false;
   
   sampT        = 1;
-  defSampT     = 1;
+  defSampT     = 0;
   tk_type      = None;
   
   cylmass      = 0.0;
@@ -415,7 +415,7 @@ void EmpCylSL::reset(int numr, int lmax, int mmax, int nord,
 
   // Initialize data and storage
   sampT = 1;
-  defSampT = 1;
+  defSampT = 0;
 
   cylmass = 0.0;
   cylmass1.resize(nthrds);
@@ -1923,7 +1923,7 @@ void EmpCylSL::setup_accumulation(int mlevel)
       }
     }
 
-    if (PCAVAR and sampT>0) {
+    if (PCAVAR and sampT>1) {
       for (int nth=0; nth<nthrds; nth++) {
 	for (unsigned T=0; T<sampT; T++) {
 	  numbT1[nth][T] = 0;
@@ -1946,7 +1946,7 @@ void EmpCylSL::setup_accumulation(int mlevel)
     if (m>0) accum_sin[m].setZero();
   }
 
-  if ( (PCAVAR or PCAEOF) and mlevel==0 and sampT>0) {
+  if ( (covar or PCAVAR or PCAEOF) and mlevel==0 and sampT>0) {
 
     for (int nth=0; nth<nthrds; nth++) {
 
@@ -1954,13 +1954,27 @@ void EmpCylSL::setup_accumulation(int mlevel)
 	for (auto & v : tvar[nth]) v.setZero();
       }
 
+      if (covar or PCAVAR) {
+	for (unsigned T=0; T<sampT; T++) {
+	  numbT1[nth][T] = 0;
+	  massT1[nth][T] = 0.0;
+	}
+      }
+
+      if (covar) {
+	for (unsigned T=0; T<sampT; T++) {
+	  for (int mm=0; mm<=MMAX; mm++) {
+	    MV[nth][T][mm].setZero();
+	    VC[nth][T][mm].setZero();
+	  }
+	}
+      }
+
       if (PCAVAR) {
 	minSNR = std::numeric_limits<double>::max();
 	maxSNR = 0.0;
 
 	for (unsigned T=0; T<sampT; T++) {
-	  numbT1[nth][T] = 0;
-	  massT1[nth][T] = 0.0;
 	  for (int mm=0; mm<=MMAX; mm++) {
 	    covV[nth][T][mm].setZero();
 	    covM[nth][T][mm].setZero();
@@ -2026,7 +2040,7 @@ void EmpCylSL::init_pca()
     for (int nth=0; nth<nthrds;nth++) {
       if (PCAEOF) {
 	tvar[nth].resize(MMAX + 1);
-	for (auto & v : tvar[nth]) v.resize(rank3, rank3);
+	for (auto & v : tvar[nth]) v.resize(NORDER, NORDER);
       }
 
       if (PCAVAR) {
@@ -2042,8 +2056,8 @@ void EmpCylSL::init_pca()
 	  covV[nth][T].resize(MMAX+1);
 	  covM[nth][T].resize(MMAX+1);
 	  for (int mm=0; mm<=MMAX; mm++) {
-	    covV[nth][T][mm].resize(rank3);
-	    covM[nth][T][mm].resize(rank3, rank3);
+	    covV[nth][T][mm].resize(NORDER);
+	    covM[nth][T][mm].resize(NORDER, NORDER);
 	  }
 	}
       }
@@ -2051,6 +2065,48 @@ void EmpCylSL::init_pca()
   }
 }
 
+void EmpCylSL::init_covar()
+{
+  if (covar) {
+
+    if (defSampT) sampT = defSampT;
+    else          sampT = floor(sqrt(nbodstot));
+
+    pthread_mutex_init(&used_lock, NULL);
+
+    MV.resize(nthrds);
+    VC.resize(nthrds);
+
+    // These are shared with PCA
+    numbT1.resize(nthrds);
+    massT1.resize(nthrds);
+    numbT .resize(sampT, 0);
+    massT .resize(sampT, 0);
+
+    for (int nth=0; nth<nthrds;nth++) {
+
+      MV[nth].resize(sampT);
+      VC[nth].resize(sampT);
+
+      for (unsigned T=0; T<sampT; T++) {
+
+	MV[nth][T].resize(MMAX + 1);
+	VC[nth][T].resize(MMAX + 1);
+
+	for (unsigned mm=0; mm<=MMAX; mm++) {
+	  MV[nth][T][mm].resize(NORDER, NORDER);
+	  VC[nth][T][mm].resize(NORDER);
+	  MV[nth][T][mm].setZero();
+	  VC[nth][T][mm].setZero();
+	}
+      }
+
+      numbT1[nth].resize(sampT, 0);
+      massT1[nth].resize(sampT, 0);
+    }
+  }
+}
+
 void EmpCylSL::setup_table()
 {
   // Create storage for EOF tables
@@ -3991,7 +4047,7 @@ void EmpCylSL::accumulate(double r, double z, double phi, double mass,
   double norm = -4.0*M_PI;
   
   unsigned whch;
-  if (compute and PCAVAR) {
+  if (compute and (covar or PCAVAR)) {
     whch = seq % sampT;
     pthread_mutex_lock(&used_lock);
     numbT1[id][whch] += 1;
@@ -4045,9 +4101,17 @@ void EmpCylSL::accumulate(double r, double z, double phi, double mass,
       }
     }
 
+    if (compute and covar) {
+      Eigen::VectorXcd vec = std::complex<double>(mcos, msin) *
+	vc[id].row(mm).transpose() * norm;
+
+      VC[id][whch][mm] += mass * vec;
+      MV[id][whch][mm] += mass * (vec * vec.adjoint());
+    }
+
     cylmass1[id] += mass;
   }
-
+  // End: M loop
 }
 
 
@@ -4137,38 +4201,69 @@ void EmpCylSL::make_coefficients(unsigned M0, bool compute)
 	      tvar[0][mm](nn, oo) += tvar[nth][mm](nn, oo);
       }
 	
-      for (unsigned T=0; T<sampT; T++) {
-	numbT1[0][T] += numbT1[nth][T];
-	massT1[0][T] += massT1[nth][T];
+      if (covar or PCAVAR) {
 
-	for (int mm=0; mm<=MMAX; mm++) {
-	  for (int nn=0; nn<rank3; nn++) {
-	    covV(0, T, mm)[nn] += covV(nth, T, mm)[nn];
-	    for (int oo=0; oo<rank3; oo++) {
-	      covM(0, T, mm)(nn, oo) += covM(nth, T, mm)(nn, oo);
+	for (unsigned T=0; T<sampT; T++) {
+	  numbT1[0][T] += numbT1[nth][T];
+	  massT1[0][T] += massT1[nth][T];
+
+	  for (int mm=0; mm<=MMAX; mm++) {
+
+	    if (covar) {
+	      VC[0][T][mm] += VC[nth][T][mm];
+	      MV[0][T][mm] += MV[nth][T][mm];
+	    }
+
+	    if (PCAVAR) {
+	      for (int nn=0; nn<rank3; nn++) {
+		covV(0, T, mm)[nn] += covV(nth, T, mm)[nn];
+		for (int oo=0; oo<rank3; oo++) {
+		  covM(0, T, mm)(nn, oo) += covM(nth, T, mm)(nn, oo);
+		}
+	      }
 	    }
+
 	  }
+	  // END: M loop
 	}
+	// END: T loop
 
-      } // T loop
+      }
+      // END: PCAEOF/covar stanza
       
     } // Thread loop
 
 
     // Mass used to compute variance in each partition
     //
-    if (PCAVAR) {
+    if (covar or PCAVAR) {
       if (use_mpi) {
+	
 	MPI_Allreduce ( &numbT1[0][0], &numbT[0], sampT,
 			MPI_UNSIGNED, MPI_SUM, MPI_COMM_WORLD);
 	MPI_Allreduce ( &massT1[0][0], &massT[0], sampT,
 			MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
-      } else {
+      }
+      else {
 	numbT = numbT1[0];
 	massT = massT1[0];
       }
     }
 
+    // Complex coariance components
+    //
+    if (covar and use_mpi) {
+
+      for (unsigned T=0; T<sampT; T++) {
+	for (int mm=0; mm<=MMAX; mm++) {
+	  MPI_Allreduce ( MPI_IN_PLACE, MV[0][T][mm].data(), MV[0][T][mm].size(),
+			  MPI_DOUBLE_COMPLEX, MPI_SUM, MPI_COMM_WORLD);
+	  MPI_Allreduce ( MPI_IN_PLACE, VC[0][T][mm].data(), VC[0][T][mm].size(),
+			  MPI_DOUBLE_COMPLEX, MPI_SUM, MPI_COMM_WORLD);
+	}
+      }
+    }
+
     // Test variance
     //
     if (PCAEOF) {
@@ -4259,6 +4354,7 @@ void EmpCylSL::make_coefficients(bool compute)
     MPIin2 .resize(rank3*rank3*(MMAX+1));
     MPIout2.resize(rank3*rank3*(MMAX+1));
   }
+
 				// Sum up over threads
 				// 
   for (unsigned M=0; M<=multistep; M++) {
@@ -4269,13 +4365,24 @@ void EmpCylSL::make_coefficients(bool compute)
 
       howmany1[M][0] += howmany1[M][nth];
 
-      if (compute and PCAVAR) {
+      if (compute and (covar or PCAVAR)) {
+
 	for (unsigned T=0; T<sampT; T++) {
 	  numbT1[0][T] += numbT1[nth][T];
 	  massT1[0][T] += massT1[nth][T];
+
+	  if (covar) {
+	    for (int mm=0; mm<=MMAX; mm++) {
+	      VC[0][T][mm] += VC[nth][T][mm];
+	      MV[0][T][mm] += MV[nth][T][mm];
+	    }
+	  }
+	  // END: covar stanza
 	}
+	// END: T loop
       }
-      
+
+
       for (int mm=0; mm<=MMAX; mm++) {
 	for (int nn=0; nn<rank3; nn++) {
 	  cosN(M)[0][mm][nn] += cosN(M)[nth][mm][nn];
@@ -4296,12 +4403,13 @@ void EmpCylSL::make_coefficients(bool compute)
 	  sinN(M)[0][mm][nn] += sinN(M)[nth][mm][nn];
 	}
       }
-
     }
+    // END: thread loop
   }
-
-				// Begin distribution loop
-
+  // END: multistep loop
+  
+  // Begin distribution loop
+  //
   howmany. resize(multistep+1);	// Resize and zero, if empty
   howmany1.resize(multistep+1);
   for (unsigned M=0; M<=multistep; M++) {
@@ -4338,37 +4446,50 @@ void EmpCylSL::make_coefficients(bool compute)
   }
   
 
-  if (compute and PCAVAR) {
+  if (compute and (covar or PCAVAR)) {
 
     for (unsigned T=0; T<sampT; T++) {
       for (int mm=0; mm<=MMAX; mm++) {
-	for (int nn=0; nn<rank3; nn++) {
-	  MPIin[mm*rank3 + nn] = covV(0, T, mm)[nn];
-	  for (int oo=0; oo<rank3; oo++) {
-	    MPIin2[mm*rank3*rank3 + nn*rank3 + oo] = covM(0, T, mm)(nn, oo);
-	  }
+
+	if (covar and use_mpi) {
+	  MPI_Allreduce ( MPI_IN_PLACE, MV[0][T][mm].data(), MV[0][T][mm].size(),
+			  MPI_DOUBLE_COMPLEX, MPI_SUM, MPI_COMM_WORLD);
+	  MPI_Allreduce ( MPI_IN_PLACE, VC[0][T][mm].data(), VC[0][T][mm].size(),
+			  MPI_DOUBLE_COMPLEX, MPI_SUM, MPI_COMM_WORLD);
 	}
-      }
 
-      if (use_mpi) {
-	MPI_Allreduce ( MPIin.data(), MPIout.data(), rank3*(MMAX+1),
-			MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+	if (PCAVAR) {
+	  for (int nn=0; nn<rank3; nn++) {
+	    MPIin[mm*rank3 + nn] = covV(0, T, mm)[nn];
+	    for (int oo=0; oo<rank3; oo++) {
+	      MPIin2[mm*rank3*rank3 + nn*rank3 + oo] = covM(0, T, mm)(nn, oo);
+	    }
+	  }
+
+	  if (use_mpi) {
+	    MPI_Allreduce ( MPIin.data(), MPIout.data(), rank3*(MMAX+1),
+			    MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
 	
-	MPI_Allreduce ( MPIin2.data(), MPIout2.data(), rank3*rank3*(MMAX+1),
-			MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
-      } else {
-	MPIout  = MPIin;
-	MPIout2 = MPIin2;
-      }
+	    MPI_Allreduce ( MPIin2.data(), MPIout2.data(), rank3*rank3*(MMAX+1),
+			    MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+	  } else {
+	    MPIout  = MPIin;
+	    MPIout2 = MPIin2;
+	  }
 
-      for (int mm=0; mm<=MMAX; mm++) {
-	for (int nn=0; nn<rank3; nn++) {
-	  covV(0, T, mm)[nn] = MPIout[mm*rank3 + nn];
-	  for (int oo=0; oo<rank3; oo++) {
-	    covM(0, T, mm)(nn, oo) = MPIout2[mm*rank3*rank3 + nn*rank3 + oo];
+	  for (int mm=0; mm<=MMAX; mm++) {
+	    for (int nn=0; nn<rank3; nn++) {
+	      covV(0, T, mm)[nn] = MPIout[mm*rank3 + nn];
+	      for (int oo=0; oo<rank3; oo++) {
+		covM(0, T, mm)(nn, oo) = MPIout2[mm*rank3*rank3 + nn*rank3 + oo];
+	      }
+	    }
 	  }
 	}
+	// END: PCAVAR stanza
+
       }
+      // M loop
 
     } // T loop
 
@@ -4397,10 +4518,9 @@ void EmpCylSL::make_coefficients(bool compute)
 	  accum_sin[mm][nn] = MPIout[mm*rank3 + nn];
   }
   
-  if (compute and PCAVAR) {
+  if (compute and (compute or PCAVAR)) {
 				// Mass used to compute variance in
 				// each partition
-
     if (use_mpi) {
       MPI_Allreduce ( &numbT1[0][0], &numbT[0], sampT,
 		      MPI_UNSIGNED, MPI_SUM, MPI_COMM_WORLD);
@@ -4412,7 +4532,6 @@ void EmpCylSL::make_coefficients(bool compute)
       massT = massT1[0];
     }
     
-    
   } // END: 'compute' stanza
   
   std::fill(coefs_made.begin(), coefs_made.end(), true);
@@ -4794,6 +4913,15 @@ void EmpCylSL::pca_hall(bool compute, bool subsamp)
       if (PCAEOF) 
 	for (auto & v : tvar[nth]) v.setZero();
 
+      if (covar) {
+	for (unsigned T=0; T<sampT; T++) {
+	  for (unsigned mm=0; mm<=MMAX; mm++) {
+	    MV[nth][T][mm].setZero();
+	    VC[nth][T][mm].setZero();
+	  }
+	}
+      }
+
       if (PCAVAR) {
 	for (unsigned T=0; T<sampT; T++) {
 	  numbT1[nth][T] = 0;
@@ -4812,6 +4940,46 @@ void EmpCylSL::pca_hall(bool compute, bool subsamp)
 }
 
 
+/** Return a vector of tuples of basis functions and the
+    covariance matrix for subsamples of particles */
+std::vector<std::vector<EmpCylSL::CoefCovarType>>
+EmpCylSL::getCoefCovariance()
+{
+  std::vector<std::vector<EmpCylSL::CoefCovarType>> ret;
+
+  if (covar) {
+    ret.resize(sampT);
+    for (unsigned T=0; T<sampT; T++) {
+      ret[T].resize(MMAX+1);
+      for (int M=0; M<=MMAX; M++)  {
+	std::get<0>(ret[T][M]) = VC[0][T][M];
+	std::get<1>(ret[T][M]) = MV[0][T][M];
+      }
+    }
+  }
+
+  return ret;
+}
+
+std::tuple<Eigen::VectorXi, Eigen::VectorXd>
+EmpCylSL::getCovarSamples()
+{
+  std::tuple<Eigen::VectorXi, Eigen::VectorXd> ret;
+
+  if (covar) {
+    std::get<0>(ret).resize(sampT);
+    std::get<1>(ret).resize(sampT);
+    std::get<0>(ret).setZero();
+    std::get<1>(ret).setZero();
+    for (unsigned T=0; T<sampT; T++) {
+      std::get<0>(ret)(T) = numbT[T];
+      std::get<1>(ret)(T) = massT[T];
+    }
+  }
+  
+  return ret;
+}
+
 void EmpCylSL::get_trimmed
 (double snr,
  std::vector<Eigen::VectorXd>& ac_cos, std::vector<Eigen::VectorXd>& ac_sin,
diff --git a/exputil/ParticleReader.cc b/exputil/ParticleReader.cc
index 52a33dd0c..9a177bdd2 100644
--- a/exputil/ParticleReader.cc
+++ b/exputil/ParticleReader.cc
@@ -294,6 +294,18 @@ namespace PR {
   
   const Particle* GadgetNative::firstParticle()
   {
+    if (done) {
+      if (myid==0 and _verbose)
+      std::cout << "---- ParticleReader::GadgetNative: "
+		<< "restarting reader" << std::endl;
+      done = false;
+      curfile = _files.begin();
+      if (not nextFile()) {	// Try opening
+	std::cerr << "GadgetNative: no files found; logic error?"
+		  << std::endl;
+      }
+    }
+
     pcount = 0;
     
     return &particles[pcount++];
@@ -305,6 +317,7 @@ namespace PR {
       return &particles[pcount++];
     } else {
       if (nextFile()) return firstParticle();
+      done = true;
       return 0;
     }
   }
@@ -650,6 +663,17 @@ namespace PR {
   
   const Particle* GadgetHDF5::firstParticle()
   {
+    if (done) {
+      if (myid==0 and _verbose)
+      std::cout << "---- ParticleReader::GadgetHDF5: "
+		<< "restarting reader" << std::endl;
+      done = false;
+      curfile = _files.begin();
+      if (not nextFile()) {
+	std::cerr << "GadgetHDF5: no files found; logic error?" << std::endl;
+      }
+    }
+
     pcount = 0;
     
     return & particles[pcount++];
@@ -661,7 +685,8 @@ namespace PR {
       return & particles[pcount++];
     } else {
       if (nextFile()) return firstParticle();
-      else return 0;
+      done = true;
+      return 0;
     }
   }
    
@@ -705,7 +730,7 @@ namespace PR {
   }
   
 
-  PSPhdf5::PSPhdf5(const std::vector<std::string>& files, bool verbose) : ParticleReader()
+  PSPhdf5::PSPhdf5(const std::vector<std::string>& files, bool verbose) : PSP(verbose)
   {
     _files   = files;
     _verbose = verbose;
@@ -1221,6 +1246,17 @@ namespace PR {
 
   const Particle* PSPhdf5::firstParticle()
   {
+    if (done) {
+      if (myid==0 and _verbose)
+	std::cout << "---- ParticleReader::PSPhdf5: "
+		  << "restarting reader" << std::endl;
+      done = false;
+      curfile = _files.begin();
+      if (not nextFile()) {
+	std::cerr << "PSPhdf5: no files found; logic error?" << std::endl;
+      }
+    }
+
     pcount = 0;
     
     return & particles[pcount++];
@@ -1232,7 +1268,8 @@ namespace PR {
       return & particles[pcount++];
     } else {
       if (nextFile()) return firstParticle();
-      else return 0;
+      done = true;
+      return 0;
     }
   }
   
diff --git a/include/EmpCylSL.H b/include/EmpCylSL.H
index 7bc412bbe..ebb2af947 100644
--- a/include/EmpCylSL.H
+++ b/include/EmpCylSL.H
@@ -94,6 +94,13 @@ protected:
   //! EOF variance computation
   using VarMat = std::vector< std::vector< std::vector< std::vector<double> > > >;
   VarMat SC, SS, SCe, SCo, SSe, SSo;
+
+  //! Complex covariance implemenation
+  using CVarMat = std::vector< std::vector< std::vector< Eigen::MatrixXcd> > >;
+  using CVarVec = std::vector< std::vector< std::vector< Eigen::VectorXcd> > >;
+  CVarMat MV;
+  CVarVec VC;
+  bool covar = false;
   //@}
 
   std::vector<Eigen::MatrixXd> var, varE, varO;
@@ -210,7 +217,7 @@ protected:
   //! Coefficient magic number
   inline static const unsigned int cmagic = 0xc0a57a3;
 
-  std::vector<short> coefs_made;
+  std::vector<bool> coefs_made;
   bool eof_made;
 
   SphModTblPtr make_sl();
@@ -630,6 +637,15 @@ public:
   //! Initialize PCA work space
   void init_pca();
 
+  //! Initialize covariance work space
+  void init_covar();
+
+  //! Turn on/off covariance computation
+  void set_covar(bool tf) {
+    covar = tf;
+    if (tf) init_covar();
+  }
+
   //! Necessary member function currently unused (change design?)
   void determine_coefficients() {};
 
@@ -930,6 +946,15 @@ public:
   //! Check orthogonality for basis (pyEXP style)
   std::vector<Eigen::MatrixXd> orthoCheck();
 
+  //! Return a vector of tuples of basis coefficients and the
+  //! covariance matrix for subsamples of particles
+  using CoefCovarType = std::tuple<Eigen::VectorXcd, Eigen::MatrixXcd>;
+  std::vector<std::vector<CoefCovarType>> getCoefCovariance();
+
+  //! Covariance statistics
+  std::tuple<Eigen::VectorXi, Eigen::VectorXd> getCovarSamples();
+  
+
 #if HAVE_LIBCUDA==1
   cudaMappingConstants getCudaMappingConstants();
 
diff --git a/include/ParticleReader.H b/include/ParticleReader.H
index b96986446..b9b02fad6 100644
--- a/include/ParticleReader.H
+++ b/include/ParticleReader.H
@@ -35,6 +35,7 @@ namespace PR
     static std::vector<std::string> readerTypes;
     int numprocs, myid;
     bool use_mpi;
+    bool done = false;
     
     //! Check for a list of directories or files and not both
     static bool parseDirectoryList(const std::vector<std::string>& files);
@@ -234,90 +235,6 @@ namespace PR
     
   };
   
-  class PSPhdf5 : public ParticleReader
-  {
-  protected:
-
-    //@{
-    //! Parameters
-    double time;
-    std::vector<Particle> particles;
-    bool _verbose;
-    std::vector<std::string> _files;
-
-    std::vector<double> mass;
-    std::vector<int> Niattrib, Ndattrib;
-    std::vector<unsigned long> npart, nptot;
-    unsigned long totalCount;
-    bool real4, gadget4;
-    int ntypes, nfiles;
-    //@}
-    
-    //! Current file
-    std::vector<std::string>::iterator curfile;
-
-    //! Components
-    std::vector<std::string> comps;
-
-    //@{
-    //! Current component and index
-    std::string curcomp;
-    int curindx;
-    //@}
-
-    unsigned pcount;
-    template <typename Scalar> void read_and_load();
-    template <typename Scalar> void read_and_load_gadget4();
-    template <typename Scalar> void read_and_load_psp();
-
-    void getInfo();
-    void packParticle();
-    bool nextFile();
-
-  public:
-    
-    //! Constructor
-    PSPhdf5(const std::vector<std::string>& file, bool verbose=false);
-    
-    //! Select a particular particle type and reset the iterator
-    virtual void SelectType(const std::string& type) {
-      auto it = std::find(comps.begin(), comps.end(), type);
-      if (it != comps.end()) {
-	curcomp = type;
-	curindx = it - comps.begin();
-      }
-      else {
-	std::cerr << "PSPhdf5 error: could not find particle component <" << type << ">" << std::endl;
-	std::cerr << "Available particle components are:";
-	for (auto s : comps) std::cerr << " " << s;
-	std::cerr << std::endl;
-	throw std::runtime_error("PSPhdf5: non-existent component");
-      }
-
-      curfile = _files.begin();	// Set to first file and open
-      nextFile();
-    }
-    
-    //! Number of particles in the chosen type
-    virtual unsigned long CurrentNumber() { return totalCount; }
-    
-    //! Return list of particle types
-    virtual std::vector<std::string> GetTypes() { return comps; }
-    
-    //! Get current time
-    virtual double CurrentTime() { return time; }
-    
-    //! Reset to beginning of particles for this component
-    virtual const Particle* firstParticle();
-    
-    //! Get the next particle
-    virtual const Particle* nextParticle();
-    
-    //! Get file cound per snapshot
-    virtual int NumFiles() { return nfiles; }
-
-  };
-  
   class PSPstanza 
   {
   public:
@@ -646,6 +563,90 @@ namespace PR
     //@}
   };
 
+  class PSPhdf5 : public PSP
+  {
+  protected:
+
+    //@{
+    //! Parameters
+    double time;
+    std::vector<Particle> particles;
+    bool _verbose;
+    std::vector<std::string> _files;
+
+    std::vector<double> mass;
+    std::vector<int> Niattrib, Ndattrib;
+    std::vector<unsigned long> npart, nptot;
+    unsigned long totalCount;
+    bool real4, gadget4;
+    int ntypes, nfiles;
+    //@}
+    
+    //! Current file
+    std::vector<std::string>::iterator curfile;
+
+    //! Components
+    std::vector<std::string> comps;
+
+    //@{
+    //! Current component and index
+    std::string curcomp;
+    int curindx;
+    //@}
+
+    unsigned pcount;
+    template <typename Scalar> void read_and_load();
+    template <typename Scalar> void read_and_load_gadget4();
+    template <typename Scalar> void read_and_load_psp();
+
+    void getInfo();
+    void packParticle();
+    bool nextFile();
+
+  public:
+    
+    //! Constructor
+    PSPhdf5(const std::vector<std::string>& file, bool verbose=false);
+    
+    //! Select a particular particle type and reset the iterator
+    virtual void SelectType(const std::string& type) {
+      auto it = std::find(comps.begin(), comps.end(), type);
+      if (it != comps.end()) {
+	curcomp = type;
+	curindx = it - comps.begin();
+      }
+      else {
+	std::cerr << "PSPhdf5 error: could not find particle component <" << type << ">" << std::endl;
+	std::cerr << "Available particle components are:";
+	for (auto s : comps) std::cerr << " " << s;
+	std::cerr << std::endl;
+	throw std::runtime_error("PSPhdf5: non-existent component");
+      }
+
+      curfile = _files.begin();	// Set to first file and open
+      nextFile();
+    }
+    
+    //! Number of particles in the chosen type
+    virtual unsigned long CurrentNumber() { return nptot[curindx]; }
+    
+    //! Return list of particle types
+    virtual std::vector<std::string> GetTypes() { return comps; }
+    
+    //! Get current time
+    virtual double CurrentTime() { return time; }
+    
+    //! Reset to beginning of particles for this component
+    virtual const Particle* firstParticle();
+    
+    //! Get the next particle
+    virtual const Particle* nextParticle();
+    
+    //! Get file cound per snapshot
+    virtual int NumFiles() { return nfiles; }
+
+  };
+  
   typedef std::shared_ptr<ParticleReader> PRptr;
 }
   
diff --git a/include/tipsy.H b/include/tipsy.H
index 01be7e3df..b598fcc7e 100644
--- a/include/tipsy.H
+++ b/include/tipsy.H
@@ -1,6 +1,7 @@
 #ifndef _tipsy_H
 #define _tipsy_H
 
+#include <cstdint>
 #include <map>
 
 // For MPI awareness
diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc
index 7c217d9c1..10300df08 100644
--- a/pyEXP/BasisWrappers.cc
+++ b/pyEXP/BasisWrappers.cc
@@ -51,12 +51,12 @@ void BasisFactoryClasses(py::module &m)
         a phase-space snapshot.  This is a specialized version of FieldBasis.
 
     Each of these bases take a YAML configuration file as input. These parameter
-    lists are as subset of and have the same structure as thosed used by EXP.
+    lists are as subset of and have the same structure as those used by EXP.
     The factory and the individual constructors will check the parameters keys
     and warn of mismatches for safety.  See the EXP documentation and the pyEXP
     examples for more detail.  The first four bases are the most often used bi-
     orthogonal basis types used for computing the potential and forces from
-    density distributions.  Other biorthgonal bases in EXP but not in pyEXP
+    density distributions.  Other biorthogonal bases in EXP but not in pyEXP
     include those for cubic and slab geometries and other special-purpose bases
     such as the Hernquist, Clutton-Brock sphere and two-dimensional disk basis.
     These will be made available in a future release if there is demand.  Note
@@ -102,7 +102,7 @@ void BasisFactoryClasses(py::module &m)
 
        def myFunctor(m, pos, vel, index):
           ret = False  # Default return value
-          # some caculation with scalar mass, pos array, vel array and
+          # some calculation with scalar mass, pos array, vel array and
           # integer index that sets ret to True if desired . . . 
           return ret
 
@@ -128,7 +128,7 @@ void BasisFactoryClasses(py::module &m)
     phase-space field functor to cylindrical or spherical velocities
     based on the 'dof' parameter.  More on 'dof' below.
 
-    Scalablility
+    Scalability
     ------------
     createFromArray() is a convenience method allows you to transform
     coordinates and preprocess phase space using your own methods and
@@ -198,7 +198,7 @@ void BasisFactoryClasses(py::module &m)
     Coefs data base and installs the interpolated coefficients for the
     current time in the basis instance.  The SingleTimeFunc interpolates on
     the Coefs data base for a single fixed time and sets the interpolated
-    coefficients once at the beginning of the integration.  This implementes
+    coefficients once at the beginning of the integration.  This implements
     a fixed potential model.  AccelFunc can be inherited by a native Python
     class and the evalcoefs() may be implemented in Python and passed to
     IntegrateOrbits in the same way as a native C++ class.
@@ -282,8 +282,8 @@ void BasisFactoryClasses(py::module &m)
       PYBIND11_OVERRIDE(FCReturn, Basis, getFieldsCoefs, x, y, z, coefs);
     }
     
-    void accumulate(double x, double y, double z, double mass) override {
-      PYBIND11_OVERRIDE_PURE(void, Basis, accumulate, mass, x, y, z, mass);
+    void accumulate(double x, double y, double z, double mass, unsigned long int indx) override {
+      PYBIND11_OVERRIDE_PURE(void, Basis, accumulate, x, y, z, mass, indx);
     }
 
     void reset_coefs(void) override {
@@ -363,9 +363,10 @@ void BasisFactoryClasses(py::module &m)
     }
 
     void accumulate(double m, double x, double y, double z,
-		    double u, double v, double w) override
+		    double u, double v, double w,
+		    unsigned long int indx) override
     {
-      PYBIND11_OVERRIDE(void, FieldBasis, accumulate, m, x, y, z, u, v, w);
+      PYBIND11_OVERRIDE(void, FieldBasis, accumulate, m, x, y, z, u, v, w, indx);
     }
     
     void reset_coefs(void) override {
@@ -422,8 +423,8 @@ void BasisFactoryClasses(py::module &m)
     // Inherit the constructors
     using BasisClasses::BiorthBasis::BiorthBasis;
 
-    void accumulate(double x, double y, double z, double mass) override {
-      PYBIND11_OVERRIDE_PURE(void, BiorthBasis, accumulate, x, y, z, mass);
+    void accumulate(double x, double y, double z, double mass, unsigned long int indx) override {
+      PYBIND11_OVERRIDE_PURE(void, BiorthBasis, accumulate, x, y, z, mass, indx);
     }
 
     void reset_coefs(void) override {
@@ -438,6 +439,23 @@ void BasisFactoryClasses(py::module &m)
       PYBIND11_OVERRIDE_PURE(void, BiorthBasis, set_coefs, coefs);
     }
 
+    using CCelement = std::tuple<Eigen::VectorXcd, Eigen::MatrixXcd>;
+    using CCreturn  = std::vector<std::vector<CCelement>>;
+
+    CCreturn getCoefCovariance(void) override {
+      PYBIND11_OVERRIDE(CCreturn, BiorthBasis, getCoefCovariance,);
+    }
+
+    using Selement = std::tuple<Eigen::VectorXi, Eigen::VectorXd>;
+
+    Selement getCovarSamples() override {
+      PYBIND11_OVERRIDE(Selement, BiorthBasis, getCovarSamples,);
+    }
+
+    void enableCoefCovariance(bool pcavar, int nsamples) override {
+      PYBIND11_OVERRIDE(void, BiorthBasis, enableCoefCovariance, pcavar, nsamples);
+    }
+
   };
 
   class PySpherical : public Spherical
@@ -487,8 +505,8 @@ void BasisFactoryClasses(py::module &m)
       PYBIND11_OVERRIDE(std::vector<double>, Spherical, getFields, x, y, z);
     }
 
-    void accumulate(double x, double y, double z, double mass) override {
-      PYBIND11_OVERRIDE(void, Spherical, accumulate, x, y, z, mass);
+    void accumulate(double x, double y, double z, double mass, unsigned long int indx) override {
+      PYBIND11_OVERRIDE(void, Spherical, accumulate, x, y, z, mass, indx);
     }
 
     void reset_coefs(void) override {
@@ -503,6 +521,10 @@ void BasisFactoryClasses(py::module &m)
       PYBIND11_OVERRIDE_PURE(std::vector<Eigen::MatrixXd>, Spherical, orthoCheck, knots);
     }
 
+    void enableCoefCovariance(bool pcavar, int nsamples) override {
+      PYBIND11_OVERRIDE(void, Spherical, enableCoefCovariance, pcavar, nsamples);
+    }
+
   };
 
 
@@ -542,8 +564,8 @@ void BasisFactoryClasses(py::module &m)
       PYBIND11_OVERRIDE(std::vector<double>, Cylindrical, getFields, x, y, z);
     }
 
-    void accumulate(double x, double y, double z, double mass) override {
-      PYBIND11_OVERRIDE(void, Cylindrical, accumulate, x, y, z, mass);
+    void accumulate(double x, double y, double z, double mass, unsigned long int indx) override {
+      PYBIND11_OVERRIDE(void, Cylindrical, accumulate, x, y, z, mass, indx);
     }
 
     void reset_coefs(void) override {
@@ -554,6 +576,10 @@ void BasisFactoryClasses(py::module &m)
       PYBIND11_OVERRIDE(void, Cylindrical, make_coefs,);
     }
 
+    void enableCoefCovariance(bool pcavar, int nsamples) override {
+      PYBIND11_OVERRIDE(void, Cylindrical, enableCoefCovariance, pcavar, nsamples);
+    }
+
   };
 
 
@@ -611,9 +637,9 @@ void BasisFactoryClasses(py::module &m)
       PYBIND11_OVERRIDE(std::vector<double>, FlatDisk, getFields, x, y, z);
     }
 
-    void accumulate(double x, double y, double z, double mass) override
+    void accumulate(double x, double y, double z, double mass, unsigned long int indx) override
     {
-      PYBIND11_OVERRIDE(void, FlatDisk, accumulate, x, y, z, mass);
+      PYBIND11_OVERRIDE(void, FlatDisk, accumulate, x, y, z, mass, indx);
     }
 
     void reset_coefs(void) override
@@ -682,9 +708,9 @@ void BasisFactoryClasses(py::module &m)
       PYBIND11_OVERRIDE(std::vector<double>, CBDisk, getFields, x, y, z);
     }
 
-    void accumulate(double x, double y, double z, double mass) override
+    void accumulate(double x, double y, double z, double mass, unsigned long int indx) override
     {
-      PYBIND11_OVERRIDE(void, CBDisk, accumulate, x, y, z, mass);
+      PYBIND11_OVERRIDE(void, CBDisk, accumulate, x, y, z, mass, indx);
     }
 
     void reset_coefs(void) override
@@ -756,9 +782,9 @@ void BasisFactoryClasses(py::module &m)
       PYBIND11_OVERRIDE(std::vector<double>, Slab, getFields, x, y, z);
     }
 
-    void accumulate(double x, double y, double z, double mass) override
+    void accumulate(double x, double y, double z, double mass, unsigned long int indx) override
     {
-      PYBIND11_OVERRIDE(void, Slab, accumulate, x, y, z, mass);
+      PYBIND11_OVERRIDE(void, Slab, accumulate, x, y, z, mass, indx);
     }
 
     void reset_coefs(void) override
@@ -830,9 +856,9 @@ void BasisFactoryClasses(py::module &m)
       PYBIND11_OVERRIDE(std::vector<double>, Cube, getFields, x, y, z);
     }
 
-    void accumulate(double x, double y, double z, double mass) override
+    void accumulate(double x, double y, double z, double mass, unsigned long int indx) override
     {
-      PYBIND11_OVERRIDE(void, Cube, accumulate, x, y, z, mass);
+      PYBIND11_OVERRIDE(void, Cube, accumulate, x, y, z, mass, indx);
     }
 
     void reset_coefs(void) override
@@ -845,6 +871,23 @@ void BasisFactoryClasses(py::module &m)
       PYBIND11_OVERRIDE(void, Cube, make_coefs,);
     }
 
+    using CCelement = std::tuple<Eigen::VectorXcd, Eigen::MatrixXcd>;
+    using CCreturn  = std::vector<std::vector<CCelement>>;
+
+    CCreturn getCoefCovariance(void) override {
+      PYBIND11_OVERRIDE(CCreturn, Cube, getCoefCovariance,);
+    }
+
+    using Selement = std::tuple<Eigen::VectorXi, Eigen::VectorXd>;
+
+    Selement getCovarSamples() override {
+      PYBIND11_OVERRIDE(Selement, Cube, getCovarSamples,);
+    }
+
+    void enableCoefCovariance(bool pcavar, int nsamples) override {
+      PYBIND11_OVERRIDE(void, Cube, enableCoefCovariance, pcavar, nsamples);
+    }
+
   };
 
 
@@ -1091,7 +1134,7 @@ void BasisFactoryClasses(py::module &m)
          See also
          --------
          setNonInertial : set non-inertial data from an Orient file
-         setNonInertialAccel : set the non-inertial acceration
+         setNonInertialAccel : set the non-inertial acceleration
          )",
 	 py::arg("N"), py::arg("times"), py::arg("pos")
          )
@@ -1118,7 +1161,7 @@ void BasisFactoryClasses(py::module &m)
          See also
          --------
          setNonInertial : set non-inertial data from a time series of values
-         setNonInertialAccel : set the non-inertial acceration
+         setNonInertialAccel : set the non-inertial acceleration
          )",
 	 py::arg("N"), py::arg("orient")
          )
@@ -1319,6 +1362,146 @@ void BasisFactoryClasses(py::module &m)
          makeFromArray : create coefficients contributions
          )",
 	 py::arg("mass"), py::arg("pos"))
+    .def("getCoefCovariance", &BasisClasses::BiorthBasis::getCoefCovariance,
+         R"(
+         Return the coefficient vectors and covariance matrices for each partition of the
+         accumulated particles. The number of partitions is set by the configuration
+         parameter 'sampT' (default: 100).
+
+         The first dimension are the time samples.  The second dimension is the angular
+         index. Each element is a tuple of the coefficient vector and covariance.
+         The values are complex128 and represents the full amplitude and phase information.
+
+         Returns
+         -------
+         arrays: list[list[tuple[np.ndarray, np.ndarray]]]
+             Each element is a tuple (coef_vector, coef_covariance_matrix),
+             where coef_vector and coef_covariance_matrix are numpy.ndarray.
+             Each coef_covariance_matrix is of shape (nmax, nmax). All values are
+             complex128.
+
+         Shape and Indexing
+         ------------------
+         - The first list index is the number of time samples.
+         - The second list index is the angular elements.  For spherical bases, all (l, m) pairs
+           are in triangular index order with l in [0,...,lmax] and m in [0,...,l] for a total of
+           (lmax+1)*(lmax+2)/2 entries.  For cylindrical bases, there are (mmax+1) harmonic
+           entries for each value m in [0,...,mmax].
+         - Each covariance matrix is of shape (nmax, nmax), where nmax is the number of basis
+           functions
+
+         See also
+         --------
+         getCovarSamples     : get counts and mass in each partition
+         writeCoefCovariance : write the coefficient vectors, covariance
+                               matrices, and metadata to an HDF5 file.
+         )")
+    .def("getCovarSamples", &BasisClasses::BiorthBasis::getCovarSamples,
+         R"(
+         Return a vector of counts and mass used to compute the partitioned 
+         vectors and covariance.  The length of both returned vectors equals 'sampT'
+
+         Parameters
+         ----------
+         None
+
+         Returns
+         -------
+         arrays: tuple(ndarray, ndarray)
+              a tuple of arrays containing the counts and mass in each
+              partitioned sample
+
+         See also
+         --------
+         getCoefCovariance   : get the coefficient vectors and covariance matrices
+                               for the partitioned phase space.
+         writeCoefCovariance : write the coefficient vectors, covariance
+                               matrices, and metadata to an HDF5 file.
+         )")
+    .def("writeCoefCovariance", &BasisClasses::BiorthBasis::writeCoefCovariance,
+         R"(
+         Write the partitioned coefficient vectors and covariance matrices
+         to a specified HDF5 file.  The number of partitions is set by the
+         configuration parameter 'sampT' and defaults to 100.
+
+         On first call, the file is created, metadata is written, and the
+         coefficient vectors and covariance matrices are stored.  On subsequent
+         calls, the file is updated with new covariance datasets.
+
+         The file will be called 'coefcovar.<compname>.<runtag>.h5'.
+
+         Parameters
+         ----------
+         compname : str
+                    the component/basis name segment used in the output HDF5 filename
+         runtag   : str
+                    the run identifier tag
+         time     : float
+                    the snapshot time
+
+         Returns
+         -------
+         None
+
+         See also
+         --------
+         getCoefCovariance : get the coefficient vectors and covariance matrices
+                             for the partitioned phase space.
+         getCovarSamples   : get counts and mass in each partition
+         )", py::arg("compname"), py::arg("runtag"), py::arg("time")=0.0)
+    .def("setCovarFloatType", &BasisClasses::BiorthBasis::setCovarFloatType,
+         R"(
+         Set the floating point type used for covariance storage in HDF5. Set to
+         true for float4.  Type float8 is the default (false).  This does not affect
+         the return type in getCoefCovariance().
+
+         Parameters
+         ----------
+         floatType : bool
+                 Use 'float32' rather than 'float64' if true
+
+         Returns
+         -------
+         None
+         )")
+    .def("enableCoefCovariance", &BasisClasses::BiorthBasis::enableCoefCovariance,
+	 R"(
+         Enable or disable the coefficient covariance computation and set the
+         default number of partitions to use for the covariance computation.
+
+         Parameters
+         ----------
+         pcavar   : bool
+                    enable (true) or disable (false) the covariance computation
+         nsamples : int
+                    number of time partitions to use for covariance computation
+
+         Returns
+         -------
+         None
+	)", py::arg("pcavar"), py::arg("nsamples")=100)
+    .def("setCovarH5Compress", &BasisClasses::BiorthBasis::setCovarH5Compress,
+         R"(
+	 Set the HDF5 compression level for covariance storage in HDF5.  The Szip
+     	 compression algorithm may also be enabled but seems to not have better
+         performance for these data than the standard Gzip compression.  Gzip with
+      	 level 5 and shuffling is enabled by default.
+
+         Parameters
+         ----------
+         compress : int
+		 HDF5 Gzip compression level 0-9 (default: 5)
+         chunkSize : int
+		 HDF5 chunk size for dataset storage (default: 1024*1024)
+         shuffle : bool
+		 Use shuffle filter if true (default: true)
+         szip : bool
+		 Use Szip compression algorithm (default: false)
+
+         Returns
+         -------
+         None
+         )", py::arg("compress")=5, py::arg("chunkSize")=1024*1024, py::arg("shuffle")=true, py::arg("azip")=false)
     .def("makeFromFunction", &BasisClasses::BiorthBasis::makeFromFunction,
 	 py::call_guard<py::gil_scoped_release>(),
 	 R"(
@@ -1375,13 +1558,13 @@ void BasisFactoryClasses(py::module &m)
 	 )
     .def("getFields", &BasisClasses::BiorthBasis::getFields,
 	 R"(
-         Return the field evaluations for a given cartesian position. The
+         Return the field evaluations for a given Cartesian position. The
          fields include density, potential, and force.  The density and
          potential evaluations are separated into full, axisymmetric and
          non-axisymmetric contributions.
 
          You can get the fields labels by using the __call__ method of the
-         basis object.  This is equilevalent to a tuple of the getFields()
+         basis object.  This is equivalent to a tuple of the getFields()
          output with a list of field labels.
 
          Parameters
@@ -1405,7 +1588,7 @@ void BasisFactoryClasses(py::module &m)
 	 py::arg("x"), py::arg("y"), py::arg("z"))
     .def("getAccel", py::overload_cast<double, double, double>(&BasisClasses::BiorthBasis::getAccel),
 	 R"(
-         Return the acceleration for a given cartesian position
+         Return the acceleration for a given Cartesian position
 
          Parameters
          ----------
@@ -1429,7 +1612,7 @@ void BasisFactoryClasses(py::module &m)
 	 py::arg("x"), py::arg("y"), py::arg("z"))
     .def("getAccel", py::overload_cast<Eigen::Ref<const Eigen::VectorXd>, Eigen::Ref<const Eigen::VectorXd>, Eigen::Ref<const Eigen::VectorXd>>(&BasisClasses::BiorthBasis::getAccel),
 	 R"(
-         Return the acceleration for a given cartesian position
+         Return the acceleration for a given Cartesian position
 
          Parameters
          ----------
@@ -1453,7 +1636,7 @@ void BasisFactoryClasses(py::module &m)
 	 py::arg("x"), py::arg("y"), py::arg("z"))
     .def("getAccel", py::overload_cast<Eigen::Ref<const RowMatrixXd>>(&BasisClasses::BiorthBasis::getAccel),
 	 R"(
-         Return the acceleration for an array of cartesian positions
+         Return the acceleration for an array of Cartesian positions
 
          Parameters
          ----------
@@ -1473,7 +1656,7 @@ void BasisFactoryClasses(py::module &m)
 	 py::arg("pos"))
     .def("getAccelArray", py::overload_cast<Eigen::Ref<const Eigen::VectorXd>, Eigen::Ref<const Eigen::VectorXd>, Eigen::Ref<const Eigen::VectorXd>>(&BasisClasses::BiorthBasis::getAccel),
 	 R"(
-         Return the acceleration for a given cartesian position
+         Return the acceleration for a given Cartesian position
 
          Parameters
          ----------
@@ -1497,7 +1680,7 @@ void BasisFactoryClasses(py::module &m)
 	 py::arg("x"), py::arg("y"), py::arg("z"))
     .def("getFieldsCoefs", &BasisClasses::BiorthBasis::getFieldsCoefs,
 	 R"(
-         Return the field evaluations for a given cartesian position
+         Return the field evaluations for a given Cartesian position
          for every frame in a coefficient set.  The field evaluations are
          produced by a call to getFields().
 
@@ -1534,7 +1717,7 @@ void BasisFactoryClasses(py::module &m)
          coordinates for SphericalSL and Bessel, cylindrical coordinates for
          Cylindrical, FlatDisk, and CBDisk, and Cartesian coordinates for the 
          Slab and Cube. This member function can be used to override the
-         default. The available coorindates are: 'spherical', 'cylindrical', 
+         default. The available coordinates are: 'spherical', 'cylindrical', 
          'cartesian'.
 
          Parameters
@@ -1559,9 +1742,9 @@ void BasisFactoryClasses(py::module &m)
          -------
          None
          )")
-    .def("accumulate", [](BasisClasses::BiorthBasis& A, double x, double y, double z, double mass)
+    .def("accumulate", [](BasisClasses::BiorthBasis& A, double x, double y, double z, double mass, unsigned long int indx)
     {
-      return A.accumulate(x, y, z, mass);
+      return A.accumulate(x, y, z, mass, indx);
     },
       R"(
          Add the contribution of a single particle to the coefficients
@@ -1576,12 +1759,14 @@ void BasisFactoryClasses(py::module &m)
              z-axis position
          mass : float
              particle mass
+         indx : unsigned long int
+             particle index for selection functor
 
          Returns
          -------
          None
         )", 
-	 py::arg("x"), py::arg("y"), py::arg("z"), py::arg("mass"))
+      py::arg("x"), py::arg("y"), py::arg("z"), py::arg("mass"), py::arg("indx")=0)
     .def("getMass",            &BasisClasses::BiorthBasis::getMass,
 	 R"(
          Return the total mass of particles contributing the current coefficient set
@@ -1601,7 +1786,7 @@ void BasisFactoryClasses(py::module &m)
          )")
     .def("make_coefs",         &BasisClasses::BiorthBasis::make_coefs,
 	 R"(
-         Create the coefficients after particle accumuluation is complete
+         Create the coefficients after particle accumulation is complete
 
          Returns
          -------
@@ -1705,7 +1890,7 @@ void BasisFactoryClasses(py::module &m)
 	   return A.orthoCheck(knots);
 	 },
 	R"(
-        Check orthgonality of basis functions by quadrature
+        Check orthogonality of basis functions by quadrature
 
         Inner-product matrix of Sturm-Liouville solutions indexed by
         harmonic order used to assess fidelity.
@@ -1739,6 +1924,7 @@ void BasisFactoryClasses(py::module &m)
       )",
       py::arg("cachefile"));
 
+
     py::class_<BasisClasses::Spherical, std::shared_ptr<BasisClasses::Spherical>, PySpherical, BasisClasses::BiorthBasis>(m, "Spherical")
       .def(py::init<const std::string&, const std::string&>(),
 	 R"(
@@ -1791,7 +1977,7 @@ void BasisFactoryClasses(py::module &m)
 	return A.orthoCheck(knots);
       },
 	R"(
-        Check orthgonality of basis functions by quadrature
+        Check orthogonality of basis functions by quadrature
 
         Inner-product matrix of Sturm-Liouville solutions indexed by
         harmonic order used to assess fidelity.
@@ -1825,6 +2011,7 @@ void BasisFactoryClasses(py::module &m)
             cache parameters
         )",
 	py::arg("cachefile"))
+
       .def_static("I", [](int l, int m)
       {
 	if (l<0) throw std::runtime_error("l must be greater than 0");
@@ -1848,9 +2035,10 @@ void BasisFactoryClasses(py::module &m)
             index array packing index
       )",
 	py::arg("l"), py::arg("m"))
+
       .def_static("invI", [](int I)
       {
-	if (I<0) std::runtime_error("I must be an interger greater than or equal to 0");
+	if (I<0) std::runtime_error("I must be an integer greater than or equal to 0");
 	int l = std::floor(0.5*(-1.0 + std::sqrt(1.0 + 8.0 * I)));
 	int m = I - int(l * (l + 1) / 2);
 	return std::tuple<int, int>(l, m);
@@ -1869,8 +2057,7 @@ void BasisFactoryClasses(py::module &m)
         -------
         (l, m) : tuple
             the harmonic indices (l, m).
-      )", py::arg("I"));
-
+        )", py::arg("I"));
   
     py::class_<BasisClasses::SphericalSL, std::shared_ptr<BasisClasses::SphericalSL>, BasisClasses::Spherical>(m, "SphericalSL")
       .def(py::init<const std::string&>(),
@@ -1923,7 +2110,7 @@ void BasisFactoryClasses(py::module &m)
 	return A.orthoCheck(knots);
       },
 	R"(
-        Check orthgonality of basis functions by quadrature
+        Check orthogonality of basis functions by quadrature
 
         Inner-product matrix of Sturm-Liouville solutions indexed by
         harmonic order used to assess fidelity.
@@ -1965,7 +2152,7 @@ void BasisFactoryClasses(py::module &m)
 	return A.orthoCheck(knots);
       },
 	R"(
-        Check orthgonality of basis functions by quadrature
+        Check orthogonality of basis functions by quadrature
 
         Inner-product matrix of Sturm-Liouville solutions indexed by
         harmonic order used to assess fidelity.
@@ -2035,7 +2222,7 @@ void BasisFactoryClasses(py::module &m)
       return A.orthoCheck();
     },
       R"(
-      Check orthgonality of basis functions by quadrature
+      Check orthogonality of basis functions by quadrature
 
       Inner-product matrix of Sturm-Liouville solutions indexed by
       harmonic order used to assess fidelity.
@@ -2123,7 +2310,7 @@ void BasisFactoryClasses(py::module &m)
       return A.orthoCheck();
     },
       R"(
-      Check orthgonality of basis functions by quadrature
+      Check orthogonality of basis functions by quadrature
 
       Inner-product matrix of Sturm-Liouville solutions indexed by
       harmonic order used to assess fidelity.
@@ -2174,7 +2361,7 @@ void BasisFactoryClasses(py::module &m)
 
          Returns
          -------
-         Cube
+         Slab
              the new instance
          )", py::arg("YAMLstring"))
     .def("getBasis", &BasisClasses::Slab::getBasis,
@@ -2221,7 +2408,7 @@ void BasisFactoryClasses(py::module &m)
       return A.orthoCheck();
     },
       R"(
-      Check orthgonality of basis functions by quadrature
+      Check orthogonality of basis functions by quadrature
 
       Inner-product matrix of indexed by flattened wave number (nx, ny, nz) where
       each of nx is in [-nmaxx, nmaxx], ny is in [-nmaxy, nmaxy] and nz is in 
@@ -2259,18 +2446,79 @@ void BasisFactoryClasses(py::module &m)
          Cube
              the new instance
          )", py::arg("YAMLstring"))
+    .def("setCovar", &BasisClasses::Cube::setCovar,
+       R"(
+       Turn on coefficient covariance matrix computation.  Set to OFF by default to
+       save disk space, computatation time, and memory.
+
+       Parameters
+       ----------
+       flag : bool
+              True to enable covariance matrix computation, False to disable it.
+
+       Returns
+       -------
+       None
+
+       Notes
+       -----
+       The covariance matrix computation can use significant memory and disk space
+       for the Cube basis especially and is disabled by default.
+       )", py::arg("flag")=true)
+    .def("index1D", &BasisClasses::Cube::index1D,
+      R"(
+      Returns a flattened 1-d index into the arrays and matrices returned by the
+      getCoefCovariance() routines from wave number indexing.
+
+      Parameters
+      ----------
+      nx : int
+           wave-number index in [0, 2*nmax_x]
+      ny : int
+           wave-number index in [0, 2*nmax_y]
+      nz : int
+           wave-number index in [0, 2*nmax_z]
+
+      Returns
+      -------
+      indx : int
+           flattened index for coefficient array and covariance matrix
+
+      Notes
+      -----
+      Wave-numbers range from 2*pi/L*[-nmax, nmax] in each dimension.  The indices
+      nx, ny, nz describe the wave-number offset by nmax.
+
+      )",
+	 py::arg("nx"), py::arg("ny"), py::arg("nz"))
+    .def("index3D", &BasisClasses::Cube::index3D,
+      R"(
+      Returns a tuple of indices for the wave-numbers (kx, ky, kz) from the flattened
+      1-d index for the arrays and matrices returned by the getCoefCovariance() routines
+
+      Parameters
+      ----------
+      index : int
+              flattened index for coefficient array and covariance matrix
+
+      Returns
+      -------
+      indices : tuple(int, int, int)
+              wave-number indices for each dimension in the ranges [0, 2*nmax-x],
+              [0, 2*nmax_y], [0, 2*nmax_z]
+      )", py::arg("index"))
     .def("orthoCheck", [](BasisClasses::Cube& A)
     {
       return A.orthoCheck();
     },
       R"(
-      Check orthgonality of basis functions by quadrature
+      Check orthogonality of basis functions by quadrature
 
       Inner-product matrix of indexed by flattened wave number (nx, ny, nz) where
       each of nx is in [-nmaxx, nmaxx], and so on for ny and nz.  Each dimension 
       has dx=2*nmaxx+1 wave numbers and similarly for dy and dz.  The index into the
       array is index=(nx+nmaxx)*dx*dy + (ny+nmaxy)*dy + (nz+nmaxz).   This is an 
-      analyic basis so the orthogonality matrix is not a check of any numerical
+      analytic basis so the orthogonality matrix is not a check of any numerical
       computation other than the quadrature itself.  It is included for completeness.
 
       Parameters
@@ -2503,7 +2751,7 @@ void BasisFactoryClasses(py::module &m)
 	return A.orthoCheck();
       },
       R"(
-        Check orthgonality of basis functions by quadrature
+        Check orthogonality of basis functions by quadrature
 
         Inner-product matrix of orthogonal functions
 
@@ -2665,4 +2913,114 @@ void BasisFactoryClasses(py::module &m)
 	py::arg("tinit"), py::arg("tfinal"), py::arg("h"),
 	py::arg("ps"), py::arg("basiscoef"), py::arg("func"),
 	py::arg("nout")=0);
+
+  py::class_<BasisClasses::CovarianceReader, std::shared_ptr<BasisClasses::CovarianceReader>>
+    (m, "CovarianceReader")
+    .def(py::init<const std::string&, int>(),
+       R"(
+       Read a covariance database from an HDF5 file
+
+       Parameters
+       ----------
+       filename : str
+                The HDF5 filename
+       stride   : int, default=1
+                The stride for reading datasets
+
+       Returns
+       -------
+       CovarianceReader
+            the new instance
+
+       )", py::arg("filename"), py::arg("stride")=1)
+    .def("Times", &BasisClasses::CovarianceReader::Times,
+       R"(
+       Get the list of evaluation times
+
+       Parameters
+       ----------
+       None
+
+       Returns
+       -------
+       times : list(float)
+            a list of evaluation times
+      )")
+  .def("getCoefCovariance", static_cast<std::vector<std::vector<std::tuple<Eigen::VectorXcd, Eigen::MatrixXcd>>>
+       (BasisClasses::CovarianceReader::*)(unsigned)>(&BasisClasses::CovarianceReader::getCoefCovariance),
+     R"(
+     Get the covariance matrices for the basis coefficients
+
+     Parameters
+     ----------
+     index : int
+             the time index
+     Returns
+     -------
+      list(list(tuple(numpy.ndarray, numpy.ndarray)))
+	    list of partitioned coefficients and their covariance matrices for
+            each subsample
+    )",
+    py::arg("index"))
+  .def("getCoefCovariance", static_cast<std::vector<std::vector<std::tuple<Eigen::VectorXcd, Eigen::MatrixXcd>>>
+       (BasisClasses::CovarianceReader::*)(double)>(&BasisClasses::CovarianceReader::getCoefCovariance),
+     R"(
+     Get the covariance matrices for the basis coefficients
+
+     Parameters
+     ----------
+     time  : float
+             the evaluation time
+
+     Returns
+     -------
+      list(list(tuple(numpy.ndarray, numpy.ndarray)))
+	    List of partitioned coefficients and their covariance matrices for
+            each subsample. The returns are complex-valued.
+    )",
+    py::arg("time"))
+   .def("getCovarSamples",
+	py::overload_cast<unsigned>(&BasisClasses::CovarianceReader::getCovarSamples),
+      R"(
+      Get sample counts for the covariance computation
+
+      Parameters
+      ----------
+      index : unsigned int
+              the time index
+
+      Returns
+      -------
+      tuple(numpy.ndarray, numpy.ndarray)
+          sample counts and masses for the covariance computation
+      )",
+      py::arg("index"))
+    .def("getCovarSamples", py::overload_cast<double>(&BasisClasses::CovarianceReader::getCovarSamples),
+      R"(
+      Get sample counts for the covariance computation
+
+      Parameters
+      ----------
+      time  : float
+              the evaluation time
+
+      Returns
+      -------
+      tuple(numpy.ndarray, numpy.ndarray)
+          sample counts and masses for the covariance computation
+      )",
+      py::arg("time"))
+   .def("basisIDname", &BasisClasses::CovarianceReader::basisIDname,
+     R"(
+     Get the basis ID name
+
+     Parameters
+     ----------
+     None
+
+     Returns
+     -------
+     BasisID : str
+               the basis ID name
+     )");
 }
diff --git a/pyEXP/CoefWrappers.cc b/pyEXP/CoefWrappers.cc
index 81ccdc266..8aa420b32 100644
--- a/pyEXP/CoefWrappers.cc
+++ b/pyEXP/CoefWrappers.cc
@@ -1053,7 +1053,7 @@ void CoefficientClasses(py::module &m) {
          type : str
              type of coefficient container
          verbose : bool
-             display verbose information.
+             display verbose information (default=False)
 
          Returns
          -------
@@ -1066,7 +1066,7 @@ void CoefficientClasses(py::module &m) {
          simulation
          )",
          py::arg("type"),
-         py::arg("verbose"))
+         py::arg("verbose")=false)
     .def("__call__",
 	 &CoefClasses::Coefs::getData,
          R"(
@@ -1380,12 +1380,12 @@ void CoefficientClasses(py::module &m) {
          Parameters
          ----------
          verbose : bool
-             display verbose information.
+             display verbose information (default=False)
 
          Returns
          -------
          SphCoefs instance
-         )")
+         )", py::arg("verbose")=false)
     .def("__call__",
 	 &CoefClasses::SphCoefs::getMatrix,
          R"(
@@ -1455,12 +1455,12 @@ void CoefficientClasses(py::module &m) {
          Parameters
          ----------
          verbose : bool
-             display verbose information.
+             display verbose information (default=False)
 
          Returns
          -------
          CylCoefs instance
-         )")
+         )", py::arg("verbose")=false)
     .def("__call__",
 	 &CoefClasses::CylCoefs::getMatrix,
          R"(
@@ -1557,12 +1557,12 @@ void CoefficientClasses(py::module &m) {
          Parameters
          ----------
          verbose : bool
-             display verbose information.
+             display verbose information (default=False)
 
          Returns
          -------
          SphFldCoefs instance
-         )")
+         )", py::arg("verbose")=false)
     .def("__call__",
 	 [](CoefClasses::SphFldCoefs& A, double time)
 	 {
@@ -1642,12 +1642,12 @@ void CoefficientClasses(py::module &m) {
          Parameters
          ----------
          verbose : bool
-             display verbose information.
+             display verbose information (default=False)
 
          Returns
          -------
          CylFldCoefs instance
-         )")
+         )", py::arg("verbose")=false)
     .def("__call__",
 	 [](CoefClasses::CylFldCoefs& A, double time)
 	 {
@@ -1721,12 +1721,12 @@ void CoefficientClasses(py::module &m) {
          Parameters
          ----------
          verbose : bool
-             display verbose information.
+             display verbose information (default=False)
 
          Returns
          -------
          SlabCoefs instance
-         )")
+         )", py::arg("verbose")=false)
     .def("__call__",
 	 &CoefClasses::SlabCoefs::getTensor,
          R"(
@@ -1815,12 +1815,12 @@ void CoefficientClasses(py::module &m) {
          Parameters
          ----------
          verbose : bool
-             display verbose information.
+             display verbose information (default=False)
 
          Returns
          -------
          CubeCoefs instance
-         )")
+         )", py::arg("verbose")=false)
     .def("__call__",
 	 &CoefClasses::CubeCoefs::getTensor,
          R"(
@@ -1908,7 +1908,7 @@ void CoefficientClasses(py::module &m) {
          Parameters
          ----------
          verbose : bool
-             display verbose information.
+             display verbose information (default=True)
 
          Returns
          -------
@@ -1936,7 +1936,7 @@ void CoefficientClasses(py::module &m) {
          type : str
              ascii table data file
          verbose : bool
-             display verbose information.
+             display verbose information (default=True)
 
          Returns
          -------
@@ -1953,7 +1953,7 @@ void CoefficientClasses(py::module &m) {
          array : ndarray
              data columns
          verbose : bool
-             display verbose information.
+             display verbose information (default=True)
 
          Returns
          -------
@@ -1978,7 +1978,7 @@ void CoefficientClasses(py::module &m) {
          Parameters
          ----------
          verbose : bool
-             display verbose information.
+             display verbose information (default=True)
 
          Returns
          -------
@@ -2006,7 +2006,7 @@ void CoefficientClasses(py::module &m) {
          type : str
              ascii table data file
          verbose : bool
-             display verbose information.
+             display verbose information (default=True)
 
          Returns
          -------
@@ -2023,7 +2023,7 @@ void CoefficientClasses(py::module &m) {
          array : ndarray
              data columns
          verbose : bool
-             display verbose information.
+             display verbose information (default=True)
 
          Returns
          -------
diff --git a/pyEXP/ParticleReaderWrappers.cc b/pyEXP/ParticleReaderWrappers.cc
index 2e01b46d5..770d01aeb 100644
--- a/pyEXP/ParticleReaderWrappers.cc
+++ b/pyEXP/ParticleReaderWrappers.cc
@@ -15,11 +15,12 @@ void ParticleReaderClasses(py::module &m) {
     "The available particle readers are:\n"
     "  1. PSPout         The monolithic EXP phase-space snapshot format\n"
     "  2. PSPspl         Like PSPout, but split into multiple file chunks\n"
-    "  3. GadgetNative   The original Gadget native format\n"
-    "  4  GadgetHDF5     The newer HDF5 Gadget format\n"
-    "  5. TipsyNative    The original Tipsy format\n"
-    "  6. TipsyXDR       The original XDR Tipsy format\n"
-    "  7. Bonsai         This is the Bonsai varient of Tipsy files\n\n"
+    "  3. PSPhdf5        The HDF5 version of the EXP phase-space format\n"
+    "  4. GadgetNative   The original Gadget native format\n"
+    "  5  GadgetHDF5     The newer HDF5 Gadget format\n"
+    "  6. TipsyNative    The original Tipsy format\n"
+    "  7. TipsyXDR       The original XDR Tipsy format\n"
+    "  8. Bonsai         This is the Bonsai varient of Tipsy files\n\n"
     "We have a helper function, getReaders, to get a list to help you\n"
     "remember.  Try: pyEXP.read.ParticleReader.getReaders()\n\n"
     "Each reader can manage snapshots split into many files by parallel,\n"
@@ -216,6 +217,38 @@ void ParticleReaderClasses(py::module &m) {
     }
   };
 
+  class PyPSPhdf5 : public PSPhdf5
+  {
+  public:
+
+    // Inherit the constructors
+    using PSPhdf5::PSPhdf5;
+
+    void SelectType(const std::string& type) override {
+      PYBIND11_OVERRIDE(void, PSPhdf5, SelectType, type);
+    }
+    
+    std::vector<std::string> GetTypes() override {
+      PYBIND11_OVERRIDE(std::vector<std::string>, PSPhdf5, GetTypes,);
+    }
+    
+    unsigned long CurrentNumber() override {
+      PYBIND11_OVERRIDE(unsigned long, PSPhdf5, CurrentNumber,);
+    }
+    
+    double CurrentTime() override {
+      PYBIND11_OVERRIDE(double, PSPhdf5, CurrentTime,);
+    }
+    
+    const Particle* firstParticle() override {
+      PYBIND11_OVERRIDE(const Particle*, PSPhdf5, firstParticle,);
+    }
+
+    const Particle* nextParticle() override {
+      PYBIND11_OVERRIDE(const Particle*, PSPhdf5, nextParticle,);
+    }
+  };
+
   class PyTipsy : public Tipsy
   {
   public:
@@ -401,8 +434,8 @@ void ParticleReaderClasses(py::module &m) {
   pr.def_static("getReaders", []()
   {
     const std::vector<std::string> formats = {
-      "PSPout", "PSPspl", "GadgetNative", "GadgetHDF5", "TipsyNative",
-      "TipsyXDR", "Bonsai"};
+      "PSPout", "PSPspl", "PSPhdf5", "GadgetNative", "GadgetHDF5",
+      "TipsyNative", "TipsyXDR", "Bonsai"};
 
     return formats;
     },
@@ -431,7 +464,7 @@ void ParticleReaderClasses(py::module &m) {
          Returns
          -------
          ParticleReader
-         )");
+         )", py::arg("files"), py::arg("verbose")=false);
 
   py::class_<GadgetNative, std::shared_ptr<GadgetNative>, PyGadgetNative, ParticleReader>(m, "GadgetNative")
     .def(py::init<const std::vector<std::string>&, bool>(),
@@ -448,7 +481,7 @@ void ParticleReaderClasses(py::module &m) {
          Returns
          -------
          ParticleReader
-         )");
+         )", py::arg("files"), py::arg("verbose")=false);
 
   py::class_<PSP, std::shared_ptr<PSP>, PyPSP, ParticleReader>(m, "PSP")
     .def(py::init<bool>(), "Base class for PSP reader")
@@ -464,7 +497,20 @@ void ParticleReaderClasses(py::module &m) {
 
   py::class_<PSPout, std::shared_ptr<PSPout>, PyPSPout, PSP>(m, "PSPout")
     .def(py::init<const std::vector<std::string>&, bool>(),
-	 "Reader for monolitic PSP format (single file written by root process)")
+	 R"(
+         Read PSP ascii monolithic format snapshot files
+
+         Parameters
+         ----------
+         files : list(str)
+             List of files with phase-space segments comprising a single snapshot
+         verbose : bool, default=False
+             Verbose, diagnostic output
+
+         Returns
+         -------
+         ParticleReader
+         )", py::arg("files"), py::arg("verbose")=false)
     .def("CurrentNumber",   &PSPout::CurrentNumber)
     .def("GetTypes",        &PSPout::GetTypes)
     .def("CurrentTime",     &PSPout::CurrentTime);
@@ -472,11 +518,46 @@ void ParticleReaderClasses(py::module &m) {
 
   py::class_<PSPspl, std::shared_ptr<PSPspl>, PyPSPspl, PSP>(m, "PSPspl")
     .def(py::init<const std::vector<std::string>&, bool>(),
-	 "Reader for split PSP files (multiple files written by each process)")
+	 R"(
+	 Reader for split PSP files (multiple files written by each process)
+
+         Parameters
+         ----------
+         files : list(str)
+             List of files with phase-space segments comprising a single snapshot
+         verbose : bool, default=False
+             Verbose, diagnostic output
+
+         Returns
+         -------
+         ParticleReader
+         )", py::arg("files"), py::arg("verbose")=false)
     .def("CurrentNumber",   &PSPspl::CurrentNumber)
     .def("GetTypes",        &PSPspl::GetTypes)
     .def("CurrentTime",     &PSPspl::CurrentTime);
 	 
+  py::class_<PSPhdf5, std::shared_ptr<PSPhdf5>, PyPSPhdf5, PSP>(m, "PSPhdf5")
+    .def(py::init<const std::vector<std::string>&, bool>(),
+	 R"(
+	 Reader for HDF5 PSP files (both single and multiple files)
+
+         Parameters
+         ----------
+         files : list(str)
+             List of files with phase-space segments comprising a single snapshot
+         verbose : bool, default=False
+             Verbose, diagnostic output
+
+         Returns
+         -------
+         ParticleReader
+         )", py::arg("files"), py::arg("verbose")=false)
+    .def("SelectType",      &PSPhdf5::SelectType)
+    .def("NumFiles",        &PSPhdf5::NumFiles)
+    .def("CurrentNumber",   &PSPhdf5::CurrentNumber)
+    .def("GetTypes",        &PSPhdf5::GetTypes)
+    .def("CurrentTime",     &PSPhdf5::CurrentTime);
+	 
   py::class_<Tipsy, std::shared_ptr<Tipsy>, PyTipsy, ParticleReader> tipsy(m, "Tipsy");
   
   py::enum_<Tipsy::TipsyType>(tipsy, "TipsyType")
@@ -485,7 +566,24 @@ void ParticleReaderClasses(py::module &m) {
     .value("bonsai", Tipsy::TipsyType::bonsai)
     .export_values();
 
-  tipsy.def(py::init<const std::string&, Tipsy::TipsyType, bool>())
+  tipsy.def(py::init<const std::string&, Tipsy::TipsyType, bool>(),
+	    R"(
+            Read Tipsy format snapshots
+
+            Parameters
+            ----------
+            file : str
+                   The Tipsy snapshot file
+            type : TipsyType
+                   The Tipsy file type (native, xdr, bonsai)
+            verbose : bool, default=False
+                   Verbose, diagnostic output
+
+            Returns
+            -------
+            ParticleReader
+
+            )", py::arg("file"), py::arg("type"), py::arg("verbose")=false)
     .def(py::init<const std::vector<std::string>&, Tipsy::TipsyType, bool>())
     .def("SelectType",      &Tipsy::SelectType)
     .def("CurrentNumber",   &Tipsy::CurrentNumber)
@@ -493,4 +591,3 @@ void ParticleReaderClasses(py::module &m) {
     .def("CurrentTime",     &Tipsy::CurrentTime);
 
 }
-
diff --git a/pyEXP/UtilWrappers.cc b/pyEXP/UtilWrappers.cc
index 2dc43a195..7970fa95f 100644
--- a/pyEXP/UtilWrappers.cc
+++ b/pyEXP/UtilWrappers.cc
@@ -6,6 +6,7 @@
 namespace py = pybind11;
 
 #include "Centering.H"
+#include "KDdensity.H"
 #include "ParticleIterator.H"
 
 void UtilityClasses(py::module &m) {
@@ -243,4 +244,93 @@ void UtilityClasses(py::module &m) {
       None
 
    )", py::arg("verbose") = false);
+
+  py::class_<Utility::KDdensity>(m, "KDdensity")
+    .def(py::init<PR::PRptr, int>(),
+	 R"(
+         Create a k-d tree density estimator for a particle set
+
+         Parameters
+         ----------
+	 reader : ParticleReader
+                  the particle-reader class instance
+         Ndens  : int, default=32
+                  number of particles per sample ball (32 is a good 
+		  compromise between accuracy and runtime; 16 is okay if 
+		  you are trying to shave off runtime. 
+
+         Returns
+	 -------
+	 KDdensity
+	          the KDdensity instance
+    )", py::arg("reader"), py::arg("Ndens")=32)
+    .def(py::init<const Eigen::VectorXd&, const KDdensity::RowMatrixXd&, int>(),
+	 R"(
+         Create a k-d tree density estimator for a particle set
+
+         Parameters
+         ----------
+	 mass   : ndarray(N)
+                  mass array
+	 pos    : ndarray(N,3)
+                  x, y, z array
+         Ndens  : int, default=32
+                  number of particles per sample ball (32 is a good 
+		  compromise between accuracy and runtime; 16 is okay if 
+		  you are trying to shave off runtime. 
+
+         Returns
+	 -------
+	 KDdensity
+	          the KDdensity instance
+    )", py::arg("mass"), py::arg("pos"), py::arg("Ndens")=32)
+    .def("getDensityAtPoint", &Utility::KDdensity::getDensityAtPoint,
+	 R"(
+         Get the density estimate at a given point
+
+         Parameters
+         ----------
+         x : float
+             the x coordinate of the evaluation point
+         y : float
+             the y coordinate of the evaluation point
+         z : float
+             the z coordinate of the evaluation point
+
+         Returns
+         -------
+         float
+                 the density estimate at the particle position
+    )", py::arg("x"), py::arg("y"), py::arg("z"))
+    .def("getDensityAtPoint",
+	 [](Utility::KDdensity& A, std::vector<double>& pos) {
+	   return A.getDensityAtPoint(pos[0], pos[1], pos[2]);
+	 },
+	 R"(
+         Get the density estimate at a given point
+
+         Parameters
+         ----------
+         pos : list(float)
+               x, y, z coordinates of the evaluation point
+
+         Returns
+         -------
+         float
+                 the density estimate at the particle position
+    )", py::arg("pos"))
+    .def("getDensityByIndex", &Utility::KDdensity::getDensityByIndex,
+	 R"(
+         Get the density estimate at a given particle index
+
+         Parameters
+         ----------
+         index : int
+                 the particle index
+
+         Returns
+         -------
+         float
+                 the density estimate at the particle position
+    )", py::arg("index"));
 }
diff --git a/src/CylEXP.cc b/src/CylEXP.cc
index 34034b6c4..5f3067d96 100644
--- a/src/CylEXP.cc
+++ b/src/CylEXP.cc
@@ -278,7 +278,7 @@ void CylEXP::compute_multistep_coefficients(unsigned mlevel)
     }
   }
 
-  coefs_made = vector<short>(multistep+1, true);
+  coefs_made = std::vector<bool>(multistep+1, true);
 }
 
 
diff --git a/utils/ICs/initial.cc b/utils/ICs/initial.cc
index 1366cfda2..fc4f43c40 100644
--- a/utils/ICs/initial.cc
+++ b/utils/ICs/initial.cc
@@ -264,6 +264,7 @@ double       HRATIO  = 1.0;
 double       DWEIGHT = 1.0;
 double       Mfac  = 1.0;
 double       HERNA = 0.10;
+bool         sech2 = true;
 
 #include "Particle.H"
 
@@ -298,8 +299,8 @@ double DiskDens(double R, double z, double phi)
     {
       double a1 = ASCALE;
       double a2 = ASCALE*ARATIO;
-      double h1 = HSCALE;
-      double h2 = HSCALE*HRATIO;
+      double h1 = sech2 ? 0.5*HSCALE : HSCALE;
+      double h2 = sech2 ? 0.5*HSCALE*HRATIO : HSCALE*HRATIO;
       double w1 = 1.0/(1.0+DWEIGHT);
       double w2 = DWEIGHT/(1.0+DWEIGHT);
 
@@ -315,7 +316,7 @@ double DiskDens(double R, double z, double phi)
   case DiskType::diskbulge:
     {
       double acyl = ASCALE;
-      double hcyl = HSCALE;
+      double hcyl = sech2 ? 0.5*HSCALE : HSCALE;
       double f = cosh(z/HSCALE);
       double rr = pow(pow(R, 2) + pow(z,2), 0.5);
       double w1 = Mfac;
@@ -335,8 +336,9 @@ double DiskDens(double R, double z, double phi)
   case DiskType::exponential:
   default:
     {
-      double f = cosh(z/HSCALE);
-      ans = exp(-R/ASCALE)/(4.0*M_PI*ASCALE*ASCALE*HSCALE*f*f);
+      double h = sech2 ? 0.5*HSCALE : HSCALE;
+      double f = cosh(z/h);
+      ans = exp(-R/ASCALE)/(4.0*M_PI*ASCALE*ASCALE*h*f*f);
     }
     break;
   }
@@ -399,7 +401,7 @@ main(int ac, char **av)
   bool         const_height, images, multi, basis, zeropos, zerovel;
   bool         report, ignore, evolved, diskmodel;
   int          nhalo, ndisk, ngas, ngparam;
-  std::string  hbods, dbods, gbods, suffix, centerfile, halofile1, halofile2;
+  std::string  hbods, dbods, gbods, outtag, runtag, centerfile, halofile1, halofile2;
   std::string  cachefile, config, gentype, dtype, dmodel, mtype, ctype;
 
   const std::string mesg("Generates a Monte Carlo realization of a halo with an\n embedded disk using Jeans' equations\n");
@@ -409,6 +411,12 @@ main(int ac, char **av)
   options.add_options()
     ("h,help", "Print this help message")
     ("T,template", "Write template options file with current and all default values")
+    ("sech2", "Use sech^2 vertical profile for disk density",
+     cxxopts::value<bool>()->default_value("true"))
+    ("mtype", "EmpCylSL spherical model type (one of: Exponential, Gaussian, Plummer, Power)",
+     cxxopts::value<string>(mtype)->default_value("Exponential"))
+    ("ppow", "Power-law index for EmpCylSL Power spherical model",
+     cxxopts::value<double>(PPower)->default_value("2.0"))
     ("c,config", "Parameter configuration file",
      cxxopts::value<string>(config))
     ("deproject", "The EmpCylSL deprojection from specified disk model (EXP or MN)",
@@ -599,10 +607,10 @@ main(int ac, char **av)
      cxxopts::value<double>(Tmin)->default_value("500.0"))
     ("centerfile", "File containing phase-space center",
      cxxopts::value<std::string>(centerfile))
-    ("runtag", "Prefix for output files",
-     cxxopts::value<std::string>(runtag)->default_value("gendisk"))
-    ("suffix", "Suffix for output files (none by default)",
-     cxxopts::value<std::string>(suffix))
+    ("suffix", "Prefix for output files",
+     cxxopts::value<std::string>(outtag)->default_value("gendisk"))
+    ("runtag", "Suffix for output files (none by default)",
+     cxxopts::value<std::string>(runtag))
     ("threads", "Number of threads to run",
      cxxopts::value<int>(nthrds)->default_value("1"))
     ("allow", "Allow multimass algorithm to generature negative masses for testing")
@@ -799,10 +807,10 @@ main(int ac, char **av)
 
   int n_particlesH, n_particlesD, n_particlesG;
 
-  if (suffix.size()>0) {
-    hbods = hbods + "." + suffix;
-    dbods = dbods + "." + suffix;
-    gbods = gbods + "." + suffix;
+  if (runtag.size()>0) {
+    hbods = hbods + "." + runtag;
+    dbods = dbods + "." + runtag;
+    gbods = gbods + "." + runtag;
   }
 
   // Divvy up the particles by core
@@ -866,7 +874,7 @@ main(int ac, char **av)
   if (vm.count("itmax"))  DiskHalo::ITMAX    = itmax;
   if (vm.count("allow"))  DiskHalo::ALLOW    = true;
   if (vm.count("nomono")) DiskHalo::use_mono = false;
-  if (suffix.size())      DiskHalo::RUNTAG   = suffix;
+  if (runtag.size())      DiskHalo::RUNTAG   = runtag;
 
   AddDisk::use_mpi      = true;
   AddDisk::Rmin         = RMIN;
@@ -978,7 +986,7 @@ main(int ac, char **av)
     // Basis orthgonality check
     //
     if (vm.count("ortho")) {
-      std::ofstream out(runtag + ".ortho_check");
+      std::ofstream out(outtag + ".ortho_check");
       expandd->ortho_check(out);
     }
   }
@@ -1179,7 +1187,7 @@ main(int ac, char **av)
 	std::cout << "Dumping coefficients halo . . . " << std::flush;
 	ostringstream sout;
 	sout << "halo_coefs.";
-	if (suffix.size()>0) sout << suffix;
+	if (outtag.size()>0) sout << outtag;
 	else                 sout << "dump";
 	ofstream out(sout.str());
 	if (out) expandh->dump_coefs(out, false);
@@ -1189,7 +1197,7 @@ main(int ac, char **av)
 	std::cout << "Dumping a probe through the halo . . . " << std::flush;
 	ostringstream sout;
 	sout << "halo_probe.";
-	if (suffix.size()>0) sout << suffix;
+	if (outtag.size()>0) sout << outtag;
 	else                 sout << "dump";
 	ofstream out(sout.str());
 	if (out) {
@@ -1279,7 +1287,7 @@ main(int ac, char **av)
 	std::cout << "Dumping coefficients . . . " << std::flush;
 	ostringstream sout;
 	sout << "disk_coefs.";
-	if (suffix.size()>0) sout << suffix;
+	if (outtag.size()>0) sout << outtag;
 	else                 sout << "dump";
 	ofstream out(sout.str());
 	if (out) expandd->dump_coefs(out);
@@ -1325,11 +1333,11 @@ main(int ac, char **av)
 
     if (ndisk) {
       int nout = 200;
-      string dumpstr = runtag + ".dump";
+      string dumpstr = outtag + ".dump";
       expandd->dump_basis(dumpstr, 0);
-      expandd->dump_images(runtag, 5.0*scale_length, 5.0*scale_height,
+      expandd->dump_images(outtag, 5.0*scale_length, 5.0*scale_height,
 			   nout, nout, false);
-      expandd->dump_images_basis(runtag, 5.0*scale_length, 5.0*scale_height,
+      expandd->dump_images_basis(outtag, 5.0*scale_length, 5.0*scale_height,
 				 nout, nout, false, 0, MMAX, 0, NMAXD-1);
     }