diff --git a/expui/BasisFactory.H b/expui/BasisFactory.H
index e270fae1e..6f3351497 100644
--- a/expui/BasisFactory.H
+++ b/expui/BasisFactory.H
@@ -196,15 +196,25 @@ namespace BasisClasses
     operator()(double x1, double x2, double x3,
 	       const Coord ctype=Coord::Spherical);
     
-    //! Evaluate fields at a point
+    //! Evaluate fields at a point in the expansion-origin frame
     virtual std::vector<double> getFields(double x, double y, double z);
+
+    //! Evaluate fields at a point in the coefficient-origin frame
+    virtual std::vector<double> getFieldsOrigin(double x, double y, double z);
     
-    //! Evaluate fields at a point for all coefficients sets
+    //! Evaluate fields at a point for all coefficients sets in the
+    //! expansion-origin frame
     virtual std::tuple<std::map<std::string, Eigen::VectorXd>,
 		       Eigen::VectorXd> getFieldsCoefs
     (double x, double y, double z, std::shared_ptr<CoefClasses::Coefs> coefs);
+
+    //! Evaluate fields at a point for all coefficients sets in the
+    //! coefficient-origin frame
+    virtual std::tuple<std::map<std::string, Eigen::VectorXd>,
+		       Eigen::VectorXd> getFieldsCoefsOrigin
+    (double x, double y, double z, std::shared_ptr<CoefClasses::Coefs> coefs);
     
-    //! Evaluate fields at a point, and provide field lables
+    //! Evaluate fields at a point, and provide field labels
     virtual std::tuple<std::vector<double>, std::vector<std::string>>
     evaluate(double x, double y, double z)
     { return {getFields(x, y, z), getFieldLabels(coordinates)}; }
diff --git a/expui/BasisFactory.cc b/expui/BasisFactory.cc
index aba544cc3..0830656b6 100644
--- a/expui/BasisFactory.cc
+++ b/expui/BasisFactory.cc
@@ -229,7 +229,13 @@ namespace BasisClasses
   }
     
   std::vector<double> Basis::getFields(double x, double y, double z)
+  { return crt_eval(x, y, z); }
+
+  std::vector<double> Basis::getFieldsOrigin(double x, double y, double z)
   {
+    x -= coefctr(0);
+    y -= coefctr(1);
+    z -= coefctr(2);
     return crt_eval(x, y, z);
   }
     
@@ -237,7 +243,7 @@ namespace BasisClasses
   Basis::getFieldsCoefs
   (double x, double y, double z, std::shared_ptr<CoefClasses::Coefs> coefs)
   {
-    // Python dictonary for return
+    // Python dictionary for return
     std::map<std::string, Eigen::VectorXd> ret;
 
     // Times for the coefficients
@@ -249,9 +255,46 @@ namespace BasisClasses
 
     // Make the return dictionary of arrays
     for (int i=0; i<times.size(); i++) {
+      // Load the coefficients for the current time
       set_coefs(coefs->getCoefStruct(times[i]));
+
       // The field evaluation
       auto v = crt_eval(x, y, z); 
+
+      // Pack the fields into the dictionary
+      for (int j=0; j<fields.size(); j++) ret[fields[j]][i] = v[j];
+    }
+
+    // An attempt at an efficient return type for the time array
+    Eigen::VectorXd T =
+      Eigen::Map<Eigen::VectorXd, Eigen::Unaligned>(times.data(), times.size());
+
+    // Return the dictionary and the time array
+    return {ret, T};
+  }
+
+  std::tuple<std::map<std::string, Eigen::VectorXd>, Eigen::VectorXd>
+  Basis::getFieldsCoefsOrigin
+  (double x, double y, double z, std::shared_ptr<CoefClasses::Coefs> coefs)
+  {
+    // Python dictionary for return
+    std::map<std::string, Eigen::VectorXd> ret;
+
+    // Times for the coefficients
+    auto times  = coefs->Times();
+
+    // Initialize the dictionary/map
+    auto fields = getFieldLabels(coordinates);
+    for (auto s : fields) ret[s].resize(times.size());
+
+    // Make the return dictionary of arrays
+    for (int i=0; i<times.size(); i++) {
+      // Load the coefficients for the current time
+      set_coefs(coefs->getCoefStruct(times[i]));
+
+      // The field evaluation
+      auto v = getFieldsOrigin(x, y, z);
+
       // Pack the fields into the dictionary
       for (int j=0; j<fields.size(); j++) ret[fields[j]][i] = v[j];
     }
diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc
index b7723992f..e9c4e4536 100644
--- a/expui/BiorthBasis.cc
+++ b/expui/BiorthBasis.cc
@@ -818,6 +818,11 @@ namespace BasisClasses
   void Spherical::computeAccel(double x, double y, double z,
 			       Eigen::Ref<Eigen::Vector3d> acc)
   {
+    // Shift to center
+    x -= coefctr(0);
+    y -= coefctr(1);
+    z -= coefctr(2);
+
     // Get polar coordinates
     double R2    = x*x + y*y;
     double r2    = R2  + z*z;
@@ -1792,8 +1797,12 @@ namespace BasisClasses
     
     tdens = sl->accumulated_dens_eval(R, z, phi, tdens0);
 
-    double tpotx = tpotR*x/R - tpotp*y/R ;
-    double tpoty = tpotR*y/R + tpotp*x/R ;
+    double tpotx = 0.0;
+    double tpoty = 0.0;
+    if (R > 0.0) {
+      tpotx = tpotR*x/R - tpotp*y/R;
+      tpoty = tpotR*y/R + tpotp*x/R;
+    }
 
     return
       {tdens0, tdens - tdens0, tdens,
@@ -1804,6 +1813,11 @@ namespace BasisClasses
   void Cylindrical::computeAccel(double x, double y, double z,
 				 Eigen::Ref<Eigen::Vector3d> acc)
   {
+    // Shift to center
+    x -= coefctr(0);
+    y -= coefctr(1);
+    z -= coefctr(2);
+
     double R = sqrt(x*x + y*y);
     double phi = atan2(y, x);
 
@@ -1813,8 +1827,12 @@ namespace BasisClasses
     
     tdens = sl->accumulated_dens_eval(R, z, phi, tdens0);
 
-    double tpotx = tpotR*x/R - tpotp*y/R ;
-    double tpoty = tpotR*y/R + tpotp*x/R ;
+    double tpotx = 0.0;
+    double tpoty = 0.0;
+    if (R > 0.0) {
+      tpotx = tpotR*x/R - tpotp*y/R;
+      tpoty = tpotR*y/R + tpotp*x/R;
+    }
 
     // Apply G to forces on return
     acc << tpotx*G, tpoty*G, tpotz*G;
@@ -2481,6 +2499,11 @@ namespace BasisClasses
   void FlatDisk::computeAccel(double x, double y, double z, 
 			      Eigen::Ref<Eigen::Vector3d> acc)
   {
+    // Shift to center
+    x -= coefctr(0);
+    y -= coefctr(1);
+    z -= coefctr(2);
+
     // Get thread id
     int tid = omp_get_thread_num();
 
@@ -2565,8 +2588,12 @@ namespace BasisClasses
     zpot *= -G;
     ppot *= -G;
 
-    double potx = rpot*x/R - ppot*y/R;
-    double poty = rpot*y/R + ppot*x/R;
+    double potx = 0.0;
+    double poty = 0.0;
+    if (R > 0.0) {
+      potx = rpot*x/R - ppot*y/R;
+      poty = rpot*y/R + ppot*x/R;
+    }
 
     acc << potx, poty, zpot;
   }
@@ -3270,6 +3297,11 @@ namespace BasisClasses
   void CBDisk::computeAccel(double x, double y, double z,
 			    Eigen::Ref<Eigen::Vector3d> acc)
   {
+    // Shift to center
+    x -= coefctr(0);
+    y -= coefctr(1);
+    z -= coefctr(2);
+
     // Get thread id
     int tid = omp_get_thread_num();
 
@@ -3330,8 +3362,12 @@ namespace BasisClasses
     rpot *= -G;
     ppot *= -G;
 
-    double potx = rpot*x/R - ppot*y/R;
-    double poty = rpot*y/R + ppot*x/R;
+    double potx = 0.0;
+    double poty = 0.0;
+    if (R > 0.0) {
+      potx = rpot*x/R - ppot*y/R;
+      poty = rpot*y/R + ppot*x/R;
+    }
 
     acc << potx, poty, zpot;
   }
@@ -3756,6 +3792,11 @@ namespace BasisClasses
   void Slab::computeAccel(double x, double y, double z,
 			  Eigen::Ref<Eigen::Vector3d> acc)
   {
+    // Shift to center
+    x -= coefctr(0);
+    y -= coefctr(1);
+    z -= coefctr(2);
+
     // Loop indices
     //
     int ix, iy, iz;
@@ -4328,6 +4369,11 @@ namespace BasisClasses
   void Cube::computeAccel(double x, double y, double z,
 			  Eigen::Ref<Eigen::Vector3d> acc)
   {
+    // Shift to center
+    x -= coefctr(0);
+    y -= coefctr(1);
+    z -= coefctr(2);
+
     // Get thread id
     int tid = omp_get_thread_num();
 
diff --git a/expui/CMakeLists.txt b/expui/CMakeLists.txt
index f25d0d7b7..15d7e4509 100644
--- a/expui/CMakeLists.txt
+++ b/expui/CMakeLists.txt
@@ -97,11 +97,15 @@ target_sources(expui PUBLIC FILE_SET HEADERS
 	${CMAKE_CURRENT_SOURCE_DIR}/../include/DiskWithHalo.H
 	${CMAKE_CURRENT_SOURCE_DIR}/../include/EmpCyl2d.H
 	${CMAKE_CURRENT_SOURCE_DIR}/../include/EXPmath.H
+	# ${CMAKE_CURRENT_SOURCE_DIR}/../include/EXPversion.H
 	${CMAKE_CURRENT_SOURCE_DIR}/../include/libvars.H
 	${CMAKE_CURRENT_SOURCE_DIR}/../include/Timer.H
 	${CMAKE_CURRENT_SOURCE_DIR}/../include/coef.H
 	${CMAKE_CURRENT_SOURCE_DIR}/../include/Covariance.H
 	${CMAKE_CURRENT_SOURCE_DIR}/../include/ExpDeproj.H
+	${CMAKE_CURRENT_SOURCE_DIR}/../include/cudaUtil.cuH
+	${CMAKE_CURRENT_SOURCE_DIR}/../include/cudaParticle.cuH
+	${CMAKE_CURRENT_SOURCE_DIR}/../include/cudaMappingConstants.cuH	
 )
 
 install(TARGETS expui FILE_SET HEADERS DESTINATION include/EXP)
diff --git a/include/BiorthCube.H b/include/BiorthCube.H
index a89ebb87d..808ed7f49 100644
--- a/include/BiorthCube.H
+++ b/include/BiorthCube.H
@@ -17,8 +17,9 @@
 #include <yaml-cpp/yaml.h>
 
 #if HAVE_LIBCUDA==1
-#include <cudaUtil.cuH>
-#include <cudaMappingConstants.cuH>
+#include "cudaUtil.cuH"
+#include "cudaParticle.cuH"
+#include "cudaMappingConstants.cuH"
 #endif
 
 // For reading and writing cache file
@@ -28,11 +29,6 @@
 #include <highfive/H5DataSpace.hpp>
 #include <highfive/H5Attribute.hpp>
 
-#if HAVE_LIBCUDA==1
-#include <cudaParticle.cuH>
-#include <cudaMappingConstants.cuH>
-#endif
-
 #include "EmpCyl2d.H"
 
 //!! BiorthCube grid class
diff --git a/include/BiorthCyl.H b/include/BiorthCyl.H
index 687433586..226120441 100644
--- a/include/BiorthCyl.H
+++ b/include/BiorthCyl.H
@@ -17,8 +17,9 @@
 #include <yaml-cpp/yaml.h>
 
 #if HAVE_LIBCUDA==1
-#include <cudaUtil.cuH>
-#include <cudaMappingConstants.cuH>
+#include "cudaUtil.cuH"
+#include "cudaParticle.cuH"
+#include "cudaMappingConstants.cuH"
 #endif
 
 // For reading and writing cache file
@@ -28,11 +29,6 @@
 #include <highfive/H5DataSpace.hpp>
 #include <highfive/H5Attribute.hpp>
 
-#if HAVE_LIBCUDA==1
-#include <cudaParticle.cuH>
-#include <cudaMappingConstants.cuH>
-#endif
-
 #include "EmpCyl2d.H"
 
 //!! BiorthCyl grid class
diff --git a/include/EmpCylSL.H b/include/EmpCylSL.H
index 9811a0059..1bd993e10 100644
--- a/include/EmpCylSL.H
+++ b/include/EmpCylSL.H
@@ -24,8 +24,8 @@
 #include "coef.H"
 
 #if HAVE_LIBCUDA==1
-#include <cudaParticle.cuH>
-#include <cudaMappingConstants.cuH>
+#include "cudaParticle.cuH"
+#include "cudaMappingConstants.cuH"
 #endif
 
 #include "libvars.H"
diff --git a/include/SLGridMP2.H b/include/SLGridMP2.H
index da209403c..dd46e28cf 100644
--- a/include/SLGridMP2.H
+++ b/include/SLGridMP2.H
@@ -19,8 +19,8 @@
 using namespace __EXP__;
 
 #if HAVE_LIBCUDA==1
-#include <cudaUtil.cuH>
-#include <cudaMappingConstants.cuH>
+#include "cudaUtil.cuH"
+#include "cudaMappingConstants.cuH"
 #endif
 
 
diff --git a/include/cudaMappingConstants.cuH b/include/cudaMappingConstants.cuH
index 067d299b6..3e845436d 100644
--- a/include/cudaMappingConstants.cuH
+++ b/include/cudaMappingConstants.cuH
@@ -1,9 +1,7 @@
-// -*- C++ -*-
-
 #ifndef CUDA_MAPPING_CONSTANTS_H
 #define CUDA_MAPPING_CONSTANTS_H
 
-#include <cudaUtil.cuH>
+#include "cudaUtil.cuH"
 
 struct cudaMappingConstants
 {
@@ -23,3 +21,5 @@ struct cudaMappingConstants
 };
 
 #endif
+
+// -*- C++ -*-
diff --git a/include/cudaParticle.cuH b/include/cudaParticle.cuH
index 5147b01dd..ccc2102b1 100644
--- a/include/cudaParticle.cuH
+++ b/include/cudaParticle.cuH
@@ -1,5 +1,3 @@
-// -*- C++ -*-
-
 #ifndef PARTICLE_CUH
 #define PARTICLE_CUH
 
@@ -16,11 +14,9 @@
 #include <thrust/host_vector.h>
 #include <thrust/sort.h>
 
-#include <config_exp.h>
-
-#include <cudaUtil.cuH>
-
-#include <Particle.H>
+#include "config_exp.h"
+#include "cudaUtil.cuH"
+#include "Particle.H"
 
 //! Simplified particle structure for use in CUDA kernel code
 struct cudaParticle
@@ -125,3 +121,5 @@ struct cuPartToChange
 };
 
 #endif
+
+// -*- C++ -*-
diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc
index 1429c6249..6da222c1a 100644
--- a/pyEXP/BasisWrappers.cc
+++ b/pyEXP/BasisWrappers.cc
@@ -562,7 +562,11 @@ void BasisFactoryClasses(py::module &m)
     using Spherical::Spherical;
 
     std::vector<double> getFields(double x, double y, double z) override {
-      PYBIND11_OVERRIDE(std::vector<double>, Spherical, getFields, x, y, z);
+      py::function override = py::get_override(static_cast<const Spherical*>(this), "getFields");
+      if (override) {
+        return py::cast<std::vector<double>>(override(x, y, z));
+      }
+      return Spherical::getFields(x, y, z);
     }
 
     void accumulate(double x, double y, double z, double mass, unsigned long int indx) override {
@@ -1522,7 +1526,7 @@ void BasisFactoryClasses(py::module &m)
          potential evaluations are separated into full, axisymmetric and
          non-axisymmetric contributions.
 
-         You can get the fields labels by using the __call__ method of the
+         You can get the field labels by using the __call__ method of the
          basis object.  This is equivalent to a tuple of the getFields()
          output with a list of field labels.
 
@@ -1541,13 +1545,38 @@ void BasisFactoryClasses(py::module &m)
 
          See also
          --------
+         getFieldsOrigin: get fields in coefficient-origin frame
          getFieldsCoefs : get fields for each coefficient set
          __call__       : same as getFields() but provides field labels in a tuple
          )",
 	 py::arg("x"), py::arg("y"), py::arg("z"))
+    .def("getFieldsOrigin", &BasisClasses::BiorthBasis::getFieldsOrigin,
+	 R"(
+         Return the field evaluations for a given Cartesian position in
+         the frame defined by the current coefficients.
+
+         Parameters
+         ----------
+         x : float
+             x-axis position
+         y : float
+             y-axis position
+         z : float
+             z-axis position
+
+         Returns
+         -------
+         fields: numpy.ndarray
+
+         See also
+         --------
+         getFields      : get fields in expansion-origin frame
+         getFieldsCoefs : get fields for each coefficient set
+         )",
+	 py::arg("x"), py::arg("y"), py::arg("z"))
     .def("getAccel", py::overload_cast<double, double, double>(&BasisClasses::BiorthBasis::getAccel),
 	 R"(
-         Return the acceleration for a given Cartesian position
+         Return the acceleration for a given Cartesian position in the frame defined by the coefficients.
 
          Parameters
          ----------
@@ -1571,7 +1600,7 @@ void BasisFactoryClasses(py::module &m)
 	 py::arg("x"), py::arg("y"), py::arg("z"))
     .def("getAccel", py::overload_cast<Eigen::Ref<const Eigen::VectorXd>, Eigen::Ref<const Eigen::VectorXd>, Eigen::Ref<const Eigen::VectorXd>>(&BasisClasses::BiorthBasis::getAccel),
 	 R"(
-         Return the acceleration for a given Cartesian position
+         Return the acceleration for a given Cartesian position in the frame defined by the coefficients.
 
          Parameters
          ----------
@@ -1595,7 +1624,7 @@ void BasisFactoryClasses(py::module &m)
 	 py::arg("x"), py::arg("y"), py::arg("z"))
     .def("getAccel", py::overload_cast<Eigen::Ref<const RowMatrixXd>>(&BasisClasses::BiorthBasis::getAccel),
 	 R"(
-         Return the acceleration for an array of Cartesian positions
+         Return the acceleration for an array of Cartesian positions in the frame defined by the coefficients.
 
          Parameters
          ----------
@@ -1615,7 +1644,7 @@ void BasisFactoryClasses(py::module &m)
 	 py::arg("pos"))
     .def("getAccelArray", py::overload_cast<Eigen::Ref<const Eigen::VectorXd>, Eigen::Ref<const Eigen::VectorXd>, Eigen::Ref<const Eigen::VectorXd>>(&BasisClasses::BiorthBasis::getAccel),
 	 R"(
-         Return the acceleration for a given Cartesian position
+         Return the acceleration for a given Cartesian position in the frame defined by the coefficients.
 
          Parameters
          ----------
@@ -1665,8 +1694,36 @@ void BasisFactoryClasses(py::module &m)
 
          See also
          --------
-         getFields  : get fields for the currently assigned coefficients
-         __call__   : same getFields() but provides field labels in a tuple
+         getFields        : get fields for the currently assigned coefficients
+         getFieldsOrigin  : get fields in coefficient-origin frame
+         __call__         : same getFields() but provides field labels in a tuple
+         )",
+	 py::arg("x"), py::arg("y"), py::arg("z"), py::arg("coefs"))
+    .def("getFieldsCoefsOrigin", &BasisClasses::BiorthBasis::getFieldsCoefsOrigin,
+	 R"(
+         Return the field evaluations for a given Cartesian position
+         for every frame in a coefficient set in the frame defined by
+         each coefficient structure.
+
+         Parameters
+         ----------
+         x : float
+             x-axis position
+         y : float
+             y-axis position
+         z : float
+             z-axis position
+         coefs: CoefClasses::Coefs
+             the coefficient set
+
+         Returns
+         -------
+         tuple of a dictionary of fields of array values, and an
+             array of evaluation times
+
+         See also
+         --------
+         getFieldsCoefs : get fields in expansion-origin frame
          )",
 	 py::arg("x"), py::arg("y"), py::arg("z"), py::arg("coefs"))
     .def("setFieldType",       &BasisClasses::BiorthBasis::setFieldType,