Update basis normalization

CMakeLists.txt

@@ -33,6 +33,9 @@ endif()
 # Required compiler features
 add_compile_options(-D_REENTRANT)
 
+# Bake in library paths (esp. useful for HPC sites with modules)
+set(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE)
+
 # Check and enforce that we are a git repository.  Necessary for
 # submodules to work correctly.
 if(EXISTS "${PROJECT_SOURCE_DIR}/.git")

expui/BiorthBess.cc

@@ -5,7 +5,7 @@
 #include <EXPException.H>
 #include <interp.H>
 #include <gaussQ.H>
-
+#include <cassert>
 
 BiorthBess::BiorthBess(double rmax, int lmax, int nmax, int RNUM) :
   rmax(rmax), lmax(lmax), nmax(nmax), RNUM(RNUM)

exputil/EmpCyl2d.cc

@@ -873,9 +873,15 @@ void EmpCyl2d::create_tables()
       if (m==0) xgrid[i] = r;
 
       for (int n=0; n<nmaxfid; n++) {
+<<<<<<< Updated upstream
 	pot(n) = basis->potl(m, n, r) / sqrt(basis->norm(n, m)*0.5);
 	den(n) = basis->dens(m, n, r) / sqrt(basis->norm(n, m)*2.0);
 	dph(n) = basis->dpot(m, n, r) / sqrt(basis->norm(n, m)*0.5);
+=======
+	pot(n) = basis->potl(m, n, r) / sqrt(basis->norm(n, m)*0.5/M_PI);
+	den(n) = basis->dens(m, n, r) / sqrt(basis->norm(n, m)*2.0*M_PI);
+	dph(n) = basis->dpot(m, n, r) / sqrt(basis->norm(n, m)*0.5/M_PI);
+>>>>>>> Stashed changes
       }
 
       pot = U.transpose() * pot;
@@ -1203,7 +1209,11 @@ double EmpCyl2d::get_potl(double r, int M, int N)
   checkMN(M, N, "get_potl");
 
   if (basis_test) {
+<<<<<<< Updated upstream
     return basis->potl(M, N, r)/sqrt(basis->norm(N, M)*0.5);
+=======
+    return basis->potl(M, N, r)/sqrt(basis->norm(N, M)*0.5/M_PI);
+>>>>>>> Stashed changes
   }
 
   int lo, hi;
@@ -1218,7 +1228,11 @@ double EmpCyl2d::get_dens(double r, int M, int N)
   checkMN(M, N, "get_dens");
 
   if (basis_test) {
+<<<<<<< Updated upstream
     return basis->dens(M, N, r)/sqrt(basis->norm(N, M)*2.0);
+=======
+    return basis->dens(M, N, r)/sqrt(basis->norm(N, M)*2.0*M_PI);
+>>>>>>> Stashed changes
   }
 
   int lo, hi;
@@ -1233,7 +1247,11 @@ double EmpCyl2d::get_dpot(double r, int M, int N)
   checkMN(M, N, "get_dpot");
 
   if (basis_test) {
+<<<<<<< Updated upstream
     return basis->dpot(M, N, r)/sqrt(basis->norm(N, M)*0.5);
+=======
+    return basis->dpot(M, N, r)/sqrt(basis->norm(N, M)*2.0*M_PI);
+>>>>>>> Stashed changes
   }
 
   int lo, hi;
@@ -1307,7 +1325,11 @@ void EmpCyl2d::checkCoefs()
 
       for (int j=0; j<nmax; j++) {
 	coefs(j) += fac * disk->dens(rr) * get_potl(rr, 0, j);
+<<<<<<< Updated upstream
 	coef0(j) += fac * disk->dens(rr) * basis->potl(0, j, rr) / sqrt(basis->norm(j, 0)*0.5);
+=======
+	coef0(j) += fac * disk->dens(rr) * basis->potl(0, j, rr) / sqrt(basis->norm(j, 0)*2.0*M_PI);
+>>>>>>> Stashed changes
       }
   }
 
@@ -1331,8 +1353,13 @@ void EmpCyl2d::checkCoefs()
     for (int j=0; j<nmax; j++) {
       uu += coefs(j) * get_potl(rr, 0, j);
       vv += coefs(j) * get_dens(rr, 0, j);
+<<<<<<< Updated upstream
       yy += coef0(j) * basis->potl(0, j, rr) / sqrt(basis->norm(j, 0)*0.5);
       zz += coef0(j) * basis->dens(0, j, rr) / sqrt(basis->norm(j, 0)*2.0);
+=======
+      yy += coef0(j) * basis->potl(0, j, rr) / sqrt(basis->norm(j, 0)*2.0*M_PI);
+      zz += coef0(j) * basis->dens(0, j, rr) / sqrt(basis->norm(j, 0)*0.5/M_PI);
+>>>>>>> Stashed changes
     }
 
     std::cout << std::setw(16) << rr

pyEXP/TensorToArray.H

@@ -27,6 +27,38 @@ py::array_t<T> make_ndarray3(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)