From 76cdbd551708481e4a4805e547b82637bbc1d5f4 Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Wed, 7 Aug 2024 17:10:57 -0400
Subject: [PATCH 01/52] Implementation of a point set evaluator for
 FieldGenerator

---
 expui/FieldGenerator.H  |  23 +++++-
 expui/FieldGenerator.cc | 178 ++++++++++++++++++++++++++++++++++++++++
 pyEXP/FieldWrappers.cc  |  86 ++++++++++++++++---
 3 files changed, 273 insertions(+), 14 deletions(-)

diff --git a/expui/FieldGenerator.H b/expui/FieldGenerator.H
index 3b4873ade..4ab573aef 100644
--- a/expui/FieldGenerator.H
+++ b/expui/FieldGenerator.H
@@ -20,6 +20,7 @@ namespace Field
     
     std::vector<double> times, pmin, pmax;
     std::vector<int>    grid;
+    Eigen::MatrixXd     mesh;
     
     //! Sanity check time vector with coefficient DB
     void check_times(CoefClasses::CoefsPtr coefs);
@@ -35,12 +36,32 @@ namespace Field
 
   public:
     
-    //! Constructor
+    //! Constructor for a rectangular grid
     FieldGenerator(const std::vector<double> &time,
 		   const std::vector<double> &pmin,
 		   const std::vector<double> &pmax,
 		   const std::vector<int>    &grid);
 
+    //! Constructor for an arbitrary point mesh.  The mesh should be
+    //! an Nx3 array
+    FieldGenerator(const std::vector<double> &time,
+		   const Eigen::MatrixXd     &mesh);
+
+    /** Get field quantities at user defined points
+
+	For example:
+	.
+	.
+	// Generate the fields for all coefficients in 'coefs'
+	auto db = points(basis, coefs);
+
+	// Get fields evaluated at Time=3.14 and for all points in
+	// your supplied mesh for density ("dens")
+	Eigen::MatrixXf points = db["3.14"]["dens"];
+    */
+    std::map<double, std::map<std::string, Eigen::VectorXf>>
+    points(BasisClasses::BasisPtr basis, CoefClasses::CoefsPtr coefs);
+
     /** Get a field slices as a map in time and type
 
 	For example:
diff --git a/expui/FieldGenerator.cc b/expui/FieldGenerator.cc
index 5e00e1c82..73a1595f7 100644
--- a/expui/FieldGenerator.cc
+++ b/expui/FieldGenerator.cc
@@ -41,6 +41,33 @@ namespace Field
     }
   }
   
+  FieldGenerator::FieldGenerator(const std::vector<double> &time,
+				 const Eigen::MatrixXd     &mesh) :
+    
+    times(time), mesh(mesh)
+  {
+    // Sanity check on mesh
+    //
+    if (mesh.cols() != 3)
+      throw std::runtime_error("FieldGenerator: bad mesh specification. The mesh must be an Nx3 array where the columns are Cartesian points");
+
+    // Check whether MPI is initialized
+    //
+    int flag;
+    MPI_Initialized(&flag);
+    if (flag) use_mpi = true;
+    else      use_mpi = false;
+
+    
+    // Fall back sanity (works for me but this needs to be fixed
+    // generally)
+    //
+    if (use_mpi) {
+      MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
+      MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+    }
+  }
+  
   void FieldGenerator::check_times(CoefClasses::CoefsPtr coefs)
   {
     std::vector<double> ctimes = coefs->Times();
@@ -867,5 +894,156 @@ namespace Field
   }
   // END histogram1d
 
+  std::map<double, std::map<std::string, Eigen::VectorXf>>
+  FieldGenerator::points(BasisClasses::BasisPtr basis,
+			 CoefClasses::CoefsPtr coefs)
+  {
+    // Sanity check on mesh
+    //
+    if (mesh.size() == 0)
+      throw std::runtime_error("FieldGenerator::points: bad mesh specification.  Did you call the mesh constructor?");
+
+    // Set midplane evaluation parameters
+    //
+    basis->setMidplane(midplane);
+    basis->setColumnHeight(colheight);
+
+    // Check
+    //
+    check_times(coefs);
+
+    std::map<double, std::map<std::string, Eigen::VectorXf>> ret;
+
+    // Now get the desired coordinate type
+    //
+    auto ctype = basis->coordinates;
+
+    // Field labels (force field labels added below)
+    //
+    auto labels = basis->getFieldLabels(ctype);
+
+    int ncnt = 0;		// Process counter for MPI
+
+    std::map<std::string, Eigen::VectorXf> frame;
+    for (auto label : labels) {
+      frame[label].resize(mesh.rows());
+    }	
+
+    for (auto T : times) {
+
+      if (ncnt++ % numprocs > 0) continue;
+
+      if (not coefs->getCoefStruct(T)) {
+	std::cout << "Could not find time=" << T << ", continuing" << std::endl;
+	continue;
+      }
+
+      basis->set_coefs(coefs->getCoefStruct(T));
+
+#pragma omp parallel for
+      for (int k=0; k<mesh.rows(); k++) {
+
+	// Cartesian to spherical for all_eval
+	//
+	double x = mesh(k, 0);
+	double y = mesh(k, 1);
+	double z = mesh(k, 2);
+
+	// Coordinate values
+	double r, costh, phi, R;
+
+	// Return values
+	double p0, p1, d0, d1, f1, f2, f3;
+	std::vector<double> v;
+
+	if (ctype == BasisClasses::Basis::Coord::Spherical) {
+	  r     = sqrt(x*x + y*y + z*z) + 1.0e-18;
+	  costh = z/r;
+	  phi   = atan2(y, x);
+	  v = (*basis)(r, costh, phi, ctype);
+	} else if (ctype == BasisClasses::Basis::Coord::Cylindrical) {
+	  R     = sqrt(x*x + y*y) + 1.0e-18;
+	  phi   = atan2(y, x);
+	  v = (*basis)(R, z, phi, ctype);
+	} else {
+	  v = (*basis)(x, y, z, BasisClasses::Basis::Coord::Cartesian);
+	}
+	
+	// Pack the frame structure
+	//
+	for (int n=0; n<labels.size(); n++)
+	  frame[labels[n]](k) = v[n];
+      }
+
+      ret[T] = frame;
+    }
+
+    if (use_mpi) {
+      
+      std::vector<char> bf(9);
+
+      for (int n=1; n<numprocs; n++) {
+
+	if (myid==n) {
+	  int sz = ret.size();
+	  MPI_Send(&sz, 1, MPI_INT, 0, 102, MPI_COMM_WORLD);
+	  for (auto & v : ret) {
+	    MPI_Send(&v.first, 1, MPI_DOUBLE, 0, 103, MPI_COMM_WORLD);
+	    int fsz = v.second.size();
+	    MPI_Send(&fsz, 1, MPI_INT, 0, 104, MPI_COMM_WORLD);
+	    
+	    for (auto & f : v.second) {
+	      MPI_Send(f.first.c_str(), f.first.length(), MPI_CHAR, 0, 105,
+		       MPI_COMM_WORLD);
+	      
+	      MPI_Send(f.second.data(), f.second.size(),
+		       MPI_FLOAT, 0, 106, MPI_COMM_WORLD);
+	    }
+	  }
+	}
+
+	if (myid==0) {
+	  MPI_Status status;
+	  int sz, fsz, l;
+	  double T;
+
+	  MPI_Recv(&sz, 1, MPI_INT, n, 102, MPI_COMM_WORLD, &status);
+	  for (int i=0; i<sz; i++) {
+	    MPI_Recv(&T, 1, MPI_DOUBLE, n, 103, MPI_COMM_WORLD, &status);
+	    MPI_Recv(&fsz, 1, MPI_INT, n, 104, MPI_COMM_WORLD, &status);
+	    
+	    for (int j=0; j<fsz; j++) {
+	      // Get the field name
+	      //
+	      MPI_Probe(n, 105, MPI_COMM_WORLD, &status);
+	      MPI_Get_count(&status, MPI_CHAR, &l);
+	      MPI_Recv(bf.data(), l, MPI_CHAR, n, 105, MPI_COMM_WORLD, &status);
+	      std::string s(bf.data(), l);
+	      
+	      // Sanity check
+	      //
+	      if (frame.find(s) == frame.end()) {
+		std::cerr << "Error finding <" << s << "> in field map"
+			  << std::endl;
+	      }
+
+	      // Get the data
+	      //
+	      MPI_Recv(frame[s].data(), frame[s].size(), MPI_FLOAT, n, 106,
+		       MPI_COMM_WORLD, &status);
+	    }
+	    
+	    ret[T] = frame;
+	  }
+	}
+      }
+    }
+
+    // Toggle off midplane evaluation
+    basis->setMidplane(false);
+
+    return ret;
+  }
+  
 }
 // END namespace Field
diff --git a/pyEXP/FieldWrappers.cc b/pyEXP/FieldWrappers.cc
index fb79bef5c..995f91358 100644
--- a/pyEXP/FieldWrappers.cc
+++ b/pyEXP/FieldWrappers.cc
@@ -24,22 +24,24 @@ void FieldGeneratorClasses(py::module &m) {
     a list of upper bounds, and a list of knots per dimension.  These
     lists all have rank 3 for (x, y, z).  For a two-dimensional surface,
     one of the knot array values must be zero.  The member functions
-    lines, slices and volumes, called with the basis and coefficient
-    objects, return a numpy.ndarray containing the field evaluations.
-    Each of these functions returns a dictionary of times to a dictionary
-    of field names to numpy.ndarrays at each time.  There are also members
-    which will write these generated fields to files. The linear probe
-    members, 'lines' and 'file_lines', evaluate 'num' field points along
-    a user-specified segment between the 3d points 'beg' and 'end'.  See
-    help(pyEXP.basis) and help(pyEXP.coefs) for info on the basis and
-    coefficient objects.
+    lines, slices, an arbitrary point-set, and volumes, called with the
+    basis and coefficient objects, return a numpy.ndarray containing the
+    field evaluations.  Each of these functions returns a dictionary of
+    times to a dictionary of field names to numpy.ndarrays at each time.
+    There are also members which will write these generated fields to files.
+    The linear probe members, 'lines' and 'file_lines', evaluate 'num'
+    field points along a user-specified segment between the 3d points 'beg'
+    and 'end'.  See help(pyEXP.basis) and help(pyEXP.coefs) for info on
+    the basis and coefficient objects.
 
     Data packing
     ------------
     All slices and volumes are returned as numpy.ndarrays in row-major
     order.  That is, the first index is x, the second is y, and the
     third is z.  The ranks are specified by the 'gridsize' array with
-    (nx, ny, nz) as input to the FieldGenerator constructor.
+    (nx, ny, nz) as input to the FieldGenerator constructor.  A point
+    mesh is rows of (x, y, z) Cartesian coordinates.  The output field
+    values returned by points is in the same order as the input array.
 
     Coordinate systems
     ------------------
@@ -116,6 +118,24 @@ void FieldGeneratorClasses(py::module &m) {
 	py::arg("times"), py::arg("lower"), py::arg("upper"),
 	py::arg("gridsize"));
 
+  f.def(py::init<const std::vector<double>, const Eigen::MatrixXd>(),
+	R"(
+        Create fields for given times and and provided point set
+
+        Parameters
+        ----------
+        times : list(float,...)
+            list of evaluation times
+        mesh : ndarray of Nx3 floats
+            point set (x, y, z) for field evaluation
+
+        Returns
+        -------
+        FieldGenerator
+            new object
+        )",
+	py::arg("times"), py::arg("mesh"));
+
   f.def("setMidplane", &Field::FieldGenerator::setMidplane,
 	R"(
         Set the field generator to generate midplane fields
@@ -160,6 +180,38 @@ void FieldGeneratorClasses(py::module &m) {
 
         See also
         --------
+        points : generate fields at an array of mesh points
+        lines : generate fields along a line given by its end points
+        volumes : generate fields in volume given by the initializtion grid
+       )", py::arg("basis"), py::arg("coefs"));
+  
+  f.def("points", &Field::FieldGenerator::points,
+	R"(
+        Return a dictionary of arrays (1d numpy arrays) indexed by time
+        and field type corresponding to the mesh points
+
+        Parameters
+        ----------
+        basis : Basis
+            basis instance of any geometry; geometry will be deduced by the generator
+        coefs : Coefs
+            coefficient container instance
+
+        Returns
+        -------
+        dict({time: {field-name: numpy.ndarray})
+            dictionary of times, field names, and data arrays
+
+        Notes
+        -----
+        Each data array in the dictionary is a 1-dimensional array with the
+        same order as the mesh points in the constructor
+        routines.
+
+        See also
+        --------
+        slices : generate fields in a surface slice given by the
+                 initializtion grid
         lines : generate fields along a line given by its end points
         volumes : generate fields in volume given by the initializtion grid
        )", py::arg("basis"), py::arg("coefs"));
@@ -196,7 +248,9 @@ void FieldGeneratorClasses(py::module &m) {
 
         See also
         --------
-        slices : generate fields in a surface slice given by the initializtion grid
+        points : generate fields at an array of mesh points
+        slices : generate fields in a surface slice given by the
+                 initializtion grid
         volumes : generate fields in volume given by the initializtion grid
         )",
 	py::arg("basis"), py::arg("coefs"),
@@ -281,7 +335,9 @@ void FieldGeneratorClasses(py::module &m) {
 
         See also
         --------
-        slices : generate fields in a surface slice given by the initializtion grid
+        points : generate fields at an array of mesh points
+        slices : generate fields in a surface slice given by the
+                 initializtion grid
         volumes : generate fields in volume given by the initializtion grid
         lines : generate fields along a line given by its end points
         )", 
@@ -314,7 +370,9 @@ void FieldGeneratorClasses(py::module &m) {
 
         See also
         --------
-        slices : generate fields in a surface slice given by the initializtion grid
+        points : generate fields at an array of mesh points
+        slices : generate fields in a surface slice given by the
+                 initializtion grid
         volumes : generate fields in volume given by the initializtion grid
         lines : generate fields along a line given by its end points
         )",
@@ -357,6 +415,8 @@ void FieldGeneratorClasses(py::module &m) {
 
         See also
         --------
+        points : generate fields at an array of mesh points
+
         lines : generate fields along a line given by its end points
         slices : generate fields in a slice given by the initializtion grid
         )");

From 9739781997d9736a0953004906bc146c593884c8 Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Thu, 26 Sep 2024 18:24:13 -0400
Subject: [PATCH 02/52] Use Eigen::Ref<> to pass arrays to prevent pybind11
 mapping failure in pyEXP

---
 expui/Coefficients.H  | 22 ++++++++++------------
 expui/Coefficients.cc | 18 +++++++++---------
 pyEXP/CoefWrappers.cc | 12 ++++++------
 3 files changed, 25 insertions(+), 27 deletions(-)

diff --git a/expui/Coefficients.H b/expui/Coefficients.H
index d6a121815..34a9c3a89 100644
--- a/expui/Coefficients.H
+++ b/expui/Coefficients.H
@@ -111,7 +111,7 @@ namespace CoefClasses
     virtual Eigen::VectorXcd& getData(double time) = 0;
 
     //! Set coefficient store at given time to the provided matrix
-    virtual void setData(double time, const Eigen::VectorXcd& data) = 0;
+    virtual void setData(double time, Eigen::Ref<Eigen::VectorXcd> data) = 0;
 
     //! Interpolate coefficient matrix at given time
     std::tuple<Eigen::VectorXcd&, bool> interpolate(double time);
@@ -273,10 +273,10 @@ namespace CoefClasses
     Eigen::MatrixXcd& getMatrix(double time);
 
     //! Set coefficient matrix at given time
-    virtual void setData(double time, const Eigen::VectorXcd& arr);
+    virtual void setData(double time, Eigen::Ref<Eigen::VectorXcd> arr);
 
     //! Natural data setter for pyEXP
-    void setMatrix(double time, const Eigen::MatrixXcd& mat);
+    void setMatrix(double time, Eigen::Ref<Eigen::MatrixXcd> mat);
 
     //! Interpolate coefficient matrix at given time
     std::tuple<Eigen::MatrixXcd&, bool> interpolate(double time)
@@ -411,10 +411,10 @@ namespace CoefClasses
     Eigen::MatrixXcd& getMatrix(double time);
     
     //! Set coefficient matrix at given time
-    virtual void setData(double time, const Eigen::VectorXcd& arr);
+    virtual void setData(double time, Eigen::Ref<Eigen::VectorXcd> arr);
     
     //! For pyEXP
-    void setMatrix(double time, const Eigen::MatrixXcd& arr);
+    void setMatrix(double time, Eigen::Ref<Eigen::MatrixXcd> arr);
     
     //! Interpolate coefficient matrix at given time
     std::tuple<Eigen::MatrixXcd&, bool> interpolate(double time)
@@ -542,7 +542,7 @@ namespace CoefClasses
     virtual Eigen3d& getTensor(double time);
     
     //! Set coefficient matrix at given time
-    virtual void setData(double time, const Eigen::VectorXcd& dat);
+    virtual void setData(double time, Eigen::Ref<Eigen::VectorXcd> dat);
 
     //! Native version
     virtual void setTensor(double time, const Eigen3d& dat);
@@ -673,7 +673,7 @@ namespace CoefClasses
     virtual Eigen3d& getTensor(double time);
     
     //! Set coefficient matrix at given time
-    virtual void setData(double time, const Eigen::VectorXcd& dat);
+    virtual void setData(double time, Eigen::Ref<Eigen::VectorXcd> dat);
 
     //! Native version
     virtual void setTensor(double time, const Eigen3d& dat);
@@ -810,9 +810,7 @@ namespace CoefClasses
     /** Set coefficient matrix at given time.  This is for pybind11,
 	since the operator() will not allow lvalue assignment, it
 	seems. */
-    virtual void setData(double time, const Eigen::VectorXcd& arr);
-
-    void setArray(double time, const Eigen::VectorXcd& arr) { setData(time, arr); }
+    virtual void setData(double time, Eigen::Ref<Eigen::VectorXcd> arr);
 
     //! Get coefficient structure at a given time
     virtual std::shared_ptr<CoefStruct> getCoefStruct(double time)
@@ -914,7 +912,7 @@ namespace CoefClasses
     SphFldStruct::coefType& getMatrix(double time);
 
     //! Set coefficient matrix at given time
-    virtual void setData(double time, const Eigen::VectorXcd& arr);
+    virtual void setData(double time, Eigen::Ref<Eigen::VectorXcd> arr);
 
     //! Natural data setter for pyEXP
     void setMatrix(double time, const SphFldStruct::coefType& mat);
@@ -1041,7 +1039,7 @@ namespace CoefClasses
     CylFldStruct::coefType & getMatrix(double time);
 
     //! Set coefficient matrix at given time
-    virtual void setData(double time, const Eigen::VectorXcd& arr);
+    virtual void setData(double time, Eigen::Ref<Eigen::VectorXcd> arr);
 
     //! Natural data setter for pyEXP
     void setMatrix(double time, const CylFldStruct::coefType& mat);
diff --git a/expui/Coefficients.cc b/expui/Coefficients.cc
index bccd2675f..dfd58e5b8 100644
--- a/expui/Coefficients.cc
+++ b/expui/Coefficients.cc
@@ -469,7 +469,7 @@ namespace CoefClasses
     return mat;
   }
   
-  void SphCoefs::setData(double time, const Eigen::VectorXcd& dat)
+  void SphCoefs::setData(double time, Eigen::Ref<Eigen::VectorXcd> dat)
   {
     auto it = coefs.find(roundTime(time));
 
@@ -484,7 +484,7 @@ namespace CoefClasses
     }
   }
   
-  void SphCoefs::setMatrix(double time, const Eigen::MatrixXcd& dat)
+  void SphCoefs::setMatrix(double time, Eigen::Ref<Eigen::MatrixXcd> dat)
   {
     auto it = coefs.find(roundTime(time));
 
@@ -891,7 +891,7 @@ namespace CoefClasses
     return mat;
   }
   
-  void CylCoefs::setData(double time, const Eigen::VectorXcd& dat)
+  void CylCoefs::setData(double time, Eigen::Ref<Eigen::VectorXcd> dat)
   {
     auto it = coefs.find(roundTime(time));
 
@@ -906,7 +906,7 @@ namespace CoefClasses
     }
   }
 
-  void CylCoefs::setMatrix(double time, const Eigen::MatrixXcd& dat)
+  void CylCoefs::setMatrix(double time, Eigen::Ref<Eigen::MatrixXcd> dat)
   {
     auto it = coefs.find(roundTime(time));
 
@@ -1247,7 +1247,7 @@ namespace CoefClasses
     return arr;
   }
 
-  void SlabCoefs::setData(double time, const Eigen::VectorXcd& dat)
+  void SlabCoefs::setData(double time, Eigen::Ref<Eigen::VectorXcd> dat)
   {
     auto it = coefs.find(roundTime(time));
 
@@ -1599,7 +1599,7 @@ namespace CoefClasses
     return arr;
   }
 
-  void CubeCoefs::setData(double time, const Eigen::VectorXcd& dat)
+  void CubeCoefs::setData(double time, Eigen::Ref<Eigen::VectorXcd> dat)
   {
     auto it = coefs.find(roundTime(time));
 
@@ -1968,7 +1968,7 @@ namespace CoefClasses
     return arr;
   }
   
-  void TableData::setData(double time, const Eigen::VectorXcd& dat)
+  void TableData::setData(double time, Eigen::Ref<Eigen::VectorXcd> dat)
   {
     auto it = coefs.find(roundTime(time));
 
@@ -2458,7 +2458,7 @@ namespace CoefClasses
     return *mat;
   }
   
-  void SphFldCoefs::setData(double time, const Eigen::VectorXcd& dat)
+  void SphFldCoefs::setData(double time, Eigen::Ref<Eigen::VectorXcd> dat)
   {
     auto it = coefs.find(roundTime(time));
 
@@ -2897,7 +2897,7 @@ namespace CoefClasses
     return *mat;
   }
   
-  void CylFldCoefs::setData(double time, const Eigen::VectorXcd& dat)
+  void CylFldCoefs::setData(double time, Eigen::Ref<Eigen::VectorXcd> dat)
   {
     auto it = coefs.find(roundTime(time));
 
diff --git a/pyEXP/CoefWrappers.cc b/pyEXP/CoefWrappers.cc
index cce752970..d170cfd33 100644
--- a/pyEXP/CoefWrappers.cc
+++ b/pyEXP/CoefWrappers.cc
@@ -155,7 +155,7 @@ void CoefficientClasses(py::module &m) {
       PYBIND11_OVERRIDE_PURE(Eigen::VectorXcd&, Coefs, getData, time);
     }
 
-    void setData(double time, const Eigen::VectorXcd& array) override {
+    void setData(double time, Eigen::Ref<Eigen::VectorXcd> array) override {
       PYBIND11_OVERRIDE_PURE(void, Coefs, setData, time, array);
     }
 
@@ -231,7 +231,7 @@ void CoefficientClasses(py::module &m) {
       PYBIND11_OVERRIDE(Eigen::VectorXcd&, SphCoefs, getData, time);
     }
 
-    void setData(double time, const Eigen::VectorXcd& array) override {
+    void setData(double time, Eigen::Ref<Eigen::VectorXcd> array) override {
       PYBIND11_OVERRIDE(void, SphCoefs, setData, time, array);
     }
 
@@ -310,7 +310,7 @@ void CoefficientClasses(py::module &m) {
       PYBIND11_OVERRIDE(Eigen::VectorXcd&, CylCoefs, getData, time);
     }
 
-    void setData(double time, const Eigen::VectorXcd& array) override {
+    void setData(double time, Eigen::Ref<Eigen::VectorXcd> array) override {
       PYBIND11_OVERRIDE(void, CylCoefs, setData, time, array);
     }
 
@@ -388,7 +388,7 @@ void CoefficientClasses(py::module &m) {
       PYBIND11_OVERRIDE(Eigen::VectorXcd&, SlabCoefs, getData, time);
     }
 
-    void setData(double time, const Eigen::VectorXcd& array) override {
+    void setData(double time, Eigen::Ref<Eigen::VectorXcd> array) override {
       PYBIND11_OVERRIDE(void, SlabCoefs, setData, time, array);
     }
 
@@ -467,7 +467,7 @@ void CoefficientClasses(py::module &m) {
       PYBIND11_OVERRIDE(Eigen::VectorXcd&, CubeCoefs, getData, time);
     }
 
-    void setData(double time, const Eigen::VectorXcd& array) override {
+    void setData(double time, Eigen::Ref<Eigen::VectorXcd> array) override {
       PYBIND11_OVERRIDE(void, CubeCoefs, setData, time, array);
     }
 
@@ -546,7 +546,7 @@ void CoefficientClasses(py::module &m) {
       PYBIND11_OVERRIDE(Eigen::VectorXcd&, TableData, getData, time);
     }
 
-    void setData(double time, const Eigen::VectorXcd& array) override {
+    void setData(double time, Eigen::Ref<Eigen::VectorXcd> array) override {
       PYBIND11_OVERRIDE(void, TableData, setData, time, array);
     }
 

From 12176dc9bf17c23297b9fc4a1a1aa41c414ac3bf Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Sat, 28 Sep 2024 11:45:58 -0400
Subject: [PATCH 03/52] Change names of Eigen::Tensor-->ndarray generators for
 clarity

---
 pyEXP/BasisWrappers.cc |  2 +-
 pyEXP/FieldWrappers.cc | 11 ++++---
 pyEXP/TensorToArray.H  | 68 ++++++++++++++++++++++--------------------
 3 files changed, 43 insertions(+), 38 deletions(-)

diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc
index df9e8346f..e0b72ab25 100644
--- a/pyEXP/BasisWrappers.cc
+++ b/pyEXP/BasisWrappers.cc
@@ -2014,7 +2014,7 @@ void BasisFactoryClasses(py::module &m)
 	  std::tie(T, O) =
 	    BasisClasses::IntegrateOrbits(tinit, tfinal, h, ps, bfe, F, stride);
 
-	  py::array_t<float> ret = make_ndarray<float>(O);
+	  py::array_t<float> ret = make_ndarray3<float>(O);
 	  return std::tuple<Eigen::VectorXd, py::array_t<float>>(T, ret);
 	},
 	R"(
diff --git a/pyEXP/FieldWrappers.cc b/pyEXP/FieldWrappers.cc
index fb79bef5c..5e0c1dd4e 100644
--- a/pyEXP/FieldWrappers.cc
+++ b/pyEXP/FieldWrappers.cc
@@ -196,7 +196,8 @@ void FieldGeneratorClasses(py::module &m) {
 
         See also
         --------
-        slices : generate fields in a surface slice given by the initializtion grid
+        slices : generate fields in a surface slice given by the
+                 initializtion grid
         volumes : generate fields in volume given by the initializtion grid
         )",
 	py::arg("basis"), py::arg("coefs"),
@@ -281,7 +282,8 @@ void FieldGeneratorClasses(py::module &m) {
 
         See also
         --------
-        slices : generate fields in a surface slice given by the initializtion grid
+        slices : generate fields in a surface slice given by the
+                 initializtion grid
         volumes : generate fields in volume given by the initializtion grid
         lines : generate fields along a line given by its end points
         )", 
@@ -314,7 +316,8 @@ void FieldGeneratorClasses(py::module &m) {
 
         See also
         --------
-        slices : generate fields in a surface slice given by the initializtion grid
+        slices : generate fields in a surface slice given by the
+                 initializtion grid
         volumes : generate fields in volume given by the initializtion grid
         lines : generate fields along a line given by its end points
         )",
@@ -329,7 +332,7 @@ void FieldGeneratorClasses(py::module &m) {
     auto vols = A.volumes(basis, coefs);
     for (auto & v : vols) {
       for (auto & u : v.second) {
-	ret[v.first][u.first] = make_ndarray<float>(u.second);
+	ret[v.first][u.first] = make_ndarray3<float>(u.second);
       }
     }
 
diff --git a/pyEXP/TensorToArray.H b/pyEXP/TensorToArray.H
index f55b1bc34..765148eac 100644
--- a/pyEXP/TensorToArray.H
+++ b/pyEXP/TensorToArray.H
@@ -3,7 +3,7 @@
 
 //! Helper function that maps the Eigen::Tensor<T, 3> into an numpy.ndarray
 template <typename T>
-py::array_t<T> make_ndarray(Eigen::Tensor<T, 3>& mat)
+py::array_t<T> make_ndarray3(Eigen::Tensor<T, 3>& mat)
 {
   // Get the tensor dimenions
   auto dims = mat.dimensions();
@@ -27,38 +27,6 @@ py::array_t<T> make_ndarray(Eigen::Tensor<T, 3>& mat)
      );
 }
 
-template <typename T>
-Eigen::Tensor<T, 3> make_tensor3(py::array_t<T> array)
-{
-  // Request a buffer descriptor from Python
-  py::buffer_info buffer_info = array.request();
-
-  // Get the array dimenions
-  T *data = static_cast<T *>(buffer_info.ptr);
-  std::vector<ssize_t> shape = buffer_info.shape;
-
-  // Check rank
-  if (shape.size() != 3) {
-    std::ostringstream sout;
-    sout << "make_tensor3: tensor rank must be 3, found "
-	 << shape.size();
-    throw std::runtime_error(sout.str());
-  }
-
-  // Build result tensor with col-major ordering
-  Eigen::Tensor<T, 3> tensor(shape[0], shape[1], shape[2]);
-  for (int i=0, l=0; i < shape[0]; i++) {
-    for (int j=0; j < shape[1]; j++) {
-      for (int k=0; k < shape[2]; k++) {
-	tensor(i, j, k) = data[l++];
-      }
-    }
-  }
-
-  return tensor;
-}
-
-
 //! Helper function that maps the Eigen::Tensor<T, 4> into an numpy.ndarray
 template <typename T>
 py::array_t<T> make_ndarray4(Eigen::Tensor<T, 4>& mat)
@@ -85,6 +53,40 @@ py::array_t<T> make_ndarray4(Eigen::Tensor<T, 4>& mat)
      );
 }
 
+template <typename T>
+Eigen::Tensor<T, 3> make_tensor3(py::array_t<T>& in)
+{
+  // Request a buffer descriptor from Python
+  py::buffer_info buffer_info = in.request();
+
+  // Extract the data 
+  T *data = static_cast<T *>(buffer_info.ptr);
+
+  // Get the data shape
+  std::vector<ssize_t> shape = buffer_info.shape;
+
+  // Check rank
+  if (shape.size() != 3) {
+    std::ostringstream sout;
+    sout << "make_tensor3: tensor rank must be 3, found "
+	 << shape.size();
+    throw std::runtime_error(sout.str());
+  }
+
+  // Reorder the data to satisfy the col-major Eigen::Tensor ordering
+  //
+  Eigen::Tensor<T, 3> tensor(shape[0], shape[1], shape[2]);
+  for (int i=0, c=0; i < shape[0]; i++) {
+    for (int j=0; j < shape[1]; j++) {
+      for (int k=0; k < shape[2]; k++, c++) {
+	tensor(i, j, k) = data[c];
+      }
+    }
+  }
+
+  return tensor;
+}
+
 template <typename T>
 Eigen::Tensor<T, 4> make_tensor4(py::array_t<T> array)
 {

From 01516c3556bd12f07c8364839bd1bfd0170f900c Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Sat, 28 Sep 2024 11:46:12 -0400
Subject: [PATCH 04/52] Updates for setMatrix() method for field coefficient
 types which need to set tensors

---
 expui/Coefficients.H  | 31 +++++++++++++++-------------
 expui/Coefficients.cc | 44 ++++++++++++++++++++--------------------
 pyEXP/CoefWrappers.cc | 47 +++++++++++++++++++++----------------------
 3 files changed, 62 insertions(+), 60 deletions(-)

diff --git a/expui/Coefficients.H b/expui/Coefficients.H
index 34a9c3a89..7718c7b7b 100644
--- a/expui/Coefficients.H
+++ b/expui/Coefficients.H
@@ -111,7 +111,7 @@ namespace CoefClasses
     virtual Eigen::VectorXcd& getData(double time) = 0;
 
     //! Set coefficient store at given time to the provided matrix
-    virtual void setData(double time, Eigen::Ref<Eigen::VectorXcd> data) = 0;
+    virtual void setData(double time, Eigen::VectorXcd& data) = 0;
 
     //! Interpolate coefficient matrix at given time
     std::tuple<Eigen::VectorXcd&, bool> interpolate(double time);
@@ -273,10 +273,10 @@ namespace CoefClasses
     Eigen::MatrixXcd& getMatrix(double time);
 
     //! Set coefficient matrix at given time
-    virtual void setData(double time, Eigen::Ref<Eigen::VectorXcd> arr);
+    virtual void setData(double time, Eigen::VectorXcd& arr);
 
     //! Natural data setter for pyEXP
-    void setMatrix(double time, Eigen::Ref<Eigen::MatrixXcd> mat);
+    void setMatrix(double time, Eigen::MatrixXcd& mat);
 
     //! Interpolate coefficient matrix at given time
     std::tuple<Eigen::MatrixXcd&, bool> interpolate(double time)
@@ -411,10 +411,10 @@ namespace CoefClasses
     Eigen::MatrixXcd& getMatrix(double time);
     
     //! Set coefficient matrix at given time
-    virtual void setData(double time, Eigen::Ref<Eigen::VectorXcd> arr);
+    virtual void setData(double time, Eigen::VectorXcd& arr);
     
     //! For pyEXP
-    void setMatrix(double time, Eigen::Ref<Eigen::MatrixXcd> arr);
+    void setMatrix(double time, Eigen::MatrixXcd& arr);
     
     //! Interpolate coefficient matrix at given time
     std::tuple<Eigen::MatrixXcd&, bool> interpolate(double time)
@@ -542,7 +542,7 @@ namespace CoefClasses
     virtual Eigen3d& getTensor(double time);
     
     //! Set coefficient matrix at given time
-    virtual void setData(double time, Eigen::Ref<Eigen::VectorXcd> dat);
+    virtual void setData(double time, Eigen::VectorXcd& dat);
 
     //! Native version
     virtual void setTensor(double time, const Eigen3d& dat);
@@ -673,7 +673,7 @@ namespace CoefClasses
     virtual Eigen3d& getTensor(double time);
     
     //! Set coefficient matrix at given time
-    virtual void setData(double time, Eigen::Ref<Eigen::VectorXcd> dat);
+    virtual void setData(double time, Eigen::VectorXcd& dat);
 
     //! Native version
     virtual void setTensor(double time, const Eigen3d& dat);
@@ -810,7 +810,10 @@ namespace CoefClasses
     /** Set coefficient matrix at given time.  This is for pybind11,
 	since the operator() will not allow lvalue assignment, it
 	seems. */
-    virtual void setData(double time, Eigen::Ref<Eigen::VectorXcd> arr);
+    virtual void setData(double time, Eigen::VectorXcd& arr);
+
+    void setArray(double time, Eigen::VectorXcd& arr)
+    { setData(time, arr); }
 
     //! Get coefficient structure at a given time
     virtual std::shared_ptr<CoefStruct> getCoefStruct(double time)
@@ -909,13 +912,13 @@ namespace CoefClasses
     virtual Eigen::VectorXcd& getData(double time);
     
     //! Operator for pyEXP
-    SphFldStruct::coefType& getMatrix(double time);
+    SphFldStruct::dataType getMatrix(double time);
 
     //! Set coefficient matrix at given time
-    virtual void setData(double time, Eigen::Ref<Eigen::VectorXcd> arr);
+    virtual void setData(double time, Eigen::VectorXcd& arr);
 
     //! Natural data setter for pyEXP
-    void setMatrix(double time, const SphFldStruct::coefType& mat);
+    void setMatrix(double time, SphFldStruct::dataType& mat);
 
     //! Interpolate coefficient matrix at given time
     std::tuple<SphFldStruct::coefType&, bool> interpolate(double time)
@@ -1036,13 +1039,13 @@ namespace CoefClasses
     virtual Eigen::VectorXcd& getData(double time);
     
     //! Operator for pyEXP
-    CylFldStruct::coefType & getMatrix(double time);
+    CylFldStruct::dataType getMatrix(double time);
 
     //! Set coefficient matrix at given time
-    virtual void setData(double time, Eigen::Ref<Eigen::VectorXcd> arr);
+    virtual void setData(double time, Eigen::VectorXcd& arr);
 
     //! Natural data setter for pyEXP
-    void setMatrix(double time, const CylFldStruct::coefType& mat);
+    void setMatrix(double time, CylFldStruct::dataType& mat);
 
     //! Interpolate coefficient matrix at given time
     std::tuple<CylFldStruct::coefType&, bool> interpolate(double time)
diff --git a/expui/Coefficients.cc b/expui/Coefficients.cc
index dfd58e5b8..fe55b0efe 100644
--- a/expui/Coefficients.cc
+++ b/expui/Coefficients.cc
@@ -469,13 +469,13 @@ namespace CoefClasses
     return mat;
   }
   
-  void SphCoefs::setData(double time, Eigen::Ref<Eigen::VectorXcd> dat)
+  void SphCoefs::setData(double time, Eigen::VectorXcd& dat)
   {
     auto it = coefs.find(roundTime(time));
 
     if (it == coefs.end()) {
       std::ostringstream str;
-      str << "SphCoefs::setMatrix: requested time=" << time << " not found";
+      str << "SphCoefs::setData: requested time=" << time << " not found";
       throw std::runtime_error(str.str());
     } else {
       it->second->store = dat;
@@ -484,7 +484,7 @@ namespace CoefClasses
     }
   }
   
-  void SphCoefs::setMatrix(double time, Eigen::Ref<Eigen::MatrixXcd> dat)
+  void SphCoefs::setMatrix(double time, Eigen::MatrixXcd& dat)
   {
     auto it = coefs.find(roundTime(time));
 
@@ -891,13 +891,13 @@ namespace CoefClasses
     return mat;
   }
   
-  void CylCoefs::setData(double time, Eigen::Ref<Eigen::VectorXcd> dat)
+  void CylCoefs::setData(double time, Eigen::VectorXcd& dat)
   {
     auto it = coefs.find(roundTime(time));
 
     if (it == coefs.end()) {
       std::ostringstream str;
-      str << "CylCoefs::setMatrix: requested time=" << time << " not found";
+      str << "CylCoefs::setData: requested time=" << time << " not found";
       throw std::runtime_error(str.str());
     } else {
       it->second->store = dat;
@@ -906,7 +906,7 @@ namespace CoefClasses
     }
   }
 
-  void CylCoefs::setMatrix(double time, Eigen::Ref<Eigen::MatrixXcd> dat)
+  void CylCoefs::setMatrix(double time, Eigen::MatrixXcd& dat)
   {
     auto it = coefs.find(roundTime(time));
 
@@ -1247,13 +1247,13 @@ namespace CoefClasses
     return arr;
   }
 
-  void SlabCoefs::setData(double time, Eigen::Ref<Eigen::VectorXcd> dat)
+  void SlabCoefs::setData(double time, Eigen::VectorXcd& dat)
   {
     auto it = coefs.find(roundTime(time));
 
     if (it == coefs.end()) {
       std::ostringstream str;
-      str << "CylCoefs::setMatrix: requested time=" << time << " not found";
+      str << "CylCoefs::setData: requested time=" << time << " not found";
       throw std::runtime_error(str.str());
     } else {
       it->second->store = dat;
@@ -1268,7 +1268,7 @@ namespace CoefClasses
 
     if (it == coefs.end()) {
       std::ostringstream str;
-      str << "CylCoefs::setMatrix: requested time=" << time << " not found";
+      str << "SlabCoefs::setTensor: requested time=" << time << " not found";
       throw std::runtime_error(str.str());
     } else {
       it->second->allocate();	// Assign storage for the flattened tensor
@@ -1599,13 +1599,13 @@ namespace CoefClasses
     return arr;
   }
 
-  void CubeCoefs::setData(double time, Eigen::Ref<Eigen::VectorXcd> dat)
+  void CubeCoefs::setData(double time, Eigen::VectorXcd& dat)
   {
     auto it = coefs.find(roundTime(time));
 
     if (it == coefs.end()) {
       std::ostringstream str;
-      str << "CylCoefs::setMatrix: requested time=" << time << " not found";
+      str << "CubeCoefs::setData: requested time=" << time << " not found";
       throw std::runtime_error(str.str());
     } else {
       it->second->store = dat;
@@ -1620,7 +1620,7 @@ namespace CoefClasses
 
     if (it == coefs.end()) {
       std::ostringstream str;
-      str << "CylCoefs::setMatrix: requested time=" << time << " not found";
+      str << "CubeCoefs::setTensor: requested time=" << time << " not found";
       throw std::runtime_error(str.str());
     } else {
       it->second->allocate();	// Assign storage for the flattened tensor
@@ -1968,13 +1968,13 @@ namespace CoefClasses
     return arr;
   }
   
-  void TableData::setData(double time, Eigen::Ref<Eigen::VectorXcd> dat)
+  void TableData::setData(double time, Eigen::VectorXcd& dat)
   {
     auto it = coefs.find(roundTime(time));
 
     if (it == coefs.end()) {
       std::ostringstream str;
-      str << "TableData::setMatrix: requested time=" << time << " not found";
+      str << "TableData::setData: requested time=" << time << " not found";
       throw std::runtime_error(str.str());
     } else {
       it->second->store = dat;
@@ -2442,7 +2442,7 @@ namespace CoefClasses
     return arr;
   }
   
-  SphFldStruct::coefType & SphFldCoefs::getMatrix(double time)
+  SphFldStruct::dataType SphFldCoefs::getMatrix(double time)
   {
     auto it = coefs.find(roundTime(time));
 
@@ -2458,13 +2458,13 @@ namespace CoefClasses
     return *mat;
   }
   
-  void SphFldCoefs::setData(double time, Eigen::Ref<Eigen::VectorXcd> dat)
+  void SphFldCoefs::setData(double time, Eigen::VectorXcd& dat)
   {
     auto it = coefs.find(roundTime(time));
 
     if (it == coefs.end()) {
       std::ostringstream str;
-      str << "SphVelCoefs::setMatrix: requested time=" << time << " not found";
+      str << "SphFldCoefs::setData: requested time=" << time << " not found";
       throw std::runtime_error(str.str());
     } else {
       it->second->store = dat;
@@ -2473,13 +2473,13 @@ namespace CoefClasses
     }
   }
   
-  void SphFldCoefs::setMatrix(double time, const SphFldStruct::coefType& dat)
+  void SphFldCoefs::setMatrix(double time, SphFldStruct::dataType& dat)
   {
     auto it = coefs.find(roundTime(time));
 
     if (it == coefs.end()) {
       std::ostringstream str;
-      str << "SphVelCoefs::setMatrix: requested time=" << time << " not found";
+      str << "SphFldCoefs::setMatrix: requested time=" << time << " not found";
       throw std::runtime_error(str.str());
     } else {
       it->second->allocate();
@@ -2882,7 +2882,7 @@ namespace CoefClasses
     return arr;
   }
   
-  CylFldStruct::coefType & CylFldCoefs::getMatrix(double time)
+  CylFldStruct::dataType CylFldCoefs::getMatrix(double time)
   {
     auto it = coefs.find(roundTime(time));
 
@@ -2897,7 +2897,7 @@ namespace CoefClasses
     return *mat;
   }
   
-  void CylFldCoefs::setData(double time, Eigen::Ref<Eigen::VectorXcd> dat)
+  void CylFldCoefs::setData(double time, Eigen::VectorXcd& dat)
   {
     auto it = coefs.find(roundTime(time));
 
@@ -2912,7 +2912,7 @@ namespace CoefClasses
     }
   }
   
-  void CylFldCoefs::setMatrix(double time, const CylFldStruct::coefType& dat)
+  void CylFldCoefs::setMatrix(double time, CylFldStruct::dataType& dat)
   {
     auto it = coefs.find(roundTime(time));
 
diff --git a/pyEXP/CoefWrappers.cc b/pyEXP/CoefWrappers.cc
index d170cfd33..aa58603b9 100644
--- a/pyEXP/CoefWrappers.cc
+++ b/pyEXP/CoefWrappers.cc
@@ -155,7 +155,7 @@ void CoefficientClasses(py::module &m) {
       PYBIND11_OVERRIDE_PURE(Eigen::VectorXcd&, Coefs, getData, time);
     }
 
-    void setData(double time, Eigen::Ref<Eigen::VectorXcd> array) override {
+    void setData(double time, Eigen::VectorXcd& array) override {
       PYBIND11_OVERRIDE_PURE(void, Coefs, setData, time, array);
     }
 
@@ -231,7 +231,7 @@ void CoefficientClasses(py::module &m) {
       PYBIND11_OVERRIDE(Eigen::VectorXcd&, SphCoefs, getData, time);
     }
 
-    void setData(double time, Eigen::Ref<Eigen::VectorXcd> array) override {
+    void setData(double time, Eigen::VectorXcd& array) override {
       PYBIND11_OVERRIDE(void, SphCoefs, setData, time, array);
     }
 
@@ -310,7 +310,7 @@ void CoefficientClasses(py::module &m) {
       PYBIND11_OVERRIDE(Eigen::VectorXcd&, CylCoefs, getData, time);
     }
 
-    void setData(double time, Eigen::Ref<Eigen::VectorXcd> array) override {
+    void setData(double time, Eigen::VectorXcd& array) override {
       PYBIND11_OVERRIDE(void, CylCoefs, setData, time, array);
     }
 
@@ -388,7 +388,7 @@ void CoefficientClasses(py::module &m) {
       PYBIND11_OVERRIDE(Eigen::VectorXcd&, SlabCoefs, getData, time);
     }
 
-    void setData(double time, Eigen::Ref<Eigen::VectorXcd> array) override {
+    void setData(double time, Eigen::VectorXcd& array) override {
       PYBIND11_OVERRIDE(void, SlabCoefs, setData, time, array);
     }
 
@@ -467,7 +467,7 @@ void CoefficientClasses(py::module &m) {
       PYBIND11_OVERRIDE(Eigen::VectorXcd&, CubeCoefs, getData, time);
     }
 
-    void setData(double time, Eigen::Ref<Eigen::VectorXcd> array) override {
+    void setData(double time, Eigen::VectorXcd& array) override {
       PYBIND11_OVERRIDE(void, CubeCoefs, setData, time, array);
     }
 
@@ -546,7 +546,7 @@ void CoefficientClasses(py::module &m) {
       PYBIND11_OVERRIDE(Eigen::VectorXcd&, TableData, getData, time);
     }
 
-    void setData(double time, Eigen::Ref<Eigen::VectorXcd> array) override {
+    void setData(double time, Eigen::VectorXcd& array) override {
       PYBIND11_OVERRIDE(void, TableData, setData, time, array);
     }
 
@@ -1185,7 +1185,7 @@ void CoefficientClasses(py::module &m) {
 	 [](CoefClasses::SphCoefs& A)
 	 {
 	   auto M = A.getAllCoefs(); // Need a copy here
-	   py::array_t<std::complex<double>> ret = make_ndarray<std::complex<double>>(M);
+	   py::array_t<std::complex<double>> ret = make_ndarray3<std::complex<double>>(M);
 	   return ret;
 	 },
 	 R"(
@@ -1260,7 +1260,7 @@ void CoefficientClasses(py::module &m) {
 	 [](CoefClasses::CylCoefs& A)
 	 {
 	   auto M = A.getAllCoefs(); // Need a copy here
-	   py::array_t<std::complex<double>> ret = make_ndarray<std::complex<double>>(M);
+	   py::array_t<std::complex<double>> ret = make_ndarray3<std::complex<double>>(M);
 	   return ret;
 	 },
 	 R"(
@@ -1321,11 +1321,11 @@ void CoefficientClasses(py::module &m) {
          SphFldCoefs instance
          )")
     .def("__call__",
-	 [](CoefClasses::SphFldCoefs& A, double t)
+	 [](CoefClasses::SphFldCoefs& A, double time)
 	 {
-	   Eigen::Tensor<std::complex<double>, 3> M = A.getMatrix(t);
-	   py::array_t<std::complex<double>> ret = make_ndarray<std::complex<double>>(M);
-	   return ret;
+	   // Need a copy here
+	   auto M = A.getMatrix(time);
+	   return make_ndarray3<std::complex<double>>(M);
 	 },
          R"(
          Return the coefficient tensor for the desired time.
@@ -1347,11 +1347,11 @@ void CoefficientClasses(py::module &m) {
          )",
          py::arg("time"))
     .def("setMatrix",
-	 [](CoefClasses::SphFldCoefs& A, double t,
-	    py::array_t<std::complex<double>> array)
+	 [](CoefClasses::SphFldCoefs& A, double time,
+	    py::array_t<std::complex<double>> mat)
 	 {
-	   auto M = make_tensor3<std::complex<double>>(array);
-	   A.setMatrix(t, M);
+	   auto M = make_tensor3<std::complex<double>>(mat);
+	   A.setMatrix(time, M);
 	 },
          R"(
          Enter and/or rewrite the coefficient tensor at the provided time
@@ -1406,11 +1406,10 @@ void CoefficientClasses(py::module &m) {
          CylFldCoefs instance
          )")
     .def("__call__",
-	 [](CoefClasses::CylFldCoefs& A, double t)
+	 [](CoefClasses::CylFldCoefs& A, double time)
 	 {
-	   Eigen::Tensor<std::complex<double>, 3> M = A.getMatrix(t);
-	   py::array_t<std::complex<double>> ret = make_ndarray<std::complex<double>>(M);
-	   return ret;
+	   auto M = A.getMatrix(time); // Need a copy here
+	   return make_ndarray3<std::complex<double>>(M);
 	 },
          R"(
          Return the coefficient tensor for the desired time.
@@ -1432,11 +1431,11 @@ void CoefficientClasses(py::module &m) {
          )",
          py::arg("time"))
     .def("setMatrix",
-	 [](CoefClasses::CylFldCoefs& A, double t,
-	    py::array_t<std::complex<double>> array)
+	 [](CoefClasses::CylFldCoefs& A, double time,
+	    py::array_t<std::complex<double>> mat)
 	 {
-	   auto M = make_tensor3<std::complex<double>>(array);
-	   A.setMatrix(t, M);
+	   auto M = make_tensor3<std::complex<double>>(mat);
+	   A.setMatrix(time, M);
 	 },
          R"(
          Enter and/or rewrite the coefficient tensor at the provided time

From 14a6ef5e3915d92e330e3e3e7011990e1a9e94b6 Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Sat, 12 Oct 2024 10:58:31 -0400
Subject: [PATCH 05/52] Added CTests for 'setMatrix()' and 'createFromArray'
 [no ci]

---
 tests/CMakeLists.txt      | 21 ++++++++++--
 tests/Halo/changeCoefs.py | 67 +++++++++++++++++++++++++++++++++++++++
 tests/Halo/createCoefs.py | 46 +++++++++++++++++++++++++++
 3 files changed, 131 insertions(+), 3 deletions(-)
 create mode 100644 tests/Halo/changeCoefs.py
 create mode 100644 tests/Halo/createCoefs.py

diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 9aee83328..0d73608f2 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -73,10 +73,10 @@ if(ENABLE_NBODY)
     ${PYTHON_EXECUTABLE} check.py
     WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo)
 
-  set_tests_properties(expNbodyCheck2TW PROPERTIES DEPENDS makeNbodyTest)
+  set_tests_properties(expNbodyCheck2TW PROPERTIES DEPENDS expNbodyTest)
 
   # This adds a coefficient read test using pyEXP only if
-  # makeNbodyTest is run and pyEXP has been built
+  # expNbodyTest is run and pyEXP has been built
   if(ENABLE_PYEXP)
      # Read coefficient file with pyEXP
      add_test(NAME pyEXPCoefReadTest
@@ -84,8 +84,23 @@ if(ENABLE_NBODY)
               PYTHONPATH=${CMAKE_BINARY_DIR}/pyEXP:$ENV{PYTHONPATH}
               ${PYTHON_EXECUTABLE} readCoefs.py
               WORKING_DIRECTORY "${CMAKE_CURRENT_LIST_DIR}/Halo")
-     set_tests_properties(pyEXPCoefReadTest PROPERTIES DEPENDS makeNbodyTest)
+     set_tests_properties(pyEXPCoefReadTest PROPERTIES DEPENDS expNbodyTest)
      set_tests_properties(pyEXPCoefReadTest PROPERTIES LABELS "long")
+
+     add_test(NAME pyEXPCoefMatrixTest
+              COMMAND ${CMAKE_COMMAND} -E env
+              PYTHONPATH=${CMAKE_BINARY_DIR}/pyEXP:$ENV{PYTHONPATH}
+              ${PYTHON_EXECUTABLE} changeCoefs.py
+              WORKING_DIRECTORY "${CMAKE_CURRENT_LIST_DIR}/Halo")
+     set_tests_properties(pyEXPCoefMatrixTest PROPERTIES DEPENDS expNbodyTest)
+     set_tests_properties(pyEXPCoefMatrixTest PROPERTIES LABELS "long")
+
+     add_test(NAME pyEXPCoefCreateTest
+              COMMAND ${CMAKE_COMMAND} -E env
+              PYTHONPATH=${CMAKE_BINARY_DIR}/pyEXP:$ENV{PYTHONPATH}
+              ${PYTHON_EXECUTABLE} createCoefs.py
+              WORKING_DIRECTORY "${CMAKE_CURRENT_LIST_DIR}/Halo")
+     set_tests_properties(pyEXPCoefCreateTest PROPERTIES LABELS "quit")
   endif()
 
   # A separate test to remove the generated files if they all exist;
diff --git a/tests/Halo/changeCoefs.py b/tests/Halo/changeCoefs.py
new file mode 100644
index 000000000..d9ebc92e3
--- /dev/null
+++ b/tests/Halo/changeCoefs.py
@@ -0,0 +1,67 @@
+import os
+import pyEXP
+import numpy as np
+
+# Read the coefs from the file created by the expNbodyTest
+#
+coefs = pyEXP.coefs.Coefs.factory('outcoef.halo.run0')
+
+print("Got coefs for name=", coefs.getName())
+
+
+times = coefs.Times()
+data  = coefs.getAllCoefs()
+
+print("The data sizes and coefficient orders are:")
+print("M orders, N orders, Times:", data.shape)
+
+
+# Make the halo basis config
+config="""
+---
+id : sphereSL
+parameters :
+  numr: 4000
+  rmin: 0.0001
+  rmax: 1.95
+  Lmax: 2
+  nmax: 10
+  rmapping : 0.0667
+  self_consistent: true
+  modelname: SLGridSph.model
+  cachename: SLGridSph.cache.run0
+...
+"""
+
+# Construct the basis instances
+#
+basis = pyEXP.basis.Basis.factory(config)
+
+# Let's zero all of the odd order data
+#
+for k in range(data.shape[0]):
+    l, m = pyEXP.basis.SphericalSL.invI(k)
+    if l%2!=0 or m%2!=0: data[k,:,:] *= 0.0
+
+# Reset the coefficient data
+#
+for i in range(data.shape[2]):
+    coefs.setMatrix(times[i], data[:,:,i])
+
+# Check that odd coefficients are really zero
+data1 = coefs.getAllCoefs()
+minZero =  1.0e30
+maxZero = -1.0e30
+for k in range(data1.shape[0]):
+    l, m = pyEXP.basis.SphericalSL.invI(k)
+    if l%2!=0 or m%2!=0:
+        for v in data1[k,:,:].flatten():
+           minZero = min([minZero, abs(v)])
+           maxZero = max([maxZero, abs(v)])
+
+print('Zero test: min, max: {:13.7e}, {:13.7e}'.format(minZero, maxZero))
+
+if minZero < -1.0e-18 or maxZero > 1.0e-18:
+    exit(1)
+else:
+    exit(0)
diff --git a/tests/Halo/createCoefs.py b/tests/Halo/createCoefs.py
new file mode 100644
index 000000000..126f6ad1e
--- /dev/null
+++ b/tests/Halo/createCoefs.py
@@ -0,0 +1,46 @@
+import os
+import pyEXP
+import numpy as np
+
+# Make the halo basis config
+config="""
+---
+id : sphereSL
+parameters :
+  numr: 4000
+  rmin: 0.0001
+  rmax: 1.95
+  Lmax: 2
+  nmax: 10
+  rmapping : 0.0667
+  self_consistent: true
+  modelname: SLGridSph.model
+  cachename: SLGridSph.cache.run0
+...
+"""
+
+# Construct the basis instances
+#
+basis = pyEXP.basis.Basis.factory(config)
+
+# Open the file and read the array.  The file was generated using
+# psp2ascii.
+#
+bodyfile = 'new.bods'
+if not os.path.exists(bodyfile):
+    print('Body file <{}> does not exist'.format(bodyfile))
+    exit(1)
+
+data = np.loadtxt(bodyfile, skiprows=1, usecols=(1, 2, 3, 4))
+print(data.shape)
+
+# Call the basis to generate coefficients
+#
+coef = basis.createFromArray(data[:,0], data[:,1:4], time=3.0)
+
+# Add the coefficient structure coefficient structure
+coefs = pyEXP.coefs.SphCoefs(True)
+coefs.add(coef)
+print("Times:", coefs.Times())
+
+exit(0)

From 1a65ec688c74aced01764fcfec8a09573aa06b60 Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Sat, 12 Oct 2024 11:42:54 -0400
Subject: [PATCH 06/52] Typo in test label

---
 tests/CMakeLists.txt | 10 ++++------
 1 file changed, 4 insertions(+), 6 deletions(-)

diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 0d73608f2..80fe838fe 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -73,7 +73,7 @@ if(ENABLE_NBODY)
     ${PYTHON_EXECUTABLE} check.py
     WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo)
 
-  set_tests_properties(expNbodyCheck2TW PROPERTIES DEPENDS expNbodyTest)
+  set_tests_properties(expNbodyCheck2TW PROPERTIES DEPENDS expNbodyTest LABELS "long")
 
   # This adds a coefficient read test using pyEXP only if
   # expNbodyTest is run and pyEXP has been built
@@ -84,23 +84,21 @@ if(ENABLE_NBODY)
               PYTHONPATH=${CMAKE_BINARY_DIR}/pyEXP:$ENV{PYTHONPATH}
               ${PYTHON_EXECUTABLE} readCoefs.py
               WORKING_DIRECTORY "${CMAKE_CURRENT_LIST_DIR}/Halo")
-     set_tests_properties(pyEXPCoefReadTest PROPERTIES DEPENDS expNbodyTest)
-     set_tests_properties(pyEXPCoefReadTest PROPERTIES LABELS "long")
+     set_tests_properties(pyEXPCoefReadTest PROPERTIES DEPENDS expNbodyTest LABELS "long")
 
      add_test(NAME pyEXPCoefMatrixTest
               COMMAND ${CMAKE_COMMAND} -E env
               PYTHONPATH=${CMAKE_BINARY_DIR}/pyEXP:$ENV{PYTHONPATH}
               ${PYTHON_EXECUTABLE} changeCoefs.py
               WORKING_DIRECTORY "${CMAKE_CURRENT_LIST_DIR}/Halo")
-     set_tests_properties(pyEXPCoefMatrixTest PROPERTIES DEPENDS expNbodyTest)
-     set_tests_properties(pyEXPCoefMatrixTest PROPERTIES LABELS "long")
+     set_tests_properties(pyEXPCoefMatrixTest PROPERTIES DEPENDS expNbodyTest LABELS "long")
 
      add_test(NAME pyEXPCoefCreateTest
               COMMAND ${CMAKE_COMMAND} -E env
               PYTHONPATH=${CMAKE_BINARY_DIR}/pyEXP:$ENV{PYTHONPATH}
               ${PYTHON_EXECUTABLE} createCoefs.py
               WORKING_DIRECTORY "${CMAKE_CURRENT_LIST_DIR}/Halo")
-     set_tests_properties(pyEXPCoefCreateTest PROPERTIES LABELS "quit")
+     set_tests_properties(pyEXPCoefCreateTest PROPERTIES LABELS "quick")
   endif()
 
   # A separate test to remove the generated files if they all exist;

From eb4bfd37dd6b43162e6521a08e162f66ad238a02 Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Sat, 12 Oct 2024 12:08:55 -0400
Subject: [PATCH 07/52] Use random data rather than depend on a body file

---
 tests/Halo/createCoefs.py | 49 ++++++++++++++++++++++++++++++---------
 1 file changed, 38 insertions(+), 11 deletions(-)

diff --git a/tests/Halo/createCoefs.py b/tests/Halo/createCoefs.py
index 126f6ad1e..20ffbfa6a 100644
--- a/tests/Halo/createCoefs.py
+++ b/tests/Halo/createCoefs.py
@@ -1,5 +1,6 @@
 import os
 import pyEXP
+import random
 import numpy as np
 
 # Make the halo basis config
@@ -23,24 +24,50 @@
 #
 basis = pyEXP.basis.Basis.factory(config)
 
-# Open the file and read the array.  The file was generated using
-# psp2ascii.
+# Make some fake body data
 #
-bodyfile = 'new.bods'
-if not os.path.exists(bodyfile):
-    print('Body file <{}> does not exist'.format(bodyfile))
-    exit(1)
-
-data = np.loadtxt(bodyfile, skiprows=1, usecols=(1, 2, 3, 4))
-print(data.shape)
 
 # Call the basis to generate coefficients
 #
-coef = basis.createFromArray(data[:,0], data[:,1:4], time=3.0)
+mass = []
+xpos = []
+ypos = []
+zpos = []
+
+for i in range(0,100):
+    mass.append(0.001)
+    xpos.append(random.random()*2.0 - 1.0)
+    ypos.append(random.random()*2.0 - 1.0)
+    zpos.append(random.random()*2.0 - 1.0)
+
+print("---- createFromArray usings lists")
+coef1 = basis.createFromArray(mass, [xpos, ypos, zpos], time=3.0)
+
+mass  = np.array(mass)
+data  = np.array([xpos, ypos, zpos])
+
+print("---- createFromArray usings numpy arrays converted from lists")
+coef2 = basis.createFromArray(mass, data, time=3.1)
+
+mass = np.ones(1000) * 1.0e6
+xpos = np.random.normal(0.0, 1.0, 1000)
+ypos = np.random.normal(0.0, 1.0, 1000)
+zpos = np.random.normal(0.0, 1.0, 1000)
+
+print("---- createFromArray using a list of numpy arrays")
+coef3 = basis.createFromArray(mass, [xpos, ypos, zpos], time=3.2)
+
+data  = np.array([xpos, ypos, zpos])
+
+print("---- createFromArray usings pure numpy arrays")
+coef4 = basis.createFromArray(mass, data, time=3.3)
 
 # Add the coefficient structure coefficient structure
 coefs = pyEXP.coefs.SphCoefs(True)
-coefs.add(coef)
+coefs.add(coef1)
+coefs.add(coef2)
+coefs.add(coef3)
+coefs.add(coef4)
 print("Times:", coefs.Times())
 
 exit(0)

From b81da9a199f6e4d21bce96f1293cb597e99aff89 Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Sat, 12 Oct 2024 12:35:19 -0400
Subject: [PATCH 08/52] Get some verbose output from github runner

---
 tests/CMakeLists.txt | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 80fe838fe..560f9e8b5 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -1,5 +1,7 @@
 include(CTest)
 
+set(CTEST_OUTPUT_ON_FAILURE ON)
+
 # Is EXP configured for pyEXP? If yes, run pyEXP tests...
 #
 if(ENABLE_PYEXP)

From 88d133871771b9ae8609266c834cd95757f48998 Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Sat, 12 Oct 2024 13:15:10 -0400
Subject: [PATCH 09/52] Make ctest log the output on failure

---
 .github/workflows/build.yml | 2 +-
 tests/Halo/createCoefs.py   | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
index b9ddde223..d6b08ab84 100644
--- a/.github/workflows/build.yml
+++ b/.github/workflows/build.yml
@@ -54,7 +54,7 @@ jobs:
 
       - name: CTest Quick
         working-directory: ./build
-        run: ctest -L quick
+        run: ctest --output-on-failure -L quick
         
       #- name: CTest Long
       #working-directory: ./build
diff --git a/tests/Halo/createCoefs.py b/tests/Halo/createCoefs.py
index 20ffbfa6a..5ce71f481 100644
--- a/tests/Halo/createCoefs.py
+++ b/tests/Halo/createCoefs.py
@@ -70,4 +70,4 @@
 coefs.add(coef4)
 print("Times:", coefs.Times())
 
-exit(0)
+exit(1)

From 41db6184222c7916baa83936d9299f158f80103b Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Sat, 12 Oct 2024 13:32:28 -0400
Subject: [PATCH 10/52] Remove intentional failure from test script

---
 tests/Halo/createCoefs.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/Halo/createCoefs.py b/tests/Halo/createCoefs.py
index 5ce71f481..20ffbfa6a 100644
--- a/tests/Halo/createCoefs.py
+++ b/tests/Halo/createCoefs.py
@@ -70,4 +70,4 @@
 coefs.add(coef4)
 print("Times:", coefs.Times())
 
-exit(1)
+exit(0)

From b9377955d1c8aea79fe667471b3de560f26c22ac Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Sat, 12 Oct 2024 18:46:36 -0400
Subject: [PATCH 11/52] Add another print statement

---
 tests/Halo/createCoefs.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/tests/Halo/createCoefs.py b/tests/Halo/createCoefs.py
index 20ffbfa6a..839a101a4 100644
--- a/tests/Halo/createCoefs.py
+++ b/tests/Halo/createCoefs.py
@@ -34,7 +34,8 @@
 ypos = []
 zpos = []
 
-for i in range(0,100):
+print("---- creating list data")
+for i in range(0, 100):
     mass.append(0.001)
     xpos.append(random.random()*2.0 - 1.0)
     ypos.append(random.random()*2.0 - 1.0)

From dd78429f39c516ebe80e46038e4bb2a23e4172d5 Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Sat, 12 Oct 2024 19:40:31 -0400
Subject: [PATCH 12/52] Add another print statement

---
 tests/Halo/createCoefs.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/tests/Halo/createCoefs.py b/tests/Halo/createCoefs.py
index 839a101a4..b85b57555 100644
--- a/tests/Halo/createCoefs.py
+++ b/tests/Halo/createCoefs.py
@@ -14,16 +14,19 @@
   Lmax: 2
   nmax: 10
   rmapping : 0.0667
-  self_consistent: true
   modelname: SLGridSph.model
   cachename: SLGridSph.cache.run0
 ...
 """
 
+print("---- about to create basis")
+
 # Construct the basis instances
 #
 basis = pyEXP.basis.Basis.factory(config)
 
+print("---- created basis")
+
 # Make some fake body data
 #
 

From e9600e46b7116d94542203070ee815fbefdd95cf Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Sat, 12 Oct 2024 20:39:29 -0400
Subject: [PATCH 13/52] Exit after basis creation; sanity check

---
 tests/Halo/createCoefs.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/tests/Halo/createCoefs.py b/tests/Halo/createCoefs.py
index b85b57555..085055b2c 100644
--- a/tests/Halo/createCoefs.py
+++ b/tests/Halo/createCoefs.py
@@ -27,6 +27,8 @@
 
 print("---- created basis")
 
+exit(0)
+
 # Make some fake body data
 #
 

From eaa6521a1d0245258855924de0424c4d32676727 Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Sat, 12 Oct 2024 21:04:52 -0400
Subject: [PATCH 14/52] Exit after data creation; sanity check

---
 tests/Halo/createCoefs.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/tests/Halo/createCoefs.py b/tests/Halo/createCoefs.py
index 085055b2c..33d7dcbdd 100644
--- a/tests/Halo/createCoefs.py
+++ b/tests/Halo/createCoefs.py
@@ -27,8 +27,6 @@
 
 print("---- created basis")
 
-exit(0)
-
 # Make some fake body data
 #
 
@@ -46,6 +44,8 @@
     ypos.append(random.random()*2.0 - 1.0)
     zpos.append(random.random()*2.0 - 1.0)
 
+exit(0)
+
 print("---- createFromArray usings lists")
 coef1 = basis.createFromArray(mass, [xpos, ypos, zpos], time=3.0)
 

From 6432520cebbcc96cf20c63c132ca990633c1f62b Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Sat, 12 Oct 2024 21:19:36 -0400
Subject: [PATCH 15/52] Exit after first coef creation; sanity check

---
 tests/Halo/createCoefs.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/tests/Halo/createCoefs.py b/tests/Halo/createCoefs.py
index 33d7dcbdd..64e688fdb 100644
--- a/tests/Halo/createCoefs.py
+++ b/tests/Halo/createCoefs.py
@@ -44,11 +44,11 @@
     ypos.append(random.random()*2.0 - 1.0)
     zpos.append(random.random()*2.0 - 1.0)
 
-exit(0)
-
 print("---- createFromArray usings lists")
 coef1 = basis.createFromArray(mass, [xpos, ypos, zpos], time=3.0)
 
+exit(0)
+
 mass  = np.array(mass)
 data  = np.array([xpos, ypos, zpos])
 

From 8d332a20f66f1f50d481166d1187edad0478ee5f Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Sat, 12 Oct 2024 21:30:47 -0400
Subject: [PATCH 16/52] Exit after third coef creation; sanity check

---
 tests/Halo/createCoefs.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/tests/Halo/createCoefs.py b/tests/Halo/createCoefs.py
index 64e688fdb..ce9aa5bb9 100644
--- a/tests/Halo/createCoefs.py
+++ b/tests/Halo/createCoefs.py
@@ -47,8 +47,6 @@
 print("---- createFromArray usings lists")
 coef1 = basis.createFromArray(mass, [xpos, ypos, zpos], time=3.0)
 
-exit(0)
-
 mass  = np.array(mass)
 data  = np.array([xpos, ypos, zpos])
 
@@ -63,6 +61,8 @@
 print("---- createFromArray using a list of numpy arrays")
 coef3 = basis.createFromArray(mass, [xpos, ypos, zpos], time=3.2)
 
+exit(0)
+
 data  = np.array([xpos, ypos, zpos])
 
 print("---- createFromArray usings pure numpy arrays")

From 2369691aaef6f5f03715f7ab4e303d41353f9f56 Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Sat, 12 Oct 2024 21:44:08 -0400
Subject: [PATCH 17/52] Exit after second coef creation; sanity check

---
 tests/Halo/createCoefs.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/tests/Halo/createCoefs.py b/tests/Halo/createCoefs.py
index ce9aa5bb9..72bbdd7f5 100644
--- a/tests/Halo/createCoefs.py
+++ b/tests/Halo/createCoefs.py
@@ -53,6 +53,8 @@
 print("---- createFromArray usings numpy arrays converted from lists")
 coef2 = basis.createFromArray(mass, data, time=3.1)
 
+exit(0)
+
 mass = np.ones(1000) * 1.0e6
 xpos = np.random.normal(0.0, 1.0, 1000)
 ypos = np.random.normal(0.0, 1.0, 1000)
@@ -61,8 +63,6 @@
 print("---- createFromArray using a list of numpy arrays")
 coef3 = basis.createFromArray(mass, [xpos, ypos, zpos], time=3.2)
 
-exit(0)
-
 data  = np.array([xpos, ypos, zpos])
 
 print("---- createFromArray usings pure numpy arrays")

From 747e6a5d58334afe4a4fa52a7e2b83b211059e23 Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Sat, 12 Oct 2024 23:14:16 -0400
Subject: [PATCH 18/52] Comment out second coef creation; sanity check

---
 tests/Halo/createCoefs.py | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/tests/Halo/createCoefs.py b/tests/Halo/createCoefs.py
index 72bbdd7f5..d7dd4bf34 100644
--- a/tests/Halo/createCoefs.py
+++ b/tests/Halo/createCoefs.py
@@ -47,13 +47,13 @@
 print("---- createFromArray usings lists")
 coef1 = basis.createFromArray(mass, [xpos, ypos, zpos], time=3.0)
 
-mass  = np.array(mass)
-data  = np.array([xpos, ypos, zpos])
+# mass  = np.array(mass)
+# data  = np.array([xpos, ypos, zpos])
 
-print("---- createFromArray usings numpy arrays converted from lists")
-coef2 = basis.createFromArray(mass, data, time=3.1)
+# print("---- createFromArray usings numpy arrays converted from lists")
+# coef2 = basis.createFromArray(mass, data, time=3.1)
 
-exit(0)
+# exit(0)
 
 mass = np.ones(1000) * 1.0e6
 xpos = np.random.normal(0.0, 1.0, 1000)
@@ -71,7 +71,7 @@
 # Add the coefficient structure coefficient structure
 coefs = pyEXP.coefs.SphCoefs(True)
 coefs.add(coef1)
-coefs.add(coef2)
+# coefs.add(coef2)
 coefs.add(coef3)
 coefs.add(coef4)
 print("Times:", coefs.Times())

From 10434365e4edd7142ce294d6b44e965420e05c35 Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Sat, 12 Oct 2024 23:30:46 -0400
Subject: [PATCH 19/52] Comment out all but first coef creation; sanity check

---
 tests/Halo/createCoefs.py | 22 +++++++++++-----------
 1 file changed, 11 insertions(+), 11 deletions(-)

diff --git a/tests/Halo/createCoefs.py b/tests/Halo/createCoefs.py
index d7dd4bf34..255d6c7b2 100644
--- a/tests/Halo/createCoefs.py
+++ b/tests/Halo/createCoefs.py
@@ -55,25 +55,25 @@
 
 # exit(0)
 
-mass = np.ones(1000) * 1.0e6
-xpos = np.random.normal(0.0, 1.0, 1000)
-ypos = np.random.normal(0.0, 1.0, 1000)
-zpos = np.random.normal(0.0, 1.0, 1000)
+# mass = np.ones(1000) * 1.0e6
+# xpos = np.random.normal(0.0, 1.0, 1000)
+# ypos = np.random.normal(0.0, 1.0, 1000)
+# zpos = np.random.normal(0.0, 1.0, 1000)
 
-print("---- createFromArray using a list of numpy arrays")
-coef3 = basis.createFromArray(mass, [xpos, ypos, zpos], time=3.2)
+# print("---- createFromArray using a list of numpy arrays")
+# coef3 = basis.createFromArray(mass, [xpos, ypos, zpos], time=3.2)
 
-data  = np.array([xpos, ypos, zpos])
+# data  = np.array([xpos, ypos, zpos])
 
-print("---- createFromArray usings pure numpy arrays")
-coef4 = basis.createFromArray(mass, data, time=3.3)
+# print("---- createFromArray usings pure numpy arrays")
+# coef4 = basis.createFromArray(mass, data, time=3.3)
 
 # Add the coefficient structure coefficient structure
 coefs = pyEXP.coefs.SphCoefs(True)
 coefs.add(coef1)
 # coefs.add(coef2)
-coefs.add(coef3)
-coefs.add(coef4)
+# coefs.add(coef3)
+# coefs.add(coef4)
 print("Times:", coefs.Times())
 
 exit(0)

From 278e1fee97a91f9099795d3971d62c2284afbc37 Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Sat, 12 Oct 2024 23:36:45 -0400
Subject: [PATCH 20/52] Restore original test script

---
 tests/Halo/createCoefs.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/tests/Halo/createCoefs.py b/tests/Halo/createCoefs.py
index 64e688fdb..b85b57555 100644
--- a/tests/Halo/createCoefs.py
+++ b/tests/Halo/createCoefs.py
@@ -47,8 +47,6 @@
 print("---- createFromArray usings lists")
 coef1 = basis.createFromArray(mass, [xpos, ypos, zpos], time=3.0)
 
-exit(0)
-
 mass  = np.array(mass)
 data  = np.array([xpos, ypos, zpos])
 

From 21045a6ab67b6ced084c338df5fb4804e85a478b Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Sun, 13 Oct 2024 10:34:38 -0400
Subject: [PATCH 21/52] Label creatCoefs test as 'long' to prevent CI failure

---
 tests/CMakeLists.txt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 560f9e8b5..9f315ce7c 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -100,7 +100,7 @@ if(ENABLE_NBODY)
               PYTHONPATH=${CMAKE_BINARY_DIR}/pyEXP:$ENV{PYTHONPATH}
               ${PYTHON_EXECUTABLE} createCoefs.py
               WORKING_DIRECTORY "${CMAKE_CURRENT_LIST_DIR}/Halo")
-     set_tests_properties(pyEXPCoefCreateTest PROPERTIES LABELS "quick")
+     set_tests_properties(pyEXPCoefCreateTest PROPERTIES LABELS "long")
   endif()
 
   # A separate test to remove the generated files if they all exist;

From 4c272053083475fd26d5a910805b8de449379ce3 Mon Sep 17 00:00:00 2001
From: Martin Weinberg <mdw@umass.edu>
Date: Sun, 13 Oct 2024 11:37:29 -0400
Subject: [PATCH 22/52] Update createCoefs.py

One more try just to make sure that this fails on all `createFromArrry` statements.
---
 tests/Halo/createCoefs.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/tests/Halo/createCoefs.py b/tests/Halo/createCoefs.py
index 773a9a0c6..f52fd313e 100644
--- a/tests/Halo/createCoefs.py
+++ b/tests/Halo/createCoefs.py
@@ -47,6 +47,8 @@
 print("---- createFromArray usings lists")
 coef1 = basis.createFromArray(mass, [xpos, ypos, zpos], time=3.0)
 
+exit(0) # TEST
+
 mass  = np.array(mass)
 data  = np.array([xpos, ypos, zpos])
 

From 6853072329bde47ce701b5b407702e8748d9ab9f Mon Sep 17 00:00:00 2001
From: Martin Weinberg <mdw@umass.edu>
Date: Sun, 13 Oct 2024 12:19:03 -0400
Subject: [PATCH 23/52] Update createCoefs.py

More bisection...
---
 tests/Halo/createCoefs.py | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/tests/Halo/createCoefs.py b/tests/Halo/createCoefs.py
index f52fd313e..ad6d9f97b 100644
--- a/tests/Halo/createCoefs.py
+++ b/tests/Halo/createCoefs.py
@@ -47,7 +47,11 @@
 print("---- createFromArray usings lists")
 coef1 = basis.createFromArray(mass, [xpos, ypos, zpos], time=3.0)
 
-exit(0) # TEST
+coefs = pyEXP.coefs.SphCoefs(True)
+coefs.add(coef1)
+
+print("Times:", coefs.Times())
+exit(0) # TEST END
 
 mass  = np.array(mass)
 data  = np.array([xpos, ypos, zpos])

From ca80b092d89945bb854bfb8894601bdb14056860 Mon Sep 17 00:00:00 2001
From: Martin Weinberg <mdw@umass.edu>
Date: Sun, 13 Oct 2024 12:54:51 -0400
Subject: [PATCH 24/52] Update createCoefs.py

Another bisection, renaming arrays.
---
 tests/Halo/createCoefs.py | 32 +++++++++++++++++---------------
 1 file changed, 17 insertions(+), 15 deletions(-)

diff --git a/tests/Halo/createCoefs.py b/tests/Halo/createCoefs.py
index ad6d9f97b..3889a828c 100644
--- a/tests/Halo/createCoefs.py
+++ b/tests/Halo/createCoefs.py
@@ -32,32 +32,34 @@
 
 # Call the basis to generate coefficients
 #
-mass = []
-xpos = []
-ypos = []
-zpos = []
+mass1 = []
+xpos1 = []
+ypos1 = []
+zpos1 = []
 
 print("---- creating list data")
 for i in range(0, 100):
-    mass.append(0.001)
-    xpos.append(random.random()*2.0 - 1.0)
-    ypos.append(random.random()*2.0 - 1.0)
-    zpos.append(random.random()*2.0 - 1.0)
+    mass1.append(0.001)
+    xpos1.append(random.random()*2.0 - 1.0)
+    ypos1.append(random.random()*2.0 - 1.0)
+    zpos1.append(random.random()*2.0 - 1.0)
 
 print("---- createFromArray usings lists")
-coef1 = basis.createFromArray(mass, [xpos, ypos, zpos], time=3.0)
+coef1 = basis.createFromArray(mass1, [xpos1, ypos1, zpos1], time=3.0)
 
 coefs = pyEXP.coefs.SphCoefs(True)
 coefs.add(coef1)
 
-print("Times:", coefs.Times())
-exit(0) # TEST END
-
-mass  = np.array(mass)
-data  = np.array([xpos, ypos, zpos])
+mass2 = np.array(mass1)
+data2 = np.array([xpos1, ypos1, zpos1])
 
 print("---- createFromArray usings numpy arrays converted from lists")
-coef2 = basis.createFromArray(mass, data, time=3.1)
+coef2 = basis.createFromArray(mass2, data2, time=3.1)
+
+coefs.add(coef2)
+
+print("Times:", coefs.Times())
+exit(0) # TEST END
 
 mass = np.ones(1000) * 1.0e6
 xpos = np.random.normal(0.0, 1.0, 1000)

From 7035069931ee16a89aa76c4d624727ba84bc2704 Mon Sep 17 00:00:00 2001
From: Martin Weinberg <mdw@umass.edu>
Date: Sun, 13 Oct 2024 13:12:25 -0400
Subject: [PATCH 25/52] Update createCoefs.py

More updates and bisection...
---
 tests/Halo/createCoefs.py | 31 ++++++++++++++++---------------
 1 file changed, 16 insertions(+), 15 deletions(-)

diff --git a/tests/Halo/createCoefs.py b/tests/Halo/createCoefs.py
index 3889a828c..d4c702ceb 100644
--- a/tests/Halo/createCoefs.py
+++ b/tests/Halo/createCoefs.py
@@ -58,29 +58,30 @@
 
 coefs.add(coef2)
 
-print("Times:", coefs.Times())
-exit(0) # TEST END
-
-mass = np.ones(1000) * 1.0e6
-xpos = np.random.normal(0.0, 1.0, 1000)
-ypos = np.random.normal(0.0, 1.0, 1000)
-zpos = np.random.normal(0.0, 1.0, 1000)
+mass3 = np.ones(1000) * 1.0e6
+xpos3 = np.random.normal(0.0, 1.0, 1000)
+ypos3 = np.random.normal(0.0, 1.0, 1000)
+zpos3 = np.random.normal(0.0, 1.0, 1000)
 
 print("---- createFromArray using a list of numpy arrays")
-coef3 = basis.createFromArray(mass, [xpos, ypos, zpos], time=3.2)
+coef3 = basis.createFromArray(mass3, [xpos3, ypos3, zpos3], time=3.2)
+
+coefs.add(coef3)
 
-data  = np.array([xpos, ypos, zpos])
+data4  = np.array([xpos3, ypos3, zpos3])
 
 print("---- createFromArray usings pure numpy arrays")
-coef4 = basis.createFromArray(mass, data, time=3.3)
+coef4 = basis.createFromArray(mass3, data4, time=3.3)
+
+coefs.add(coef4)
 
 # Add the coefficient structure coefficient structure
 #
-coefs = pyEXP.coefs.SphCoefs(True)
-coefs.add(coef1)
-coefs.add(coef2)
-coefs.add(coef3)
-coefs.add(coef4)
+# coefs = pyEXP.coefs.SphCoefs(True)
+# coefs.add(coef1)
+# coefs.add(coef2)
+# coefs.add(coef3)
+# coefs.add(coef4)
 print("Times:", coefs.Times())
 
 exit(0)

From f2932d20e5dd7e6c08c104f4575920f138265df4 Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Sun, 13 Oct 2024 13:33:00 -0400
Subject: [PATCH 26/52] Spruce up for final commit (I hope)

---
 tests/Halo/createCoefs.py | 18 ++++++++++--------
 1 file changed, 10 insertions(+), 8 deletions(-)

diff --git a/tests/Halo/createCoefs.py b/tests/Halo/createCoefs.py
index d4c702ceb..786c2c98e 100644
--- a/tests/Halo/createCoefs.py
+++ b/tests/Halo/createCoefs.py
@@ -27,8 +27,11 @@
 
 print("---- created basis")
 
-# Make some fake body data
+# Create a coefficient structure
 #
+coefs = pyEXP.coefs.SphCoefs(True)
+
+print("---- created coefficients")
 
 # Call the basis to generate coefficients
 #
@@ -47,9 +50,10 @@
 print("---- createFromArray usings lists")
 coef1 = basis.createFromArray(mass1, [xpos1, ypos1, zpos1], time=3.0)
 
-coefs = pyEXP.coefs.SphCoefs(True)
 coefs.add(coef1)
 
+print("---- creating array data from list data")
+
 mass2 = np.array(mass1)
 data2 = np.array([xpos1, ypos1, zpos1])
 
@@ -58,6 +62,8 @@
 
 coefs.add(coef2)
 
+print("---- creating array data from list of numpy arrays")
+
 mass3 = np.ones(1000) * 1.0e6
 xpos3 = np.random.normal(0.0, 1.0, 1000)
 ypos3 = np.random.normal(0.0, 1.0, 1000)
@@ -69,19 +75,15 @@
 coefs.add(coef3)
 
 data4  = np.array([xpos3, ypos3, zpos3])
+print("---- position data shape is:", data4.shape)
 
 print("---- createFromArray usings pure numpy arrays")
 coef4 = basis.createFromArray(mass3, data4, time=3.3)
 
 coefs.add(coef4)
 
-# Add the coefficient structure coefficient structure
+# Check the coefficient structure for the 4 added times
 #
-# coefs = pyEXP.coefs.SphCoefs(True)
-# coefs.add(coef1)
-# coefs.add(coef2)
-# coefs.add(coef3)
-# coefs.add(coef4)
 print("Times:", coefs.Times())
 
 exit(0)

From 8a3d555e23c3d9496cf36ae5ff982478d605c2c0 Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Sun, 13 Oct 2024 15:03:20 -0400
Subject: [PATCH 27/52] Change to pass arrays by copy rather than reference

---
 expui/BasisFactory.H      |  4 ++--
 expui/BasisFactory.cc     |  2 +-
 expui/BiorthBasis.H       |  4 ++--
 expui/BiorthBasis.cc      |  4 ++--
 expui/Coefficients.H      | 30 +++++++++++++-------------
 expui/Coefficients.cc     | 26 +++++++++++------------
 expui/FieldBasis.H        |  2 +-
 expui/FieldBasis.cc       |  2 +-
 pyEXP/BasisWrappers.cc    |  6 +++---
 pyEXP/CoefWrappers.cc     | 12 +++++------
 tests/Halo/createCoefs.py | 44 +++++++++++++++++++++------------------
 11 files changed, 70 insertions(+), 66 deletions(-)

diff --git a/expui/BasisFactory.H b/expui/BasisFactory.H
index 089887934..d27e8e686 100644
--- a/expui/BasisFactory.H
+++ b/expui/BasisFactory.H
@@ -205,13 +205,13 @@ namespace BasisClasses
 
     //! Accumulate coefficient contributions from arrays
     virtual void
-    addFromArray(Eigen::VectorXd& m, RowMatrixXd& p, bool roundrobin,
+    addFromArray(Eigen::VectorXd m, RowMatrixXd p, bool roundrobin,
 		 bool posvelrows) = 0;
 
     //! Generate coeffients from an array and optional center location
     //! for the expansion
     CoefClasses::CoefStrPtr createFromArray
-    (Eigen::VectorXd& m, RowMatrixXd& p, double time=0.0,
+    (Eigen::VectorXd m, RowMatrixXd p, double time=0.0,
      std::vector<double> center={0.0, 0.0, 0.0},
      bool roundrobin=true, bool posvelrows=false);
 
diff --git a/expui/BasisFactory.cc b/expui/BasisFactory.cc
index 5b574a8e2..6569f0ee4 100644
--- a/expui/BasisFactory.cc
+++ b/expui/BasisFactory.cc
@@ -272,7 +272,7 @@ namespace BasisClasses
   // Generate coefficients from a phase-space table
   //
   CoefClasses::CoefStrPtr Basis::createFromArray
-  (Eigen::VectorXd& m, RowMatrixXd& p, double time, std::vector<double> ctr,
+  (Eigen::VectorXd m, RowMatrixXd p, double time, std::vector<double> ctr,
    bool roundrobin, bool posvelrows)
   {
     initFromArray(ctr);
diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H
index a4af27c1a..518025c73 100644
--- a/expui/BiorthBasis.H
+++ b/expui/BiorthBasis.H
@@ -97,7 +97,7 @@ namespace BasisClasses
     //! Generate coeffients from an array and optional center location
     //! for the expansion
     CoefClasses::CoefStrPtr createFromArray
-    (Eigen::VectorXd& m, RowMatrixXd& p, double time=0.0,
+    (Eigen::VectorXd m, RowMatrixXd p, double time=0.0,
      std::vector<double> center={0.0, 0.0, 0.0},
      bool roundrobin=true, bool posvelrows=false);
     
@@ -113,7 +113,7 @@ namespace BasisClasses
     //! Initialize accumulating coefficients from arrays.  This is
     //! called once to initialize the accumulation.
     void addFromArray
-    (Eigen::VectorXd& m, RowMatrixXd& p,
+    (Eigen::VectorXd m, RowMatrixXd p,
      bool roundrobin=true, bool posvelrows=false);
     
     //! Create and the coefficients from the array accumulation with the
diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc
index 9625e4f2d..553d7106b 100644
--- a/expui/BiorthBasis.cc
+++ b/expui/BiorthBasis.cc
@@ -2847,7 +2847,7 @@ namespace BasisClasses
   }
 
   // Accumulate coefficient contributions from arrays
-  void BiorthBasis::addFromArray(Eigen::VectorXd& m, RowMatrixXd& p,
+  void BiorthBasis::addFromArray(Eigen::VectorXd m, RowMatrixXd p,
 				 bool RoundRobin, bool PosVelRows)
   {
     // Sanity check: is coefficient instance created?  This is not
@@ -2957,7 +2957,7 @@ namespace BasisClasses
   // Generate coefficients from a phase-space table
   //
   CoefClasses::CoefStrPtr BiorthBasis::createFromArray
-  (Eigen::VectorXd& m, RowMatrixXd& p, double time, std::vector<double> ctr,
+  (Eigen::VectorXd m, RowMatrixXd p, double time, std::vector<double> ctr,
    bool RoundRobin, bool PosVelRows)
   {
     initFromArray(ctr);
diff --git a/expui/Coefficients.H b/expui/Coefficients.H
index 7718c7b7b..46664237f 100644
--- a/expui/Coefficients.H
+++ b/expui/Coefficients.H
@@ -111,7 +111,7 @@ namespace CoefClasses
     virtual Eigen::VectorXcd& getData(double time) = 0;
 
     //! Set coefficient store at given time to the provided matrix
-    virtual void setData(double time, Eigen::VectorXcd& data) = 0;
+    virtual void setData(double time, Eigen::VectorXcd data) = 0;
 
     //! Interpolate coefficient matrix at given time
     std::tuple<Eigen::VectorXcd&, bool> interpolate(double time);
@@ -273,10 +273,10 @@ namespace CoefClasses
     Eigen::MatrixXcd& getMatrix(double time);
 
     //! Set coefficient matrix at given time
-    virtual void setData(double time, Eigen::VectorXcd& arr);
+    virtual void setData(double time, Eigen::VectorXcd arr);
 
     //! Natural data setter for pyEXP
-    void setMatrix(double time, Eigen::MatrixXcd& mat);
+    void setMatrix(double time, Eigen::MatrixXcd mat);
 
     //! Interpolate coefficient matrix at given time
     std::tuple<Eigen::MatrixXcd&, bool> interpolate(double time)
@@ -411,10 +411,10 @@ namespace CoefClasses
     Eigen::MatrixXcd& getMatrix(double time);
     
     //! Set coefficient matrix at given time
-    virtual void setData(double time, Eigen::VectorXcd& arr);
+    virtual void setData(double time, Eigen::VectorXcd arr);
     
     //! For pyEXP
-    void setMatrix(double time, Eigen::MatrixXcd& arr);
+    void setMatrix(double time, Eigen::MatrixXcd arr);
     
     //! Interpolate coefficient matrix at given time
     std::tuple<Eigen::MatrixXcd&, bool> interpolate(double time)
@@ -542,10 +542,10 @@ namespace CoefClasses
     virtual Eigen3d& getTensor(double time);
     
     //! Set coefficient matrix at given time
-    virtual void setData(double time, Eigen::VectorXcd& dat);
+    virtual void setData(double time, Eigen::VectorXcd dat);
 
     //! Native version
-    virtual void setTensor(double time, const Eigen3d& dat);
+    virtual void setTensor(double time, const Eigen3d dat);
     
     //! Interpolate coefficient tensor at given time
     std::tuple<Eigen::Tensor<std::complex<double>, 3>&, bool>
@@ -673,10 +673,10 @@ namespace CoefClasses
     virtual Eigen3d& getTensor(double time);
     
     //! Set coefficient matrix at given time
-    virtual void setData(double time, Eigen::VectorXcd& dat);
+    virtual void setData(double time, Eigen::VectorXcd dat);
 
     //! Native version
-    virtual void setTensor(double time, const Eigen3d& dat);
+    virtual void setTensor(double time, const Eigen3d dat);
     
     //! Interpolate coefficient tensor at given time
     std::tuple<Eigen::Tensor<std::complex<double>, 3>&, bool>
@@ -810,9 +810,9 @@ namespace CoefClasses
     /** Set coefficient matrix at given time.  This is for pybind11,
 	since the operator() will not allow lvalue assignment, it
 	seems. */
-    virtual void setData(double time, Eigen::VectorXcd& arr);
+    virtual void setData(double time, Eigen::VectorXcd arr);
 
-    void setArray(double time, Eigen::VectorXcd& arr)
+    void setArray(double time, Eigen::VectorXcd arr)
     { setData(time, arr); }
 
     //! Get coefficient structure at a given time
@@ -915,10 +915,10 @@ namespace CoefClasses
     SphFldStruct::dataType getMatrix(double time);
 
     //! Set coefficient matrix at given time
-    virtual void setData(double time, Eigen::VectorXcd& arr);
+    virtual void setData(double time, Eigen::VectorXcd arr);
 
     //! Natural data setter for pyEXP
-    void setMatrix(double time, SphFldStruct::dataType& mat);
+    void setMatrix(double time, SphFldStruct::dataType mat);
 
     //! Interpolate coefficient matrix at given time
     std::tuple<SphFldStruct::coefType&, bool> interpolate(double time)
@@ -1042,10 +1042,10 @@ namespace CoefClasses
     CylFldStruct::dataType getMatrix(double time);
 
     //! Set coefficient matrix at given time
-    virtual void setData(double time, Eigen::VectorXcd& arr);
+    virtual void setData(double time, Eigen::VectorXcd arr);
 
     //! Natural data setter for pyEXP
-    void setMatrix(double time, CylFldStruct::dataType& mat);
+    void setMatrix(double time, CylFldStruct::dataType mat);
 
     //! Interpolate coefficient matrix at given time
     std::tuple<CylFldStruct::coefType&, bool> interpolate(double time)
diff --git a/expui/Coefficients.cc b/expui/Coefficients.cc
index fe55b0efe..83e352a94 100644
--- a/expui/Coefficients.cc
+++ b/expui/Coefficients.cc
@@ -469,7 +469,7 @@ namespace CoefClasses
     return mat;
   }
   
-  void SphCoefs::setData(double time, Eigen::VectorXcd& dat)
+  void SphCoefs::setData(double time, Eigen::VectorXcd dat)
   {
     auto it = coefs.find(roundTime(time));
 
@@ -484,7 +484,7 @@ namespace CoefClasses
     }
   }
   
-  void SphCoefs::setMatrix(double time, Eigen::MatrixXcd& dat)
+  void SphCoefs::setMatrix(double time, Eigen::MatrixXcd dat)
   {
     auto it = coefs.find(roundTime(time));
 
@@ -891,7 +891,7 @@ namespace CoefClasses
     return mat;
   }
   
-  void CylCoefs::setData(double time, Eigen::VectorXcd& dat)
+  void CylCoefs::setData(double time, Eigen::VectorXcd dat)
   {
     auto it = coefs.find(roundTime(time));
 
@@ -906,7 +906,7 @@ namespace CoefClasses
     }
   }
 
-  void CylCoefs::setMatrix(double time, Eigen::MatrixXcd& dat)
+  void CylCoefs::setMatrix(double time, Eigen::MatrixXcd dat)
   {
     auto it = coefs.find(roundTime(time));
 
@@ -1247,7 +1247,7 @@ namespace CoefClasses
     return arr;
   }
 
-  void SlabCoefs::setData(double time, Eigen::VectorXcd& dat)
+  void SlabCoefs::setData(double time, Eigen::VectorXcd dat)
   {
     auto it = coefs.find(roundTime(time));
 
@@ -1262,7 +1262,7 @@ namespace CoefClasses
     }
   }
 
-  void SlabCoefs::setTensor(double time, const Eigen3d& dat)
+  void SlabCoefs::setTensor(double time, const Eigen3d dat)
   {
     auto it = coefs.find(roundTime(time));
 
@@ -1599,7 +1599,7 @@ namespace CoefClasses
     return arr;
   }
 
-  void CubeCoefs::setData(double time, Eigen::VectorXcd& dat)
+  void CubeCoefs::setData(double time, Eigen::VectorXcd dat)
   {
     auto it = coefs.find(roundTime(time));
 
@@ -1614,7 +1614,7 @@ namespace CoefClasses
     }
   }
 
-  void CubeCoefs::setTensor(double time, const Eigen3d& dat)
+  void CubeCoefs::setTensor(double time, const Eigen3d dat)
   {
     auto it = coefs.find(roundTime(time));
 
@@ -1968,7 +1968,7 @@ namespace CoefClasses
     return arr;
   }
   
-  void TableData::setData(double time, Eigen::VectorXcd& dat)
+  void TableData::setData(double time, Eigen::VectorXcd dat)
   {
     auto it = coefs.find(roundTime(time));
 
@@ -2458,7 +2458,7 @@ namespace CoefClasses
     return *mat;
   }
   
-  void SphFldCoefs::setData(double time, Eigen::VectorXcd& dat)
+  void SphFldCoefs::setData(double time, Eigen::VectorXcd dat)
   {
     auto it = coefs.find(roundTime(time));
 
@@ -2473,7 +2473,7 @@ namespace CoefClasses
     }
   }
   
-  void SphFldCoefs::setMatrix(double time, SphFldStruct::dataType& dat)
+  void SphFldCoefs::setMatrix(double time, SphFldStruct::dataType dat)
   {
     auto it = coefs.find(roundTime(time));
 
@@ -2897,7 +2897,7 @@ namespace CoefClasses
     return *mat;
   }
   
-  void CylFldCoefs::setData(double time, Eigen::VectorXcd& dat)
+  void CylFldCoefs::setData(double time, Eigen::VectorXcd dat)
   {
     auto it = coefs.find(roundTime(time));
 
@@ -2912,7 +2912,7 @@ namespace CoefClasses
     }
   }
   
-  void CylFldCoefs::setMatrix(double time, CylFldStruct::dataType& dat)
+  void CylFldCoefs::setMatrix(double time, CylFldStruct::dataType dat)
   {
     auto it = coefs.find(roundTime(time));
 
diff --git a/expui/FieldBasis.H b/expui/FieldBasis.H
index 25ef5dbe3..b04653e3d 100644
--- a/expui/FieldBasis.H
+++ b/expui/FieldBasis.H
@@ -146,7 +146,7 @@ namespace BasisClasses
     //! Initialize accumulating coefficients from arrays.  This is
     //! called once to initialize the accumulation.
     void addFromArray
-    (Eigen::VectorXd& m, RowMatrixXd& p, bool roundrobin=true, bool posvelrows=false);
+    (Eigen::VectorXd m, RowMatrixXd p, bool roundrobin=true, bool posvelrows=false);
     
     //! Accumulate new coefficients
     virtual void accumulate(double mass,
diff --git a/expui/FieldBasis.cc b/expui/FieldBasis.cc
index 7c85e0aee..2d25a1866 100644
--- a/expui/FieldBasis.cc
+++ b/expui/FieldBasis.cc
@@ -619,7 +619,7 @@ namespace BasisClasses
   }
 
   // Accumulate coefficient contributions from arrays
-  void FieldBasis::addFromArray(Eigen::VectorXd& m, RowMatrixXd& p,
+  void FieldBasis::addFromArray(Eigen::VectorXd m, RowMatrixXd p,
 				bool roundrobin, bool posvelrows)
   {
     // Sanity check: is coefficient instance created?  This is not
diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc
index e0b72ab25..07e2369f0 100644
--- a/pyEXP/BasisWrappers.cc
+++ b/pyEXP/BasisWrappers.cc
@@ -286,7 +286,7 @@ void BasisFactoryClasses(py::module &m)
     }
 
     virtual void
-    addFromArray(Eigen::VectorXd& m, RowMatrixXd& p, bool roundrobin, bool posvelrows) override {
+    addFromArray(Eigen::VectorXd m, RowMatrixXd p, bool roundrobin, bool posvelrows) override {
       PYBIND11_OVERRIDE_PURE(void, Basis, addFromArray, m, p, roundrobin, posvelrows);
     }
   };
@@ -973,7 +973,7 @@ void BasisFactoryClasses(py::module &m)
          )",
 	 py::arg("center") = std::vector<double>(3, 0.0))
     .def("addFromArray",
-	 [](BasisClasses::BiorthBasis& A, Eigen::VectorXd& mass, RowMatrixXd& pos)
+	 [](BasisClasses::BiorthBasis& A, Eigen::VectorXd mass, RowMatrixXd pos)
 	 {
 	   return A.addFromArray(mass, pos);
 	 },
@@ -1804,7 +1804,7 @@ void BasisFactoryClasses(py::module &m)
          )",
 	 py::arg("center") = std::vector<double>(3, 0.0))
     .def("addFromArray",
-	 [](BasisClasses::FieldBasis& A, Eigen::VectorXd& mass, RowMatrixXd& ps)
+	 [](BasisClasses::FieldBasis& A, Eigen::VectorXd mass, RowMatrixXd ps)
 	 {
 	   return A.addFromArray(mass, ps);
 	 },
diff --git a/pyEXP/CoefWrappers.cc b/pyEXP/CoefWrappers.cc
index aa58603b9..da2ce7227 100644
--- a/pyEXP/CoefWrappers.cc
+++ b/pyEXP/CoefWrappers.cc
@@ -155,7 +155,7 @@ void CoefficientClasses(py::module &m) {
       PYBIND11_OVERRIDE_PURE(Eigen::VectorXcd&, Coefs, getData, time);
     }
 
-    void setData(double time, Eigen::VectorXcd& array) override {
+    void setData(double time, Eigen::VectorXcd array) override {
       PYBIND11_OVERRIDE_PURE(void, Coefs, setData, time, array);
     }
 
@@ -231,7 +231,7 @@ void CoefficientClasses(py::module &m) {
       PYBIND11_OVERRIDE(Eigen::VectorXcd&, SphCoefs, getData, time);
     }
 
-    void setData(double time, Eigen::VectorXcd& array) override {
+    void setData(double time, Eigen::VectorXcd array) override {
       PYBIND11_OVERRIDE(void, SphCoefs, setData, time, array);
     }
 
@@ -310,7 +310,7 @@ void CoefficientClasses(py::module &m) {
       PYBIND11_OVERRIDE(Eigen::VectorXcd&, CylCoefs, getData, time);
     }
 
-    void setData(double time, Eigen::VectorXcd& array) override {
+    void setData(double time, Eigen::VectorXcd array) override {
       PYBIND11_OVERRIDE(void, CylCoefs, setData, time, array);
     }
 
@@ -388,7 +388,7 @@ void CoefficientClasses(py::module &m) {
       PYBIND11_OVERRIDE(Eigen::VectorXcd&, SlabCoefs, getData, time);
     }
 
-    void setData(double time, Eigen::VectorXcd& array) override {
+    void setData(double time, Eigen::VectorXcd array) override {
       PYBIND11_OVERRIDE(void, SlabCoefs, setData, time, array);
     }
 
@@ -467,7 +467,7 @@ void CoefficientClasses(py::module &m) {
       PYBIND11_OVERRIDE(Eigen::VectorXcd&, CubeCoefs, getData, time);
     }
 
-    void setData(double time, Eigen::VectorXcd& array) override {
+    void setData(double time, Eigen::VectorXcd array) override {
       PYBIND11_OVERRIDE(void, CubeCoefs, setData, time, array);
     }
 
@@ -546,7 +546,7 @@ void CoefficientClasses(py::module &m) {
       PYBIND11_OVERRIDE(Eigen::VectorXcd&, TableData, getData, time);
     }
 
-    void setData(double time, Eigen::VectorXcd& array) override {
+    void setData(double time, Eigen::VectorXcd array) override {
       PYBIND11_OVERRIDE(void, TableData, setData, time, array);
     }
 
diff --git a/tests/Halo/createCoefs.py b/tests/Halo/createCoefs.py
index 786c2c98e..3a761c82e 100644
--- a/tests/Halo/createCoefs.py
+++ b/tests/Halo/createCoefs.py
@@ -35,50 +35,54 @@
 
 # Call the basis to generate coefficients
 #
-mass1 = []
-xpos1 = []
-ypos1 = []
-zpos1 = []
+mass = []
+xpos = []
+ypos = []
+zpos = []
 
 print("---- creating list data")
 for i in range(0, 100):
-    mass1.append(0.001)
-    xpos1.append(random.random()*2.0 - 1.0)
-    ypos1.append(random.random()*2.0 - 1.0)
-    zpos1.append(random.random()*2.0 - 1.0)
+    mass.append(0.01)
+    xpos.append(random.random()*2.0 - 1.0)
+    ypos.append(random.random()*2.0 - 1.0)
+    zpos.append(random.random()*2.0 - 1.0)
 
 print("---- createFromArray usings lists")
-coef1 = basis.createFromArray(mass1, [xpos1, ypos1, zpos1], time=3.0)
+coef1 = basis.createFromArray(mass, [xpos, ypos, zpos], time=3.0)
 
 coefs.add(coef1)
 
 print("---- creating array data from list data")
 
-mass2 = np.array(mass1)
-data2 = np.array([xpos1, ypos1, zpos1])
+# Note: this overwrites mass and data variables.  But data is now
+# passed by rference so pybind11 should have copies and not break the
+# references though garbage collection on the Python side
+#
+mass = np.array(mass)
+data = np.array([xpos, ypos, zpos])
 
 print("---- createFromArray usings numpy arrays converted from lists")
-coef2 = basis.createFromArray(mass2, data2, time=3.1)
+coef2 = basis.createFromArray(mass, data, time=3.1)
 
 coefs.add(coef2)
 
 print("---- creating array data from list of numpy arrays")
 
-mass3 = np.ones(1000) * 1.0e6
-xpos3 = np.random.normal(0.0, 1.0, 1000)
-ypos3 = np.random.normal(0.0, 1.0, 1000)
-zpos3 = np.random.normal(0.0, 1.0, 1000)
+mass = np.ones(100) * 1.0e-02
+xpos = np.random.normal(0.0, 1.0, 100)
+ypos = np.random.normal(0.0, 1.0, 100)
+zpos = np.random.normal(0.0, 1.0, 100)
 
 print("---- createFromArray using a list of numpy arrays")
-coef3 = basis.createFromArray(mass3, [xpos3, ypos3, zpos3], time=3.2)
+coef3 = basis.createFromArray(mass, [xpos, ypos, zpos], time=3.2)
 
 coefs.add(coef3)
 
-data4  = np.array([xpos3, ypos3, zpos3])
-print("---- position data shape is:", data4.shape)
+data  = np.array([xpos, ypos, zpos])
+print("---- position data shape is:", data.shape)
 
 print("---- createFromArray usings pure numpy arrays")
-coef4 = basis.createFromArray(mass3, data4, time=3.3)
+coef4 = basis.createFromArray(mass, data, time=3.3)
 
 coefs.add(coef4)
 

From fd69366b4d79235d08c4adaf0a1d3773646ddedf Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Sun, 13 Oct 2024 15:29:16 -0400
Subject: [PATCH 28/52] Restore createArray test to 'quick'

---
 tests/CMakeLists.txt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 9f315ce7c..560f9e8b5 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -100,7 +100,7 @@ if(ENABLE_NBODY)
               PYTHONPATH=${CMAKE_BINARY_DIR}/pyEXP:$ENV{PYTHONPATH}
               ${PYTHON_EXECUTABLE} createCoefs.py
               WORKING_DIRECTORY "${CMAKE_CURRENT_LIST_DIR}/Halo")
-     set_tests_properties(pyEXPCoefCreateTest PROPERTIES LABELS "long")
+     set_tests_properties(pyEXPCoefCreateTest PROPERTIES LABELS "quick")
   endif()
 
   # A separate test to remove the generated files if they all exist;

From 34bedcec339d32ac73458369b4e84374456cb4e6 Mon Sep 17 00:00:00 2001
From: Martin Weinberg <mdw@umass.edu>
Date: Sun, 13 Oct 2024 16:32:26 -0400
Subject: [PATCH 29/52] Update createCoefs.py

Bisect again (?!)
---
 tests/Halo/createCoefs.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/tests/Halo/createCoefs.py b/tests/Halo/createCoefs.py
index 3a761c82e..48a2eb2b2 100644
--- a/tests/Halo/createCoefs.py
+++ b/tests/Halo/createCoefs.py
@@ -52,6 +52,10 @@
 
 coefs.add(coef1)
 
+print("Times:", coefs.Times())
+
+exit(0)
+
 print("---- creating array data from list data")
 
 # Note: this overwrites mass and data variables.  But data is now

From fa3740ce36c8be2e89a2cf6984d6ed1a37d99a6c Mon Sep 17 00:00:00 2001
From: Martin Weinberg <mdw@umass.edu>
Date: Sun, 13 Oct 2024 16:49:34 -0400
Subject: [PATCH 30/52] Update createCoefs.py

add 'sleep()' to try to get some diagnostic output from the runner
---
 tests/Halo/createCoefs.py | 9 +++++++--
 1 file changed, 7 insertions(+), 2 deletions(-)

diff --git a/tests/Halo/createCoefs.py b/tests/Halo/createCoefs.py
index 48a2eb2b2..851817b64 100644
--- a/tests/Halo/createCoefs.py
+++ b/tests/Halo/createCoefs.py
@@ -1,4 +1,4 @@
-import os
+import os, time
 import pyEXP
 import random
 import numpy as np
@@ -27,12 +27,16 @@
 
 print("---- created basis")
 
+time.sleep(1)
+
 # Create a coefficient structure
 #
 coefs = pyEXP.coefs.SphCoefs(True)
 
 print("---- created coefficients")
 
+time.sleep(1)
+
 # Call the basis to generate coefficients
 #
 mass = []
@@ -48,8 +52,9 @@
     zpos.append(random.random()*2.0 - 1.0)
 
 print("---- createFromArray usings lists")
+time.sleep(1)
 coef1 = basis.createFromArray(mass, [xpos, ypos, zpos], time=3.0)
-
+time.sleep(1)
 coefs.add(coef1)
 
 print("Times:", coefs.Times())

From 996054d594e494a5682d09d69443e93e0eaf4d8a Mon Sep 17 00:00:00 2001
From: Martin Weinberg <mdw@umass.edu>
Date: Sun, 13 Oct 2024 17:08:35 -0400
Subject: [PATCH 31/52] Update createCoefs.py

---
 tests/Halo/createCoefs.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/tests/Halo/createCoefs.py b/tests/Halo/createCoefs.py
index 851817b64..88d356c54 100644
--- a/tests/Halo/createCoefs.py
+++ b/tests/Halo/createCoefs.py
@@ -37,6 +37,8 @@
 
 time.sleep(1)
 
+exit(0)
+
 # Call the basis to generate coefficients
 #
 mass = []

From 0343410742b159661ad3e9572531070bb510fd56 Mon Sep 17 00:00:00 2001
From: Martin Weinberg <mdw@umass.edu>
Date: Sun, 13 Oct 2024 18:29:41 -0400
Subject: [PATCH 32/52] Update createCoefs.py

---
 tests/Halo/createCoefs.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/tests/Halo/createCoefs.py b/tests/Halo/createCoefs.py
index 88d356c54..851817b64 100644
--- a/tests/Halo/createCoefs.py
+++ b/tests/Halo/createCoefs.py
@@ -37,8 +37,6 @@
 
 time.sleep(1)
 
-exit(0)
-
 # Call the basis to generate coefficients
 #
 mass = []

From 0fda70eae42a49af5bf17b7f8bb28e61676ceda8 Mon Sep 17 00:00:00 2001
From: Martin Weinberg <mdw@umass.edu>
Date: Sun, 13 Oct 2024 23:31:16 -0400
Subject: [PATCH 33/52] Update CMakeLists.txt

Relabel `pyEXPCoefCreateTest` to `long` (again)
---
 tests/CMakeLists.txt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 560f9e8b5..9f315ce7c 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -100,7 +100,7 @@ if(ENABLE_NBODY)
               PYTHONPATH=${CMAKE_BINARY_DIR}/pyEXP:$ENV{PYTHONPATH}
               ${PYTHON_EXECUTABLE} createCoefs.py
               WORKING_DIRECTORY "${CMAKE_CURRENT_LIST_DIR}/Halo")
-     set_tests_properties(pyEXPCoefCreateTest PROPERTIES LABELS "quick")
+     set_tests_properties(pyEXPCoefCreateTest PROPERTIES LABELS "long")
   endif()
 
   # A separate test to remove the generated files if they all exist;

From 42be626843e96b1823a2a61716fbfc70bef145b0 Mon Sep 17 00:00:00 2001
From: michael-petersen <petersen.michael.s@gmail.com>
Date: Mon, 14 Oct 2024 13:25:02 +0100
Subject: [PATCH 34/52] Add extra checks to CI for runner diagnostic
 information

---
 .github/workflows/build.yml | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
index d6b08ab84..85f50c5bb 100644
--- a/.github/workflows/build.yml
+++ b/.github/workflows/build.yml
@@ -12,14 +12,14 @@ jobs:
   exp:
     strategy:
       matrix:
-        os: [ubuntu-latest]
-        cc: [gcc]
+        os: [ubuntu-latest, ubuntu-22.04]
+        cc: [gcc, mpicc]
 
     name: "Test pyEXP Build"
     runs-on: ${{ matrix.os }}
     steps:
       - name: Checkout
-        uses: actions/checkout@v3
+        uses: actions/checkout@v4
 
       - name: Install core dependencies - ubuntu
         if: runner.os == 'Linux'
@@ -28,6 +28,7 @@ jobs:
           sudo apt-get install -y build-essential libeigen3-dev libfftw3-dev libhdf5-dev libopenmpi-dev
           sudo pip install numpy
 
+
       - name: Setup submodule and build
         run: |
           git submodule update --init --recursive

From f9c0bacea875d19ce27ef1de0675eb8ce8efec21 Mon Sep 17 00:00:00 2001
From: michael-petersen <petersen.michael.s@gmail.com>
Date: Mon, 14 Oct 2024 13:39:05 +0100
Subject: [PATCH 35/52] restore coef check to quick to induce failures

---
 tests/CMakeLists.txt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 9f315ce7c..560f9e8b5 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -100,7 +100,7 @@ if(ENABLE_NBODY)
               PYTHONPATH=${CMAKE_BINARY_DIR}/pyEXP:$ENV{PYTHONPATH}
               ${PYTHON_EXECUTABLE} createCoefs.py
               WORKING_DIRECTORY "${CMAKE_CURRENT_LIST_DIR}/Halo")
-     set_tests_properties(pyEXPCoefCreateTest PROPERTIES LABELS "long")
+     set_tests_properties(pyEXPCoefCreateTest PROPERTIES LABELS "quick")
   endif()
 
   # A separate test to remove the generated files if they all exist;

From bfbfbc848bccc79d00ca62b4dfdecd541379a9dd Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Mon, 14 Oct 2024 10:11:18 -0400
Subject: [PATCH 36/52] Rollback to reference version of setMatrix() and
 createFromArray(); remove debugging statements from createCoefs.py

---
 expui/BasisFactory.H      |  4 ++--
 expui/BasisFactory.cc     |  2 +-
 expui/BiorthBasis.H       |  4 ++--
 expui/BiorthBasis.cc      |  4 ++--
 expui/Coefficients.H      | 30 +++++++++++++++---------------
 expui/Coefficients.cc     | 26 +++++++++++++-------------
 expui/FieldBasis.H        |  2 +-
 expui/FieldBasis.cc       |  2 +-
 pyEXP/BasisWrappers.cc    |  6 +++---
 pyEXP/CoefWrappers.cc     | 12 ++++++------
 tests/Halo/createCoefs.py | 12 +++---------
 11 files changed, 49 insertions(+), 55 deletions(-)

diff --git a/expui/BasisFactory.H b/expui/BasisFactory.H
index d27e8e686..089887934 100644
--- a/expui/BasisFactory.H
+++ b/expui/BasisFactory.H
@@ -205,13 +205,13 @@ namespace BasisClasses
 
     //! Accumulate coefficient contributions from arrays
     virtual void
-    addFromArray(Eigen::VectorXd m, RowMatrixXd p, bool roundrobin,
+    addFromArray(Eigen::VectorXd& m, RowMatrixXd& p, bool roundrobin,
 		 bool posvelrows) = 0;
 
     //! Generate coeffients from an array and optional center location
     //! for the expansion
     CoefClasses::CoefStrPtr createFromArray
-    (Eigen::VectorXd m, RowMatrixXd p, double time=0.0,
+    (Eigen::VectorXd& m, RowMatrixXd& p, double time=0.0,
      std::vector<double> center={0.0, 0.0, 0.0},
      bool roundrobin=true, bool posvelrows=false);
 
diff --git a/expui/BasisFactory.cc b/expui/BasisFactory.cc
index 6569f0ee4..5b574a8e2 100644
--- a/expui/BasisFactory.cc
+++ b/expui/BasisFactory.cc
@@ -272,7 +272,7 @@ namespace BasisClasses
   // Generate coefficients from a phase-space table
   //
   CoefClasses::CoefStrPtr Basis::createFromArray
-  (Eigen::VectorXd m, RowMatrixXd p, double time, std::vector<double> ctr,
+  (Eigen::VectorXd& m, RowMatrixXd& p, double time, std::vector<double> ctr,
    bool roundrobin, bool posvelrows)
   {
     initFromArray(ctr);
diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H
index 518025c73..a4af27c1a 100644
--- a/expui/BiorthBasis.H
+++ b/expui/BiorthBasis.H
@@ -97,7 +97,7 @@ namespace BasisClasses
     //! Generate coeffients from an array and optional center location
     //! for the expansion
     CoefClasses::CoefStrPtr createFromArray
-    (Eigen::VectorXd m, RowMatrixXd p, double time=0.0,
+    (Eigen::VectorXd& m, RowMatrixXd& p, double time=0.0,
      std::vector<double> center={0.0, 0.0, 0.0},
      bool roundrobin=true, bool posvelrows=false);
     
@@ -113,7 +113,7 @@ namespace BasisClasses
     //! Initialize accumulating coefficients from arrays.  This is
     //! called once to initialize the accumulation.
     void addFromArray
-    (Eigen::VectorXd m, RowMatrixXd p,
+    (Eigen::VectorXd& m, RowMatrixXd& p,
      bool roundrobin=true, bool posvelrows=false);
     
     //! Create and the coefficients from the array accumulation with the
diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc
index 553d7106b..9625e4f2d 100644
--- a/expui/BiorthBasis.cc
+++ b/expui/BiorthBasis.cc
@@ -2847,7 +2847,7 @@ namespace BasisClasses
   }
 
   // Accumulate coefficient contributions from arrays
-  void BiorthBasis::addFromArray(Eigen::VectorXd m, RowMatrixXd p,
+  void BiorthBasis::addFromArray(Eigen::VectorXd& m, RowMatrixXd& p,
 				 bool RoundRobin, bool PosVelRows)
   {
     // Sanity check: is coefficient instance created?  This is not
@@ -2957,7 +2957,7 @@ namespace BasisClasses
   // Generate coefficients from a phase-space table
   //
   CoefClasses::CoefStrPtr BiorthBasis::createFromArray
-  (Eigen::VectorXd m, RowMatrixXd p, double time, std::vector<double> ctr,
+  (Eigen::VectorXd& m, RowMatrixXd& p, double time, std::vector<double> ctr,
    bool RoundRobin, bool PosVelRows)
   {
     initFromArray(ctr);
diff --git a/expui/Coefficients.H b/expui/Coefficients.H
index 46664237f..7718c7b7b 100644
--- a/expui/Coefficients.H
+++ b/expui/Coefficients.H
@@ -111,7 +111,7 @@ namespace CoefClasses
     virtual Eigen::VectorXcd& getData(double time) = 0;
 
     //! Set coefficient store at given time to the provided matrix
-    virtual void setData(double time, Eigen::VectorXcd data) = 0;
+    virtual void setData(double time, Eigen::VectorXcd& data) = 0;
 
     //! Interpolate coefficient matrix at given time
     std::tuple<Eigen::VectorXcd&, bool> interpolate(double time);
@@ -273,10 +273,10 @@ namespace CoefClasses
     Eigen::MatrixXcd& getMatrix(double time);
 
     //! Set coefficient matrix at given time
-    virtual void setData(double time, Eigen::VectorXcd arr);
+    virtual void setData(double time, Eigen::VectorXcd& arr);
 
     //! Natural data setter for pyEXP
-    void setMatrix(double time, Eigen::MatrixXcd mat);
+    void setMatrix(double time, Eigen::MatrixXcd& mat);
 
     //! Interpolate coefficient matrix at given time
     std::tuple<Eigen::MatrixXcd&, bool> interpolate(double time)
@@ -411,10 +411,10 @@ namespace CoefClasses
     Eigen::MatrixXcd& getMatrix(double time);
     
     //! Set coefficient matrix at given time
-    virtual void setData(double time, Eigen::VectorXcd arr);
+    virtual void setData(double time, Eigen::VectorXcd& arr);
     
     //! For pyEXP
-    void setMatrix(double time, Eigen::MatrixXcd arr);
+    void setMatrix(double time, Eigen::MatrixXcd& arr);
     
     //! Interpolate coefficient matrix at given time
     std::tuple<Eigen::MatrixXcd&, bool> interpolate(double time)
@@ -542,10 +542,10 @@ namespace CoefClasses
     virtual Eigen3d& getTensor(double time);
     
     //! Set coefficient matrix at given time
-    virtual void setData(double time, Eigen::VectorXcd dat);
+    virtual void setData(double time, Eigen::VectorXcd& dat);
 
     //! Native version
-    virtual void setTensor(double time, const Eigen3d dat);
+    virtual void setTensor(double time, const Eigen3d& dat);
     
     //! Interpolate coefficient tensor at given time
     std::tuple<Eigen::Tensor<std::complex<double>, 3>&, bool>
@@ -673,10 +673,10 @@ namespace CoefClasses
     virtual Eigen3d& getTensor(double time);
     
     //! Set coefficient matrix at given time
-    virtual void setData(double time, Eigen::VectorXcd dat);
+    virtual void setData(double time, Eigen::VectorXcd& dat);
 
     //! Native version
-    virtual void setTensor(double time, const Eigen3d dat);
+    virtual void setTensor(double time, const Eigen3d& dat);
     
     //! Interpolate coefficient tensor at given time
     std::tuple<Eigen::Tensor<std::complex<double>, 3>&, bool>
@@ -810,9 +810,9 @@ namespace CoefClasses
     /** Set coefficient matrix at given time.  This is for pybind11,
 	since the operator() will not allow lvalue assignment, it
 	seems. */
-    virtual void setData(double time, Eigen::VectorXcd arr);
+    virtual void setData(double time, Eigen::VectorXcd& arr);
 
-    void setArray(double time, Eigen::VectorXcd arr)
+    void setArray(double time, Eigen::VectorXcd& arr)
     { setData(time, arr); }
 
     //! Get coefficient structure at a given time
@@ -915,10 +915,10 @@ namespace CoefClasses
     SphFldStruct::dataType getMatrix(double time);
 
     //! Set coefficient matrix at given time
-    virtual void setData(double time, Eigen::VectorXcd arr);
+    virtual void setData(double time, Eigen::VectorXcd& arr);
 
     //! Natural data setter for pyEXP
-    void setMatrix(double time, SphFldStruct::dataType mat);
+    void setMatrix(double time, SphFldStruct::dataType& mat);
 
     //! Interpolate coefficient matrix at given time
     std::tuple<SphFldStruct::coefType&, bool> interpolate(double time)
@@ -1042,10 +1042,10 @@ namespace CoefClasses
     CylFldStruct::dataType getMatrix(double time);
 
     //! Set coefficient matrix at given time
-    virtual void setData(double time, Eigen::VectorXcd arr);
+    virtual void setData(double time, Eigen::VectorXcd& arr);
 
     //! Natural data setter for pyEXP
-    void setMatrix(double time, CylFldStruct::dataType mat);
+    void setMatrix(double time, CylFldStruct::dataType& mat);
 
     //! Interpolate coefficient matrix at given time
     std::tuple<CylFldStruct::coefType&, bool> interpolate(double time)
diff --git a/expui/Coefficients.cc b/expui/Coefficients.cc
index 83e352a94..fe55b0efe 100644
--- a/expui/Coefficients.cc
+++ b/expui/Coefficients.cc
@@ -469,7 +469,7 @@ namespace CoefClasses
     return mat;
   }
   
-  void SphCoefs::setData(double time, Eigen::VectorXcd dat)
+  void SphCoefs::setData(double time, Eigen::VectorXcd& dat)
   {
     auto it = coefs.find(roundTime(time));
 
@@ -484,7 +484,7 @@ namespace CoefClasses
     }
   }
   
-  void SphCoefs::setMatrix(double time, Eigen::MatrixXcd dat)
+  void SphCoefs::setMatrix(double time, Eigen::MatrixXcd& dat)
   {
     auto it = coefs.find(roundTime(time));
 
@@ -891,7 +891,7 @@ namespace CoefClasses
     return mat;
   }
   
-  void CylCoefs::setData(double time, Eigen::VectorXcd dat)
+  void CylCoefs::setData(double time, Eigen::VectorXcd& dat)
   {
     auto it = coefs.find(roundTime(time));
 
@@ -906,7 +906,7 @@ namespace CoefClasses
     }
   }
 
-  void CylCoefs::setMatrix(double time, Eigen::MatrixXcd dat)
+  void CylCoefs::setMatrix(double time, Eigen::MatrixXcd& dat)
   {
     auto it = coefs.find(roundTime(time));
 
@@ -1247,7 +1247,7 @@ namespace CoefClasses
     return arr;
   }
 
-  void SlabCoefs::setData(double time, Eigen::VectorXcd dat)
+  void SlabCoefs::setData(double time, Eigen::VectorXcd& dat)
   {
     auto it = coefs.find(roundTime(time));
 
@@ -1262,7 +1262,7 @@ namespace CoefClasses
     }
   }
 
-  void SlabCoefs::setTensor(double time, const Eigen3d dat)
+  void SlabCoefs::setTensor(double time, const Eigen3d& dat)
   {
     auto it = coefs.find(roundTime(time));
 
@@ -1599,7 +1599,7 @@ namespace CoefClasses
     return arr;
   }
 
-  void CubeCoefs::setData(double time, Eigen::VectorXcd dat)
+  void CubeCoefs::setData(double time, Eigen::VectorXcd& dat)
   {
     auto it = coefs.find(roundTime(time));
 
@@ -1614,7 +1614,7 @@ namespace CoefClasses
     }
   }
 
-  void CubeCoefs::setTensor(double time, const Eigen3d dat)
+  void CubeCoefs::setTensor(double time, const Eigen3d& dat)
   {
     auto it = coefs.find(roundTime(time));
 
@@ -1968,7 +1968,7 @@ namespace CoefClasses
     return arr;
   }
   
-  void TableData::setData(double time, Eigen::VectorXcd dat)
+  void TableData::setData(double time, Eigen::VectorXcd& dat)
   {
     auto it = coefs.find(roundTime(time));
 
@@ -2458,7 +2458,7 @@ namespace CoefClasses
     return *mat;
   }
   
-  void SphFldCoefs::setData(double time, Eigen::VectorXcd dat)
+  void SphFldCoefs::setData(double time, Eigen::VectorXcd& dat)
   {
     auto it = coefs.find(roundTime(time));
 
@@ -2473,7 +2473,7 @@ namespace CoefClasses
     }
   }
   
-  void SphFldCoefs::setMatrix(double time, SphFldStruct::dataType dat)
+  void SphFldCoefs::setMatrix(double time, SphFldStruct::dataType& dat)
   {
     auto it = coefs.find(roundTime(time));
 
@@ -2897,7 +2897,7 @@ namespace CoefClasses
     return *mat;
   }
   
-  void CylFldCoefs::setData(double time, Eigen::VectorXcd dat)
+  void CylFldCoefs::setData(double time, Eigen::VectorXcd& dat)
   {
     auto it = coefs.find(roundTime(time));
 
@@ -2912,7 +2912,7 @@ namespace CoefClasses
     }
   }
   
-  void CylFldCoefs::setMatrix(double time, CylFldStruct::dataType dat)
+  void CylFldCoefs::setMatrix(double time, CylFldStruct::dataType& dat)
   {
     auto it = coefs.find(roundTime(time));
 
diff --git a/expui/FieldBasis.H b/expui/FieldBasis.H
index b04653e3d..25ef5dbe3 100644
--- a/expui/FieldBasis.H
+++ b/expui/FieldBasis.H
@@ -146,7 +146,7 @@ namespace BasisClasses
     //! Initialize accumulating coefficients from arrays.  This is
     //! called once to initialize the accumulation.
     void addFromArray
-    (Eigen::VectorXd m, RowMatrixXd p, bool roundrobin=true, bool posvelrows=false);
+    (Eigen::VectorXd& m, RowMatrixXd& p, bool roundrobin=true, bool posvelrows=false);
     
     //! Accumulate new coefficients
     virtual void accumulate(double mass,
diff --git a/expui/FieldBasis.cc b/expui/FieldBasis.cc
index 2d25a1866..7c85e0aee 100644
--- a/expui/FieldBasis.cc
+++ b/expui/FieldBasis.cc
@@ -619,7 +619,7 @@ namespace BasisClasses
   }
 
   // Accumulate coefficient contributions from arrays
-  void FieldBasis::addFromArray(Eigen::VectorXd m, RowMatrixXd p,
+  void FieldBasis::addFromArray(Eigen::VectorXd& m, RowMatrixXd& p,
 				bool roundrobin, bool posvelrows)
   {
     // Sanity check: is coefficient instance created?  This is not
diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc
index 07e2369f0..e0b72ab25 100644
--- a/pyEXP/BasisWrappers.cc
+++ b/pyEXP/BasisWrappers.cc
@@ -286,7 +286,7 @@ void BasisFactoryClasses(py::module &m)
     }
 
     virtual void
-    addFromArray(Eigen::VectorXd m, RowMatrixXd p, bool roundrobin, bool posvelrows) override {
+    addFromArray(Eigen::VectorXd& m, RowMatrixXd& p, bool roundrobin, bool posvelrows) override {
       PYBIND11_OVERRIDE_PURE(void, Basis, addFromArray, m, p, roundrobin, posvelrows);
     }
   };
@@ -973,7 +973,7 @@ void BasisFactoryClasses(py::module &m)
          )",
 	 py::arg("center") = std::vector<double>(3, 0.0))
     .def("addFromArray",
-	 [](BasisClasses::BiorthBasis& A, Eigen::VectorXd mass, RowMatrixXd pos)
+	 [](BasisClasses::BiorthBasis& A, Eigen::VectorXd& mass, RowMatrixXd& pos)
 	 {
 	   return A.addFromArray(mass, pos);
 	 },
@@ -1804,7 +1804,7 @@ void BasisFactoryClasses(py::module &m)
          )",
 	 py::arg("center") = std::vector<double>(3, 0.0))
     .def("addFromArray",
-	 [](BasisClasses::FieldBasis& A, Eigen::VectorXd mass, RowMatrixXd ps)
+	 [](BasisClasses::FieldBasis& A, Eigen::VectorXd& mass, RowMatrixXd& ps)
 	 {
 	   return A.addFromArray(mass, ps);
 	 },
diff --git a/pyEXP/CoefWrappers.cc b/pyEXP/CoefWrappers.cc
index da2ce7227..aa58603b9 100644
--- a/pyEXP/CoefWrappers.cc
+++ b/pyEXP/CoefWrappers.cc
@@ -155,7 +155,7 @@ void CoefficientClasses(py::module &m) {
       PYBIND11_OVERRIDE_PURE(Eigen::VectorXcd&, Coefs, getData, time);
     }
 
-    void setData(double time, Eigen::VectorXcd array) override {
+    void setData(double time, Eigen::VectorXcd& array) override {
       PYBIND11_OVERRIDE_PURE(void, Coefs, setData, time, array);
     }
 
@@ -231,7 +231,7 @@ void CoefficientClasses(py::module &m) {
       PYBIND11_OVERRIDE(Eigen::VectorXcd&, SphCoefs, getData, time);
     }
 
-    void setData(double time, Eigen::VectorXcd array) override {
+    void setData(double time, Eigen::VectorXcd& array) override {
       PYBIND11_OVERRIDE(void, SphCoefs, setData, time, array);
     }
 
@@ -310,7 +310,7 @@ void CoefficientClasses(py::module &m) {
       PYBIND11_OVERRIDE(Eigen::VectorXcd&, CylCoefs, getData, time);
     }
 
-    void setData(double time, Eigen::VectorXcd array) override {
+    void setData(double time, Eigen::VectorXcd& array) override {
       PYBIND11_OVERRIDE(void, CylCoefs, setData, time, array);
     }
 
@@ -388,7 +388,7 @@ void CoefficientClasses(py::module &m) {
       PYBIND11_OVERRIDE(Eigen::VectorXcd&, SlabCoefs, getData, time);
     }
 
-    void setData(double time, Eigen::VectorXcd array) override {
+    void setData(double time, Eigen::VectorXcd& array) override {
       PYBIND11_OVERRIDE(void, SlabCoefs, setData, time, array);
     }
 
@@ -467,7 +467,7 @@ void CoefficientClasses(py::module &m) {
       PYBIND11_OVERRIDE(Eigen::VectorXcd&, CubeCoefs, getData, time);
     }
 
-    void setData(double time, Eigen::VectorXcd array) override {
+    void setData(double time, Eigen::VectorXcd& array) override {
       PYBIND11_OVERRIDE(void, CubeCoefs, setData, time, array);
     }
 
@@ -546,7 +546,7 @@ void CoefficientClasses(py::module &m) {
       PYBIND11_OVERRIDE(Eigen::VectorXcd&, TableData, getData, time);
     }
 
-    void setData(double time, Eigen::VectorXcd array) override {
+    void setData(double time, Eigen::VectorXcd& array) override {
       PYBIND11_OVERRIDE(void, TableData, setData, time, array);
     }
 
diff --git a/tests/Halo/createCoefs.py b/tests/Halo/createCoefs.py
index 851817b64..c72b5566d 100644
--- a/tests/Halo/createCoefs.py
+++ b/tests/Halo/createCoefs.py
@@ -1,4 +1,4 @@
-import os, time
+import os
 import pyEXP
 import random
 import numpy as np
@@ -27,16 +27,12 @@
 
 print("---- created basis")
 
-time.sleep(1)
-
 # Create a coefficient structure
 #
 coefs = pyEXP.coefs.SphCoefs(True)
 
 print("---- created coefficients")
 
-time.sleep(1)
-
 # Call the basis to generate coefficients
 #
 mass = []
@@ -52,15 +48,13 @@
     zpos.append(random.random()*2.0 - 1.0)
 
 print("---- createFromArray usings lists")
-time.sleep(1)
+
 coef1 = basis.createFromArray(mass, [xpos, ypos, zpos], time=3.0)
-time.sleep(1)
+
 coefs.add(coef1)
 
 print("Times:", coefs.Times())
 
-exit(0)
-
 print("---- creating array data from list data")
 
 # Note: this overwrites mass and data variables.  But data is now

From d2866db40c628f1ff317db332472044a97d694aa Mon Sep 17 00:00:00 2001
From: Martin Weinberg <mdw@umass.edu>
Date: Mon, 14 Oct 2024 11:54:53 -0400
Subject: [PATCH 37/52] Update CMakeLists.txt

to reassign `pyEXPCoefCreateTest` to `long`
---
 tests/CMakeLists.txt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 560f9e8b5..9f315ce7c 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -100,7 +100,7 @@ if(ENABLE_NBODY)
               PYTHONPATH=${CMAKE_BINARY_DIR}/pyEXP:$ENV{PYTHONPATH}
               ${PYTHON_EXECUTABLE} createCoefs.py
               WORKING_DIRECTORY "${CMAKE_CURRENT_LIST_DIR}/Halo")
-     set_tests_properties(pyEXPCoefCreateTest PROPERTIES LABELS "quick")
+     set_tests_properties(pyEXPCoefCreateTest PROPERTIES LABELS "long")
   endif()
 
   # A separate test to remove the generated files if they all exist;

From a0077057d1b5b7f0a189645afae05d7405145b8c Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Mon, 14 Oct 2024 17:41:16 -0400
Subject: [PATCH 38/52] git submodule updated

---
 extern/pybind11 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/extern/pybind11 b/extern/pybind11
index 914c06fb2..58c382a8e 160000
--- a/extern/pybind11
+++ b/extern/pybind11
@@ -1 +1 @@
-Subproject commit 914c06fb252b6cc3727d0eedab6736e88a3fcb01
+Subproject commit 58c382a8e3d7081364d2f5c62e7f429f0412743b

From 02d09fb8cd99cb2c3f3cff06eddf574ac395a480 Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Mon, 14 Oct 2024 18:31:40 -0400
Subject: [PATCH 39/52] Relabel test as 'quick' after updating pybind11

---
 .github/workflows/build.yml | 4 ++--
 extern/HighFive             | 2 +-
 extern/yaml-cpp             | 2 +-
 tests/CMakeLists.txt        | 2 +-
 4 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
index 85f50c5bb..ee04f774a 100644
--- a/.github/workflows/build.yml
+++ b/.github/workflows/build.yml
@@ -25,7 +25,7 @@ jobs:
         if: runner.os == 'Linux'
         run: |
           sudo apt-get update
-          sudo apt-get install -y build-essential libeigen3-dev libfftw3-dev libhdf5-dev libopenmpi-dev
+          sudo apt-get install -y build-essential libeigen3-dev libfftw3-dev libhdf5-dev libopenmpi-dev cmake gdb python3-dev libpython3-dev
           sudo pip install numpy
 
 
@@ -43,7 +43,7 @@ jobs:
           cmake
           -DENABLE_NBODY=YES
           -DENABLE_PYEXP=YES
-          -DCMAKE_BUILD_TYPE=Release
+          -DCMAKE_BUILD_TYPE=Debug
           -DEigen3_DIR=/usr/include/eigen3/share/eigen3/cmake
           -DCMAKE_INSTALL_PREFIX=./install
           -Wno-dev
diff --git a/extern/HighFive b/extern/HighFive
index 12064079d..6c9624792 160000
--- a/extern/HighFive
+++ b/extern/HighFive
@@ -1 +1 @@
-Subproject commit 12064079d9533c0636310f25606571b3dbb977cf
+Subproject commit 6c9624792cada7e4cd5e6d4a8c53204242cc8455
diff --git a/extern/yaml-cpp b/extern/yaml-cpp
index 1d8ca1f35..da82fd982 160000
--- a/extern/yaml-cpp
+++ b/extern/yaml-cpp
@@ -1 +1 @@
-Subproject commit 1d8ca1f35eb3a9c9142462b28282a848e5d29a91
+Subproject commit da82fd982c260e7f335ce5acbceff24b270544d1
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 9f315ce7c..560f9e8b5 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -100,7 +100,7 @@ if(ENABLE_NBODY)
               PYTHONPATH=${CMAKE_BINARY_DIR}/pyEXP:$ENV{PYTHONPATH}
               ${PYTHON_EXECUTABLE} createCoefs.py
               WORKING_DIRECTORY "${CMAKE_CURRENT_LIST_DIR}/Halo")
-     set_tests_properties(pyEXPCoefCreateTest PROPERTIES LABELS "long")
+     set_tests_properties(pyEXPCoefCreateTest PROPERTIES LABELS "quick")
   endif()
 
   # A separate test to remove the generated files if they all exist;

From b215b496466223ad89df46b29a155f7a408a9543 Mon Sep 17 00:00:00 2001
From: Martin Weinberg <mdw@umass.edu>
Date: Mon, 14 Oct 2024 19:57:08 -0400
Subject: [PATCH 40/52] Update build.yml

Reverted from `Debug` to `Release`
---
 .github/workflows/build.yml | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
index ee04f774a..76e13a4e3 100644
--- a/.github/workflows/build.yml
+++ b/.github/workflows/build.yml
@@ -25,7 +25,7 @@ jobs:
         if: runner.os == 'Linux'
         run: |
           sudo apt-get update
-          sudo apt-get install -y build-essential libeigen3-dev libfftw3-dev libhdf5-dev libopenmpi-dev cmake gdb python3-dev libpython3-dev
+          sudo apt-get install -y build-essential libeigen3-dev libfftw3-dev libhdf5-dev libopenmpi-dev
           sudo pip install numpy
 
 
@@ -43,7 +43,7 @@ jobs:
           cmake
           -DENABLE_NBODY=YES
           -DENABLE_PYEXP=YES
-          -DCMAKE_BUILD_TYPE=Debug
+          -DCMAKE_BUILD_TYPE=Release
           -DEigen3_DIR=/usr/include/eigen3/share/eigen3/cmake
           -DCMAKE_INSTALL_PREFIX=./install
           -Wno-dev
@@ -59,4 +59,4 @@ jobs:
         
       #- name: CTest Long
       #working-directory: ./build
-      #run: ctest -L long
+      #run: ctest --output-on-failure -L long

From ac91538f7fe87d2ecf1ffc8623c25412e3e4593a Mon Sep 17 00:00:00 2001
From: Martin Weinberg <mdw@umass.edu>
Date: Tue, 15 Oct 2024 15:09:43 -0400
Subject: [PATCH 41/52] Update build.yml

test gcc and clang matrix
---
 .github/workflows/build.yml | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
index 76e13a4e3..deac202ad 100644
--- a/.github/workflows/build.yml
+++ b/.github/workflows/build.yml
@@ -13,7 +13,13 @@ jobs:
     strategy:
       matrix:
         os: [ubuntu-latest, ubuntu-22.04]
-        cc: [gcc, mpicc]
+        c_compiler: [gcc, clang]
+        cpp_compiler: [g++, clang++]
+        include:
+          - c_compiler: gcc
+            cpp_compiler: g++
+          - c_compiler: clang
+            cpp_compiler: clang++
 
     name: "Test pyEXP Build"
     runs-on: ${{ matrix.os }}

From cdcab974e7a3e9aa167160ce9bf862891a9bf068 Mon Sep 17 00:00:00 2001
From: Martin Weinberg <mdw@umass.edu>
Date: Tue, 15 Oct 2024 15:34:35 -0400
Subject: [PATCH 42/52] Update build.yml

Another try.
---
 .github/workflows/build.yml | 11 +++++------
 1 file changed, 5 insertions(+), 6 deletions(-)

diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
index deac202ad..02e491330 100644
--- a/.github/workflows/build.yml
+++ b/.github/workflows/build.yml
@@ -13,13 +13,11 @@ jobs:
     strategy:
       matrix:
         os: [ubuntu-latest, ubuntu-22.04]
-        c_compiler: [gcc, clang]
-        cpp_compiler: [g++, clang++]
         include:
-          - c_compiler: gcc
-            cpp_compiler: g++
-          - c_compiler: clang
-            cpp_compiler: clang++
+          - cc: gcc
+            cxx: g++
+          - cc: clang
+            cxx: clang++
 
     name: "Test pyEXP Build"
     runs-on: ${{ matrix.os }}
@@ -44,6 +42,7 @@ jobs:
         if: runner.os == 'Linux'
         env:
           CC: ${{ matrix.cc }}
+          CXX: ${{ matrix.cxx }}
         working-directory: ./build
         run: >-
           cmake

From 28fa507299fd3a99a0a3eda796ad60c09e6cecc3 Mon Sep 17 00:00:00 2001
From: Martin Weinberg <mdw@umass.edu>
Date: Tue, 15 Oct 2024 15:47:05 -0400
Subject: [PATCH 43/52] Update build.yml

Need OpenMP support for clang
---
 .github/workflows/build.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
index 02e491330..2b3280ea9 100644
--- a/.github/workflows/build.yml
+++ b/.github/workflows/build.yml
@@ -29,7 +29,7 @@ jobs:
         if: runner.os == 'Linux'
         run: |
           sudo apt-get update
-          sudo apt-get install -y build-essential libeigen3-dev libfftw3-dev libhdf5-dev libopenmpi-dev
+          sudo apt-get install -y build-essential libeigen3-dev libfftw3-dev libhdf5-dev libopenmpi-dev libomp-dev
           sudo pip install numpy
 
 

From ae7fa0ad2c388ba64700fd577d9a87f5ba767502 Mon Sep 17 00:00:00 2001
From: Martin Weinberg <mdw@umass.edu>
Date: Tue, 15 Oct 2024 16:05:37 -0400
Subject: [PATCH 44/52] Update build.yml

Need flags for clang++ to use stdc++ lib.
---
 .github/workflows/build.yml | 10 ++++++++--
 1 file changed, 8 insertions(+), 2 deletions(-)

diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
index 2b3280ea9..b2826df0f 100644
--- a/.github/workflows/build.yml
+++ b/.github/workflows/build.yml
@@ -13,11 +13,16 @@ jobs:
     strategy:
       matrix:
         os: [ubuntu-latest, ubuntu-22.04]
+        cc: [gcc, clang]
+        cxx: [g++, clang++]
         include:
+          - cxx: clang++
+            cxxflags: -stdlib=libstdc++
+        exclude:
           - cc: gcc
-            cxx: g++
-          - cc: clang
             cxx: clang++
+          - cc: clang
+            cxx: g++
 
     name: "Test pyEXP Build"
     runs-on: ${{ matrix.os }}
@@ -51,6 +56,7 @@ jobs:
           -DCMAKE_BUILD_TYPE=Release
           -DEigen3_DIR=/usr/include/eigen3/share/eigen3/cmake
           -DCMAKE_INSTALL_PREFIX=./install
+          -DCMAKE_CXX_FLAGS=${{ matrix.cxxflags }}
           -Wno-dev
           ..
 

From 60ffa3d7538a3a3e1db496767c9e94f0d1d78dcc Mon Sep 17 00:00:00 2001
From: Martin Weinberg <mdw@umass.edu>
Date: Tue, 15 Oct 2024 16:50:21 -0400
Subject: [PATCH 45/52] Update CMakeLists.txt

Remove hardwired use of llvm's libc++
---
 CMakeLists.txt | 1 -
 1 file changed, 1 deletion(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index dbede3fe3..f61876725 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -13,7 +13,6 @@ set(CMAKE_CXX_EXTENSIONS OFF)
 # Compiler flags.  Not all tested thoroughly...
 if(CMAKE_CXX_COMPILER_ID MATCHES "Clang")
   # using Clang
-  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -stdlib=libc++")
 elseif(CMAKE_CXX_COMPILER_ID MATCHES "GNU")
   # using GCC
 elseif(CMAKE_CXX_COMPILER_ID MATCHES "Intel")

From 7d9cf631c990ceb188f16b2ef0104c163e60c8f8 Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Tue, 15 Oct 2024 17:22:02 -0400
Subject: [PATCH 46/52] Decrease of the clang chatter; add 24.04 explicitly to
 the CI matrix

---
 .github/workflows/build.yml | 2 +-
 include/DiskWithHalo.H      | 6 +++---
 include/EXPini.H            | 3 ++-
 include/SLGridMP2.H         | 3 +++
 include/massmodel.H         | 6 +++---
 5 files changed, 12 insertions(+), 8 deletions(-)

diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
index b2826df0f..b81c922b4 100644
--- a/.github/workflows/build.yml
+++ b/.github/workflows/build.yml
@@ -12,7 +12,7 @@ jobs:
   exp:
     strategy:
       matrix:
-        os: [ubuntu-latest, ubuntu-22.04]
+        os: [ubuntu-24.04, ubuntu-22.04]
         cc: [gcc, clang]
         cxx: [g++, clang++]
         include:
diff --git a/include/DiskWithHalo.H b/include/DiskWithHalo.H
index 79e9038f7..8f96ce36d 100644
--- a/include/DiskWithHalo.H
+++ b/include/DiskWithHalo.H
@@ -66,13 +66,13 @@ public:
     dur = dur1 + dur2;
   }
   
-  double get_mass(const double x1, const double x2, const double x3)
+  double get_mass(const double x1, const double x2, const double x3) override
   { return get_mass(sqrt(x1*x1 + x2*x2 + x3*x3)); }
   
-  double get_density(const double x1, const double x2, const double x3)
+  double get_density(const double x1, const double x2, const double x3) override
   { return get_density(sqrt(x1*x1 + x2*x2 + x3*x3)); }
   
-  double get_pot(const double x1, const double x2, const double x3)
+  double get_pot(const double x1, const double x2, const double x3) override
   { return get_pot(sqrt(x1*x1 + x2*x2 + x3*x3)); }
 
   // Addiional member functions
diff --git a/include/EXPini.H b/include/EXPini.H
index b76e39ff0..79ac22e59 100644
--- a/include/EXPini.H
+++ b/include/EXPini.H
@@ -87,8 +87,9 @@ cxxopts::ParseResult LoadConfig(cxxopts::Options& options,
 {
   YAML::Node conf = YAML::LoadFile(config);
 
-  int count = conf.size()*2+1, cnt = 1;
+  const int count = conf.size()*2+1;
   char* data[count];
+  int cnt = 1;
 
   data[0] = new char [11];
   strcpy(data[0], "LoadConfig"); // Emulate the caller name
diff --git a/include/SLGridMP2.H b/include/SLGridMP2.H
index f83d98fe9..e521b90ef 100644
--- a/include/SLGridMP2.H
+++ b/include/SLGridMP2.H
@@ -314,6 +314,9 @@ private:
     //! Constructor
     CoordMap(double H) : H(H) {}
   
+    //! Destructor
+    virtual ~CoordMap() {}
+
     //! Convert from vertical to mapped coordinate
     virtual double z_to_xi  (double z) = 0;
 
diff --git a/include/massmodel.H b/include/massmodel.H
index a91af5a8e..88939322d 100644
--- a/include/massmodel.H
+++ b/include/massmodel.H
@@ -188,13 +188,13 @@ public:
   virtual double get_dpot2(const double) = 0;
   virtual void get_pot_dpot(const double, double&, double&) = 0;
   
-  double get_mass(const double x1, const double x2, const double x3)
+  double get_mass(const double x1, const double x2, const double x3) override
   { return get_mass(sqrt(x1*x1 + x2*x2 + x3*x3)); }
   
-  double get_density(const double x1, const double x2, const double x3)
+  double get_density(const double x1, const double x2, const double x3) override
   { return get_density(sqrt(x1*x1 + x2*x2 + x3*x3)); }
   
-  double get_pot(const double x1, const double x2, const double x3)
+  double get_pot(const double x1, const double x2, const double x3) override
   { return get_pot(sqrt(x1*x1 + x2*x2 + x3*x3)); }
   //@}
   

From 136fa1606349384b9a21b9991f609231b9693bc4 Mon Sep 17 00:00:00 2001
From: Martin Weinberg <mdw@umass.edu>
Date: Tue, 15 Oct 2024 17:26:10 -0400
Subject: [PATCH 47/52] Update build.yml

Use apt numpy package for 24.04
---
 .github/workflows/build.yml | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
index b81c922b4..f530765b1 100644
--- a/.github/workflows/build.yml
+++ b/.github/workflows/build.yml
@@ -34,9 +34,7 @@ jobs:
         if: runner.os == 'Linux'
         run: |
           sudo apt-get update
-          sudo apt-get install -y build-essential libeigen3-dev libfftw3-dev libhdf5-dev libopenmpi-dev libomp-dev
-          sudo pip install numpy
-
+          sudo apt-get install -y build-essential libeigen3-dev libfftw3-dev libhdf5-dev libopenmpi-dev libomp-dev python3-numpy
 
       - name: Setup submodule and build
         run: |

From 0218739a308e2243d85333b9ac5d52ccf48c8d2d Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Tue, 15 Oct 2024 20:34:44 -0400
Subject: [PATCH 48/52] This commit addresses all of the warnings and
 non-compliance issues identified by Clang and Clang++.  They are all minor. 
 The result passes 'ctest -L long' but pyEXP-example notebooks have not been
 run.

---
 expui/BiorthBasis.cc        | 10 +++++-----
 expui/CoefContainer.cc      |  6 +++---
 exputil/BarrierWrapper.cc   |  7 ++++---
 exputil/EmpCylSL.cc         |  2 +-
 exputil/GaussCore.c         | 20 ++++++++------------
 exputil/GaussCore.h         |  6 ++++--
 exputil/Hermite.c           | 13 ++++---------
 exputil/Hermite.h           |  9 ++++++---
 exputil/Jacobi.c            | 18 ++++++------------
 exputil/Jacobi.h            | 11 ++++++++---
 exputil/Laguerre.c          |  9 +++------
 exputil/Laguerre.h          |  6 ++++--
 exputil/ParticleReader.cc   | 13 +++++++------
 exputil/YamlConfig.cc       |  2 +-
 include/EXPini.H            |  4 ++--
 utils/Analysis/gas2dcyl.cc  |  6 +++---
 utils/Analysis/psphisto.cc  |  9 +++++----
 utils/ICs/SphericalSL.cc    |  2 +-
 utils/ICs/gensph.cc         | 14 +++++++-------
 utils/PhaseSpace/pspmono.cc |  6 +++---
 20 files changed, 85 insertions(+), 88 deletions(-)

diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc
index 9625e4f2d..8e319109b 100644
--- a/expui/BiorthBasis.cc
+++ b/expui/BiorthBasis.cc
@@ -368,7 +368,7 @@ namespace BasisClasses
   {
     // Sanity check on derived class type
     //
-    if (typeid(*coef) != typeid(CoefClasses::SphStruct))
+    if (typeid(coef.get()) != typeid(CoefClasses::SphStruct*))
       throw std::runtime_error("Spherical::set_coefs: you must pass a CoefClasses::SphStruct");
 
     // Sanity check on dimensionality
@@ -1431,7 +1431,7 @@ namespace BasisClasses
 
   void Cylindrical::set_coefs(CoefClasses::CoefStrPtr coef)
   {
-    if (typeid(*coef) != typeid(CoefClasses::CylStruct))
+    if (typeid(coef.get()) != typeid(CoefClasses::CylStruct*))
       throw std::runtime_error("Cylindrical::set_coefs: you must pass a CoefClasses::CylStruct");
 
     CoefClasses::CylStruct* cf = dynamic_cast<CoefClasses::CylStruct*>(coef.get());
@@ -1708,7 +1708,7 @@ namespace BasisClasses
   {
     // Sanity check on derived class type
     //
-    if (typeid(*coef) != typeid(CoefClasses::CylStruct))
+    if (typeid(coef.get()) != typeid(CoefClasses::CylStruct*))
       throw std::runtime_error("FlatDisk::set_coefs: you must pass a CoefClasses::CylStruct");
 
     // Sanity check on dimensionality
@@ -2140,7 +2140,7 @@ namespace BasisClasses
   {
     // Sanity check on derived class type
     //
-    if (typeid(*coef) != typeid(CoefClasses::SlabStruct))
+    if (typeid(coef) != typeid(CoefClasses::SlabStruct*))
       throw std::runtime_error("Slab::set_coefs: you must pass a CoefClasses::SlabStruct");
 
     // Sanity check on dimensionality
@@ -2572,7 +2572,7 @@ namespace BasisClasses
   {
     // Sanity check on derived class type
     //
-    if (typeid(*coef) != typeid(CoefClasses::CubeStruct))
+    if (typeid(coef.get()) != typeid(CoefClasses::CubeStruct*))
       throw std::runtime_error("Cube::set_coefs: you must pass a CoefClasses::CubeStruct");
 
     // Sanity check on dimensionality
diff --git a/expui/CoefContainer.cc b/expui/CoefContainer.cc
index d5deb13da..0d48b17fe 100644
--- a/expui/CoefContainer.cc
+++ b/expui/CoefContainer.cc
@@ -72,7 +72,7 @@ namespace MSSA
       unpack_cylfld();
     }
     else {
-      throw std::runtime_error(std::string("CoefDB::unpack_channels(): can not reflect coefficient type=") + typeid(*coefs.get()).name());
+      throw std::runtime_error(std::string("CoefDB::unpack_channels(): can not reflect coefficient type=") + typeid(coefs.get()).name());
     }
 
     return coefs;
@@ -91,7 +91,7 @@ namespace MSSA
     else if (dynamic_cast<CoefClasses::CylFldCoefs*>(coefs.get()))
       { } // Do nothing
     else {
-      throw std::runtime_error(std::string("CoefDB::background(): can not reflect coefficient type=") + typeid(*coefs.get()).name());
+      throw std::runtime_error(std::string("CoefDB::background(): can not reflect coefficient type=") + typeid(coefs.get()).name());
     }
   }
 
@@ -112,7 +112,7 @@ namespace MSSA
     else if (dynamic_cast<CoefClasses::CylFldCoefs*>(coefs.get()))
       pack_cylfld();
     else {
-      throw std::runtime_error(std::string("CoefDB::pack_channels(): can not reflect coefficient type=") + typeid(*coefs.get()).name());
+      throw std::runtime_error(std::string("CoefDB::pack_channels(): can not reflect coefficient type=") + typeid(coefs.get()).name());
     }
   }
 
diff --git a/exputil/BarrierWrapper.cc b/exputil/BarrierWrapper.cc
index 03755f62d..514aa4e05 100644
--- a/exputil/BarrierWrapper.cc
+++ b/exputil/BarrierWrapper.cc
@@ -101,12 +101,13 @@ void BarrierWrapper::light_operator(const string& label,
 
 				// Compare adjacent strings in the list
     if (localid==0) {
-      char tmp[cbufsz];		// Working buffer
+				// Working buffer
+      auto tmp = std::make_unique<char[]>(cbufsz);
       bool firstime   = true;
-      string one, two = strncpy(tmp, &bufferT[0], cbufsz);
+      string one, two = strncpy(tmp.get(), &bufferT[0], cbufsz);
       for (int n=1; n<commsize; n++) {
 	one = two;
-	two = strncpy(tmp, &bufferT[cbufsz*n], cbufsz);
+	two = strncpy(tmp.get(), &bufferT[cbufsz*n], cbufsz);
 	if (one.compare(two) != 0) {
 	  if (firstime) {
 	    cout << std::string(60,'-') << std::endl << left
diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc
index e69be02aa..01b941501 100644
--- a/exputil/EmpCylSL.cc
+++ b/exputil/EmpCylSL.cc
@@ -5875,7 +5875,7 @@ void EmpCylSL::dump_images(const string& OUTFILE,
   //============
   // Open files
   //============
-  std::ofstream out[Number];
+  auto out = std::make_unique<std::ofstream[]>(Number);
   for (int j=0; j<Number; j++) {
     Name = OUTFILE + Types[j] + ".eof_recon";
     out[j].open(Name.c_str());
diff --git a/exputil/GaussCore.c b/exputil/GaussCore.c
index 3d3b6523c..52acb210b 100644
--- a/exputil/GaussCore.c
+++ b/exputil/GaussCore.c
@@ -26,19 +26,18 @@
 /*
    Forward declaration for function QQp:
 */
-static int QQp();
+static int QQp(double, double*, double*);
 
 /*
    Function to test a real value for exceeding -1,
    and quit on an error if condition not met.
 */
-void GaussCheck(value)
-double value;
+void GaussCheck(double value)
 {
-    if (value <= (-1.0)) {
-	fprintf(stderr, "Gauss package: parameter out of range: %g\n", value);
-	exit(1);
-    }
+  if (value <= (-1.0)) {
+    fprintf(stderr, "Gauss package: parameter out of range: %g\n", value);
+    exit(1);
+  }
 }
 
 
@@ -54,11 +53,8 @@ static int n1;
 /*
    The workhorse.
 */
-void GaussMaster(n, alpha, beta, conflag, abscis, weight)
-int n;
-double alpha, beta;
-int conflag;
-double abscis[], weight[];
+void GaussMaster(int n, double alpha, double beta, int conflag,
+		 double abscis[], double weight[])
 {
 #define FALSE 0
 #define  TRUE 1
diff --git a/exputil/GaussCore.h b/exputil/GaussCore.h
index d8ef4fc03..61386c76a 100644
--- a/exputil/GaussCore.h
+++ b/exputil/GaussCore.h
@@ -12,8 +12,10 @@
 #define   CONFLUENT 1
 #define NOCONFLUENT 0
 
-extern void GaussMaster();
-extern void GaussCheck();
+extern void GaussMaster(int n, double alpha, double beta, int conflag,
+			double abscis[], double weight[]);
+
+extern void GaussCheck(double value);
 
 /* ARGUMENT LISTS:
 
diff --git a/exputil/Hermite.c b/exputil/Hermite.c
index b7ed38b3a..38be268a4 100644
--- a/exputil/Hermite.c
+++ b/exputil/Hermite.c
@@ -23,9 +23,8 @@
 */
 #define odd(n) ((unsigned)(n) & 01)
 
-void Odd_Hermite(n, alpha, abscis, weight, w0)
-int n;
-double alpha, abscis[], weight[], *w0;
+void Odd_Hermite(int n, double alpha, double abscis[], double weight[],
+		 double* w0)
 {
     int k;
 
@@ -41,9 +40,7 @@ double alpha, abscis[], weight[], *w0;
 };
 
 
-void Even_Hermite(n, alpha, abscis, weight)
-int n;
-double alpha, abscis[], weight[];
+void Even_Hermite(int n, double alpha, double abscis[], double weight[])
 {
     int k;
 
@@ -59,9 +56,7 @@ double alpha, abscis[], weight[];
 };
 
 
-void Hermite(n, alpha, abscis, weight)
-int n;
-double alpha, abscis[], weight[];
+void Hermite(int n, double alpha, double abscis[], double weight[])
 {
     int n1, k;
 
diff --git a/exputil/Hermite.h b/exputil/Hermite.h
index b5ff49370..76d6c078a 100644
--- a/exputil/Hermite.h
+++ b/exputil/Hermite.h
@@ -15,9 +15,12 @@
    Translation of module Hermite (export).
 */
 
-extern void Odd_Hermite();
-extern void Even_Hermite();
-extern void Hermite();
+extern void Odd_Hermite(int n, double alpha, double abscis[], double weight[],
+			double* w0);
+
+extern void Even_Hermite(int n, double alpha, double abscis[], double weight[]);
+
+extern void Hermite(int n, double alpha, double abscis[], double weight[]);
 
 /* ARGUMENT LISTS:
    
diff --git a/exputil/Jacobi.c b/exputil/Jacobi.c
index afe4557d1..a71647812 100644
--- a/exputil/Jacobi.c
+++ b/exputil/Jacobi.c
@@ -20,10 +20,7 @@
 
 #define sqr(x) ((x)*(x))
 
-void Jacobi(n, alpha, beta, abscis, weight)
-int n;
-double alpha, beta;
-double abscis[], weight[];
+void Jacobi(int n, double alpha, double beta, double abscis[], double weight[])
 {
     int k;
 
@@ -34,10 +31,8 @@ double abscis[], weight[];
 }
 
 
-void Radau_Jacobi(n, alpha, beta, abscis, weight, leftw)
-int n;
-double alpha, beta;
-double abscis[], weight[], *leftw;
+void Radau_Jacobi(int n, double alpha, double beta,
+		  double abscis[], double weight[], double *leftw)
 {
     int k;
     double temp;
@@ -62,10 +57,9 @@ double abscis[], weight[], *leftw;
 };
 
 
-void Lobatto_Jacobi(n, alpha, beta, abscis, weight, leftw, rightw)
-int n;
-double alpha, beta;
-double abscis[], weight[], *leftw, *rightw;
+void Lobatto_Jacobi(int n, double alpha, double beta,
+		    double abscis[], double weight[],
+		    double* leftw, double* rightw)
 {
     int k;
     double temp1, temp2;
diff --git a/exputil/Jacobi.h b/exputil/Jacobi.h
index d719820c0..9ae931820 100644
--- a/exputil/Jacobi.h
+++ b/exputil/Jacobi.h
@@ -15,10 +15,15 @@
    Translation of module Jacobi (export).
 */
 
+extern void Jacobi(int n, double alpha, double beta,
+		   double abscis[], double weight[]);
 
-extern void Jacobi();
-extern void Radau_Jacobi();
-extern void Lobatto_Jacobi();
+extern void Radau_Jacobi(int n, double alpha, double beta,
+			 double abscis[], double weight[], double *leftw);
+
+extern void Lobatto_Jacobi(int n, double alpha, double beta,
+			   double abscis[], double weight[],
+			   double* leftw, double* rightw);
 
 /* ARGUMENT LISTS:
    void Jacobi(n, alpha, beta, abscis, weight)
diff --git a/exputil/Laguerre.c b/exputil/Laguerre.c
index 2bc97ab9a..d73046c37 100644
--- a/exputil/Laguerre.c
+++ b/exputil/Laguerre.c
@@ -15,18 +15,15 @@
 #include "GaussCore.h"
 #include "Laguerre.h"
 
-void Laguerre(n, alpha, abscis, weight)
-int n;
-double alpha, abscis[], weight[];
+void Laguerre(int n, double alpha, double abscis[], double weight[])
 {
     GaussCheck(alpha);
     GaussMaster(n, alpha, 0.0, CONFLUENT, abscis, weight);
 };
 
 
-void Radau_Laguerre(n, alpha, abscis, weight, leftw)
-int n;
-double alpha, abscis[], weight[], *leftw;
+void Radau_Laguerre(int n, double alpha, double abscis[], double weight[],
+		    double* leftw)
 {
     int k;
     double temp;
diff --git a/exputil/Laguerre.h b/exputil/Laguerre.h
index 2afdfc948..df53af5f6 100644
--- a/exputil/Laguerre.h
+++ b/exputil/Laguerre.h
@@ -15,8 +15,10 @@
    Translation of module Laguerre (export).
 */
 
-extern void Laguerre();
-extern void Radau_Laguerre();
+extern void Laguerre(int n, double alpha, double abscis[], double weight[]);
+
+extern void Radau_Laguerre(int n, double alpha,
+			   double abscis[], double weight[], double* leftw);
 
 /* 
    ARGUMENT LISTS:
diff --git a/exputil/ParticleReader.cc b/exputil/ParticleReader.cc
index 99f09a8c4..c0934aebb 100644
--- a/exputil/ParticleReader.cc
+++ b/exputil/ParticleReader.cc
@@ -3,6 +3,7 @@
 #include <iomanip>
 #include <fstream>
 #include <sstream>
+#include <memory>
 #include <numeric>
 #include <string>
 #include <vector>
@@ -530,11 +531,11 @@ namespace PR {
 	    
 	    int rank = dataspace.getSimpleExtentNdims();
 	    
-	    hsize_t dims[rank];
+	    auto dims = std::make_unique<hsize_t[]>(rank);
 	    
-	    int ndims = dataspace.getSimpleExtentDims(dims, NULL);
+	    int ndims = dataspace.getSimpleExtentDims(dims.get(), NULL);
 	    
-	    H5::DataSpace mspace(rank, dims);
+	    H5::DataSpace mspace(rank, dims.get());
 	    
 	    std::vector<float> masses(dims[0]);
 	    dataset.read(&masses[0], H5::PredType::NATIVE_FLOAT, mspace, dataspace );
@@ -569,11 +570,11 @@ namespace PR {
 	  
 	  int rank = dataspace.getSimpleExtentNdims();
 	  
-	  hsize_t dims[rank];
+	  auto dims = std::make_unique<hsize_t[]>(rank);
 	  
-	  int ndims = dataspace.getSimpleExtentDims( dims, NULL);
+	  int ndims = dataspace.getSimpleExtentDims( dims.get(), NULL);
 	  
-	  H5::DataSpace mspace(rank, dims);
+	  H5::DataSpace mspace(rank, dims.get());
 	  
 	  std::vector<unsigned> seq(dims[0]);
 	  dataset.read(&seq[0], H5::PredType::NATIVE_UINT32, mspace, dataspace );
diff --git a/exputil/YamlConfig.cc b/exputil/YamlConfig.cc
index 8f7124de1..aa2666d22 100644
--- a/exputil/YamlConfig.cc
+++ b/exputil/YamlConfig.cc
@@ -71,7 +71,7 @@ cxxopts::ParseResult LoadConfig(cxxopts::Options& options,
   YAML::Node conf = YAML::LoadFile(config);
 
   int count = conf.size()*2+1, cnt = 1;
-  char* data[count];
+  std::vector<char*> data(count);
 
   data[0] = new char [11];
   strcpy(data[0], "LoadConfig"); // Emulate the caller name
diff --git a/include/EXPini.H b/include/EXPini.H
index 79ac22e59..0690eea98 100644
--- a/include/EXPini.H
+++ b/include/EXPini.H
@@ -88,7 +88,7 @@ cxxopts::ParseResult LoadConfig(cxxopts::Options& options,
   YAML::Node conf = YAML::LoadFile(config);
 
   const int count = conf.size()*2+1;
-  char* data[count];
+  std::vector<char*> data(count);
   int cnt = 1;
 
   data[0] = new char [11];
@@ -131,7 +131,7 @@ cxxopts::ParseResult LoadConfig(cxxopts::Options& options,
   
   auto vm = options.parse(cnt, &data[0]);
 
-  for (int i=0; i<cnt; i++) delete [] data[i];
+  for (auto & v : data) delete [] v;
 
   return vm;
 }
diff --git a/utils/Analysis/gas2dcyl.cc b/utils/Analysis/gas2dcyl.cc
index f37d5aae7..1fd2e6722 100644
--- a/utils/Analysis/gas2dcyl.cc
+++ b/utils/Analysis/gas2dcyl.cc
@@ -171,9 +171,9 @@ main(int argc, char **argv)
       char *c = const_cast<char*>(files[n].c_str());
       MPI_Bcast(c, sz+1, MPI_CHAR, 0, MPI_COMM_WORLD);
     } else {
-      char l[sz+1];
-      MPI_Bcast(&l[0], sz+1, MPI_CHAR, 0, MPI_COMM_WORLD);
-      files.push_back(l);
+      auto l = std::make_unique<char[]>(sz+1);
+      MPI_Bcast(l.get(), sz+1, MPI_CHAR, 0, MPI_COMM_WORLD);
+      files.push_back(l.get());
     }
     MPI_Barrier(MPI_COMM_WORLD);
   }
diff --git a/utils/Analysis/psphisto.cc b/utils/Analysis/psphisto.cc
index 01a57c70f..cf04001b7 100644
--- a/utils/Analysis/psphisto.cc
+++ b/utils/Analysis/psphisto.cc
@@ -34,8 +34,9 @@
 #include <iomanip>
 #include <fstream>
 #include <sstream>
-#include <cmath>
+#include <memory>
 #include <string>
+#include <cmath>
 
                                 // System libs
 #include <sys/time.h>
@@ -177,9 +178,9 @@ main(int argc, char **argv)
       char *c = const_cast<char*>(files[n].c_str());
       MPI_Bcast(c, sz+1, MPI_CHAR, 0, MPI_COMM_WORLD);
     } else {
-      char l[sz+1];
-      MPI_Bcast(&l[0], sz+1, MPI_CHAR, 0, MPI_COMM_WORLD);
-      files.push_back(l);
+      auto l = std::make_unique<char[]>(sz+1);
+      MPI_Bcast(l.get(), sz+1, MPI_CHAR, 0, MPI_COMM_WORLD);
+      files.push_back(l.get());
     }
     MPI_Barrier(MPI_COMM_WORLD);
   }
diff --git a/utils/ICs/SphericalSL.cc b/utils/ICs/SphericalSL.cc
index 04171b51c..a4390b5c2 100644
--- a/utils/ICs/SphericalSL.cc
+++ b/utils/ICs/SphericalSL.cc
@@ -245,7 +245,7 @@ void SphericalSL::compute_coefficients_single(vector<Particle> &part)
 
 void SphericalSL::compute_coefficients_thread(vector<Particle>& part)
 {
-  std::thread t[nthrds];
+  std::vector<std::thread> t(nthrds);
  
   // Launch the threads
   for (int id=0; id<nthrds; ++id) {
diff --git a/utils/ICs/gensph.cc b/utils/ICs/gensph.cc
index d921c883d..60c6d2701 100644
--- a/utils/ICs/gensph.cc
+++ b/utils/ICs/gensph.cc
@@ -355,13 +355,13 @@ main(int argc, char **argv)
 #ifdef HAVE_FFTW
     if (SMOOTH>0.0) {
 
-      double a[NUMR], b[NUMR], c[NUMR];
-      fftw_complex A[NUMR], B[NUMR], C[NUMR];
+      std::vector<double> a(NUMR), b(NUMR), c(NUMR);
+      std::vector<fftw_complex> A(NUMR), B(NUMR), C(NUMR);
       fftw_plan pa, pb, pinv;
 
-      pa   = fftw_plan_dft_r2c_1d(NUMR, a, A, FFTW_ESTIMATE);
-      pb   = fftw_plan_dft_r2c_1d(NUMR, b, B, FFTW_ESTIMATE);
-      pinv = fftw_plan_dft_c2r_1d(NUMR, C, c, FFTW_ESTIMATE);
+      pa   = fftw_plan_dft_r2c_1d(NUMR, a.data(), A.data(), FFTW_ESTIMATE);
+      pb   = fftw_plan_dft_r2c_1d(NUMR, b.data(), B.data(), FFTW_ESTIMATE);
+      pinv = fftw_plan_dft_c2r_1d(NUMR, C.data(), c.data(), FFTW_ESTIMATE);
 
       double xmin=rmin, xmax=rmax;
 
@@ -370,7 +370,7 @@ main(int argc, char **argv)
       double scale = dk / NUMR;
       int indx;
 
-      ofstream tin, tout;
+      std::ofstream tin, tout;
 
       if (myid==0) {
 	tin.open("ebar_fft.input");
@@ -414,7 +414,7 @@ main(int argc, char **argv)
       // ...
 
       double mbar=-1.0;
-      vector<double> rr(NUMR), mm(NUMR);
+      std::vector<double> rr(NUMR), mm(NUMR);
       double fac;
 
       for (int i=0; i<NUMR; i++) {
diff --git a/utils/PhaseSpace/pspmono.cc b/utils/PhaseSpace/pspmono.cc
index 417960df2..9ce2269b7 100644
--- a/utils/PhaseSpace/pspmono.cc
+++ b/utils/PhaseSpace/pspmono.cc
@@ -185,15 +185,15 @@ main(int argc, char **argv)
 
       double val;
       static int buf_size = 1024;
-      char buf[buf_size];
+      auto buf = std::make_unique<char[]>(buf_size);
       
       std::vector<double> or_time, or_c[3];
 
       while (!in.eof()) {
-	in.getline(buf, buf_size);
+	in.getline(buf.get(), buf_size);
 	if (in.eof()) break;
 	
-	istringstream ins(buf);
+	istringstream ins(buf.get());
 	ins >> val;
 	or_time.push_back(val);
 	for (int k=0; k<5; k++) ins >> val;

From 755bd905688070ab62996a1d16b14a7852252028 Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Tue, 15 Oct 2024 21:06:32 -0400
Subject: [PATCH 49/52] Update for 24.04: remove pip install in favor of apt
 install of numpy

---
 .github/workflows/build.yml | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
index b81c922b4..60deed1b0 100644
--- a/.github/workflows/build.yml
+++ b/.github/workflows/build.yml
@@ -34,8 +34,7 @@ jobs:
         if: runner.os == 'Linux'
         run: |
           sudo apt-get update
-          sudo apt-get install -y build-essential libeigen3-dev libfftw3-dev libhdf5-dev libopenmpi-dev libomp-dev
-          sudo pip install numpy
+          sudo apt-get install -y build-essential libeigen3-dev libfftw3-dev libhdf5-dev libopenmpi-dev libomp-dev python3-numpy
 
 
       - name: Setup submodule and build

From e7122075e8c404319f582f44e60240f25b4d356d Mon Sep 17 00:00:00 2001
From: michael-petersen <petersen.michael.s@gmail.com>
Date: Wed, 16 Oct 2024 10:38:29 +0100
Subject: [PATCH 50/52] custom MacOS fix for fftw_complex destructor problems

---
 utils/ICs/gensph.cc | 25 ++++++++++++++++++-------
 1 file changed, 18 insertions(+), 7 deletions(-)

diff --git a/utils/ICs/gensph.cc b/utils/ICs/gensph.cc
index 60c6d2701..15f583584 100644
--- a/utils/ICs/gensph.cc
+++ b/utils/ICs/gensph.cc
@@ -356,12 +356,17 @@ main(int argc, char **argv)
     if (SMOOTH>0.0) {
 
       std::vector<double> a(NUMR), b(NUMR), c(NUMR);
-      std::vector<fftw_complex> A(NUMR), B(NUMR), C(NUMR);
+
+      // Custom fix for MacOS compile: declare explicit vectors for complex numbers
+      // see handling below
+      std::vector<double> A(2 * NUMR), B(2 * NUMR), C(2 * NUMR);
+
       fftw_plan pa, pb, pinv;
 
-      pa   = fftw_plan_dft_r2c_1d(NUMR, a.data(), A.data(), FFTW_ESTIMATE);
-      pb   = fftw_plan_dft_r2c_1d(NUMR, b.data(), B.data(), FFTW_ESTIMATE);
-      pinv = fftw_plan_dft_c2r_1d(NUMR, C.data(), c.data(), FFTW_ESTIMATE);
+      // more MacOS fixes with explicit casts
+      pa   = fftw_plan_dft_r2c_1d(NUMR, a.data(), reinterpret_cast<fftw_complex*>(A.data()), FFTW_ESTIMATE);
+      pb   = fftw_plan_dft_r2c_1d(NUMR, b.data(), reinterpret_cast<fftw_complex*>(B.data()), FFTW_ESTIMATE);
+      pinv = fftw_plan_dft_c2r_1d(NUMR, reinterpret_cast<fftw_complex*>(C.data()), c.data(), FFTW_ESTIMATE);
 
       double xmin=rmin, xmax=rmax;
 
@@ -402,9 +407,15 @@ main(int argc, char **argv)
       fftw_execute(pa);
       fftw_execute(pb);
 
-      for (int i=0; i<NUMR; i++) {
-	C[i][0] = (A[i][0] * B[i][0] - A[i][1] * B[i][1]) * scale;
-	C[i][1] = (A[i][0] * B[i][1] + A[i][1] * B[i][0]) * scale;
+      // Perform the complex multiplication and scaling: explicit for MacOS compile
+      for (int i = 0; i < NUMR; i++) {
+          double A_real = A[2 * i];
+          double A_imag = A[2 * i + 1];
+          double B_real = B[2 * i];
+          double B_imag = B[2 * i + 1];
+
+          C[2 * i]     = (A_real * B_real - A_imag * B_imag) * scale;  // Real part of C[i]
+          C[2 * i + 1] = (A_real * B_imag + A_imag * B_real) * scale;  // Imaginary part of C[i]
       }
 
       // inverse transform to get c, the convolution of a and b;

From 693c1f90df5c3392c5e96c35754503a8e9b7a461 Mon Sep 17 00:00:00 2001
From: michael-petersen <petersen.michael.s@gmail.com>
Date: Wed, 16 Oct 2024 11:01:28 +0100
Subject: [PATCH 51/52] require FFTW

---
 CMakeLists.txt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index f61876725..db6280fbd 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -62,7 +62,7 @@ set(CMAKE_FIND_PACKAGE_SORT_DIRECTION DEC)
 
 find_package(MPI REQUIRED COMPONENTS C CXX)
 find_package(OpenMP)
-find_package(FFTW)
+find_package(FFTW REQUIRED)
 find_package(HDF5 COMPONENTS C CXX HL REQUIRED)
 find_package(TIRPC)	       # Check for alternative Sun rpc support
 find_package(Eigen3 REQUIRED)

From 455c03ffbd9f25dff10a7bb5e62953a5c61703c2 Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Wed, 16 Oct 2024 09:47:19 -0400
Subject: [PATCH 52/52] Use reinterpret cast from std::complex to handle FFTW
 complex types

---
 exputil/TransformFFT.cc | 31 +++++++++++++++++--------------
 include/TransformFFT.H  |  7 ++++---
 utils/ICs/gensph.cc     | 20 ++++++++++++++------
 3 files changed, 35 insertions(+), 23 deletions(-)

diff --git a/exputil/TransformFFT.cc b/exputil/TransformFFT.cc
index 08bfdac56..8b0bb54a7 100644
--- a/exputil/TransformFFT.cc
+++ b/exputil/TransformFFT.cc
@@ -12,9 +12,11 @@ TransformFFT::TransformFFT(double DR, std::vector<double>& Y)
   dk = 2.0*M_PI/dr/N;
 
   in = Y;
-  out = new fftw_complex [N/2+1];
+  out.resize(N/2+1);
   
-  p = fftw_plan_dft_r2c_1d(N, in.data(), out, FFTW_ESTIMATE);
+  p = fftw_plan_dft_r2c_1d(N, in.data(),
+			   reinterpret_cast<fftw_complex*>(out.data()),
+			   FFTW_ESTIMATE);
 
   fftw_execute(p);
 }
@@ -30,9 +32,11 @@ TransformFFT::TransformFFT(double DR, Eigen::VectorXd& Y)
   in.resize(N);
   for (int j=0; j<N; j++) in[j] = Y[j];
 
-  out = new fftw_complex [N/2+1];
+  out.resize(N/2+1);
   
-  p = fftw_plan_dft_r2c_1d(N, in.data(), out, FFTW_ESTIMATE);
+  p = fftw_plan_dft_r2c_1d(N, in.data(),
+			   reinterpret_cast<fftw_complex*>(out.data()),
+			   FFTW_ESTIMATE);
 
   
   fftw_execute(p);
@@ -41,7 +45,6 @@ TransformFFT::TransformFFT(double DR, Eigen::VectorXd& Y)
 
 TransformFFT::~TransformFFT()
 {
-  delete [] out;
   fftw_destroy_plan(p);
 }
 
@@ -53,15 +56,15 @@ void TransformFFT::Power(Eigen::VectorXd& F, Eigen::VectorXd& P)
   double d2 = dr*dr;
 
   F(0) = 0.0;
-  P(0) = d2 * (out[0][0]*out[0][0] + out[0][1]*out[0][1]);
+  P(0) = d2 * std::norm(out[0]);
 
   for (int j=1; j<N/2; j++) {
     F(j) = dk * j;
-    P(j) = d2 * (out[j][0]*out[j][0] + out[j][1]*out[j][1]) * 2.0;
+    P(j) = d2 * std::norm(out[j]) * 2.0;
   }
   
   F(N/2) = dk * N/2;
-  P(N/2) = d2 * (out[N/2][0]*out[N/2][0] + out[N/2][1]*out[N/2][1]);
+  P(N/2) = d2 * std::norm(out[N/2]);
 }
 
 void TransformFFT::Power(std::vector<double>& F, std::vector<double>& P)
@@ -72,15 +75,15 @@ void TransformFFT::Power(std::vector<double>& F, std::vector<double>& P)
   double d2 = dr*dr;
 
   F[0] = 0.0;
-  P[0] = d2 * ( out[0][0]*out[0][0] + out[0][1]*out[0][1] );
+  P[0] = d2 * std::norm(out[0]);
 
   for (int j=1; j<N/2; j++) {
     F[j] = dk * j;
-    P[j] = d2 * out[j][0]*out[j][0] + out[j][1]*out[j][1] * 2.0;
+    P[j] = d2 * std::norm(out[j]) * 2.0;
   }
   
   F[N/2] = dk * N/2;
-  P[N/2] = d2 * (out[N/2][0]*out[N/2][0] + out[N/2][1]*out[N/2][1]);
+  P[N/2] = d2 * std::norm(out[N/2]);
 }
 
 
@@ -91,7 +94,7 @@ void TransformFFT::Inverse(Eigen::VectorXd& F, Eigen::VectorXcd& W)
 
   for (int j=0; j<N/2+1; j++) {
     F[j] = dk*j;
-    W[j] = std::complex<double>(out[j][0], out[j][1]) * dr;
+    W[j] = out[j] * dr;
   }
   
 }
@@ -105,8 +108,8 @@ void TransformFFT::Inverse(std::vector<double>& F,
 
   for (int j=0; j<N/2+1; j++) {
     F[j] = dk*j;
-    Wr[j] = out[j][0] * dr;
-    Wi[j] = out[j][1] * dr;
+    Wr[j] = std::real(out[j]) * dr;
+    Wi[j] = std::imag(out[j]) * dr;
   }
 }
 
diff --git a/include/TransformFFT.H b/include/TransformFFT.H
index ad2e8d4f4..28086490a 100644
--- a/include/TransformFFT.H
+++ b/include/TransformFFT.H
@@ -2,11 +2,12 @@
 #ifndef _TransformFFT_h
 #define _TransformFFT_h
 
+#include <complex>
+#include <vector>
+
 #include <Eigen/Dense>
 #include <fftw3.h>
 
-#include <complex>
-#include <vector>
 
 class TransformFFT
 {
@@ -15,7 +16,7 @@ class TransformFFT
 
   fftw_plan p;
   std::vector<double> in;
-  fftw_complex* out;
+  std::vector<std::complex<double>> out;
 
  public:
   
diff --git a/utils/ICs/gensph.cc b/utils/ICs/gensph.cc
index 60c6d2701..b42b397bc 100644
--- a/utils/ICs/gensph.cc
+++ b/utils/ICs/gensph.cc
@@ -57,6 +57,7 @@ using namespace __EXP__;
 #include <config_exp.h>
 
 #ifdef HAVE_FFTW
+#include <complex>
 #include <fftw3.h>
 #endif
 
@@ -356,12 +357,20 @@ main(int argc, char **argv)
     if (SMOOTH>0.0) {
 
       std::vector<double> a(NUMR), b(NUMR), c(NUMR);
-      std::vector<fftw_complex> A(NUMR), B(NUMR), C(NUMR);
+      std::vector<std::complex<double>> A(NUMR), B(NUMR), C(NUMR);
       fftw_plan pa, pb, pinv;
 
-      pa   = fftw_plan_dft_r2c_1d(NUMR, a.data(), A.data(), FFTW_ESTIMATE);
-      pb   = fftw_plan_dft_r2c_1d(NUMR, b.data(), B.data(), FFTW_ESTIMATE);
-      pinv = fftw_plan_dft_c2r_1d(NUMR, C.data(), c.data(), FFTW_ESTIMATE);
+      pa   = fftw_plan_dft_r2c_1d(NUMR, a.data(),
+				  reinterpret_cast<fftw_complex*>(A.data()),
+				  FFTW_ESTIMATE);
+
+      pb   = fftw_plan_dft_r2c_1d(NUMR, b.data(),
+				  reinterpret_cast<fftw_complex*>(B.data()),
+				  FFTW_ESTIMATE);
+
+      pinv = fftw_plan_dft_c2r_1d(NUMR,
+				  reinterpret_cast<fftw_complex*>(C.data()),
+				  c.data(), FFTW_ESTIMATE);
 
       double xmin=rmin, xmax=rmax;
 
@@ -403,8 +412,7 @@ main(int argc, char **argv)
       fftw_execute(pb);
 
       for (int i=0; i<NUMR; i++) {
-	C[i][0] = (A[i][0] * B[i][0] - A[i][1] * B[i][1]) * scale;
-	C[i][1] = (A[i][0] * B[i][1] + A[i][1] * B[i][0]) * scale;
+	C[i] = A[i] * B[i] * scale;
       }
 
       // inverse transform to get c, the convolution of a and b;