EXP-code · The9Cat · May 14, 2026 · May 14, 2026 · May 14, 2026 · May 14, 2026
diff --git a/expui/BasisFactory.H b/expui/BasisFactory.H
@@ -197,17 +197,25 @@ namespace BasisClasses
 	       const Coord ctype=Coord::Spherical);
 
     //! Evaluate fields at a point
-    virtual std::vector<double> getFields(double x, double y, double z);
+    //! @note API change in v7.10.3: Added 'origin' parameter (default=true).
+    //!       This is an ABI-breaking change for external C++ code that
+    //!       subclasses Basis or calls these methods via vtable.
+    virtual std::vector<double> getFields(double x, double y, double z,
+					  bool origin=true);
 
     //! Evaluate fields at a point for all coefficients sets
+    //! @note API change in v7.10.3: Added 'origin' parameter (default=true).
+    //!       This is an ABI-breaking change for external C++ code that
+    //!       subclasses Basis or calls these methods via vtable.
     virtual std::tuple<std::map<std::string, Eigen::VectorXd>,
 		       Eigen::VectorXd> getFieldsCoefs
-    (double x, double y, double z, std::shared_ptr<CoefClasses::Coefs> coefs);
+    (double x, double y, double z, std::shared_ptr<CoefClasses::Coefs> coefs,
+     bool origin=true);
 
-    //! 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)}; }
+    { return {getFields(x, y, z, true), getFieldLabels(coordinates)}; }
 
     //! Retrieve the coefficients 
     virtual CoefClasses::CoefStrPtr getCoefficients()

diff --git a/expui/BasisFactory.cc b/expui/BasisFactory.cc
@@ -228,14 +228,21 @@ namespace BasisClasses
     };
   }
 
-  std::vector<double> Basis::getFields(double x, double y, double z)
+  std::vector<double> Basis::getFields(double x, double y, double z,
+				       bool origin)
   {
+    if (not origin) {
+      x -= coefctr(0);
+      y -= coefctr(1);
+      z -= coefctr(2);
+    }
     return crt_eval(x, y, z);
   }
 
   std::tuple<std::map<std::string, Eigen::VectorXd>, Eigen::VectorXd>
   Basis::getFieldsCoefs
-  (double x, double y, double z, std::shared_ptr<CoefClasses::Coefs> coefs)
+  (double x, double y, double z, std::shared_ptr<CoefClasses::Coefs> coefs,
+   bool origin)
   {
     // Python dictonary for return
     std::map<std::string, Eigen::VectorXd> ret;
@@ -249,9 +256,20 @@ 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]));
+
+      double xc = x, yc = y, zc = z;
+      // Apply centering
+      if (not origin) {
+	xc -= coefctr(0);
+	yc -= coefctr(1);
+	zc -= coefctr(2);
+      }
+
       // The field evaluation
-      auto v = crt_eval(x, y, z); 
+      auto v = crt_eval(xc, yc, zc); 
+
       // 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
@@ -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;
@@ -1804,6 +1809,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 +1823,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 +2495,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 +2584,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 +3293,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 +3358,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 +3788,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 +4365,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();
 
@@ -4781,7 +4823,7 @@ namespace BasisClasses
       for (int k=0; k<3; k++) pp(k) = ps(n, k) - ctr(k);
       pp = rot * pp;
 
-      auto v = basis->getFields(pp(0), pp(1), pp(2));
+      auto v = basis->getFields(pp(0), pp(1), pp(2), true);
 
       // First 6 fields are density and potential, followed by acceleration
       for (int k=0; k<3; k++) accel(n, k) += v[6+k] - basis->pseudo(k);

diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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++ -*-