From b307d7a51c6a69a71221863f147c1d2d457d6662 Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Wed, 8 May 2024 14:19:43 -0400
Subject: [PATCH 01/25] Updated pybind11, HighFive, and yaml-cpp to latest
 upstream versions; added a HighFive template specification for the new API
 [no ci]

---
 exputil/SLGridMP2.cc | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/exputil/SLGridMP2.cc b/exputil/SLGridMP2.cc
index c74f84115..ba0819260 100644
--- a/exputil/SLGridMP2.cc
+++ b/exputil/SLGridMP2.cc
@@ -541,8 +541,8 @@ bool SLGridSph::ReadH5Cache(void)
       
       // Table arrays will be allocated
       //
-      arrays.getDataSet("ev").read(table[l].ev);
-      arrays.getDataSet("ef").read(table[l].ef);
+      arrays.getDataSet("ev").read<Eigen::VectorXd>(table[l].ev);
+      arrays.getDataSet("ef").read<Eigen::MatrixXd>(table[l].ef);
     }
     
     if (myid==0)

From 03300c68de6ed6607c1a6d286011139e9380d193 Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Thu, 9 May 2024 12:09:32 -0400
Subject: [PATCH 02/25] Added a 'look alive' time-step marker for interactive
 use

---
 src/step.cc | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/src/step.cc b/src/step.cc
index e17db5cd1..c5d562f1d 100644
--- a/src/step.cc
+++ b/src/step.cc
@@ -340,6 +340,9 @@ void do_step(int n)
 				// Stop the total step timer
   if (step_timing) timer_tot.stop();
 
+  if (VERBOSE==2 and myid==0) 	// Time step marker
+    std::cout << << std::endl << ">>>" << this_step << "<<<" << std::endl;
+
 				// Timer output
   if (step_timing && this_step!=0 && (this_step % tskip) == 0) {
     if (myid==0) {

From 4bd74a02c231770aad6696c358e41a935b2139c1 Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Thu, 9 May 2024 13:24:42 -0400
Subject: [PATCH 03/25] Fixes for FieldGenerator (eval functions not following
 the correct field-ordering convention) [no ci]

---
 expui/BiorthBasis.cc | 57 ++++++++++++++++++++++----------------------
 src/step.cc          |  2 +-
 2 files changed, 29 insertions(+), 30 deletions(-)

diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc
index b13e39e56..68897a4cc 100644
--- a/expui/BiorthBasis.cc
+++ b/expui/BiorthBasis.cc
@@ -531,15 +531,15 @@ namespace BasisClasses
     double potlfac = 1.0/scale;
     
     return
-      {den0 * densfac,
-       den1 * densfac,
-       (den0 + den1) * densfac,
-       pot0 * potlfac,
-       pot1 * potlfac,
-       (pot0 + pot1) * densfac,
-       potr * (-potlfac)/scale,
-       pott * (-potlfac),
-       potp * (-potlfac)};
+      {den0 * densfac,		// 0
+       den1 * densfac,		// 1
+       (den0 + den1) * densfac,	// 2
+       pot0 * potlfac,		// 3
+       pot1 * potlfac,		// 4
+       (pot0 + pot1) * densfac,	// 5
+       potr * (-potlfac)/scale,	// 6
+       pott * (-potlfac),	// 7
+       potp * (-potlfac)};	// 8
     //         ^
     //         |
     // Return force not potential gradient
@@ -554,10 +554,10 @@ namespace BasisClasses
     
     auto v = sph_eval(r, costh, phi);
     
-    double potR = v[4]*sinth + v[5]*costh;
-    double potz = v[4]*costh - v[5]*sinth;
+    double potR = v[6]*sinth + v[7]*costh;
+    double potz = v[6]*costh - v[7]*sinth;
 
-    return {v[0], v[1], v[2], v[3], potR, potz, v[6]};
+    return {v[0], v[1], v[2], v[3], v[4], v[5], potR, potz, v[8]};
   }
 
 
@@ -571,12 +571,10 @@ namespace BasisClasses
 
     auto v = cyl_eval(R, z, phi);
 
-    //tdens0, tdens, tpotl0, tpotl, tpotR, tpotz, tpotp
+    double tpotx = v[6]*x/R - v[8]*y/R ;
+    double tpoty = v[6]*y/R + v[8]*x/R ;
     
-    double tpotx = v[4]*x/R - v[6]*y/R ;
-    double tpoty = v[4]*y/R + v[6]*x/R ;
-    
-    return {v[0], v[1], v[2], v[3], tpotx, tpoty, v[5]};
+    return {v[0], v[1], v[2], v[3], v[4], v[5], tpotx, tpoty, v[7]};
   }
   
 
@@ -1277,14 +1275,15 @@ namespace BasisClasses
 
     if (midplane) {
       height = sl->accumulated_midplane_eval(R, -colh*hcyl, colh*hcyl, phi);
-
       return
 	{tdens0, tdens - tdens0, tdens,
-	 tpotl0, tpotl - tpotl0, tpotl, tpotR, tpotz, tpotp, height};
+	 tpotl0, tpotl - tpotl0, tpotl,
+	 tpotR, tpotz, tpotp, height};
     } else {
       return
 	{tdens0, tdens - tdens0, tdens,
-	 tpotl0, tpotl - tpotl0, tpotl, tpotR, tpotz, tpotp};
+	 tpotl0, tpotl - tpotl0, tpotl,
+	 tpotR, tpotz, tpotp};
     }
   }
   
@@ -1741,7 +1740,7 @@ namespace BasisClasses
       rpot = -totalMass*R/(r*r2 + 10.0*std::numeric_limits<double>::min());
       zpot = -totalMass*z/(r*r2 + 10.0*std::numeric_limits<double>::min());
       
-      return {den0, den1, pot0, pot1, rpot, zpot, ppot};
+      return {den0, den1, den0+den1, pot0, pot1, pot0+pot1, rpot, zpot, ppot};
     }
 
     // Get the basis fields
@@ -1827,16 +1826,16 @@ namespace BasisClasses
     // Cylindrical coords
     //
     double sinth = sqrt(fabs(1.0 - costh*costh));
-    double R = r*sinth, z = r*costh, potR, potz;
+    double R = r*sinth, z = r*costh;
 
     auto v = cyl_eval(R, z, phi);
   
     // Spherical force element converstion
     //
-    double potr = potR*sinth + potz*costh;
-    double pott = potR*costh - potz*sinth;
+    double potr = v[6]*sinth + v[7]*costh;
+    double pott = v[6]*costh - v[7]*sinth;
 
-    return {v[0], v[1], v[2], v[3], potr, pott, v[6]};
+    return {v[0], v[1], v[2], v[3], v[4], v[5], potr, pott, v[8]};
   }
 
   std::vector<double> FlatDisk::crt_eval(double x, double y, double z)
@@ -1848,10 +1847,10 @@ namespace BasisClasses
 
     auto v = cyl_eval(R, z, phi);
 
-    double potx = v[4]*x/R - v[6]*y/R;
-    double poty = v[4]*y/R + v[6]*x/R;
+    double potx = v[4]*x/R - v[8]*y/R;
+    double poty = v[4]*y/R + v[8]*x/R;
 
-    return {v[0], v[1], v[2], v[3], potx, poty, v[5]};
+    return {v[0], v[1], v[2], v[3], v[4], v[5], potx, poty, v[7]};
   }
 
   std::vector<Eigen::MatrixXd> FlatDisk::orthoCheck()
@@ -2582,7 +2581,7 @@ namespace BasisClasses
     double frcy = -frc(1).real();
     double frcz = -frc(2).real();
 
-    return {0, den1, 0, pot1, frcx, frcy, frcz};
+    return {0, den1, den1, 0, pot1, pot1, frcx, frcy, frcz};
   }
 
   std::vector<double> Cube::cyl_eval(double R, double z, double phi)
diff --git a/src/step.cc b/src/step.cc
index c5d562f1d..a057ec247 100644
--- a/src/step.cc
+++ b/src/step.cc
@@ -341,7 +341,7 @@ void do_step(int n)
   if (step_timing) timer_tot.stop();
 
   if (VERBOSE==2 and myid==0) 	// Time step marker
-    std::cout << << std::endl << ">>>" << this_step << "<<<" << std::endl;
+    std::cout << std::endl << ">>>" << this_step << "<<<" << std::endl;
 
 				// Timer output
   if (step_timing && this_step!=0 && (this_step % tskip) == 0) {

From 3ec4c56c83d673d019aa18dbde6f3485d11e49af Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Fri, 10 May 2024 17:54:18 -0400
Subject: [PATCH 04/25] Generalized KD for carrying velocity data [no ci]

---
 include/KDtree.H | 66 ++++++++++++++++++++++++++++++++++++++++++------
 1 file changed, 58 insertions(+), 8 deletions(-)

diff --git a/include/KDtree.H b/include/KDtree.H
index 57a35bc00..40a845b32 100644
--- a/include/KDtree.H
+++ b/include/KDtree.H
@@ -40,20 +40,39 @@ public:
   point() {}
   
   //! Constructor
-  point(std::array<coordinate_type, dimensions> c, double m=1) : coords_(c), mass_(m) {}
+  point(std::array<coordinate_type, dimensions> c, double m=1,
+	unsigned long indx=0) : coords_(c), mass_(m), indx_(indx) {}
+  
+  point(std::array<coordinate_type, dimensions> c,
+	std::array<coordinate_type, dimensions> v,
+	double m=1, unsigned long indx=0) : coords_(c), vels_(v),
+					    mass_(m), indx_(indx) {}
   
   //! List constructor
-  point(std::initializer_list<coordinate_type> list, double m=1) : mass_(m)
+  point(std::initializer_list<coordinate_type> list, double m=1,
+	unsigned long indx=0) : mass_(m), indx_(indx)
+  {
+    size_t n = std::min(dimensions, list.size());
+    std::copy_n(list.begin(), n, coords_.begin());
+  }
+
+  point(std::initializer_list<coordinate_type> list,
+	std::initializer_list<coordinate_type> vlst,
+	double m=1,
+	unsigned long indx=0) : mass_(m), indx_(indx)
   {
     size_t n = std::min(dimensions, list.size());
     std::copy_n(list.begin(), n, coords_.begin());
+    std::copy_n(vlst.begin(), n, vels_.begin());
   }
 
   //! Copy constructor
   point(const point& p)
   {
     for (size_t n=0; n<dimensions; n++) coords_[n] = p.coords_[n];
+    for (size_t n=0; n<dimensions; n++) vels_  [n] = p.vels_  [n];
     mass_ = p.mass_;
+    indx_ = p.indx_;
   }
 
   /**
@@ -67,6 +86,11 @@ public:
     return coords_[index];
   }
   
+  coordinate_type vel(size_t index) const
+  {
+    return vels_[index];
+  }
+  
   /**
    * Returns the distance squared from this point to another
    * point.
@@ -85,14 +109,39 @@ public:
     return dist;
   }
 
+  /**
+   * Returns the squared velocity between this point to another
+   * point.
+   *
+   * @param pt another point
+   * @return squared speed between this point and the other point
+   */
+  double speed(const point& pt) const
+  {
+    double rvel = 0;
+    for (size_t i = 0; i < dimensions; ++i)
+      {
+	double d = vel(i) - pt.vel(i);
+	rvel += d * d;
+      }
+    return rvel;
+  }
+
   double mass() const
   {
     return mass_;
   }
 
+  unsigned long indx() const
+  {
+    return indx_;
+  }
+
 private:
   std::array<coordinate_type, dimensions> coords_;
+  std::array<coordinate_type, dimensions> vels_;
   double mass_;
+  unsigned long indx_;
 };
 
 //! For iostream printing of points
@@ -163,7 +212,7 @@ private:
     size_t n = begin + (end - begin)/2;
     std::nth_element(&nodes_[begin], &nodes_[n], &nodes_[end], node_cmp(index));
     index = (index + 1) % dimensions;
-    nodes_[n].left_ = make_tree(begin, n, index);
+    nodes_[n].left_  = make_tree(begin, n, index);
     nodes_[n].right_ = make_tree(n + 1, end, index);
     return &nodes_[n];
   }
@@ -175,19 +224,20 @@ private:
     ++visited_;
 
     double d = root->distance(point);
-    if (best_.size()==0 || d < best_.rbegin()->first) {
+    if (best_.size()<N || d < best_.rbegin()->first) {
       best_.add(d, root);
     }
 
-    // This is only correct is the test point is never in the data set . . .
+    // This is only correct if the test point is never in the data set . . .
     // if (best_.begin()->first == 0) return;
 
     double dx = root->get(index) - point.get(index);
     index = (index + 1) % dimensions;
-    nearestN(dx > 0 ? root->left_ : root->right_, point, index, N);
 
-    if (dx * dx >= best_.rbegin()->first) return;
-    nearestN(dx > 0 ? root->right_ : root->left_, point, index, N);
+    nearestN(dx > 0 ? root->left_  : root->right_, point, index, N);
+
+    if (best_.size()>=N and dx * dx >= best_.rbegin()->first) return;
+    nearestN(dx > 0 ? root->right_ : root->left_,  point, index, N);
   }
   
 public:

From 19ae54824e7d254507de2a42610bd7c3ca436e9d Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Sat, 11 May 2024 12:09:40 -0400
Subject: [PATCH 05/25] Add basis cache version - Create version string
 variable for all HighFive-generated HDF5 caches - This will effectively
 expire all existing caches - Emit a message that the cache is rebuiling for
 an API change

---
 exputil/BiorthCyl.cc | 19 +++++++++++++++++++
 exputil/EmpCyl2d.cc  | 41 ++++++++++++++++++++++++++++-------------
 exputil/EmpCylSL.cc  | 15 ++++++++++++++-
 exputil/SLGridMP2.cc | 13 +++++++++++++
 include/BiorthCyl.H  |  3 +++
 include/EmpCyl2d.H   |  3 +++
 include/EmpCylSL.H   |  3 +++
 include/SLGridMP2.H  |  3 +++
 8 files changed, 86 insertions(+), 14 deletions(-)

diff --git a/exputil/BiorthCyl.cc b/exputil/BiorthCyl.cc
index 8f87e63a3..3509cc50b 100644
--- a/exputil/BiorthCyl.cc
+++ b/exputil/BiorthCyl.cc
@@ -27,6 +27,9 @@
 #include <libvars.H>
 using namespace __EXP__;	// For reference to n-body globals
 
+// Cache version
+std::string BiorthCyl::Version = "1.0";
+
 // Constructor
 BiorthCyl::BiorthCyl(const YAML::Node& conf) : conf(conf)
 {
@@ -547,6 +550,10 @@ void BiorthCyl::WriteH5Cache()
     //
     file.createAttribute<std::string>("forceID", HighFive::DataSpace::From(forceID)).write(forceID);
       
+    // Cache version
+    //
+    file.createAttribute<std::string>("Version",  HighFive::DataSpace::From(Version)).write(Version);
+
     // Stash the basis configuration (this is not yet implemented in EXP)
     //
     std::ostringstream sout; sout << conf;
@@ -631,6 +638,18 @@ bool BiorthCyl::ReadH5Cache()
     if (not checkStr(geometry, "geometry"))  return false;
     if (not checkStr(forceID,  "forceID"))   return false;
 
+    // Version check
+    //
+    if (h5file.hasAttribute("Version")) {
+      if (not checkStr(Version, "Version"))  return false;
+    } else {
+      if (myid==0)
+	std::cout << "---- BiorthCyl::ReadH5Cache: "
+		  << "recomputing cache for HighFive API change"
+		  << std::endl;
+      return false;
+    }
+
     // Parameter check
     //
     if (not checkInt(mmax,     "mmax"))      return false;
diff --git a/exputil/EmpCyl2d.cc b/exputil/EmpCyl2d.cc
index fd3408eaa..9164af353 100644
--- a/exputil/EmpCyl2d.cc
+++ b/exputil/EmpCyl2d.cc
@@ -997,6 +997,8 @@ void EmpCyl2d::WriteH5Cache()
     //
     HighFive::File file(cache_name_2d, HighFive::File::Overwrite);
     
+
+
     // Workaround for lack of HighFive boolean support
     int ilogr = 0, icmap = 0;
     if (logr) ilogr = 1;
@@ -1008,19 +1010,20 @@ void EmpCyl2d::WriteH5Cache()
 
     // Parameters
     //
-    file.createAttribute<int>        ("mmax",   HighFive::DataSpace::From(mmax)).  write(mmax);
-    file.createAttribute<int>        ("nmaxfid",   HighFive::DataSpace::From(nmaxfid)).  write(nmaxfid);
-    file.createAttribute<int>        ("nmax", HighFive::DataSpace::From(nmax)).write(nmax);
-    file.createAttribute<int>        ("numr",   HighFive::DataSpace::From(numr)).  write(numr);
-    file.createAttribute<int>        ("knots",  HighFive::DataSpace::From(knots)). write(knots);
-    file.createAttribute<int>        ("ilogr",  HighFive::DataSpace::From(ilogr)). write(ilogr);
-    file.createAttribute<int>        ("icmap",  HighFive::DataSpace::From(icmap)). write(icmap);
-    file.createAttribute<double>     ("rmin",   HighFive::DataSpace::From(rmin)).  write(rmin);
-    file.createAttribute<double>     ("rmax",   HighFive::DataSpace::From(rmax)).  write(rmax);
-    file.createAttribute<double>     ("scale",  HighFive::DataSpace::From(scale)). write(scale);
-    file.createAttribute<std::string>("params", HighFive::DataSpace::From(params)).write(params);
-    file.createAttribute<std::string>("model",  HighFive::DataSpace::From(model)). write(model);
-    file.createAttribute<std::string>("biorth", HighFive::DataSpace::From(biorth)).write(biorth);
+    file.createAttribute<std::string>("Version",   HighFive::DataSpace::From(Version)).write(Version);
+    file.createAttribute<int>        ("mmax",      HighFive::DataSpace::From(mmax)).   write(mmax);
+    file.createAttribute<int>        ("nmaxfid",   HighFive::DataSpace::From(nmaxfid)).write(nmaxfid);
+    file.createAttribute<int>        ("nmax",      HighFive::DataSpace::From(nmax)).   write(nmax);
+    file.createAttribute<int>        ("numr",      HighFive::DataSpace::From(numr)).   write(numr);
+    file.createAttribute<int>        ("knots",     HighFive::DataSpace::From(knots)).  write(knots);
+    file.createAttribute<int>        ("ilogr",     HighFive::DataSpace::From(ilogr)).  write(ilogr);
+    file.createAttribute<int>        ("icmap",     HighFive::DataSpace::From(icmap)).  write(icmap);
+    file.createAttribute<double>     ("rmin",      HighFive::DataSpace::From(rmin)).   write(rmin);
+    file.createAttribute<double>     ("rmax",      HighFive::DataSpace::From(rmax)).   write(rmax);
+    file.createAttribute<double>     ("scale",     HighFive::DataSpace::From(scale)).  write(scale);
+    file.createAttribute<std::string>("params",    HighFive::DataSpace::From(params)). write(params);
+    file.createAttribute<std::string>("model",     HighFive::DataSpace::From(model)).  write(model);
+    file.createAttribute<std::string>("biorth",    HighFive::DataSpace::From(biorth)). write(biorth);
       
     // Arrays
     //
@@ -1080,6 +1083,18 @@ bool EmpCyl2d::ReadH5Cache()
 
     //
 
+    // Version check
+    //
+    if (file.hasAttribute("Version")) {
+      if (not checkStr(Version, "Version"))  return false;
+    } else {
+      if (myid==0)
+	std::cout << "---- EmpCyl2d::ReadH5Cache: "
+		  << "recomputing cache for HighFive API change"
+		  << std::endl;
+      return false;
+    }
+
     // Serialize the config and make a string for checking
     YAML::Emitter y; y << Params;
     std::string params(y.c_str());
diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc
index 170ec88fa..9dc4860cd 100644
--- a/exputil/EmpCylSL.cc
+++ b/exputil/EmpCylSL.cc
@@ -7040,7 +7040,8 @@ void EmpCylSL::WriteH5Cache()
     std::string model = EmpModelLabs[mtype];
 
     file.createAttribute<std::string>("geometry", HighFive::DataSpace::From(geometry)).write(geometry);
-    file.createAttribute<std::string>("forceID", HighFive::DataSpace::From(forceID)).write(forceID);
+    file.createAttribute<std::string>("forceID",  HighFive::DataSpace::From(forceID)).write(forceID);
+    file.createAttribute<std::string>("Version",  HighFive::DataSpace::From(Version)).write(Version);
       
     // Write the specific parameters
     //
@@ -7161,6 +7162,18 @@ bool EmpCylSL::ReadH5Cache()
       return false;
     };
 
+    // Version check
+    //
+    if (file.hasAttribute("Version")) {
+      if (not checkStr(Version, "Version"))  return false;
+    } else {
+      if (myid==0)
+	std::cout << "---- EmpCylSL::ReadH5Cache: "
+		  << "recomputing cache for HighFive API change"
+		  << std::endl;
+      return false;
+    }
+
     if (not checkStr(geometry, "geometry"))  return false;
     if (not checkStr(forceID,  "forceID"))   return false;
     if (not checkInt(MMAX,     "mmax"))      return false;
diff --git a/exputil/SLGridMP2.cc b/exputil/SLGridMP2.cc
index ba0819260..7db5367d7 100644
--- a/exputil/SLGridMP2.cc
+++ b/exputil/SLGridMP2.cc
@@ -506,6 +506,18 @@ bool SLGridSph::ReadH5Cache(void)
     if (not checkStr(geometry, "geometry"))  return false;
     if (not checkStr(forceID,  "forceID"))   return false;
 
+    // Version check
+    //
+    if (h5file.hasAttribute("Version")) {
+      if (not checkStr(Version, "Version"))  return false;
+    } else {
+      if (myid==0)
+	std::cout << "---- SLGridSph::ReadH5Cache: "
+		  << "recomputing cache for HighFive API change"
+		  << std::endl;
+      return false;
+    }
+
     // Parameter check
     //
     if (not checkStr(modl,     "model"))     return false;
@@ -604,6 +616,7 @@ void SLGridSph::WriteH5Cache(void)
 
     file.createAttribute<std::string>("geometry",  HighFive::DataSpace::From(geometry)).write(geometry);
     file.createAttribute<std::string>("forceID",   HighFive::DataSpace::From(forceID)).write(forceID);
+    file.createAttribute<std::string>("Version",   HighFive::DataSpace::From(Version)).write(Version);
       
     // Write parameters
     //
diff --git a/include/BiorthCyl.H b/include/BiorthCyl.H
index 5b3f775b3..b052c7979 100644
--- a/include/BiorthCyl.H
+++ b/include/BiorthCyl.H
@@ -100,6 +100,9 @@ protected:
   //! Read the HDF5 cache
   virtual bool ReadH5Cache();
 
+  //! Cache versioning
+  static std::string Version;
+
 public:
 
   //! Flag for MPI enabled (default: 0=off)
diff --git a/include/EmpCyl2d.H b/include/EmpCyl2d.H
index e062644b3..440ca3667 100644
--- a/include/EmpCyl2d.H
+++ b/include/EmpCyl2d.H
@@ -127,6 +127,9 @@ protected:
   bool ReadH5Cache();
   void WriteH5Cache();
 
+  //! Cache versioning
+  inline static const std::string Version = "1.0";
+
   //! Basis magic number
   inline static const unsigned int hmagic = 0xc0a57a1;
   
diff --git a/include/EmpCylSL.H b/include/EmpCylSL.H
index 0f36e026e..30a4a6587 100644
--- a/include/EmpCylSL.H
+++ b/include/EmpCylSL.H
@@ -231,6 +231,9 @@ protected:
   //! Read HDF5 cache
   bool ReadH5Cache();
 
+  //! Cache versioning
+  inline static const std::string Version = "1.0";
+
   //! The cache file name
   std::string cachefile;
 				// 1=write, 0=read
diff --git a/include/SLGridMP2.H b/include/SLGridMP2.H
index 7581c56ff..95a9bf6b7 100644
--- a/include/SLGridMP2.H
+++ b/include/SLGridMP2.H
@@ -91,6 +91,9 @@ private:
   //! Read HDF5 cache
   bool ReadH5Cache();
 
+  //! Cache versioning
+  inline static const std::string Version = "1.0";
+
 public:
 
   //! Flag for MPI enabled (default: 0=off)

From f7ba79b86028a1287bb64e372ae929785c7961a0 Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Sat, 11 May 2024 16:24:06 -0400
Subject: [PATCH 06/25] Missing submodule updates; missing API updates as a
 consequence

---
 expui/Coefficients.cc | 10 +++-------
 expui/expMSSA.cc      |  8 ++++----
 exputil/BiorthCyl.cc  |  2 +-
 exputil/EmpCyl2d.cc   |  8 ++++----
 exputil/EmpCylSL.cc   | 23 ++++++++++-------------
 exputil/SLGridMP2.cc  |  6 ++----
 extern/HighFive       |  2 +-
 extern/yaml-cpp       |  2 +-
 include/EmpCyl2d.H    |  6 ++----
 9 files changed, 28 insertions(+), 39 deletions(-)

diff --git a/expui/Coefficients.cc b/expui/Coefficients.cc
index 452ebecba..5385be271 100644
--- a/expui/Coefficients.cc
+++ b/expui/Coefficients.cc
@@ -116,9 +116,7 @@ namespace CoefClasses
 
       if (Time < Tmin or Time > Tmax) continue;
 
-      int ldim = (Lmax+1)*(Lmax+2)/2;
-      Eigen::MatrixXcd in(ldim, Nmax);
-      stanza.getDataSet("coefficients").read(in);
+      auto in = stanza.getDataSet("coefficients").read<Eigen::MatrixXcd>();
       
       // Pack the data into the coefficient variable
       //
@@ -286,8 +284,7 @@ namespace CoefClasses
       std::array<long int, 3> shape;
       stanza.getAttribute("shape").read(shape);
 
-      Eigen::VectorXcd in(shape[0]*shape[1]*shape[2]);
-      stanza.getDataSet("coefficients").read(in);
+      auto in = stanza.getDataSet("coefficients").read<Eigen::VectorXcd>();
       
       // Pack the data into the coefficient variable
       //
@@ -378,8 +375,7 @@ namespace CoefClasses
       std::array<long int, 3> shape;
       stanza.getAttribute("shape").read(shape);
 
-      Eigen::VectorXcd in(shape[0]*shape[1]*shape[2]);
-      stanza.getDataSet("coefficients").read(in);
+      auto in = stanza.getDataSet("coefficients").read<Eigen::VectorXcd>();
       
       // Pack the data into the coefficient variable
       //
diff --git a/expui/expMSSA.cc b/expui/expMSSA.cc
index 342f6d110..713d34ecb 100644
--- a/expui/expMSSA.cc
+++ b/expui/expMSSA.cc
@@ -1513,10 +1513,10 @@ namespace MSSA {
 
       auto analysis = h5file.getGroup("mssa_analysis");
 
-      analysis.getDataSet("Y" ).read(Y );
-      analysis.getDataSet("S" ).read(S );
-      analysis.getDataSet("U" ).read(U );
-      analysis.getDataSet("PC").read(PC);
+      Y  = analysis.getDataSet("Y" ).read<Eigen::MatrixXd>();
+      S  = analysis.getDataSet("S" ).read<Eigen::VectorXd>();
+      U  = analysis.getDataSet("U" ).read<Eigen::MatrixXd>();
+      PC = analysis.getDataSet("PC").read<Eigen::MatrixXd>();
 
       numK = numT - numW + 1;	// Recompute numK, needed for
 				// reconstruction
diff --git a/exputil/BiorthCyl.cc b/exputil/BiorthCyl.cc
index 3509cc50b..cd4777cb8 100644
--- a/exputil/BiorthCyl.cc
+++ b/exputil/BiorthCyl.cc
@@ -504,7 +504,7 @@ void BiorthCyl::WriteH5Arrays(HighFive::Group& harmonic)
       sout << n;
       auto arrays = order.createGroup(sout.str());
 
-      HighFive::DataSet ds1 = arrays.createDataSet("density",   dens  [m][n]);
+      HighFive::DataSet ds1 = arrays.createDataSet<Eigen::MatrixXd>("density",   dens  [m][n]);
       HighFive::DataSet ds2 = arrays.createDataSet("potential", pot   [m][n]);
       HighFive::DataSet ds3 = arrays.createDataSet("rforce",    rforce[m][n]);
       HighFive::DataSet ds4 = arrays.createDataSet("zforce",    zforce[m][n]);
diff --git a/exputil/EmpCyl2d.cc b/exputil/EmpCyl2d.cc
index 9164af353..af0901e4e 100644
--- a/exputil/EmpCyl2d.cc
+++ b/exputil/EmpCyl2d.cc
@@ -1146,10 +1146,10 @@ bool EmpCyl2d::ReadH5Cache()
       sout << m;
       auto harmonic = file.getGroup(sout.str());
       
-      harmonic.getDataSet("potl").read(potl_array[m]);
-      harmonic.getDataSet("dens").read(dens_array[m]);
-      harmonic.getDataSet("dpot").read(dpot_array[m]);
-      harmonic.getDataSet("rot" ).read(rot_matrix[m]);
+      potl_array[m] = harmonic.getDataSet("potl").read<Eigen::MatrixXd>();
+      dens_array[m] = harmonic.getDataSet("dens").read<Eigen::MatrixXd>();
+      dpot_array[m] = harmonic.getDataSet("dpot").read<Eigen::MatrixXd>();
+      rot_matrix[m] = harmonic.getDataSet("rot" ).read<Eigen::MatrixXd>();
     }
 
   } catch (HighFive::Exception& err) {
diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc
index 9dc4860cd..23d579c0f 100644
--- a/exputil/EmpCylSL.cc
+++ b/exputil/EmpCylSL.cc
@@ -26,11 +26,8 @@
 
 // For reading and writing HDF5 cache files
 //
-#include <highfive/H5File.hpp>
-#include <highfive/H5DataSet.hpp>
-#include <highfive/H5DataSpace.hpp>
-#include <highfive/H5Attribute.hpp>
-
+#include <highfive/highfive.hpp>
+#include <highfive/eigen.hpp>
 
 #ifdef HAVE_OMP_H
 #include <omp.h>		// For multithreading basis construction
@@ -7244,10 +7241,10 @@ bool EmpCylSL::ReadH5Cache()
 	sout << n;
 	auto order = harmonic.getGroup(sout.str());
       
-	order.getDataSet("potC")   .read(potC   [m][n]);
-	order.getDataSet("rforceC").read(rforceC[m][n]);
-	order.getDataSet("zforceC").read(zforceC[m][n]);
-	order.getDataSet("densC")  .read(densC[m][n]);
+	potC   [m][n] = order.getDataSet("potC")   .read<Eigen::MatrixXd>();
+	rforceC[m][n] = order.getDataSet("rforceC").read<Eigen::MatrixXd>();
+	zforceC[m][n] = order.getDataSet("zforceC").read<Eigen::MatrixXd>();
+	densC  [m][n] = order.getDataSet("densC")  .read<Eigen::MatrixXd>();
       }
     }
 
@@ -7267,10 +7264,10 @@ bool EmpCylSL::ReadH5Cache()
 	sout << n;
 	auto order = harmonic.getGroup(sout.str());
       
-	order.getDataSet("potS")   .read(potS   [m][n]);
-	order.getDataSet("rforceS").read(rforceS[m][n]);
-	order.getDataSet("zforceS").read(zforceS[m][n]);
-	order.getDataSet("densS")  .read(densS[m][n]);
+	potS   [m][n] = order.getDataSet("potS")   .read<Eigen::MatrixXd>();
+	rforceS[m][n] = order.getDataSet("rforceS").read<Eigen::MatrixXd>();
+	zforceS[m][n] = order.getDataSet("zforceS").read<Eigen::MatrixXd>();
+	densS  [m][n] = order.getDataSet("densS")  .read<Eigen::MatrixXd>();
       }
     }
 
diff --git a/exputil/SLGridMP2.cc b/exputil/SLGridMP2.cc
index 7db5367d7..0984fff00 100644
--- a/exputil/SLGridMP2.cc
+++ b/exputil/SLGridMP2.cc
@@ -30,10 +30,8 @@
 
 // For reading and writing cache file
 //
-#include <highfive/H5File.hpp>
-#include <highfive/H5DataSet.hpp>
-#include <highfive/H5DataSpace.hpp>
-#include <highfive/H5Attribute.hpp>
+#include <highfive/highfive.hpp>
+#include <highfive/eigen.hpp>
 
 // For fortran call
 // (This should work both 32-bit and 64-bit . . . )
diff --git a/extern/HighFive b/extern/HighFive
index a794a4d65..12064079d 160000
--- a/extern/HighFive
+++ b/extern/HighFive
@@ -1 +1 @@
-Subproject commit a794a4d651b309069fec5c628d0fb82c2be92662
+Subproject commit 12064079d9533c0636310f25606571b3dbb977cf
diff --git a/extern/yaml-cpp b/extern/yaml-cpp
index 0579ae3d9..1d8ca1f35 160000
--- a/extern/yaml-cpp
+++ b/extern/yaml-cpp
@@ -1 +1 @@
-Subproject commit 0579ae3d976091d7d664aa9d2527e0d0cff25763
+Subproject commit 1d8ca1f35eb3a9c9142462b28282a848e5d29a91
diff --git a/include/EmpCyl2d.H b/include/EmpCyl2d.H
index 440ca3667..703e8e610 100644
--- a/include/EmpCyl2d.H
+++ b/include/EmpCyl2d.H
@@ -17,10 +17,8 @@
 
 // For reading and writing cache file
 //
-#include <highfive/H5File.hpp>
-#include <highfive/H5DataSet.hpp>
-#include <highfive/H5DataSpace.hpp>
-#include <highfive/H5Attribute.hpp>
+#include <highfive/highfive.hpp>
+#include <highfive/eigen.hpp>
 
 /**
    A class that implements most of the members for an Exp force routine

From cd1de795d268ac23f4a717c2910709ef5eafd73d Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Sat, 11 May 2024 21:58:25 -0400
Subject: [PATCH 07/25] Missing updates for HighFive 3.0.0

---
 expui/Coefficients.cc | 15 +++++----------
 expui/Koopman.cc      | 24 +++++++++++-------------
 expui/expMSSA.cc      |  8 +++-----
 3 files changed, 19 insertions(+), 28 deletions(-)

diff --git a/expui/Coefficients.cc b/expui/Coefficients.cc
index 5385be271..17f642933 100644
--- a/expui/Coefficients.cc
+++ b/expui/Coefficients.cc
@@ -12,10 +12,8 @@
 #include <Eigen/Dense>
 #include <yaml-cpp/yaml.h>
 
-#include <highfive/H5File.hpp>
-#include <highfive/H5DataSet.hpp>
-#include <highfive/H5DataSpace.hpp>
-#include <highfive/H5Attribute.hpp>
+#include <highfive/highfive.hpp>
+#include <highfive/eigen.hpp>
 
 #include <Coefficients.H>
 
@@ -802,8 +800,7 @@ namespace CoefClasses
 
       if (Time < Tmin or Time > Tmax) continue;
 
-      Eigen::MatrixXcd in(Mmax+1, Nmax);
-      stanza.getDataSet("coefficients").read(in);
+      auto in = stanza.getDataSet("coefficients").read<Eigen::MatrixXcd>();
       
       // Work around for previous unitiaized data bug; enforces real data
       //
@@ -1173,8 +1170,7 @@ namespace CoefClasses
 
       if (Time < Tmin or Time > Tmax) continue;
 
-      Eigen::VectorXcd in;
-      stanza.getDataSet("coefficients").read(in);
+      auto in = stanza.getDataSet("coefficients").read<Eigen::VectorXcd>();
       
       Eigen::TensorMap<Eigen3d> dat(in.data(), 2*NmaxX+1, 2*NmaxY+1, NmaxZ);
 
@@ -1526,8 +1522,7 @@ namespace CoefClasses
 
       if (Time < Tmin or Time > Tmax) continue;
 
-      Eigen::VectorXcd in;
-      stanza.getDataSet("coefficients").read(in);
+      auto in = stanza.getDataSet("coefficients").read<Eigen::VectorXcd>();
       
       Eigen::TensorMap<Eigen3d> dat(in.data(), 2*NmaxX+1, 2*NmaxY+1, 2*NmaxZ+1);
 
diff --git a/expui/Koopman.cc b/expui/Koopman.cc
index 14f14cf94..b32f00b92 100644
--- a/expui/Koopman.cc
+++ b/expui/Koopman.cc
@@ -31,10 +31,8 @@
 
 #include <Eigen/Dense>
 
-#include <highfive/H5File.hpp>
-#include <highfive/H5DataSet.hpp>
-#include <highfive/H5DataSpace.hpp>
-#include <highfive/H5Attribute.hpp>
+#include <highfive/highfive.hpp>
+#include <highfive/eigen.hpp>
 
 /* For debugging
    #undef eigen_assert
@@ -795,15 +793,15 @@ namespace MSSA {
 
       auto analysis = h5file.getGroup("koopman_analysis");
 
-      analysis.getDataSet("Phi").read(Phi);
-      analysis.getDataSet("X0" ).read(X0 );
-      analysis.getDataSet("X1" ).read(X1 );
-      analysis.getDataSet("U"  ).read(U  );
-      analysis.getDataSet("V"  ).read(V  );
-      analysis.getDataSet("A"  ).read(A  );
-      analysis.getDataSet("L"  ).read(L  );
-      analysis.getDataSet("W"  ).read(W  );
-      analysis.getDataSet("Y"  ).read(Y  );
+      Phi = analysis.getDataSet("Phi").read<Eigen::MatrixXcd>();
+      X0  = analysis.getDataSet("X0" ).read<Eigen::MatrixXd >();
+      X1  = analysis.getDataSet("X1" ).read<Eigen::MatrixXd >();
+      U   = analysis.getDataSet("U"  ).read<Eigen::MatrixXd >();
+      V   = analysis.getDataSet("V"  ).read<Eigen::MatrixXd >();
+      A   = analysis.getDataSet("A"  ).read<Eigen::MatrixXd >();
+      L   = analysis.getDataSet("L"  ).read<Eigen::VectorXcd>();
+      W   = analysis.getDataSet("W"  ).read<Eigen::MatrixXcd>();
+      Y   = analysis.getDataSet("Y"  ).read<Eigen::MatrixXd >();
 
       computed = true;
 
diff --git a/expui/expMSSA.cc b/expui/expMSSA.cc
index 713d34ecb..0aca81d71 100644
--- a/expui/expMSSA.cc
+++ b/expui/expMSSA.cc
@@ -24,10 +24,8 @@
 
 #include <Eigen/Dense>
 
-#include <highfive/H5File.hpp>
-#include <highfive/H5DataSet.hpp>
-#include <highfive/H5DataSpace.hpp>
-#include <highfive/H5Attribute.hpp>
+#include <highfive/highfive.hpp>
+#include <highfive/eigen.hpp>
 
 /* For debugging
    #undef eigen_assert
@@ -1532,7 +1530,7 @@ namespace MSSA {
 	for (int n=0; n<keylist.size(); n++) {
 	  std::ostringstream scnt;
 	  scnt << "RC_" << n;
-	  recon.getDataSet(scnt.str()).read(RC[keylist[n]]);
+	  RC[keylist[n]] = recon.getDataSet(scnt.str()).read<Eigen::MatrixXd>();
 	}
 
 	reconstructed = true;

From 92d9492e4f3518f9f8756ae6eb1ded5df7d9da3f Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Sun, 12 May 2024 12:32:19 -0400
Subject: [PATCH 08/25] Minor code clean-up and fixed a bunch of linger
 fence-post issues (that I already fixed?)

---
 exputil/sbessz.cc |  15 ++-----
 src/Bessel.H      |   2 +-
 src/Bessel.cc     | 104 ++++++++++++++++++++--------------------------
 3 files changed, 50 insertions(+), 71 deletions(-)

diff --git a/exputil/sbessz.cc b/exputil/sbessz.cc
index d8be4065e..91abe4478 100644
--- a/exputil/sbessz.cc
+++ b/exputil/sbessz.cc
@@ -16,33 +16,26 @@
 #define STEPS 6
 #define TOL 1.0e-7
 
-static int NN;
-
-static double zbess(double z)
-{
-  return EXPmath::sph_bessel(NN, z);
-}
-
 Eigen::VectorXd sbessjz(int n, int m)
 {
   Eigen::VectorXd a(m);
 
-  auto zfunc = [n](double z) { return EXPmath::cyl_bessel_j(n, z); };
+  auto zfunc = [n](double z) { return EXPmath::sph_bessel(n, z); };
 
   double dz = M_PI/STEPS;
   double z  = 0.5+fabs((double)n);
   double zl = z, fl, f;
   for (int i=0; i<m; i++) {
-    fl = EXPmath::sph_bessel(n,z);
+    fl = EXPmath::sph_bessel(n, z);
     z += dz;
     f  = EXPmath::sph_bessel(n, z);
     while (f*fl>0) {
       zl = z;
       fl = f;
       z += dz;
-      f = EXPmath::sph_bessel(n,z);
+      f = EXPmath::sph_bessel(n, z);
     }
-    a[i] = zbrent(zbess, zl, z, TOL);
+    a[i] = zbrent(zfunc, zl, z, TOL);
     zl = z;
     fl = f;
   }
diff --git a/src/Bessel.H b/src/Bessel.H
index c768347d7..56bb00028 100644
--- a/src/Bessel.H
+++ b/src/Bessel.H
@@ -56,7 +56,7 @@ private:
     Eigen::VectorXd a;
 
     Roots(int L, int nmax) : l(L), n(nmax) {
-      a = sbessjz(l-1, n);
+      a = sbessjz(l, n);
     }
 
     ~Roots() {}
diff --git a/src/Bessel.cc b/src/Bessel.cc
index 1889b7acd..fdf41b919 100644
--- a/src/Bessel.cc
+++ b/src/Bessel.cc
@@ -24,24 +24,20 @@ Bessel::Bessel(Component* c0, const YAML::Node& conf, MixtureBasis* m) : Spheric
 void Bessel::get_dpotl(int lmax, int nmax, double r, 
 		       Eigen::MatrixXd& p, Eigen::MatrixXd& dp, int tid)
 {
-  double a,aa,aaa,b,bb,bbb;
-  int klo, khi;
-  int l, n;
+  int klo = (int)( (r-r_grid[1])/r_grid_del ) + 1;
+  if (klo < 0) klo = 0;
+  if (klo > RNUM - 2) klo = RNUM - 2;
+  int khi = klo + 1;
 
-  klo = (int)( (r-r_grid[1])/r_grid_del ) + 1;
-  if (klo < 1) klo = 1;
-  if (klo >= RNUM) klo = RNUM - 1;
-  khi = klo + 1;
-
-  a = (r_grid[khi] - r)/r_grid_del;
-  b = (r - r_grid[klo])/r_grid_del;
-
-  aa = a*(a*a-1.0)*r_grid_del*r_grid_del/6.0;
-  bb = b*(b*b-1.0)*r_grid_del*r_grid_del/6.0;
-  aaa = -(3.0*a*a - 1.0)*r_grid_del/6.0;
-  bbb =  (3.0*b*b - 1.0)*r_grid_del/6.0;
+  double a = (r_grid[khi] - r)/r_grid_del;
+  double b = (r - r_grid[klo])/r_grid_del;
+  
+  double aa = a*(a*a-1.0)*r_grid_del*r_grid_del/6.0;
+  double bb = b*(b*b-1.0)*r_grid_del*r_grid_del/6.0;
+  double aaa = -(3.0*a*a - 1.0)*r_grid_del/6.0;
+  double bbb =  (3.0*b*b - 1.0)*r_grid_del/6.0;
 
-  for (l=0; l<=lmax; l++) {
+  for (int l=0; l<=lmax; l++) {
     for (int n=0; n<nmax; n++) {
       p(l, n) = a*potl_grid[l].rw(n, klo) + b*potl_grid[l].rw(n, khi) +
 	aa*potl_grid[l].rw2(n, klo) + bb*potl_grid[l].rw2(n, khi);
@@ -49,24 +45,20 @@ void Bessel::get_dpotl(int lmax, int nmax, double r,
 	aaa*potl_grid[l].rw2(n, klo) + bbb*potl_grid[l].rw2(n, khi);
     }
   }
-
 }
 
 void Bessel::get_potl(int lmax, int nmax, double r, Eigen::MatrixXd& p, int tid)
 {
-  double a,aa,b,bb;
-  int klo, khi;
+  int klo = (int)( (r-r_grid[1])/r_grid_del ) + 1;
+  if (klo < 0) klo = 0;
+  if (klo > RNUM - 2) klo = RNUM - 2;
+  int khi = klo + 1;
 
-  klo = (int)( (r-r_grid[1])/r_grid_del ) + 1;
-  if (klo < 1) klo = 1;
-  if (klo >= RNUM) klo = RNUM - 1;
-  khi = klo + 1;
+  double a = (r_grid[khi] - r)/r_grid_del;
+  double b = (r - r_grid[klo])/r_grid_del;
 
-  a = (r_grid[khi] - r)/r_grid_del;
-  b = (r - r_grid[klo])/r_grid_del;
-
-  aa = a*(a*a-1.0)*r_grid_del*r_grid_del/6.0;
-  bb = b*(b*b-1.0)*r_grid_del*r_grid_del/6.0;
+  double aa = a*(a*a-1.0)*r_grid_del*r_grid_del/6.0;
+  double bb = b*(b*b-1.0)*r_grid_del*r_grid_del/6.0;
 
   for (int l=0; l<=lmax; l++) {
     for (int n=0; n<nmax; n++) {
@@ -78,22 +70,19 @@ void Bessel::get_potl(int lmax, int nmax, double r, Eigen::MatrixXd& p, int tid)
 
 void Bessel::get_dens(int lmax, int nmax, double r, Eigen::MatrixXd& p, int tid)
 {
-  double a, aa, b, bb;
-  int klo, khi;
-
-  klo = (int)( (r-r_grid[1])/r_grid_del ) + 1;
-  if (klo < 1) klo = 1;
-  if (klo >= RNUM) klo = RNUM - 1;
-  khi = klo + 1;
+  int klo = (int)( (r-r_grid[1])/r_grid_del ) + 1;
+  if (klo < 0) klo = 0;
+  if (klo > RNUM - 2) klo = RNUM - 2;
+  int khi = klo + 1;
 
-  a = (r_grid[khi] - r)/r_grid_del;
-  b = (r - r_grid[klo])/r_grid_del;
+  double a = (r_grid[khi] - r)/r_grid_del;
+  double b = (r - r_grid[klo])/r_grid_del;
 
-  aa = a*(a*a-1.0)*r_grid_del*r_grid_del/6.0;
-  bb = b*(b*b-1.0)*r_grid_del*r_grid_del/6.0;
+  double aa = a*(a*a-1.0)*r_grid_del*r_grid_del/6.0;
+  double bb = b*(b*b-1.0)*r_grid_del*r_grid_del/6.0;
 
   for (int l=0; l<=lmax; l++) {
-    for (int n=1; n<=nmax; n++) {
+    for (int n=0; n<nmax; n++) {
       p(l, n) = a*dens_grid[l].rw(n, klo) + b*dens_grid[l].rw(n, khi) +
 	aa*dens_grid[l].rw2(n, klo) + bb*dens_grid[l].rw2(n, khi);
     }
@@ -104,19 +93,16 @@ void Bessel::get_dens(int lmax, int nmax, double r, Eigen::MatrixXd& p, int tid)
 void Bessel::get_potl_dens(int lmax, int nmax, double r, 
 			   Eigen::MatrixXd& p, Eigen::MatrixXd& d, int tid)
 {
-  double a,aa,b,bb;
-  int klo, khi;
+  int klo = (int)( (r-r_grid[1])/r_grid_del ) + 1;
+  if (klo < 0) klo = 0;
+  if (klo > RNUM - 2) klo = RNUM - 2;
+  int khi = klo + 1;
 
-  klo = (int)( (r-r_grid[1])/r_grid_del ) + 1;
-  if (klo < 1) klo = 1;
-  if (klo >= RNUM) klo = RNUM - 1;
-  khi = klo + 1;
+  double a = (r_grid[khi] - r)/r_grid_del;
+  double b = (r - r_grid[klo])/r_grid_del;
 
-  a = (r_grid[khi] - r)/r_grid_del;
-  b = (r - r_grid[klo])/r_grid_del;
-
-  aa = a*(a*a-1.0)*r_grid_del*r_grid_del/6.0;
-  bb = b*(b*b-1.0)*r_grid_del*r_grid_del/6.0;
+  double aa = a*(a*a-1.0)*r_grid_del*r_grid_del/6.0;
+  double bb = b*(b*b-1.0)*r_grid_del*r_grid_del/6.0;
 
   for (int l=0; l<=lmax; l++) {
     for (int n=0; n<nmax; n++) {
@@ -136,7 +122,7 @@ double Bessel::dens(double r, int n)
     throw GenericError("Routine dens() called with n out of bounds", __FILE__, __LINE__, 1001, true);
 
   alpha = p->a[n];
-  return alpha*M_SQRT2/fabs(EXPmath::sph_bessel(p->l, alpha)) * pow(rmax,-2.5) *
+  return alpha*M_SQRT2/fabs(EXPmath::sph_bessel(p->l+1, alpha)) * pow(rmax,-2.5) *
     EXPmath::sph_bessel(p->l, alpha*r/rmax);
 }
 
@@ -148,7 +134,7 @@ double Bessel::potl(double r, int n)
     throw GenericError("Routine potl() called with n out of bounds", __FILE__, __LINE__, 1002, true);
 
   alpha = p->a[n];
-  return M_SQRT2/fabs(alpha*EXPmath::sph_bessel(p->l,alpha)) * pow(rmax,-0.5) *
+  return M_SQRT2/fabs(alpha*EXPmath::sph_bessel(p->l+1,alpha)) * pow(rmax,-0.5) *
     EXPmath::sph_bessel(p->l,alpha*r/rmax);
 }
 
@@ -174,9 +160,9 @@ void Bessel::make_grid(double rmin, double rmax, int lmax, int nmax)
 
     for (int n=0; n<nmax; n++) {
       r = 0.0;
-      for (int ir=1; ir<=RNUM; ir++, r+=r_grid_del) {
-	potl_grid[l].rw(n, ir) = potl(r, n+1);
-	dens_grid[l].rw(n, ir) = dens(r, n+1);
+      for (int ir=0; ir<RNUM; ir++, r+=r_grid_del) {
+	potl_grid[l].rw(n, ir) = potl(r, n);
+	dens_grid[l].rw(n, ir) = dens(r, n);
       }
       
       {
@@ -200,12 +186,12 @@ void Bessel::make_grid(double rmin, double rmax, int lmax, int nmax)
 
   // check table
 
-  for (int ir=1; ir<=RNUM; ir++) assert(!std::isnan(r_grid[ir]));
+  for (int ir=0; ir<RNUM; ir++) assert(!std::isnan(r_grid[ir]));
   for (int l=0; l<=lmax; l++) {
     assert(!std::isnan(potl_grid[l].nmax));
     assert(!std::isnan(dens_grid[l].nmax));
-    for (int n=1; n<=nmax; n++) {
-      for (int ir=1; ir<=RNUM; ir++) {
+    for (int n=0; n<nmax; n++) {
+      for (int ir=0; ir<RNUM; ir++) {
 	assert(!std::isnan(potl_grid[l].rw (n, ir)));
 	assert(!std::isnan(potl_grid[l].rw2(n, ir)));
 	assert(!std::isnan(dens_grid[l].rw (n, ir)));

From ab6d0c7df97137d801512b120f6bb187e26d413f Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Sun, 12 May 2024 13:01:34 -0400
Subject: [PATCH 09/25] Add a parsing stanza and let RNUM be a configurable
 parameter [no ci]

---
 src/Bessel.H  | 20 ++++++++++++++++----
 src/Bessel.cc | 29 ++++++++++++++++++++++++++++-
 2 files changed, 44 insertions(+), 5 deletions(-)

diff --git a/src/Bessel.H b/src/Bessel.H
index 56bb00028..1f21f5709 100644
--- a/src/Bessel.H
+++ b/src/Bessel.H
@@ -17,6 +17,8 @@ class Bessel : public SphericalBasis
 
 private:
 
+  //@{
+  //! Required members for Spherical Basis
   void get_pot_coefs(int l, double *coef, double *p, double *dp);
   void get_pot_coefs_safe(int l, double *coef, double *p, double *dp,
 			      double **potd1, double **dpot1);
@@ -26,8 +28,10 @@ private:
   void get_dens(int lmax, int nmax, double r, Eigen::MatrixXd& p, int tid);
   void get_potl_dens(int lmax, int nmax, double r, Eigen::MatrixXd& p, Eigen::MatrixXd& d,int tid);
   double get_dens(double r, int l, double *coef);
+  //@}
 
-  void initialize() {}
+  //! Initialize parameters from YAML
+  void initialize();
 
   bool firstime_coef;
   bool firstime_accel;
@@ -41,9 +45,12 @@ private:
     int nmax;
   };
   
+  //@{
+  //! Grid storage and parameters
   std::vector<RGrid> dens_grid, potl_grid;
   Eigen::VectorXd r_grid;
   double r_grid_del;
+  //@}
 
   //! Cache roots for spherical Bessel functions
   class Roots
@@ -62,21 +69,26 @@ private:
     ~Roots() {}
   };
 
+  //! Root database isntance
   std::shared_ptr<Roots> p;
 
+  //@{
+  //! Density and potential members
   double dens(double r, int n);
   double potl(double r, int n);
+  //@}
 
+  //! Make the density and potential grid
   void make_grid(double rmin, double rmax, int lmax, int nmax);
   
   //! Valid keys for YAML configurations
   static const std::set<std::string> valid_keys;
 
+  //! Number of entries in the fixed table
+  int RNUM;
+
 public:
   
-  //! Number of entries in fixed table (static variable)
-  static int RNUM;
-
   //! Constructor
   Bessel(Component* c0, const YAML::Node& conf, MixtureBasis* m=0);
 
diff --git a/src/Bessel.cc b/src/Bessel.cc
index fdf41b919..657a9d1b4 100644
--- a/src/Bessel.cc
+++ b/src/Bessel.cc
@@ -6,7 +6,34 @@
 #include <Bessel.H>
 #include <EXPmath.H>
 
-int Bessel::RNUM = 1000;
+const std::set<std::string>
+Bessel::valid_keys = {
+  "rnum"
+};
+
+void Bessel::initialize()
+{
+  // Remove matched keys
+  //
+  for (auto v : valid_keys) current_keys.erase(v);
+  
+  // Assign values from YAML
+  //
+  try {
+    if (conf["rnum"])      RNUM       = conf["rnum"].as<int>();
+    else                   RNUM       = 1000;
+  }
+  catch (YAML::Exception & error) {
+    if (myid==0) std::cout << "Error parsing parameters in Sphere: "
+			   << error.what() << std::endl
+			   << std::string(60, '-') << std::endl
+			   << "Config node"        << std::endl
+			   << std::string(60, '-') << std::endl
+			   << conf                 << std::endl
+			   << std::string(60, '-') << std::endl;
+    throw std::runtime_error("Sphere::initialze: error parsing YAML");
+  }
+}
 
 Bessel::Bessel(Component* c0, const YAML::Node& conf, MixtureBasis* m) : SphericalBasis(c0, conf, m)
 {

From 1f09674115ad661a096cd671103222d7376a51e9 Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Sun, 12 May 2024 14:00:25 -0400
Subject: [PATCH 10/25] Add the git submodule commands to the CMake workflow

---
 CMakeLists.txt | 13 +++++++++++++
 1 file changed, 13 insertions(+)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index c105459ef..598bf4ad9 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -187,6 +187,19 @@ execute_process(
   OUTPUT_STRIP_TRAILING_WHITESPACE
 )
 
+# Git submodule updates
+execute_process(
+  COMMAND git submodule update --init --recursive
+  WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}
+  RESULT_VARIABLE GIT_SUBMOD_RESULT
+)
+
+if(NOT GIT_SUBMOD_RESULT EQUAL "0")
+   message(FATAL_ERROR "git submodule update --init --recursive failed ${GIT_SUBMOD_RESULT}, please checkout submodules")
+else()
+   message(STATUS "Submodules updated successfully - good")
+endif()
+
 # Get the latest abbreviated commit hash of the working branch
 execute_process(
   COMMAND git rev-parse HEAD

From 254df3d0761ee6a04d75f102e7e9e38e37252b08 Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Sun, 12 May 2024 22:00:26 -0400
Subject: [PATCH 11/25] Preliminary implementation of Bessel basis for pyEXP

---
 expui/BasisFactory.cc  |   3 +
 expui/BiorthBasis.H    | 226 ++++++++++++++++++++-------
 expui/BiorthBasis.cc   | 246 +++++++++++++++++++++---------
 expui/CMakeLists.txt   |   2 +-
 pyEXP/BasisWrappers.cc | 339 ++++++++++++++++++++++++++++-------------
 src/Bessel.cc          |   8 +-
 6 files changed, 587 insertions(+), 237 deletions(-)

diff --git a/expui/BasisFactory.cc b/expui/BasisFactory.cc
index 599a0eafb..5b574a8e2 100644
--- a/expui/BasisFactory.cc
+++ b/expui/BasisFactory.cc
@@ -175,6 +175,9 @@ namespace BasisClasses
       if ( !name.compare("sphereSL") ) {
 	basis = std::make_shared<SphericalSL>(conf);
       }
+      else if ( !name.compare("bessel") ) {
+	basis = std::make_shared<Bessel>(conf);
+      }
       else if ( !name.compare("cylinder") ) {
 	basis = std::make_shared<Cylindrical>(conf);
       }
diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H
index ce88f02f7..a4af27c1a 100644
--- a/expui/BiorthBasis.H
+++ b/expui/BiorthBasis.H
@@ -10,6 +10,7 @@
 
 #include <ParticleReader.H>
 #include <Coefficients.H>
+#include <BiorthBess.H>
 #include <BasisFactory.H>
 #include <BiorthCube.H>
 #include <SLGridMP2.H>
@@ -60,7 +61,6 @@ namespace BasisClasses
     virtual std::vector<double>
     sph_eval(double r, double costh, double phi) = 0;
 
-
     //! Evaluate fields in cylindrical coordinates in centered coordinate system
     virtual std::vector<double>
     cyl_eval(double r, double costh, double phi) = 0;
@@ -164,10 +164,10 @@ namespace BasisClasses
   };
   
   /**
-     Uses SLGridSph basis to evaluate expansion coeffients and provide
-     potential and density basis fields
+     An abstract spherical basis to evaluate expansion coeffients and
+     provide potential and density basis fields
   */
-  class SphericalSL : public BiorthBasis
+  class Spherical : public BiorthBasis
   {
     
   public:
@@ -175,45 +175,13 @@ namespace BasisClasses
     using BasisMap   = std::map<std::string, Eigen::VectorXd>;
     using BasisArray = std::vector<std::vector<BasisMap>>;
 
-  private:
+  protected:
     
     //! Helper for constructor
     void initialize();
 
-    std::shared_ptr<SLGridSph> sl;
-    std::shared_ptr<SphericalModelTable> mod;
-    
-    std::string model_file;
-    int lmax, nmax, cmap, numr;
-    double rmin, rmax, rmap;
-    
-    bool NO_L0, NO_L1, EVEN_L, EVEN_M, M0_only;
-    
-    std::vector<Eigen::MatrixXd> potd, dpot, dpt2, dend;
-    std::vector<Eigen::MatrixXd> legs, dlegs, d2legs;
-    
-    Eigen::MatrixXd factorial;
-    Eigen::MatrixXd expcoef;
-    double scale;
-    int N1, N2;
-    int used;
-    
-    using matT = std::vector<Eigen::MatrixXd>;
-    using vecT = std::vector<Eigen::VectorXd>;
-    
-    double totalMass;
-    int npart;
-    
-    Eigen::VectorXd work;
-    
-    //! For coefficient writing
-    typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>
-    EigenColMajor;
-    
-  protected:
-
-    //! Load coefficients into the new CoefStruct
-    virtual void load_coefs(CoefClasses::CoefStrPtr coefs, double time);
+    //! Load coefficients for a particular time
+    virtual void load_coefs(CoefClasses::CoefStrPtr coef, double time);
 
     //! Set coefficients
     virtual void set_coefs(CoefClasses::CoefStrPtr coefs);
@@ -222,7 +190,7 @@ namespace BasisClasses
     static const std::set<std::string> valid_keys;
 
     //! Return readable class name
-    virtual const std::string classname() { return "SphericalSL";}
+    virtual const std::string classname() { return "Spherical";}
 
     //! Subspace index
     virtual const std::string harmonic() { return "l";}
@@ -239,16 +207,60 @@ namespace BasisClasses
     virtual std::vector<double>
     cyl_eval(double R, double z, double phi);
 
+    //@{
+    //! Required basis members
+
+    //! Get potential
+    virtual void get_pot(Eigen::MatrixXd& tab, double x) = 0;
+
+    //! Get density
+    virtual void get_dens(Eigen::MatrixXd& tab, double x) = 0;
+
+    //! Get force
+    virtual void get_force(Eigen::MatrixXd& tab, double x) = 0;
+
+    //@}
+
+    //@{
+    //! Internal parameters and storage
+    int lmax, nmax, cmap, numr;
+    double rmin, rmax, rmap;
+    
+    bool NO_L0, NO_L1, EVEN_L, EVEN_M, M0_only;
+
+    std::vector<Eigen::MatrixXd> potd, dpot, dpt2, dend;
+    std::vector<Eigen::MatrixXd> legs, dlegs, d2legs;
+
+    Eigen::MatrixXd factorial;
+    Eigen::MatrixXd expcoef;
+    double scale;
+    int N1, N2;
+    int used;
+
+    using matT = std::vector<Eigen::MatrixXd>;
+    using vecT = std::vector<Eigen::VectorXd>;
+
+    double totalMass;
+    int npart;
+    
+    Eigen::VectorXd work;
+    
+    //! For coefficient writing
+    typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>
+    EigenColMajor;
+
+    //@}
+
   public:
     
     //! Constructor from YAML node
-    SphericalSL(const YAML::Node& conf);
+    Spherical(const YAML::Node& conf, const std::string& forceID);
     
     //! Constructor from YAML string
-    SphericalSL(const std::string& confstr);
+    Spherical(const std::string& confstr, const std::string& forceID);
     
     //! Destructor
-    virtual ~SphericalSL(void) {}
+    virtual ~Spherical(void) {}
     
     //! Print and return the cache parameters
     static std::map<std::string, std::string>
@@ -276,29 +288,134 @@ namespace BasisClasses
     int getNmax() { return nmax; }
 
     //! Return potential-density pair of a vector of a vector of 1d
-    //! basis-function grids for SphericalSL, logarithmically spaced
-    //! between [logxmin, logxmax] (base 10).
-    BasisArray getBasis
-    (double logxmin=-3.0, double logxmax=0.5, int numgrid=2000);
+    //! basis-function grids for Spherical [rmin, rmax]
+    virtual BasisArray getBasis
+    (double rmin=0.0, double rax=1.0, int numgrid=2000);
 
     //! Compute the orthogonality of the basis by returning inner
     //! produce matrices
-    std::vector<Eigen::MatrixXd> orthoCheck(int knots=40)
-    {
-      return sl->orthoCheck(knots);
-    }
+    virtual std::vector<Eigen::MatrixXd> orthoCheck(int knots) = 0;
 
     //! Biorthogonality sanity check
     bool orthoTest(int knots=100)
     {
-      auto [ret, worst, lworst] = orthoCompute(sl->orthoCheck(knots));
+      auto [ret, worst, lworst] = orthoCompute(orthoCheck(knots));
       // For the CTest log
-      std::cout << "SphericalSL::orthoTest: worst=" << worst << std::endl;
+      std::cout << "Spherical::orthoTest: worst=" << worst << std::endl;
       return ret;
     }
 
   };
+
+  /**
+     Uses SLGridSph basis to evaluate expansion coeffients and provide
+     potential and density basis fields
+  */
+  class SphericalSL : public Spherical
+  {
+    
+  protected:
+    
+    //! Helper for constructor
+    void initialize();
+
+    static const std::set<std::string> valid_keys;
+
+    std::shared_ptr<SLGridSph> sl;
+    std::shared_ptr<SphericalModelTable> mod;
+    
+    std::string model_file;
+    
+    //! Return readable class name
+    const std::string classname() { return "SphericalSL";}
+
+    // Get potential
+    void get_pot(Eigen::MatrixXd& tab, double r)
+    { sl->get_pot(tab, r); }
+
+    // Get density
+    void get_dens(Eigen::MatrixXd& tab, double r)
+    { sl->get_dens(tab, r); }
+
+    // Get force
+    void get_force(Eigen::MatrixXd& tab, double r)
+    { sl->get_force(tab, r); }
+
+  public:
+    
+    //! Constructor from YAML node
+    SphericalSL(const YAML::Node& conf);
+    
+    //! Constructor from YAML string
+    SphericalSL(const std::string& confstr);
+    
+    //! Destructor
+    virtual ~SphericalSL(void) {}
+    
+    //! Return potential-density pair of a vector of a vector of 1d
+    //! basis-function grids for SphericalSL, logarithmically spaced
+    //! between [logxmin, logxmax] (base 10).
+    BasisArray getBasis(double logxmin=-3.0, double logxmax=0.5, int numgrid=2000);
+
+    //! Compute the orthogonality of the basis by returning inner
+    //! produce matrices
+    std::vector<Eigen::MatrixXd> orthoCheck(int knots=40)
+    { return sl->orthoCheck(knots); }
+  };
+  
+  /**
+     Uses Bessel basis to evaluate expansion coeffients and provide
+     potential and density basis fields
+  */
+  class Bessel : public Spherical
+  {
+    
+  protected:
+    
+    //! Helper for constructor
+    void initialize();
+
+    static const std::set<std::string> valid_keys;
+
+    //! Return readable class name
+    const std::string classname() { return "Bessel";}
+
+    //! Grid size for Bessel function table
+    int rnum = 2000;
+
+    //! Biorthgonal Bessel function generator
+    std::shared_ptr<BiorthBess> bess;
+    
+    // Get potential
+    void get_pot(Eigen::MatrixXd& tab, double r)
+    { bess->get_potl(r, tab); }
+
+    // Get density
+    void get_dens(Eigen::MatrixXd& tab, double r)
+    { bess->get_dens(r, tab); }
+
+    // Get force
+    void get_force(Eigen::MatrixXd& tab, double r)
+    { bess->get_dpotl(r, tab); }
+
+  public:
+    
+    //! Constructor from YAML node
+    Bessel(const YAML::Node& conf);
+    
+    //! Constructor from YAML string
+    Bessel(const std::string& confstr);
+    
+    //! Destructor
+    virtual ~Bessel(void) {}
+
+    //! Compute the orthogonality of the basis by returning inner
+    //! produce matrices
+    std::vector<Eigen::MatrixXd> orthoCheck(int knots=40)
+    { return bess->orthoCheck(knots); }
+  };
   
+
   /**
      Uses the BiorthCyl basis to evaluate expansion coeffients and
      provide potential and density basis fields
@@ -408,8 +525,7 @@ namespace BasisClasses
     //! Return a vector of a vector of 1d basis-function grids for
     //! FlatDisk, logarithmically spaced between [logxmin, logxmax]
     //! (base 10).
-    BasisArray getBasis
-    (double logxmin=-3.0, double logxmax=0.5, int numgrid=2000);
+    BasisArray getBasis(double logxmin=-3.0, double logxmax=0.5, int numgrid=2000);
 
     //! Compute the orthogonality of the basis by returning inner
     //! produce matrices
diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc
index 68897a4cc..9625e4f2d 100644
--- a/expui/BiorthBasis.cc
+++ b/expui/BiorthBasis.cc
@@ -14,7 +14,7 @@
 namespace BasisClasses
 {
   const std::set<std::string>
-  SphericalSL::valid_keys = {
+  Spherical::valid_keys = {
     "rmapping",
     "cmap",
     "Lmax",
@@ -32,7 +32,6 @@ namespace BasisClasses
     "tkcum",
     "tk_type",
     "nmax",
-    "modelname",
     "scale",
     "rmin",
     "rmax",
@@ -59,7 +58,9 @@ namespace BasisClasses
     "logr",
     "plummer",
     "self_consistent",
-    "cachename"
+    "cachename",
+    "modelname",
+    "rnum"
   };
 
   std::vector<std::string> BiorthBasis::getFieldLabels(const Coord ctype)
@@ -90,19 +91,39 @@ namespace BasisClasses
     return labels;
   }
 
-  SphericalSL::SphericalSL(const YAML::Node& CONF) :
-    BiorthBasis(CONF, "sphereSL")
+  Spherical::Spherical(const YAML::Node& CONF, const std::string& forceID) :
+    BiorthBasis(CONF, forceID)
   {
     initialize();
   }
 
-  SphericalSL::SphericalSL(const std::string& confstr) :
-    BiorthBasis(confstr, "sphereSL")
+  Spherical::Spherical(const std::string& confstr, const std::string& forceID) :
+    BiorthBasis(confstr, forceID)
   {
     initialize();
   }
 
-  void SphericalSL::initialize()
+  SphericalSL::SphericalSL(const YAML::Node& CONF) : Spherical(CONF, "sphereSL")
+  {
+    initialize();
+  }
+
+  SphericalSL::SphericalSL(const std::string& confstr) : Spherical(confstr, "sphereSL")
+  {
+    initialize();
+  }
+
+  Bessel::Bessel(const YAML::Node& CONF) :  Spherical(CONF, "Bessel")
+  {
+    initialize();
+  }
+
+  Bessel::Bessel(const std::string& confstr) : Spherical(confstr, "Bessel")
+  {
+    initialize();
+  }
+
+  void Spherical::initialize()
   {
 
     // Assign some defaults
@@ -110,27 +131,17 @@ namespace BasisClasses
     cmap       = 1;
     lmax       = 6;
     nmax       = 18;
-    model_file = "SLGridSph.model";
     
     // Check for unmatched keys
     //
     auto unmatched = YamlCheck(conf, valid_keys);
     if (unmatched.size())
-      throw YamlConfigError("Basis::Basis::SphericalSL", "parameter", unmatched, __FILE__, __LINE__);
+      throw YamlConfigError("Basis::Basis::Spherical", "parameter", unmatched, __FILE__, __LINE__);
     
-    // Default cachename, empty by default
-    //
-    std::string cachename;
-
-    // Assign values from YAML
-    //
-    double rmap = 1.0;
-
     try {
       if (conf["cmap"])      cmap       = conf["cmap"].as<int>();
       if (conf["Lmax"])      lmax       = conf["Lmax"].as<int>();
       if (conf["nmax"])      nmax       = conf["nmax"].as<int>();
-      if (conf["modelname"]) model_file = conf["modelname"].as<std::string>();
       
       if (conf["rmapping"]) 
 	rmap = conf["rmapping"].as<double>();
@@ -168,7 +179,6 @@ 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["cachename"]) cachename = conf["cachename"].as<std::string>();
     } 
     catch (YAML::Exception & error) {
       if (myid==0) std::cout << "Error parsing parameter stanza for <"
@@ -178,43 +188,9 @@ namespace BasisClasses
 			     << conf                 << std::endl
 			     << std::string(60, '-') << std::endl;
       
-      throw std::runtime_error("SphericalSL: error parsing YAML");
-    }
-    
-    // Check for non-null cache file name.  This must be specified
-    // to prevent recomputation and unexpected behavior.
-    //
-    if (cachename.size() == 0) {
-      throw std::runtime_error
-	("SphericalSL requires a specified cachename in your YAML config\n"
-	 "for consistency with previous invocations and existing coefficient\n"
-	 "sets.  Please add explicitly add 'cachename: name' to your config\n"
-	 "with new 'name' for creating a basis or an existing 'name' for\n"
-	 "reading a previously generated basis cache\n");
+      throw std::runtime_error("Spherical: error parsing YAML");
     }
 
-    // Set MPI flag in SLGridSph from MPI_Initialized
-    SLGridSph::mpi = use_mpi ? 1 : 0;
-    
-    // Instantiate to get min/max radius from the model
-    mod = std::make_shared<SphericalModelTable>(model_file);
-    
-    // Set rmin to a sane value if not specified
-    if (not conf["rmin"] or rmin < mod->get_min_radius()) 
-      rmin = mod->get_min_radius();
-
-    // Set rmax to a sane value if not specified
-    if (not conf["rmax"] or rmax > mod->get_max_radius()) 
-      rmax = mod->get_max_radius()*0.99;
-    
-    // Finally, make the Sturm-Lioville basis...
-    sl = std::make_shared<SLGridSph>
-      (model_file, lmax, nmax, numr, rmin, rmax, true, cmap, rmap,
-       0, 1, cachename);
-    
-    // Test basis for consistency
-    orthoTest(200);
-
     // Number of possible threads
     int nthrds = omp_get_max_threads();
     
@@ -258,7 +234,94 @@ namespace BasisClasses
     coordinates = Coord::Spherical;
   }
   
-  void SphericalSL::reset_coefs(void)
+  void SphericalSL::initialize()
+  {
+
+    // Assign some defaults
+    //
+    model_file = "SLGridSph.model";
+    
+    // Default cachename, empty by default
+    //
+    std::string cachename;
+
+    try {
+      if (conf["modelname"]) model_file = conf["modelname"].as<std::string>();
+      if (conf["cachename"]) cachename  = conf["cachename"].as<std::string>();
+    } 
+    catch (YAML::Exception & error) {
+      if (myid==0) std::cout << "Error parsing parameter stanza for <"
+			     << name << ">: "
+			     << error.what() << std::endl
+			     << std::string(60, '-') << std::endl
+			     << conf                 << std::endl
+			     << std::string(60, '-') << std::endl;
+      
+      throw std::runtime_error("SphericalSL: error parsing YAML");
+    }
+    
+    // Check for non-null cache file name.  This must be specified
+    // to prevent recomputation and unexpected behavior.
+    //
+    if (cachename.size() == 0) {
+      throw std::runtime_error
+	("SphericalSL requires a specified cachename in your YAML config\n"
+	 "for consistency with previous invocations and existing coefficient\n"
+	 "sets.  Please add explicitly add 'cachename: name' to your config\n"
+	 "with new 'name' for creating a basis or an existing 'name' for\n"
+	 "reading a previously generated basis cache\n");
+    }
+
+    // Set MPI flag in SLGridSph from MPI_Initialized
+    SLGridSph::mpi = use_mpi ? 1 : 0;
+    
+    // Instantiate to get min/max radius from the model
+    mod = std::make_shared<SphericalModelTable>(model_file);
+    
+    // Set rmin to a sane value if not specified
+    if (not conf["rmin"] or rmin < mod->get_min_radius()) 
+      rmin = mod->get_min_radius();
+
+    // Set rmax to a sane value if not specified
+    if (not conf["rmax"] or rmax > mod->get_max_radius()) 
+      rmax = mod->get_max_radius()*0.99;
+    
+    // Finally, make the Sturm-Lioville basis...
+    sl = std::make_shared<SLGridSph>
+      (model_file, lmax, nmax, numr, rmin, rmax, true, cmap, rmap,
+       0, 1, cachename);
+    
+    // Test basis for consistency
+    orthoTest(200);
+  }
+  
+  void Bessel::initialize()
+  {
+    try {
+      if (conf["rnum"])
+	rnum = conf["rnum"].as<int>();
+      else
+	rnum = 2000;
+    } 
+    catch (YAML::Exception & error) {
+      if (myid==0) std::cout << "Error parsing parameter stanza for <"
+			     << name << ">: "
+			     << error.what() << std::endl
+			     << std::string(60, '-') << std::endl
+			     << conf                 << std::endl
+			     << std::string(60, '-') << std::endl;
+      
+      throw std::runtime_error("SphericalSL: error parsing YAML");
+    }
+    
+    // Finally, make the Sturm-Lioville basis...
+    bess = std::make_shared<BiorthBess>(rmax, lmax, nmax, rnum);
+
+    // Test basis for consistency
+    orthoTest(200);
+  }
+  
+  void Spherical::reset_coefs(void)
   {
     if (expcoef.rows()>0 && expcoef.cols()>0) expcoef.setZero();
     totalMass = 0.0;
@@ -266,7 +329,7 @@ namespace BasisClasses
   }
   
   
-  void SphericalSL::load_coefs(CoefClasses::CoefStrPtr coef, double time)
+  void Spherical::load_coefs(CoefClasses::CoefStrPtr coef, double time)
   {
     CoefClasses::SphStruct* cf = dynamic_cast<CoefClasses::SphStruct*>(coef.get());
 
@@ -301,12 +364,12 @@ namespace BasisClasses
     }
   }
 
-  void SphericalSL::set_coefs(CoefClasses::CoefStrPtr coef)
+  void Spherical::set_coefs(CoefClasses::CoefStrPtr coef)
   {
     // Sanity check on derived class type
     //
     if (typeid(*coef) != typeid(CoefClasses::SphStruct))
-      throw std::runtime_error("SphericalSL::set_coefs: you must pass a CoefClasses::SphStruct");
+      throw std::runtime_error("Spherical::set_coefs: you must pass a CoefClasses::SphStruct");
 
     // Sanity check on dimensionality
     //
@@ -318,7 +381,7 @@ namespace BasisClasses
       int rexp = (lmax+1)*(lmax+2)/2;
       if (rows != rexp or cols != nmax) {
 	std::ostringstream sout;
-	sout << "SphericalSL::set_coefs: the basis has (lmax, nmax)=("
+	sout << "Spherical::set_coefs: the basis has (lmax, nmax)=("
 	     << lmax << ", " << nmax
 	     << ") and the dimensions must be (rows, cols)=("
 	     << rexp << ", " << nmax
@@ -357,7 +420,7 @@ namespace BasisClasses
       coefctr = {0.0, 0.0, 0.0};
   }
 
-  void SphericalSL::accumulate(double x, double y, double z, double mass)
+  void Spherical::accumulate(double x, double y, double z, double mass)
   {
     double fac, fac1, fac2, fac4;
     double norm = -4.0*M_PI;
@@ -381,7 +444,7 @@ namespace BasisClasses
     used++;
     totalMass += mass;
     
-    sl->get_pot(potd[tid], rs);
+    get_pot(potd[tid], rs);
     
     legendre_R(lmax, costh, legs[tid]);
     
@@ -420,7 +483,7 @@ namespace BasisClasses
     
   }
   
-  void SphericalSL::make_coefs()
+  void Spherical::make_coefs()
   {
     if (use_mpi) {
       
@@ -437,7 +500,7 @@ namespace BasisClasses
   }
   
   std::vector<double>
-  SphericalSL::sph_eval(double r, double costh, double phi)
+  Spherical::sph_eval(double r, double costh, double phi)
   {
     // Get thread id
     int tid = omp_get_thread_num();
@@ -447,9 +510,9 @@ namespace BasisClasses
     
     fac1 = factorial(0, 0);
     
-    sl->get_dens (dend[tid], r/scale);
-    sl->get_pot  (potd[tid], r/scale);
-    sl->get_force(dpot[tid], r/scale);
+    get_dens (dend[tid], r/scale);
+    get_pot  (potd[tid], r/scale);
+    get_force(dpot[tid], r/scale);
     
     legendre_R(lmax, costh, legs[tid], dlegs[tid]);
     
@@ -547,7 +610,7 @@ namespace BasisClasses
 
 
   std::vector<double>
-  SphericalSL::cyl_eval(double R, double z, double phi)
+  Spherical::cyl_eval(double R, double z, double phi)
   {
     double r = sqrt(R*R + z*z) + 1.0e-18;
     double costh = z/r, sinth = R/r;
@@ -563,7 +626,7 @@ namespace BasisClasses
 
   // Evaluate in cartesian coordinates
   std::vector<double>
-  SphericalSL::crt_eval
+  Spherical::crt_eval
   (double x, double y, double z)
   {
     double R = sqrt(x*x + y*y);
@@ -578,7 +641,7 @@ namespace BasisClasses
   }
   
 
-  SphericalSL::BasisArray SphericalSL::getBasis
+  Spherical::BasisArray SphericalSL::getBasis
   (double logxmin, double logxmax, int numgrid)
   {
     // Assing return storage
@@ -599,9 +662,45 @@ namespace BasisClasses
     Eigen::MatrixXd tabpot, tabden, tabfrc;
 
     for (int i=0; i<numgrid; i++) {
-      sl->get_pot  (tabpot, pow(10.0, logxmin + dx*i));
-      sl->get_dens (tabden, pow(10.0, logxmin + dx*i));
-      sl->get_force(tabfrc, pow(10.0, logxmin + dx*i));
+      get_pot  (tabpot, pow(10.0, logxmin + dx*i));
+      get_dens (tabden, pow(10.0, logxmin + dx*i));
+      get_force(tabfrc, pow(10.0, logxmin + dx*i));
+      for (int l=0; l<=lmax; l++) {
+	for (int n=0; n<nmax;n++){
+	  ret[l][n]["potential"](i) = tabpot(l, n);
+	  ret[l][n]["density"  ](i) = tabden(l, n);
+	  ret[l][n]["rforce"   ](i) = tabfrc(l, n);
+	}
+      }
+    }
+    
+    return ret;
+  }
+
+  Spherical::BasisArray Spherical::getBasis
+  (double rmin, double rmax, int numgrid)
+  {
+    // Assing return storage
+    BasisArray ret (lmax+1);
+    for (auto & v : ret) {
+      v.resize(nmax);
+      for (auto & u : v) {
+	u["potential"].resize(numgrid); // Potential
+	u["density"  ].resize(numgrid); // Density
+	u["rforce"   ].resize(numgrid); // Radial force
+      }
+    }
+
+    // Radial grid spacing
+    double dr = (rmax - rmin)/(numgrid-1);
+
+    // Basis storage
+    Eigen::MatrixXd tabpot, tabden, tabfrc;
+
+    for (int i=0; i<numgrid; i++) {
+      get_pot  (tabpot, rmin + dr*i);
+      get_dens (tabden, rmin + dr*i);
+      get_force(tabfrc, rmin + dr*i);
       for (int l=0; l<=lmax; l++) {
 	for (int n=0; n<nmax;n++){
 	  ret[l][n]["potential"](i) = tabpot(l, n);
@@ -614,6 +713,7 @@ namespace BasisClasses
     return ret;
   }
 
+
 #define MINEPS 1.0e-10
   
   void BiorthBasis::legendre_R(int lmax, double x, Eigen::MatrixXd& p)
diff --git a/expui/CMakeLists.txt b/expui/CMakeLists.txt
index 8193c278c..8690817f7 100644
--- a/expui/CMakeLists.txt
+++ b/expui/CMakeLists.txt
@@ -34,7 +34,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)
+  Koopman.cc BiorthBess.cc)
 add_library(expui ${expui_SOURCES})
 set_target_properties(expui PROPERTIES OUTPUT_NAME expui)
 target_include_directories(expui PUBLIC ${common_INCLUDE})
diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc
index 2eb64e749..df9e8346f 100644
--- a/pyEXP/BasisWrappers.cc
+++ b/pyEXP/BasisWrappers.cc
@@ -20,29 +20,32 @@ void BasisFactoryClasses(py::module &m)
     This module provides a factory class that will create biorthogonal
     bases from input YAML configuration files.  Each basis can then be
     used to compute coefficients, provide field quantities such as
-    forces and, together with the FieldGenerator, surfaces and fields for
-    visualization.
+    forces and, together with the FieldGenerator, surfaces and fields
+    for visualization.
 
-    Seven bases are currently implemented:
+    Eight bases are currently implemented:
 
      1. SphericalSL, the Sturm-Liouiville spherical basis;
 
-     2. Cylindrical, created created by computing empirical orthogonal functions
+     2. Bessel, the classic spherical biorthogonal constructed from
+        the eigenfunctions of the spherical Laplacian;
+
+     3. Cylindrical, created created by computing empirical orthogonal functions
         over a densely sampled SphericalSL basis;
 
-     3. FlatDisk, an EOF rotation of the finite Bessel basis; and
+     4. FlatDisk, an EOF rotation of the finite Bessel basis; and
 
-     4. Slab, a biorthogonal basis for a slab geometry with a finite
+     5. Slab, a biorthogonal basis for a slab geometry with a finite
         finite vertical extent.  The basis is constructed from direct
         solution of the Sturm-Liouville equation.
 
-     5. Cube, a periodic cube basis whose functions are the Cartesian
+     6. Cube, a periodic cube basis whose functions are the Cartesian
         eigenfunctions of the Cartesian Laplacian: sines and cosines.
 
-     6. FieldBasis, for computing user-provided quantities from a
+     7. FieldBasis, for computing user-provided quantities from a
         phase-space snapshot.
 
-     7. VelocityBasis, for computing the mean field velocity fields from
+     8. VelocityBasis, for computing the mean field velocity fields from
         a phase-space snapshot.  This is a specialized version of FieldBasis.
 
     Each of these bases take a YAML configuration file as input. These parameter
@@ -284,7 +287,7 @@ void BasisFactoryClasses(py::module &m)
 
     virtual void
     addFromArray(Eigen::VectorXd& m, RowMatrixXd& p, bool roundrobin, bool posvelrows) override {
-PYBIND11_OVERRIDE_PURE(void, Basis, addFromArray, m, p, roundrobin, posvelrows);
+      PYBIND11_OVERRIDE_PURE(void, Basis, addFromArray, m, p, roundrobin, posvelrows);
     }
   };
 
@@ -409,49 +412,67 @@ PYBIND11_OVERRIDE_PURE(void, Basis, addFromArray, m, p, roundrobin, posvelrows);
 
   };
 
-  class PySphericalSL : public SphericalSL
+  class PySpherical : public Spherical
   {
   protected:
 
     void load_coefs(CoefClasses::CoefStrPtr coefs, double time) override {
-      PYBIND11_OVERRIDE(void, SphericalSL, load_coefs, coefs, time);
+      PYBIND11_OVERRIDE(void, Spherical, load_coefs, coefs, time);
     }
     
     void set_coefs(CoefClasses::CoefStrPtr coefs) override {
-      PYBIND11_OVERRIDE(void, SphericalSL, set_coefs, coefs);
+      PYBIND11_OVERRIDE(void, Spherical, set_coefs, coefs);
     }
 
     const std::string classname() override {
-      PYBIND11_OVERRIDE(std::string, SphericalSL, classname);
+      PYBIND11_OVERRIDE(std::string, Spherical, classname);
     }
 
     const std::string harmonic() override {
-      PYBIND11_OVERRIDE(std::string, SphericalSL, harmonic);
+      PYBIND11_OVERRIDE(std::string, Spherical, harmonic);
+    }
+
+    void get_pot(Eigen::MatrixXd& tab, double x) override {
+      PYBIND11_OVERRIDE_PURE(void, Spherical, get_pot, tab, x);
+    }
+    
+    void get_dens(Eigen::MatrixXd& tab, double x) override {
+      PYBIND11_OVERRIDE_PURE(void, Spherical, get_dens, tab, x);
     }
 
+    void get_force(Eigen::MatrixXd& tab, double x) override {
+      PYBIND11_OVERRIDE_PURE(void, Spherical, get_force, tab, x);
+    }
+
+
   public:
 
     // Inherit the constructors
-    using SphericalSL::SphericalSL;
+    using Spherical::Spherical;
 
     std::vector<double> getFields(double x, double y, double z) override {
-      PYBIND11_OVERRIDE(std::vector<double>, SphericalSL, getFields, x, y, z);
+      PYBIND11_OVERRIDE(std::vector<double>, Spherical, getFields, x, y, z);
     }
 
     void accumulate(double x, double y, double z, double mass) override {
-      PYBIND11_OVERRIDE(void, SphericalSL, accumulate, x, y, z, mass);
+      PYBIND11_OVERRIDE(void, Spherical, accumulate, x, y, z, mass);
     }
 
     void reset_coefs(void) override {
-      PYBIND11_OVERRIDE(void, SphericalSL, reset_coefs,);
+      PYBIND11_OVERRIDE(void, Spherical, reset_coefs,);
     }
 
     void make_coefs(void) override {
-      PYBIND11_OVERRIDE(void, SphericalSL, make_coefs,);
+      PYBIND11_OVERRIDE(void, Spherical, make_coefs,);
+    }
+
+    std::vector<Eigen::MatrixXd> orthoCheck(int knots) override {
+      PYBIND11_OVERRIDE_PURE(std::vector<Eigen::MatrixXd>, Spherical, orthoCheck, knots);
     }
 
   };
 
+
   class PyCylindrical : public Cylindrical
   {
   protected:
@@ -809,7 +830,7 @@ PYBIND11_OVERRIDE_PURE(void, Basis, addFromArray, m, p, roundrobin, posvelrows);
 
          Returns
          -------
-         CoefStructure
+         CoefStruct
              the coefficient structure created from the particles
 
          See also
@@ -1042,7 +1063,7 @@ PYBIND11_OVERRIDE_PURE(void, Basis, addFromArray, m, p, roundrobin, posvelrows);
          R"(
          Set the coordinate system for force evaluations.  The natural 
          coordinates for the basis class are the default; spherical
-         coordinates for SphericalSL, cylindrical coordinates for
+         coordinates for SphericalSL and Bessel, cylindrical coordinates for
          Cylindrical and FlatDisk, and Cartesian coordinates for the Slab
          and Cube. This member function can be used to override the default.  
          The available coorindates are: 'spherical', 'cylindrical', 
@@ -1131,37 +1152,135 @@ PYBIND11_OVERRIDE_PURE(void, Basis, addFromArray, m, p, roundrobin, posvelrows);
          None
          )", py::arg("coefs"));
 
-    py::class_<BasisClasses::SphericalSL, std::shared_ptr<BasisClasses::SphericalSL>, PySphericalSL, BasisClasses::BiorthBasis>(m, "SphericalSL")
-      .def(py::init<const std::string&>(),
+  py::class_<BasisClasses::Cylindrical, std::shared_ptr<BasisClasses::Cylindrical>, PyCylindrical, BasisClasses::BiorthBasis>(m, "Cylindrical")
+    .def(py::init<const std::string&>(),
 	 R"(
-         Create a spherical Sturm-Liouville basis
+	 Create a cylindrical EOF basis
 
          Parameters
          ----------
          YAMLstring : str
-             The YAML configuration for the spherical basis
+             The YAML configuration for the cylindrical basis
 
          Returns
          -------
-         SphericalSL
+         Cylindrical
               the new instance
          )", py::arg("YAMLstring"))
+    .def("getBasis", &BasisClasses::Cylindrical::getBasis,
+	 R"(
 
-      .def("getBasis", &BasisClasses::SphericalSL::getBasis,
+         Evaluate basis on grid for visualization
+
+         Evaluate the potential-density basis functions on a linearly spaced
+	 2d-grid for inspection.  The structure is a two-grid of dimension
+	 lmax by nmax each pointing to a dictionary of 2-d arrays ('potential',
+	 'density', 'rforce', 'zforce') of dimension numr X numz.
+
+         Parameters
+         ----------
+         xmin : float, default=0.0
+              minimum value in mapped radius
+         xmax : float, default=1.0
+              maximum value in mapped radius
+         numr : int, default=40
+              number of linearly-space evaluation points in radius
+         zmin : float, default=-0.1
+              minimum value in vertical height
+         zmax : float, default=0.1
+              maximum value in vertical height
+         numz : int, default=40
+              number of linearly-space evaluation points in height
+         linear : bool, default=True
+              use linear spacing
+
+         Returns
+         -------
+         list(list(dict))
+             dictionaries of basis functions as lists indexed by m, n
+         )",
+	 py::arg("xmin")=0.0,
+	 py::arg("xmax")=1.0,
+	 py::arg("numr")=40,
+	 py::arg("zmin")=-0.1,
+	 py::arg("zmax")=0.1,
+	 py::arg("numz")=40,
+	 py::arg("linear")=true)
+    // The following member needs to be a lambda capture because
+    // orthoCheck is not in the base class and needs to have different
+    // parameters depending on the basis type.  Here, the quadrature
+    // is determined by the scale of the meridional grid.
+    .def("orthoCheck", [](BasisClasses::Cylindrical& A)
+	 {
+	   return A.orthoCheck();
+	 },
+	R"(
+        Check orthgonality of basis functions by quadrature
+
+        Inner-product matrix of Sturm-Liouville solutions indexed by
+        harmonic order used to assess fidelity.
+
+        Parameters
+        ----------
+        knots : int
+            Number of quadrature knots
+
+        Returns
+        -------
+        list(numpy.ndarray)
+	    list of numpy.ndarrays from [0, ... , Mmax]
+        )")
+    .def_static("cacheInfo", [](std::string cachefile)
+    {
+      return BasisClasses::Cylindrical::cacheInfo(cachefile);
+    },
+      R"(
+      Report the parameters in a basis cache file and return a dictionary
+
+      Parameters
+      ----------
+      cachefile : str
+          name of cache file
+
+      Returns
+      -------
+      dict({tag: value},...)
+          cache parameters
+      )",
+      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"(
+         Create a spherical basis
+
+         Parameters
+         ----------
+         YAMLstring : str
+             The YAML configuration for the spherical basis
+         ForceID : str
+             The string identifier for this force type
+
+         Returns
+         -------
+         Spherical
+              the new instance
+         )", py::arg("YAMLstring"), py::arg("ForceID"))
+      .def("getBasis", &BasisClasses::Spherical::getBasis,
 	   R"(
            Get basis functions
 
-	   Evaluate the potential-density basis functions on a logarithmically
+	   Evaluate the potential-density basis functions on a linearly
 	   spaced grid for inspection. The structure is a two-grid of dimension
 	   lmax by nmax each pointing to a dictionary of 1-d arrays ('potential',
 	   'density', 'rforce') of dimension numr.
 
            Parameters
            ----------
-           logxmin : float, default=-3.0
-                minimum mapped radius in log10 units
-           logxmax : float, default=0.5
-                maximum mapped radius in log10 units
+           rmin : float, default=0.0
+                minimum radius
+           rmax : float, default=1.0
+                maximum radius
            numr : int, default=400
                 number of equally spaced output points
 
@@ -1170,14 +1289,14 @@ PYBIND11_OVERRIDE_PURE(void, Basis, addFromArray, m, p, roundrobin, posvelrows);
            list(list(dict))
                dictionaries of basis functions as lists indexed by l, n
            )",
-	   py::arg("logxmin")=-3.0,
-	   py::arg("logxmax")=0.5,
+	   py::arg("rmin")=0.0,
+	   py::arg("rmax")=1.0,
 	   py::arg("numr")=400)
       // The following member needs to be a lambda capture because
       // orthoCheck is not in the base class and needs to have
       // different parameters depending on the basis type.  Here the
       // user can and will often need to specify a quadrature value.
-      .def("orthoCheck", [](BasisClasses::SphericalSL& A, int knots)
+      .def("orthoCheck", [](BasisClasses::Spherical& A, int knots)
       {
 	return A.orthoCheck(knots);
       },
@@ -1200,7 +1319,7 @@ PYBIND11_OVERRIDE_PURE(void, Basis, addFromArray, m, p, roundrobin, posvelrows);
 	py::arg("knots")=40)
       .def_static("cacheInfo", [](std::string cachefile)
       {
-	return BasisClasses::SphericalSL::cacheInfo(cachefile);
+	return BasisClasses::Spherical::cacheInfo(cachefile);
       },
 	R"(
         Report the parameters in a basis cache file and return a dictionary
@@ -1261,68 +1380,98 @@ PYBIND11_OVERRIDE_PURE(void, Basis, addFromArray, m, p, roundrobin, posvelrows);
       )", py::arg("I"));
 
   
-  py::class_<BasisClasses::Cylindrical, std::shared_ptr<BasisClasses::Cylindrical>, PyCylindrical, BasisClasses::BiorthBasis>(m, "Cylindrical")
-    .def(py::init<const std::string&>(),
+    py::class_<BasisClasses::SphericalSL, std::shared_ptr<BasisClasses::SphericalSL>, BasisClasses::Spherical>(m, "SphericalSL")
+      .def(py::init<const std::string&>(),
 	 R"(
-	 Create a cylindrical EOF basis
+         Create a spherical Sturm-Liouville basis
 
          Parameters
          ----------
          YAMLstring : str
-             The YAML configuration for the cylindrical basis
+             The YAML configuration for the spherical basis
 
          Returns
          -------
-         Cylindrical
+         SphericalSL
               the new instance
          )", py::arg("YAMLstring"))
-    .def("getBasis", &BasisClasses::Cylindrical::getBasis,
-	 R"(
 
-         Evaluate basis on grid for visualization
+      .def("getBasis", &BasisClasses::SphericalSL::getBasis,
+	   R"(
+           Get basis functions
 
-         Evaluate the potential-density basis functions on a linearly spaced
-	 2d-grid for inspection.  The structure is a two-grid of dimension
-	 lmax by nmax each pointing to a dictionary of 2-d arrays ('potential',
-	 'density', 'rforce', 'zforce') of dimension numr X numz.
+	   Evaluate the potential-density basis functions on a logarithmically
+	   spaced grid for inspection. The structure is a two-grid of dimension
+	   lmax by nmax each pointing to a dictionary of 1-d arrays ('potential',
+	   'density', 'rforce') of dimension numr.
+
+           Parameters
+           ----------
+           logxmin : float, default=-3.0
+                minimum mapped radius in log10 units
+           logxmax : float, default=0.5
+                maximum mapped radius in log10 units
+           numr : int, default=400
+                number of equally spaced output points
+
+           Returns
+           -------
+           list(list(dict))
+               dictionaries of basis functions as lists indexed by l, n
+           )",
+	   py::arg("logxmin")=-3.0,
+	   py::arg("logxmax")=0.5,
+	   py::arg("numr")=400)
+      // The following member needs to be a lambda capture because
+      // orthoCheck is not in the base class and needs to have
+      // different parameters depending on the basis type.  Here the
+      // user can and will often need to specify a quadrature value.
+      .def("orthoCheck", [](BasisClasses::SphericalSL& A, int knots)
+      {
+	return A.orthoCheck(knots);
+      },
+	R"(
+        Check orthgonality of basis functions by quadrature
+
+        Inner-product matrix of Sturm-Liouville solutions indexed by
+        harmonic order used to assess fidelity.
+
+        Parameters
+        ----------
+        knots : int, default=40
+            Number of quadrature knots
+
+        Returns
+        -------
+        list(numpy.ndarray)
+	    list of numpy.ndarrays from [0, ... , Lmax]
+        )",
+	py::arg("knots")=40);
+
+  
+    py::class_<BasisClasses::Bessel, std::shared_ptr<BasisClasses::Bessel>, BasisClasses::Spherical>(m, "Bessel")
+      .def(py::init<const std::string&>(),
+	 R"(
+         Create a spherical Bessel-function basis
 
          Parameters
          ----------
-         xmin : float, default=0.0
-              minimum value in mapped radius
-         xmax : float, default=1.0
-              maximum value in mapped radius
-         numr : int, default=40
-              number of linearly-space evaluation points in radius
-         zmin : float, default=-0.1
-              minimum value in vertical height
-         zmax : float, default=0.1
-              maximum value in vertical height
-         numz : int, default=40
-              number of linearly-space evaluation points in height
-         linear : bool, default=True
-              use linear spacing
+         YAMLstring : str
+             The YAML configuration for the spherical Bessel-function basis
 
          Returns
          -------
-         list(list(dict))
-             dictionaries of basis functions as lists indexed by m, n
-         )",
-	 py::arg("xmin")=0.0,
-	 py::arg("xmax")=1.0,
-	 py::arg("numr")=40,
-	 py::arg("zmin")=-0.1,
-	 py::arg("zmax")=0.1,
-	 py::arg("numz")=40,
-	 py::arg("linear")=true)
-    // The following member needs to be a lambda capture because
-    // orthoCheck is not in the base class and needs to have different
-    // parameters depending on the basis type.  Here, the quadrature
-    // is determined by the scale of the meridional grid.
-    .def("orthoCheck", [](BasisClasses::Cylindrical& A)
-	 {
-	   return A.orthoCheck();
-	 },
+         Bessel
+              the new instance
+         )", py::arg("YAMLstring"))
+      // The following member needs to be a lambda capture because
+      // orthoCheck is not in the base class and needs to have
+      // different parameters depending on the basis type.  Here the
+      // user can and will often need to specify a quadrature value.
+      .def("orthoCheck", [](BasisClasses::Bessel& A, int knots)
+      {
+	return A.orthoCheck(knots);
+      },
 	R"(
         Check orthgonality of basis functions by quadrature
 
@@ -1331,32 +1480,15 @@ PYBIND11_OVERRIDE_PURE(void, Basis, addFromArray, m, p, roundrobin, posvelrows);
 
         Parameters
         ----------
-        knots : int
+        knots : int, default=40
             Number of quadrature knots
 
         Returns
         -------
         list(numpy.ndarray)
-	    list of numpy.ndarrays from [0, ... , Mmax]
-        )")
-    .def_static("cacheInfo", [](std::string cachefile)
-    {
-      return BasisClasses::Cylindrical::cacheInfo(cachefile);
-    },
-      R"(
-      Report the parameters in a basis cache file and return a dictionary
-
-      Parameters
-      ----------
-      cachefile : str
-          name of cache file
-
-      Returns
-      -------
-      dict({tag: value},...)
-          cache parameters
-      )",
-      py::arg("cachefile"));
+	    list of numpy.ndarrays from [0, ... , Lmax]
+        )",
+	py::arg("knots")=40);
 
   py::class_<BasisClasses::FlatDisk, std::shared_ptr<BasisClasses::FlatDisk>, PyFlatDisk, BasisClasses::BiorthBasis>(m, "FlatDisk")
     .def(py::init<const std::string&>(),
@@ -1714,7 +1846,7 @@ PYBIND11_OVERRIDE_PURE(void, Basis, addFromArray, m, p, roundrobin, posvelrows);
 
          Returns
          -------
-         CoefStructure
+         CoefStruct
              the coefficient structure created from the particles
 
          See also
@@ -1920,5 +2052,4 @@ PYBIND11_OVERRIDE_PURE(void, Basis, addFromArray, m, p, roundrobin, posvelrows);
 	py::arg("tinit"), py::arg("tfinal"), py::arg("h"),
 	py::arg("ps"), py::arg("basiscoef"), py::arg("func"),
 	py::arg("nout")=std::numeric_limits<int>::max());
-
 }
diff --git a/src/Bessel.cc b/src/Bessel.cc
index 657a9d1b4..eedf4592f 100644
--- a/src/Bessel.cc
+++ b/src/Bessel.cc
@@ -51,7 +51,7 @@ Bessel::Bessel(Component* c0, const YAML::Node& conf, MixtureBasis* m) : Spheric
 void Bessel::get_dpotl(int lmax, int nmax, double r, 
 		       Eigen::MatrixXd& p, Eigen::MatrixXd& dp, int tid)
 {
-  int klo = (int)( (r-r_grid[1])/r_grid_del ) + 1;
+  int klo = (int)( (r-r_grid[0])/r_grid_del );
   if (klo < 0) klo = 0;
   if (klo > RNUM - 2) klo = RNUM - 2;
   int khi = klo + 1;
@@ -76,7 +76,7 @@ void Bessel::get_dpotl(int lmax, int nmax, double r,
 
 void Bessel::get_potl(int lmax, int nmax, double r, Eigen::MatrixXd& p, int tid)
 {
-  int klo = (int)( (r-r_grid[1])/r_grid_del ) + 1;
+  int klo = (int)( (r-r_grid[0])/r_grid_del );
   if (klo < 0) klo = 0;
   if (klo > RNUM - 2) klo = RNUM - 2;
   int khi = klo + 1;
@@ -97,7 +97,7 @@ void Bessel::get_potl(int lmax, int nmax, double r, Eigen::MatrixXd& p, int tid)
 
 void Bessel::get_dens(int lmax, int nmax, double r, Eigen::MatrixXd& p, int tid)
 {
-  int klo = (int)( (r-r_grid[1])/r_grid_del ) + 1;
+  int klo = (int)( (r-r_grid[0])/r_grid_del );
   if (klo < 0) klo = 0;
   if (klo > RNUM - 2) klo = RNUM - 2;
   int khi = klo + 1;
@@ -120,7 +120,7 @@ void Bessel::get_dens(int lmax, int nmax, double r, Eigen::MatrixXd& p, int tid)
 void Bessel::get_potl_dens(int lmax, int nmax, double r, 
 			   Eigen::MatrixXd& p, Eigen::MatrixXd& d, int tid)
 {
-  int klo = (int)( (r-r_grid[1])/r_grid_del ) + 1;
+  int klo = (int)( (r-r_grid[0])/r_grid_del );
   if (klo < 0) klo = 0;
   if (klo > RNUM - 2) klo = RNUM - 2;
   int khi = klo + 1;

From 10c01fe7cafe1678a0876bf6d12f15c6237eb5fa Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Sun, 12 May 2024 22:23:18 -0400
Subject: [PATCH 12/25] Check in missing Bessel implementation

---
 expui/BiorthBess.H  |  93 ++++++++++++++++
 expui/BiorthBess.cc | 251 ++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 344 insertions(+)
 create mode 100644 expui/BiorthBess.H
 create mode 100644 expui/BiorthBess.cc

diff --git a/expui/BiorthBess.H b/expui/BiorthBess.H
new file mode 100644
index 000000000..ed0dcb57e
--- /dev/null
+++ b/expui/BiorthBess.H
@@ -0,0 +1,93 @@
+#ifndef _BiorthBess_H_
+#define _BiorthBess_H_
+
+#include <vector>
+#include <memory>
+#include <cmath>
+
+#include <Eigen/Eigen>
+
+Eigen::VectorXd sbessjz(int n, int m);
+
+class BiorthBess
+{
+private:
+
+  //! Grid to hold tabulated basis
+  class RGrid
+  {
+  public:
+    Eigen::MatrixXd rw;
+    Eigen::MatrixXd rw2;
+    int nmax;
+  };
+  
+  //@{
+  //! Grid storage and parameters
+  std::vector<RGrid> dens_grid, potl_grid;
+  Eigen::VectorXd r_grid;
+  double r_grid_del;
+  //@}
+
+  //! Cache roots for spherical Bessel functions
+  class Roots
+  {
+  public:
+
+    int l;
+    int n;
+
+    Eigen::VectorXd a;
+
+    Roots(int L, int nmax) : l(L), n(nmax) {
+      a = sbessjz(l, n);
+    }
+
+    ~Roots() {}
+  };
+
+  //! Root database isntance
+  std::shared_ptr<Roots> p;
+
+  //@{
+  //! Density and potential members
+  double dens(double r, int n);
+  double potl(double r, int n);
+  //@}
+
+  //@{
+  //! Parameters
+  double rmin=0.0, rmax;
+  int lmax, nmax, RNUM=2000;
+  //@}
+
+  //! Make the density and potential grid
+  void make_grid();
+  
+public:
+
+  //! Constructor
+  BiorthBess(double rmax, int lmax, int nmax, int rnum);
+
+  //! Destructor
+  virtual ~BiorthBess() {}
+	
+  //! Potential evaluation
+  void get_potl(double r, Eigen::MatrixXd& p);
+
+  //! Density evaluation
+  void get_dens(double r, Eigen::MatrixXd& p);
+
+  //! Force evaluation
+  void get_dpotl(double r, Eigen::MatrixXd& p);
+
+  //! Potential and density evaluation
+  void get_potl_dens(double r, Eigen::MatrixXd& p, Eigen::MatrixXd& d);
+
+  //! Compute the orthogonality of the basis by returning inner
+  //! produce matrices
+  std::vector<Eigen::MatrixXd> orthoCheck(int knots=40);
+};
+
+#endif
+
diff --git a/expui/BiorthBess.cc b/expui/BiorthBess.cc
new file mode 100644
index 000000000..6114da63c
--- /dev/null
+++ b/expui/BiorthBess.cc
@@ -0,0 +1,251 @@
+#include <BiorthBess.H>
+#include <EXPmath.H>
+#include <EXPException.H>
+#include <interp.H>
+#include <gaussQ.H>
+
+
+BiorthBess::BiorthBess(double rmax, int lmax, int nmax, int RNUM) :
+  rmax(rmax), lmax(lmax), nmax(nmax), RNUM(RNUM)
+{
+  // Initialize radial grids
+  make_grid();
+}
+
+// Get potential functions by from table
+void BiorthBess::get_dpotl(double r, Eigen::MatrixXd& p)
+{
+  int klo = (int)( (r-r_grid[0])/r_grid_del );
+  if (klo < 0) klo = 0;
+  if (klo > RNUM - 2) klo = RNUM - 2;
+  int khi = klo + 1;
+
+  double a = (r_grid[khi] - r)/r_grid_del;
+  double b = (r - r_grid[klo])/r_grid_del;
+  
+  double aa = a*(a*a-1.0)*r_grid_del*r_grid_del/6.0;
+  double bb = b*(b*b-1.0)*r_grid_del*r_grid_del/6.0;
+  double aaa = -(3.0*a*a - 1.0)*r_grid_del/6.0;
+  double bbb =  (3.0*b*b - 1.0)*r_grid_del/6.0;
+
+  p.resize(lmax+1, nmax);
+
+  for (int l=0; l<=lmax; l++) {
+    for (int n=0; n<nmax; n++) {
+      p(l, n) = (-potl_grid[l].rw(n, klo)+potl_grid[l].rw(n, khi))/r_grid_del+
+	aaa*potl_grid[l].rw2(n, klo) + bbb*potl_grid[l].rw2(n, khi);
+    }
+  }
+}
+
+void BiorthBess::get_potl(double r, Eigen::MatrixXd& p)
+{
+  int klo = (int)( (r-r_grid[0])/r_grid_del );
+  if (klo < 0) klo = 0;
+  if (klo > RNUM - 2) klo = RNUM - 2;
+  int khi = klo + 1;
+
+  double a = (r_grid[khi] - r)/r_grid_del;
+  double b = (r - r_grid[klo])/r_grid_del;
+
+  double aa = a*(a*a-1.0)*r_grid_del*r_grid_del/6.0;
+  double bb = b*(b*b-1.0)*r_grid_del*r_grid_del/6.0;
+
+  p.resize(lmax+1, nmax);
+
+  for (int l=0; l<=lmax; l++) {
+    for (int n=0; n<nmax; n++) {
+      p(l, n) = a*potl_grid[l].rw(n, klo) + b*potl_grid[l].rw(n, khi) +
+	aa*potl_grid[l].rw2(n, klo) + bb*potl_grid[l].rw2(n, khi);
+    }
+  }
+}
+
+void BiorthBess::get_dens(double r, Eigen::MatrixXd& p)
+{
+  int klo = (int)( (r-r_grid[0])/r_grid_del );
+  if (klo < 0) klo = 0;
+  if (klo > RNUM - 2) klo = RNUM - 2;
+  int khi = klo + 1;
+
+  double a = (r_grid[khi] - r)/r_grid_del;
+  double b = (r - r_grid[klo])/r_grid_del;
+
+  double aa = a*(a*a-1.0)*r_grid_del*r_grid_del/6.0;
+  double bb = b*(b*b-1.0)*r_grid_del*r_grid_del/6.0;
+
+  p.resize(lmax+1, nmax);
+
+  for (int l=0; l<=lmax; l++) {
+    for (int n=0; n<nmax; n++) {
+      p(l, n) = a*dens_grid[l].rw(n, klo) + b*dens_grid[l].rw(n, khi) +
+	aa*dens_grid[l].rw2(n, klo) + bb*dens_grid[l].rw2(n, khi);
+    }
+  }
+
+  p *= -1.0;
+}
+
+
+void BiorthBess::get_potl_dens(double r, Eigen::MatrixXd& p, Eigen::MatrixXd& d)
+{
+  int klo = (int)( (r-r_grid[0])/r_grid_del );
+  if (klo < 0) klo = 0;
+  if (klo > RNUM - 2) klo = RNUM - 2;
+  int khi = klo + 1;
+
+  double a = (r_grid[khi] - r)/r_grid_del;
+  double b = (r - r_grid[klo])/r_grid_del;
+
+  double aa = a*(a*a-1.0)*r_grid_del*r_grid_del/6.0;
+  double bb = b*(b*b-1.0)*r_grid_del*r_grid_del/6.0;
+
+  p.resize(lmax+1, nmax);
+  d.resize(lmax+1, nmax);
+
+  for (int l=0; l<=lmax; l++) {
+    for (int n=0; n<nmax; n++) {
+      p(l, n) = a*potl_grid[l].rw(n, klo) + b*potl_grid[l].rw(n, khi) +
+	aa*potl_grid[l].rw2(n, klo) + bb*potl_grid[l].rw2(n, khi);
+      d(l, n) = a*dens_grid[l].rw(n, klo) + b*dens_grid[l].rw(n, khi) +
+	aa*dens_grid[l].rw2(n, klo) + bb*dens_grid[l].rw2(n, khi);
+    }
+  }
+
+  d *= -1.0;
+}
+
+double BiorthBess::dens(double r, int n)
+{
+  double alpha;
+
+  if (n>p->n)
+    throw GenericError("Routine dens() called with n out of bounds", __FILE__, __LINE__, 1001, true);
+
+  alpha = p->a[n];
+  return alpha*M_SQRT2/fabs(EXPmath::sph_bessel(p->l+1, alpha)) * pow(rmax,-2.5) *
+    EXPmath::sph_bessel(p->l, alpha*r/rmax);
+}
+
+double BiorthBess::potl(double r, int n)
+{
+  double alpha;
+
+  if (n>p->n)
+    throw GenericError("Routine potl() called with n out of bounds", __FILE__, __LINE__, 1002, true);
+
+  alpha = p->a[n];
+  return M_SQRT2/fabs(alpha*EXPmath::sph_bessel(p->l+1,alpha)) * pow(rmax,-0.5) *
+    EXPmath::sph_bessel(p->l,alpha*r/rmax);
+}
+
+void BiorthBess::make_grid()
+{
+  
+  potl_grid.resize(lmax+1);
+  dens_grid.resize(lmax+1);
+  
+  r_grid.resize(RNUM);
+
+  r_grid_del = rmax/(double)(RNUM-1);
+  double r = 0.0;
+  for (int ir=0; ir<RNUM; ir++, r+=r_grid_del) r_grid[ir] = r;
+
+  for (int l=0; l<=lmax; l++) {
+    potl_grid[l].rw .resize(nmax, RNUM);
+    potl_grid[l].rw2.resize(nmax, RNUM);
+    dens_grid[l].rw .resize(nmax, RNUM);
+    dens_grid[l].rw2.resize(nmax, RNUM);
+
+    p = std::make_shared<Roots>(l, nmax);
+
+    for (int n=0; n<nmax; n++) {
+      r = 0.0;
+      for (int ir=0; ir<RNUM; ir++, r+=r_grid_del) {
+	potl_grid[l].rw(n, ir) = potl(r, n);
+	dens_grid[l].rw(n, ir) = dens(r, n);
+      }
+      
+      {
+	Eigen::VectorXd X(potl_grid[l].rw.row(n));
+	Eigen::VectorXd Y(X.size());
+	Spline(r_grid, X, 1.0e30, 1.0e30, Y);
+	potl_grid[l].rw2.row(n) = Y;
+      }
+
+      {
+	Eigen::VectorXd X(dens_grid[l].rw.row(n));
+	Eigen::VectorXd Y(X.size());
+	Spline(r_grid, X, 1.0e30, 1.0e30, Y);
+	dens_grid[l].rw2.row(n) = Y;
+      }
+    }
+
+    potl_grid[l].nmax = nmax;
+    dens_grid[l].nmax = nmax;
+  }
+
+  // check table
+
+  for (int ir=0; ir<RNUM; ir++) assert(!std::isnan(r_grid[ir]));
+  for (int l=0; l<=lmax; l++) {
+    assert(!std::isnan(potl_grid[l].nmax));
+    assert(!std::isnan(dens_grid[l].nmax));
+    for (int n=0; n<nmax; n++) {
+      for (int ir=0; ir<RNUM; ir++) {
+	assert(!std::isnan(potl_grid[l].rw (n, ir)));
+	assert(!std::isnan(potl_grid[l].rw2(n, ir)));
+	assert(!std::isnan(dens_grid[l].rw (n, ir)));
+	assert(!std::isnan(dens_grid[l].rw2(n, ir)));
+      }
+    }
+  }
+  
+}
+
+std::vector<Eigen::MatrixXd> BiorthBess::orthoCheck(int num)
+{
+  // Allocate storage
+  //
+  std::vector<Eigen::MatrixXd> one(lmax+1);
+  for (auto & u : one) {
+    u.resize(nmax, nmax);
+    u.setZero();
+  }
+
+  // Number of knots
+  //
+  LegeQuad wk(num);
+  
+  // Radial range
+  //
+  double dr = rmax - rmin;
+
+  Eigen::MatrixXd p(lmax+1, nmax), d(lmax+1, nmax);
+
+
+  // Biorthogonal integral loop
+  //
+  for (int i=0; i<num; i++) {
+
+    double r = rmin + dr*wk.knot(i);
+    double w = dr*wk.weight(i) * r * r;
+
+    // Evaluate basis at radius r
+    //
+    get_potl(r, p);
+    get_dens(r, d);
+
+    // Contribution to the integrand
+    //
+    for (int L=0, cnt=0; L<=lmax; L++) {
+      for (int n1=0; n1<nmax; n1++) {
+	for (int n2=0; n2<nmax; n2++, cnt++) {
+	  one[L](n1, n2) += w * p(L, n1) * d(L, n2);
+	}
+      }
+    }
+  }
+
+  return one;
+}

From 6d876d861cdaedeeb2bc2787da92035b657d6713 Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Mon, 13 May 2024 07:34:59 -0400
Subject: [PATCH 13/25] Revert the density sign change because it breaks the
 orthoTest [no ci]

---
 expui/BiorthBess.cc | 4 ----
 1 file changed, 4 deletions(-)

diff --git a/expui/BiorthBess.cc b/expui/BiorthBess.cc
index 6114da63c..82a0c7593 100644
--- a/expui/BiorthBess.cc
+++ b/expui/BiorthBess.cc
@@ -82,8 +82,6 @@ void BiorthBess::get_dens(double r, Eigen::MatrixXd& p)
 	aa*dens_grid[l].rw2(n, klo) + bb*dens_grid[l].rw2(n, khi);
     }
   }
-
-  p *= -1.0;
 }
 
 
@@ -111,8 +109,6 @@ void BiorthBess::get_potl_dens(double r, Eigen::MatrixXd& p, Eigen::MatrixXd& d)
 	aa*dens_grid[l].rw2(n, klo) + bb*dens_grid[l].rw2(n, khi);
     }
   }
-
-  d *= -1.0;
 }
 
 double BiorthBess::dens(double r, int n)

From 2c79f85d14d67b3f1c218bda21047e334ba69b9e Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Tue, 14 May 2024 22:39:43 -0400
Subject: [PATCH 14/25] Enforce a minimum small integer step into sorted
 particle array to choose inner and outer bin edges to prevent empty bins at
 small and large radii [no ci]

---
 src/Sphere.H  |  1 +
 src/Sphere.cc | 17 +++++++++++++++--
 2 files changed, 16 insertions(+), 2 deletions(-)

diff --git a/src/Sphere.H b/src/Sphere.H
index 1178b996b..0ec8e9c7a 100644
--- a/src/Sphere.H
+++ b/src/Sphere.H
@@ -92,6 +92,7 @@ private:
   double tnext, dtime;
   int    numr;
   int    nums;
+  int    noff;
   int    cmap;
   int    diverge;
   double dfac;
diff --git a/src/Sphere.cc b/src/Sphere.cc
index df11809e9..6c3fda20b 100644
--- a/src/Sphere.cc
+++ b/src/Sphere.cc
@@ -14,6 +14,7 @@ Sphere::valid_keys = {
   "rmapping",
   "numr",
   "nums",
+  "noff",
   "cmap",
   "diverge",
   "dfac",
@@ -32,6 +33,7 @@ Sphere::Sphere(Component* c0, const YAML::Node& conf, MixtureBasis* m) :
   rmap = 0.067*rmax;
   numr = 2000;
   nums = 2000;
+  noff = 32;
   cmap = 1;
   diverge = 0;
   dfac = 1.0;
@@ -108,6 +110,7 @@ void Sphere::initialize()
     if (conf["rmapping"])  rmap       = conf["rmapping"].as<double>();
     if (conf["numr"])      numr       = conf["numr"].as<int>();
     if (conf["nums"])      nums       = conf["nums"].as<int>();
+    if (conf["noff"])      noff       = conf["noff"].as<int>();
     if (conf["cmap"])      cmap       = conf["cmap"].as<int>();
     if (conf["diverge"])   diverge    = conf["diverge"].as<int>();
     if (conf["dfac"])      dfac       = conf["dfac"].as<double>();
@@ -251,8 +254,18 @@ void Sphere::make_model_bin()
   // Make mass array
   //
   int    numR = std::min<int>(sqrt(Ntot), nums);
-  double Rmin = std::max<double>(RM.begin()->first, rmin);
-  double Rmax = RM.rbegin()->first;
+
+  // Compute offsets into mass array for bin edges
+  //
+  int   nstep = std::max<int>(noff, std::floor(0.5*Ntot/numR));
+
+  auto   ibeg = RM.begin();
+  auto   iend = RM.rbegin();
+  std::advance(ibeg, nstep);
+  std::advance(iend, nstep);
+
+  double Rmin = std::max<double>(ibeg->first, rmin);
+  double Rmax = iend->first;
   bool   logR = false;
 
   if (logr and Rmin>0.0) {

From 80cfca0ec9efb04dc54096dfe6e44cf431ece686 Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Thu, 16 May 2024 13:41:19 -0400
Subject: [PATCH 15/25] Add runtag to model-build diagnostics [no ci]

---
 src/Sphere.cc | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/src/Sphere.cc b/src/Sphere.cc
index 6c3fda20b..cf1dc514e 100644
--- a/src/Sphere.cc
+++ b/src/Sphere.cc
@@ -316,7 +316,8 @@ void Sphere::make_model_bin()
   //
   if (myid==0) {
     static int cnt = 0;
-    std::ostringstream sout; sout << "SphereMassModel." << cnt++;
+    std::ostringstream sout;
+    sout << runtag << ".SphereMassModel." << cnt++;
     std::ofstream dbg(sout.str());
     if (dbg.good()) {
       for (int i=0; i<numR; i++)
@@ -337,6 +338,11 @@ void Sphere::make_model_bin()
   // Regenerate Sturm-Liouville grid
   //
   std::string cachename = outdir  + cache_file;
+  
+  if (logR) {
+    Rmin = exp(Rmin);
+    Rmax = exp(Rmax);
+  }
   ortho = std::make_shared<SLGridSph>(mod, Lmax, nmax, numR, Rmin, Rmax, false, 1, 1.0, cachename);
 
   // Test for basis consistency (will generate an exception if maximum

From 3fb2284760093c6921944511ef9380c57bf61ed9 Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Thu, 16 May 2024 13:42:28 -0400
Subject: [PATCH 16/25] Add additional diagnostic output; fix use_mpi logic to
 prevent MPI_Finalize calls [no ci]

---
 utils/SL/slcheck.cc | 68 +++++++++++++++++++++++----------------------
 1 file changed, 35 insertions(+), 33 deletions(-)

diff --git a/utils/SL/slcheck.cc b/utils/SL/slcheck.cc
index 4a28596e8..f182b8d8e 100644
--- a/utils/SL/slcheck.cc
+++ b/utils/SL/slcheck.cc
@@ -15,8 +15,8 @@
 
 int main(int argc, char** argv)
 {
-  bool use_mpi, use_logr;
-  double rmin, rmax, rs, dfac;
+  bool use_mpi, use_logr, diag;
+  double rmin, rmax, rmapping, dfac;
   int numr, cmap, diverge, Lmax, nmax;
   string filename, cachefile;
 
@@ -27,43 +27,38 @@ int main(int argc, char** argv)
   cxxopts::Options options(argv[0], "Check the consistency a spherical SL basis");
 
   options.add_options()
-   ("h,help", "Print this help message")
-   ("mpi", "using parallel computation",
+    ("h,help", "Print this help message")
+    ("d,diag", "Turn on Sledge diagnostics for SL computation",
+     cxxopts::value<bool>(diag)->default_value("false"))
+    ("mpi", "using parallel computation",
      cxxopts::value<bool>(use_mpi)->default_value("false"))
-   ("logr", "logarithmic spacing for orthogonality check",
+    ("logr", "logarithmic spacing for orthogonality check",
      cxxopts::value<bool>(use_logr)->default_value("false"))
-   ("cmap", "coordinates in SphereSL: use mapped (1) or linear(0) coordinates",
+    ("cmap", "coordinates in SphereSL: use mapped (1) or linear(0) coordinates",
      cxxopts::value<int>(cmap)->default_value("1"))
-   ("Lmax", "maximum number of angular harmonics in the expansion",
+    ("Lmax", "maximum number of angular harmonics in the expansion",
      cxxopts::value<int>(Lmax)->default_value("2"))
-   ("nmax", "maximum number of radial harmonics in the expansion",
+    ("nmax", "maximum number of radial harmonics in the expansion",
      cxxopts::value<int>(nmax)->default_value("10"))
-   ("numr", "radial knots for the SL grid",
+    ("numr", "radial knots for the SL grid",
      cxxopts::value<int>(numr)->default_value("1000"))
-   ("rmin", "minimum radius for the SL grid",
+    ("rmin", "minimum radius for the SL grid",
      cxxopts::value<double>(rmin)->default_value("-1.0"))
-   ("rmax", "maximum radius for the SL grid",
+    ("rmax", "maximum radius for the SL grid",
      cxxopts::value<double>(rmax)->default_value("-1.0"))
-   ("rs", "cmap scale factor",
-     cxxopts::value<double>(rs)->default_value("0.067"))
-   ("diverge", "cusp divergence for spherical model",
+    ("rmapping", "cmap scale factor",
+     cxxopts::value<double>(rmapping)->default_value("0.067"))
+    ("diverge", "cusp divergence for spherical model",
      cxxopts::value<int>(diverge)->default_value("0"))
-   ("dfac", "cusp divergence exponent for spherical model",
+    ("dfac", "cusp divergence exponent for spherical model",
      cxxopts::value<double>(dfac)->default_value("1.0"))
-   ("filename", "model file",
+    ("filename", "model file",
      cxxopts::value<string>(filename)->default_value("SLGridSph.model"))
-   ("cache", "cache file",
+    ("cache", "cache file",
      cxxopts::value<string>(cachefile)->default_value(".slgrid_sph_cache"))
     ;
 
   
-  //===================
-  // MPI preliminaries 
-  //===================
-  if (use_mpi) {
-    local_init_mpi(argc, argv);
-  }
-
   //===================
   // Parse options
   //===================
@@ -78,6 +73,13 @@ int main(int argc, char** argv)
     return 2;
   }
 
+  //===================
+  // MPI preliminaries 
+  //===================
+  if (use_mpi) {
+    local_init_mpi(argc, argv);
+  }
+
   // Print help message and exit
   //
   if (vm.count("help")) {
@@ -109,14 +111,14 @@ int main(int argc, char** argv)
   //
   try {
     ortho = std::make_shared<SLGridSph>(filename, Lmax, nmax, numr, rmin, rmax, 
-					   true, cmap, rs, 0, 1.0, cachefile, false);
-    //                                     ^               ^                  ^
-    //                                     |               |                  |
-    // Use cache file----------------------+               |                  |
-    //                                                     |                  |
-    // Cusp extrapolation----------------------------------+                  |
-    //                                                                        |
-    // Turn on diagnostic output in SL creation-------------------------------+
+					   true, cmap, rmapping, 0, 1.0, cachefile, diag);
+    //                                     ^                     ^                  ^
+    //                                     |                     |                  |
+    // Use cache file----------------------+                     |                  |
+    //                                                           |                  |
+    // Cusp extrapolation----------------------------------------+                  |
+    //                                                                              |
+    // Turn on diagnostic output in SL creation-------------------------------------+
   }
   catch (const EXPException& err) {
     if (myid==0) {
@@ -131,7 +133,7 @@ int main(int argc, char** argv)
   }
 
   if (bad) {
-    MPI_Finalize();
+    if (use_mpi) MPI_Finalize();
     exit(0);
   }
 				// Do what?

From a9ef74b9a68ddb741e2cb98ef78908660d114d2a Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Thu, 16 May 2024 13:42:50 -0400
Subject: [PATCH 17/25] Fix fencepost error in grid generation assignment [no
 ci]

---
 exputil/realize_model.cc | 1 +
 1 file changed, 1 insertion(+)

diff --git a/exputil/realize_model.cc b/exputil/realize_model.cc
index ab3b0264f..0adfc7b0a 100644
--- a/exputil/realize_model.cc
+++ b/exputil/realize_model.cc
@@ -1093,6 +1093,7 @@ Eigen::VectorXd SphericalModelMulti::gen_point(int& ierr)
       ibeg[n] = dN*n;
       iend[n] = dN*(n+1);
     }
+    iend[numprocs-1] = gen_N;
 
     std::vector<double> gen_emax(gen_N, 0.0), gen_vmax(gen_N, 0.0);
 

From 1fcc289605f890bf87473eed194c6418dff06faa Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Fri, 17 May 2024 17:52:10 -0400
Subject: [PATCH 18/25] Fix grid asymmetry error; very minor [no ci]

---
 expui/FieldGenerator.cc | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/expui/FieldGenerator.cc b/expui/FieldGenerator.cc
index 3f1631a0b..5e00e1c82 100644
--- a/expui/FieldGenerator.cc
+++ b/expui/FieldGenerator.cc
@@ -93,7 +93,7 @@ namespace Field
     // Compute the probe length
     //
     std::vector<double> dd(3);
-    for (int k=0; k<3; k++) dd[k] = (end[k] - beg[k])/num;
+    for (int k=0; k<3; k++) dd[k] = (end[k] - beg[k])/(num-1);
     double dlen = sqrt(dd[0]*dd[0] + dd[1]*dd[1] + dd[2]*dd[2]);
 
     for (int icnt=0; icnt<times.size(); icnt++) {

From 60a0140e8bc718ddd527a4efe21e6437ad3c5ebf Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Mon, 20 May 2024 14:50:31 -0400
Subject: [PATCH 19/25] Fixed mistake in exception that prevented message from
 showing; fixed mistake in default value for Coef::EvenOddPower() [no ci]

---
 expui/Coefficients.cc | 18 +++++++++++-------
 pyEXP/CoefWrappers.cc |  2 +-
 2 files changed, 12 insertions(+), 8 deletions(-)

diff --git a/expui/Coefficients.cc b/expui/Coefficients.cc
index 17f642933..78c3f2827 100644
--- a/expui/Coefficients.cc
+++ b/expui/Coefficients.cc
@@ -1085,17 +1085,21 @@ namespace CoefClasses
 	  node = YAML::Load(getYAML());
 	}
 	catch (const std::runtime_error& error) {
-	  std::cout << "CylCoefs::EvenOddPower: found a problem while loading "
-		    << "the YAML config" << std::endl;
-	  throw;
+	  throw std::runtime_error
+	    (
+	     "CylCoefs::EvenOddPower: found a problem while loading the "
+	     "YAML config"
+	     );
 	}
 
 	if (node["ncylodd"]) nodd = node["ncylodd"].as<int>();
 	else {
-	  std::cout << "CylCoefs::EvenOddPower: ncylodd is not in the YAML "
-		    << "config stanza.  Please specify this explicitly as "
-		    << "the first argument to EvenOddPower()" << std::endl;
-	  throw;
+	  throw std::runtime_error
+	    (
+	    "CylCoefs::EvenOddPower: ncylodd is not in the YAML config "
+	    "stanza.  Please specify this explicitly as the first argument "
+	    "to EvenOddPower()"
+	     );
 	}
       }
 
diff --git a/pyEXP/CoefWrappers.cc b/pyEXP/CoefWrappers.cc
index 477718c79..b1ef1d91e 100644
--- a/pyEXP/CoefWrappers.cc
+++ b/pyEXP/CoefWrappers.cc
@@ -1302,7 +1302,7 @@ void CoefficientClasses(py::module &m) {
          configuration. If in doubt, use the default.
          )",
 	 py::arg("nodd")=-1, py::arg("min")=0,
-	 py::arg("max") = std::numeric_limits<double>::max());
+	 py::arg("max")=std::numeric_limits<int>::max());
 
 
   py::class_<CoefClasses::SphFldCoefs, std::shared_ptr<CoefClasses::SphFldCoefs>, CoefClasses::Coefs>

From be9aaee58b9c467a4a5c6bcca58e08b2598ae093 Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Wed, 22 May 2024 12:20:37 -0400
Subject: [PATCH 20/25] Add a version tag to HDF5 coefficient sets and use that
 to toggle backward compatibility with the old (buggy) HighFive Eigen wrappers

---
 expui/Coefficients.H  |  9 +++++++++
 expui/Coefficients.cc | 36 ++++++++++++++++++++++++++++++++++++
 pyEXP/CoefWrappers.cc | 13 +++++++++++++
 3 files changed, 58 insertions(+)

diff --git a/expui/Coefficients.H b/expui/Coefficients.H
index 6cc6c0df8..737881916 100644
--- a/expui/Coefficients.H
+++ b/expui/Coefficients.H
@@ -54,6 +54,9 @@ namespace CoefClasses
     //! Verbose debugging output
     bool verbose;
     
+    //! Backward compatibility flag for HighFive
+    static bool H5BackCompat;
+
     //! Time vector
     std::vector<double> times;
 
@@ -64,6 +67,9 @@ namespace CoefClasses
     //! Get the YAML config for the basis (to be implemented by EXP)
     virtual std::string getYAML() = 0;
     
+    //! Coefficient file versioning
+    inline static const std::string Version = "1.0";
+
     //! Write parameter attributes (needed for derived classes)
     virtual void WriteH5Params(HighFive::File& file) = 0;
     
@@ -195,6 +201,9 @@ namespace CoefClasses
     //! Set maximum grid interpolation offset
     void setDeltaT(double dT) { deltaT = dT; }
 
+    //! Override backward compatibility for HighFive
+    static void setNewH5() { H5BackCompat = false; }
+
     class CoefsError : public std::runtime_error
     {
     public:
diff --git a/expui/Coefficients.cc b/expui/Coefficients.cc
index 78c3f2827..227125751 100644
--- a/expui/Coefficients.cc
+++ b/expui/Coefficients.cc
@@ -20,6 +20,8 @@
 namespace CoefClasses
 {
   
+  bool Coefs::H5BackCompat = true;
+
   void Coefs::copyfields(std::shared_ptr<Coefs> p)
   {
     // These variables will copy data, not pointers
@@ -91,6 +93,9 @@ namespace CoefClasses
     file.getAttribute("geometry").read(geometry);
     file.getAttribute("forceID" ).read(forceID );
     
+    bool H5back = true;
+    if (file.hasAttribute("Version")) H5back = false;
+
     // Open the snapshot group
     //
     auto snaps = file.getGroup("snapshots");
@@ -115,6 +120,18 @@ namespace CoefClasses
       if (Time < Tmin or Time > Tmax) continue;
 
       auto in = stanza.getDataSet("coefficients").read<Eigen::MatrixXcd>();
+
+      if (H5back and H5BackCompat) {
+
+	auto in2 = stanza.getDataSet("coefficients").read<Eigen::MatrixXcd>();
+	in2.transposeInPlace();
+    
+	for (size_t c=0, n=0; c<in.cols(); c++) {
+	  for (size_t r=0, n=0; r<in.rows(); r++) {
+	    in(r, c) = in2.data()[n++];
+	  }
+	}
+      }
       
       // Pack the data into the coefficient variable
       //
@@ -777,6 +794,9 @@ namespace CoefClasses
     file.getAttribute("config" ).read(config);
     file.getDataSet  ("count"  ).read(count );
     
+    bool H5back = true;
+    if (file.hasAttribute("Version")) H5back = false;
+
     // Open the snapshot group
     //
     auto snaps = file.getGroup("snapshots");
@@ -801,6 +821,18 @@ namespace CoefClasses
       if (Time < Tmin or Time > Tmax) continue;
 
       auto in = stanza.getDataSet("coefficients").read<Eigen::MatrixXcd>();
+
+      if (H5back and H5BackCompat) {
+
+	auto in2 = stanza.getDataSet("coefficients").read<Eigen::MatrixXcd>();
+	in2.transposeInPlace();
+    
+	for (size_t c=0, n=0; c<in.cols(); c++) {
+	  for (size_t r=0; r<in.rows(); r++) {
+	    in(r, c) = in2.data()[n++];
+	  }
+	}
+      }
       
       // Work around for previous unitiaized data bug; enforces real data
       //
@@ -2253,6 +2285,10 @@ namespace CoefClasses
 			  HighFive::File::ReadWrite |
 			  HighFive::File::Create);
       
+      // Write the Version string
+      //
+      file.createAttribute<std::string>("Version", HighFive::DataSpace::From(Version)).write(Version);
+
       // We write the coefficient file geometry
       //
       file.createAttribute<std::string>("geometry", HighFive::DataSpace::From(geometry)).write(geometry);
diff --git a/pyEXP/CoefWrappers.cc b/pyEXP/CoefWrappers.cc
index b1ef1d91e..f5736a3f7 100644
--- a/pyEXP/CoefWrappers.cc
+++ b/pyEXP/CoefWrappers.cc
@@ -867,6 +867,19 @@ void CoefficientClasses(py::module &m) {
          are found at the requested time
          )",
          py::arg("time"))
+    .def("newH5",
+	 &CoefClasses::Coefs::setNewH5,
+	 R"(
+         Override backwards compatibility flag for HighFive API change
+
+         Parameters
+         ----------
+         None
+
+         Returns
+         -------
+         None
+         )")
     .def("setData",
          &CoefClasses::Coefs::setData,
          R"(

From 5a6e3ebc59cf80c4adb3dd54b1169dbbc57b5c89 Mon Sep 17 00:00:00 2001
From: Martin Weinberg <martin.weinberg@comcast.net>
Date: Wed, 22 May 2024 17:42:35 -0400
Subject: [PATCH 21/25] Update expui/Coefficients.cc

Co-authored-by: Michael Petersen <petersen.michael.s@gmail.com>
---
 expui/Coefficients.cc | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/expui/Coefficients.cc b/expui/Coefficients.cc
index 227125751..1c73c5207 100644
--- a/expui/Coefficients.cc
+++ b/expui/Coefficients.cc
@@ -821,7 +821,7 @@ namespace CoefClasses
       if (Time < Tmin or Time > Tmax) continue;
 
       auto in = stanza.getDataSet("coefficients").read<Eigen::MatrixXcd>();
-
+     // If an older version of the coefficients and backwards compatibility is desired, re-order the coefficients to match the cache.
       if (H5back and H5BackCompat) {
 
 	auto in2 = stanza.getDataSet("coefficients").read<Eigen::MatrixXcd>();

From 17ea22ce7495f0ee68032e92eb7ce12929b181ce Mon Sep 17 00:00:00 2001
From: Martin Weinberg <martin.weinberg@comcast.net>
Date: Wed, 22 May 2024 17:42:52 -0400
Subject: [PATCH 22/25] Update expui/Coefficients.cc

Co-authored-by: Michael Petersen <petersen.michael.s@gmail.com>
---
 expui/Coefficients.cc | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/expui/Coefficients.cc b/expui/Coefficients.cc
index 1c73c5207..158e9f8bb 100644
--- a/expui/Coefficients.cc
+++ b/expui/Coefficients.cc
@@ -795,7 +795,7 @@ namespace CoefClasses
     file.getDataSet  ("count"  ).read(count );
     
     bool H5back = true;
-    if (file.hasAttribute("Version")) H5back = false;
+    if (file.hasAttribute("CoefficientOutputVersion")) H5back = false;
 
     // Open the snapshot group
     //

From c48a11a76e3846a28191adb1fa019cd460684240 Mon Sep 17 00:00:00 2001
From: Martin Weinberg <martin.weinberg@comcast.net>
Date: Wed, 22 May 2024 17:43:09 -0400
Subject: [PATCH 23/25] Update expui/Coefficients.cc

Co-authored-by: Michael Petersen <petersen.michael.s@gmail.com>
---
 expui/Coefficients.cc | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/expui/Coefficients.cc b/expui/Coefficients.cc
index 158e9f8bb..59a2321a0 100644
--- a/expui/Coefficients.cc
+++ b/expui/Coefficients.cc
@@ -2287,7 +2287,7 @@ namespace CoefClasses
       
       // Write the Version string
       //
-      file.createAttribute<std::string>("Version", HighFive::DataSpace::From(Version)).write(Version);
+      file.createAttribute<std::string>("CoefficientOutputVersion", HighFive::DataSpace::From(CoefficientOutputVersion)).write(CoefficientOutputVersion);
 
       // We write the coefficient file geometry
       //

From a46897a8fe34efa253c834ac5075965510f38ddc Mon Sep 17 00:00:00 2001
From: Martin Weinberg <martin.weinberg@comcast.net>
Date: Wed, 22 May 2024 17:43:16 -0400
Subject: [PATCH 24/25] Update expui/Coefficients.H

Co-authored-by: Michael Petersen <petersen.michael.s@gmail.com>
---
 expui/Coefficients.H | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/expui/Coefficients.H b/expui/Coefficients.H
index 737881916..8615c22aa 100644
--- a/expui/Coefficients.H
+++ b/expui/Coefficients.H
@@ -68,7 +68,7 @@ namespace CoefClasses
     virtual std::string getYAML() = 0;
     
     //! Coefficient file versioning
-    inline static const std::string Version = "1.0";
+    inline static const std::string CoefficientOutputVersion = "1.0";
 
     //! Write parameter attributes (needed for derived classes)
     virtual void WriteH5Params(HighFive::File& file) = 0;

From aef851a552718ef5233076294b7e8853cca0b453 Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Wed, 22 May 2024 17:56:52 -0400
Subject: [PATCH 25/25] Some fixes for Mike's suggested code changes

---
 expui/Coefficients.H  |  6 ------
 expui/Coefficients.cc | 22 ++++++++++++++++------
 pyEXP/CoefWrappers.cc | 13 -------------
 3 files changed, 16 insertions(+), 25 deletions(-)

diff --git a/expui/Coefficients.H b/expui/Coefficients.H
index 8615c22aa..d6a121815 100644
--- a/expui/Coefficients.H
+++ b/expui/Coefficients.H
@@ -54,9 +54,6 @@ namespace CoefClasses
     //! Verbose debugging output
     bool verbose;
     
-    //! Backward compatibility flag for HighFive
-    static bool H5BackCompat;
-
     //! Time vector
     std::vector<double> times;
 
@@ -201,9 +198,6 @@ namespace CoefClasses
     //! Set maximum grid interpolation offset
     void setDeltaT(double dT) { deltaT = dT; }
 
-    //! Override backward compatibility for HighFive
-    static void setNewH5() { H5BackCompat = false; }
-
     class CoefsError : public std::runtime_error
     {
     public:
diff --git a/expui/Coefficients.cc b/expui/Coefficients.cc
index 59a2321a0..230291ec7 100644
--- a/expui/Coefficients.cc
+++ b/expui/Coefficients.cc
@@ -20,8 +20,6 @@
 namespace CoefClasses
 {
   
-  bool Coefs::H5BackCompat = true;
-
   void Coefs::copyfields(std::shared_ptr<Coefs> p)
   {
     // These variables will copy data, not pointers
@@ -93,8 +91,11 @@ namespace CoefClasses
     file.getAttribute("geometry").read(geometry);
     file.getAttribute("forceID" ).read(forceID );
     
+    // Look for Coef output version to toggle backward compatibility
+    // with legacy storage order
+    //
     bool H5back = true;
-    if (file.hasAttribute("Version")) H5back = false;
+    if (file.hasAttribute("CoefficientOutputVersion")) H5back = false;
 
     // Open the snapshot group
     //
@@ -121,7 +122,10 @@ namespace CoefClasses
 
       auto in = stanza.getDataSet("coefficients").read<Eigen::MatrixXcd>();
 
-      if (H5back and H5BackCompat) {
+      // If we have a legacy set of coefficients, re-order the
+      // coefficients to match the new HighFive/Eigen ordering
+      //
+      if (H5back) {
 
 	auto in2 = stanza.getDataSet("coefficients").read<Eigen::MatrixXcd>();
 	in2.transposeInPlace();
@@ -794,6 +798,9 @@ namespace CoefClasses
     file.getAttribute("config" ).read(config);
     file.getDataSet  ("count"  ).read(count );
     
+    // Look for Coef output version to toggle backward compatibility
+    // with legacy storage order
+    //
     bool H5back = true;
     if (file.hasAttribute("CoefficientOutputVersion")) H5back = false;
 
@@ -821,8 +828,11 @@ namespace CoefClasses
       if (Time < Tmin or Time > Tmax) continue;
 
       auto in = stanza.getDataSet("coefficients").read<Eigen::MatrixXcd>();
-     // If an older version of the coefficients and backwards compatibility is desired, re-order the coefficients to match the cache.
-      if (H5back and H5BackCompat) {
+
+      // If we have a legacy set of coefficients, re-order the
+      // coefficients to match the new HighFive/Eigen ordering
+      //
+      if (H5back) {
 
 	auto in2 = stanza.getDataSet("coefficients").read<Eigen::MatrixXcd>();
 	in2.transposeInPlace();
diff --git a/pyEXP/CoefWrappers.cc b/pyEXP/CoefWrappers.cc
index f5736a3f7..b1ef1d91e 100644
--- a/pyEXP/CoefWrappers.cc
+++ b/pyEXP/CoefWrappers.cc
@@ -867,19 +867,6 @@ void CoefficientClasses(py::module &m) {
          are found at the requested time
          )",
          py::arg("time"))
-    .def("newH5",
-	 &CoefClasses::Coefs::setNewH5,
-	 R"(
-         Override backwards compatibility flag for HighFive API change
-
-         Parameters
-         ----------
-         None
-
-         Returns
-         -------
-         None
-         )")
     .def("setData",
          &CoefClasses::Coefs::setData,
          R"(