diff --git a/CMakeLists.txt b/CMakeLists.txt
index 2c95fb0e8..6d9c87c39 100644
--- a/CMakeLists.txt
+++ b/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")
diff --git a/expui/BiorthBess.cc b/expui/BiorthBess.cc
index ce35f0124..0b16c8ce0 100644
--- a/expui/BiorthBess.cc
+++ b/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)
diff --git a/exputil/EmpCyl2d.cc b/exputil/EmpCyl2d.cc
index 26a5488cc..1eaf0635e 100644
--- a/exputil/EmpCyl2d.cc
+++ b/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
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/pybind11 b/extern/pybind11
index 914c06fb2..58c382a8e 160000
--- a/extern/pybind11
+++ b/extern/pybind11
@@ -1 +1 @@
-Subproject commit 914c06fb252b6cc3727d0eedab6736e88a3fcb01
+Subproject commit 58c382a8e3d7081364d2f5c62e7f429f0412743b
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/pyEXP/TensorToArray.H b/pyEXP/TensorToArray.H
index 765148eac..d2ba60209 100644
--- a/pyEXP/TensorToArray.H
+++ b/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)