diff --git a/.gitignore b/.gitignore
index 319d66928..3687bf902 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,5 +1,6 @@
 
 build/*
+_codeql_build_dir/
 
 Makefile
 *.lo
@@ -25,3 +26,4 @@ src/exp
 src/user/CylindricalDisk.cc
 src/user/EllipsoidForce.cc
 src/user/SLSphere.cc
+_codeql_build_dir
diff --git a/_codeql_detected_source_root b/_codeql_detected_source_root
new file mode 120000
index 000000000..945c9b46d
--- /dev/null
+++ b/_codeql_detected_source_root
@@ -0,0 +1 @@
+.
\ No newline at end of file
diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H
index f40875985..65cd8aee1 100644
--- a/expui/BiorthBasis.H
+++ b/expui/BiorthBasis.H
@@ -968,7 +968,7 @@ namespace BasisClasses
     int lmaxfid, nmaxfid, mmax, mlim, nmax;
     int ncylodd, ncylnx, ncylny, ncylr, cmap, cmapR, cmapZ, vflag;
     int rnum, pnum, tnum;
-    double rcylmin, rcylmax, acyl, hcyl;
+    double rcylmin, rcylmax, acyl, hcyl, bias;
     bool expcond, logarithmic, density, EVEN_M, sech2 = true;
     
     std::vector<Eigen::MatrixXd> potd, dpot, dpt2, dend;
diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc
index d8574afcb..e718d9f2f 100644
--- a/expui/BiorthBasis.cc
+++ b/expui/BiorthBasis.cc
@@ -1204,6 +1204,7 @@ namespace BasisClasses
     "rcylmin",
     "rcylmax",
     "acyl",
+    "bias",
     "hcyl",
     "sech2",
     "snr",
@@ -1391,6 +1392,7 @@ namespace BasisClasses
     rcylmin     = 0.001;
     rcylmax     = 20.0;
     acyl        = 0.01;
+    bias        = 1.0;
     hcyl        = 0.002;
     nmax        = 18;
     mmax        = 6;
@@ -1415,7 +1417,7 @@ namespace BasisClasses
     EVEN_M      = false;
     cmapR       = 1;
     cmapZ       = 1;
-    mtype       = "exponential";
+    mtype       = "ExpSphere";
     dtype       = "exponential";
     vflag       = 0;
     
@@ -1469,6 +1471,7 @@ namespace BasisClasses
       if (conf["rcylmax"   ])    rcylmax  = conf["rcylmax"   ].as<double>();
       
       if (conf["acyl"      ])       acyl  = conf["acyl"      ].as<double>();
+      if (conf["bias"      ])       bias  = conf["bias"      ].as<double>();
       if (conf["hcyl"      ])       hcyl  = conf["hcyl"      ].as<double>();
       if (conf["sech2"     ])      sech2  = conf["sech2"     ].as<bool>();
       if (conf["lmaxfid"   ])    lmaxfid  = conf["lmaxfid"   ].as<int>();
@@ -1555,6 +1558,17 @@ namespace BasisClasses
     pnum = std::max<int>(1,  pnum);
     tnum = std::max<int>(10, tnum);
     
+    // Validate bias parameter
+    //
+    if (!std::isfinite(bias)) {
+      throw std::runtime_error("Cylindrical: 'bias' parameter must be finite");
+    }
+    if (bias <= 0.0) {
+      std::ostringstream sout;
+      sout << "Cylindrical: 'bias' parameter must be positive, got " << bias;
+      throw std::runtime_error(sout.str());
+    }
+    
     EmpCylSL::RMIN        = rcylmin;
     EmpCylSL::RMAX        = rcylmax;
     EmpCylSL::NUMX        = ncylnx;
@@ -1565,6 +1579,36 @@ namespace BasisClasses
     EmpCylSL::logarithmic = logarithmic;
     EmpCylSL::VFLAG       = vflag;
 
+    // Set deprojected model type
+    //
+    // Convert mtype string to lower case
+    std::transform(mtype.begin(), mtype.end(), mtype.begin(),
+		   [](unsigned char c){ return std::tolower(c); });
+      
+    // Set EmpCylSL mtype.  This is the spherical function used to
+    // generate the EOF basis.  If "deproject" is set, this will be
+    // overriden in EmpCylSL.
+    //
+    EmpCylSL::mtype = EmpCylSL::ExpSphere; // Default
+    if (mtype.compare("exponential")==0)
+      EmpCylSL::mtype = EmpCylSL::Exponential;
+    else if (mtype.compare("expsphere")==0)
+      EmpCylSL::mtype = EmpCylSL::ExpSphere;
+    else if (mtype.compare("gaussian")==0)
+      EmpCylSL::mtype = EmpCylSL::Gaussian;
+    else if (mtype.compare("plummer")==0)
+      EmpCylSL::mtype = EmpCylSL::Plummer;
+    else if (mtype.compare("power")==0) {
+      EmpCylSL::mtype = EmpCylSL::Power;
+      EmpCylSL::PPOW  = ppow;
+    } else {
+      if (myid==0) std::cout << "No EmpCylSL EmpModel named <"
+			     << mtype << ">, valid types are: "
+			     << "Exponential, ExpSphere, Gaussian, Plummer, Power "
+			     << "(not case sensitive)" << std::endl;
+      throw std::runtime_error("Cylindrical:initialize: EmpCylSL bad parameter");
+    }
+
     // Check for non-null cache file name.  This must be specified
     // to prevent recomputation and unexpected behavior.
     //
@@ -1580,7 +1624,7 @@ namespace BasisClasses
     // Make the empirical orthogonal basis instance
     //
     sl = std::make_shared<EmpCylSL>
-      (nmaxfid, lmaxfid, mmax, nmax, acyl, hcyl, ncylodd, cachename);
+      (nmaxfid, lmaxfid, mmax, nmax, bias*acyl, hcyl, ncylodd, cachename);
     
     // Set azimuthal harmonic order restriction?
     //
@@ -1612,33 +1656,6 @@ namespace BasisClasses
 		  << "----               remove 'ignore' from your YAML configuration." << std::endl;
       }
 
-      // Convert mtype string to lower case
-      //
-      std::transform(mtype.begin(), mtype.end(), mtype.begin(),
-		     [](unsigned char c){ return std::tolower(c); });
-      
-      // Set EmpCylSL mtype.  This is the spherical function used to
-      // generate the EOF basis.  If "deproject" is set, this will be
-      // overriden in EmpCylSL.
-      //
-      EmpCylSL::mtype = EmpCylSL::Exponential;
-      if (mtype.compare("exponential")==0)
-	EmpCylSL::mtype = EmpCylSL::Exponential;
-      else if (mtype.compare("gaussian")==0)
-	EmpCylSL::mtype = EmpCylSL::Gaussian;
-      else if (mtype.compare("plummer")==0)
-	EmpCylSL::mtype = EmpCylSL::Plummer;
-      else if (mtype.compare("power")==0) {
-	EmpCylSL::mtype = EmpCylSL::Power;
-	EmpCylSL::PPOW  = ppow;
-      } else {
-	if (myid==0) std::cout << "No EmpCylSL EmpModel named <"
-			       << mtype << ">, valid types are: "
-			       << "Exponential, Gaussian, Plummer, Power "
-			       << "(not case sensitive)" << std::endl;
-	throw std::runtime_error("Cylindrical:initialize: EmpCylSL bad parameter");
-      }
-
       // Convert dtype string to lower case
       //
       std::transform(dtype.begin(), dtype.end(), dtype.begin(),
@@ -5124,6 +5141,7 @@ namespace BasisClasses
     file.createAttribute<double>("rcylmin", HighFive::DataSpace::From(rcylmin)).write(rcylmin);
     file.createAttribute<double>("rcylmax", HighFive::DataSpace::From(rcylmax)).write(rcylmax);
     file.createAttribute<double>("acyl", HighFive::DataSpace::From(acyl)).write(acyl);
+    file.createAttribute<double>("bias", HighFive::DataSpace::From(bias)).write(bias);
     file.createAttribute<double>("hcyl", HighFive::DataSpace::From(hcyl)).write(hcyl);
   }
   
diff --git a/exputil/CMakeLists.txt b/exputil/CMakeLists.txt
index 0a029c2ec..1124b101e 100644
--- a/exputil/CMakeLists.txt
+++ b/exputil/CMakeLists.txt
@@ -11,7 +11,7 @@ set(QUAD_SRC qadapt.cc gauleg.cc qadapt2d.cc gint2.cc rombe2d.cc
 set(UTIL_SRC nrutil.cc elemfunc.cc euler.cc euler_slater.cc
 	     rotmatrix.cc wordSplit.cc FileUtils.cc BarrierWrapper.cc
 	     stack.cc localmpi.cc TableGrid.cc writePVD.cc libvars.cc
-	     TransformFFT.cc QDHT.cc YamlCheck.cc # Hankel.cc
+	     TransformFFT.cc QDHT.cc YamlCheck.cc ExpDeproj.cc # Hankel.cc
 	     parseVersionString.cc EXPmath.cc laguerre_polynomial.cpp
 	     YamlConfig.cc orthoTest.cc OrthoFunction.cc VtkGrid.cc
 	     Sutils.cc fpetrap.cc)
diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc
index a89afa68e..f474ce3cc 100644
--- a/exputil/EmpCylSL.cc
+++ b/exputil/EmpCylSL.cc
@@ -16,6 +16,7 @@
 #include <thread>
 #include "exp_thread.h"
 #include "EXPException.H"
+#include "ExpDeproj.H"
 #include "exputils.H"
 
 #include <Eigen/Eigenvalues>
@@ -68,10 +69,11 @@ double   EmpCylSL::PPOW            = 4.0;
 bool     EmpCylSL::NewCoefs        = true;
  
 
-EmpCylSL::EmpModel EmpCylSL::mtype = Exponential;
+EmpCylSL::EmpModel EmpCylSL::mtype = ExpSphere;
 
 std::map<EmpCylSL::EmpModel, std::string> EmpCylSL::EmpModelLabs =
   { {Exponential, "Exponential"},
+    {ExpSphere,   "ExpSphere"  },
     {Gaussian,    "Gaussian"   },
     {Plummer,     "Plummer"    },
     {Power,       "Power"      },
@@ -551,6 +553,9 @@ double EmpCylSL::massR(double R)
   case Exponential:
     ans = 1.0 - (1.0 + R)*exp(-R); 
     break;
+  case ExpSphere:
+    ans = expdeproj.mass(R);
+    break;
   case Gaussian:
     arg = 0.5*R*R;
     ans = 1.0 - exp(-arg);
@@ -588,6 +593,9 @@ double EmpCylSL::densR(double R)
   case Exponential:
     ans = exp(-R)/(4.0*M_PI*R);
     break;
+  case ExpSphere:
+    ans = expdeproj.density(R);
+    break;
   case Gaussian:
     arg = 0.5*R*R;
     ans = exp(-arg)/(4.0*M_PI*R);
@@ -2398,6 +2406,15 @@ void EmpCylSL::generate_eof(int numr, int nump, int numt,
 
   int cntr = 0;			// Loop counter for spreading load to nodes
   
+  // Debug info output
+  //
+  if (VFLAG & 8 && myid==0) {
+    std::cout << "---- EmpCylSL: Generating EOF with"
+	      << " Rmin=" << xi_to_r(XMIN)
+      	      << " Rmax=" << xi_to_r(XMAX)
+	      << " numt=" << numt << std::endl;
+  }
+
   // *** Radial quadrature loop
   //
   for (int qr=0; qr<numr; qr++) { 
@@ -7516,6 +7533,7 @@ bool EmpCylSL::ReadH5Cache()
 
     if (not checkStr(geometry, "geometry"))  return false;
     if (not checkStr(forceID,  "forceID"))   return false;
+    if (not checkStr(model,    "model"))     return false;
     if (not checkInt(MMAX,     "mmax"))      return false;
     if (not checkInt(NUMX,     "numx"))      return false;
     if (not checkInt(NUMY,     "numy"))      return false;
diff --git a/exputil/ExpDeproj.cc b/exputil/ExpDeproj.cc
new file mode 100644
index 000000000..183711d6d
--- /dev/null
+++ b/exputil/ExpDeproj.cc
@@ -0,0 +1,85 @@
+#include "ExpDeproj.H"
+
+#include <cmath>
+#include <stdexcept>
+#include <iostream>
+#include <fstream>
+#include <cassert>
+#include <vector>
+#include <algorithm>
+#include <numeric>
+#include <functional>
+
+void ExpDeproj::initialize()
+{
+  if (ngrid < 2) {
+    throw std::invalid_argument("n must be at least 2");
+  }
+
+  rv.resize(ngrid);
+  mv.resize(ngrid);
+
+  std::vector<double> dv(ngrid);
+
+  double log_rmin = std::log(rmin);
+  double log_rmax = std::log(rmax);
+  double dlogr = (log_rmax - log_rmin)/static_cast<double>(ngrid - 1);
+
+  for (int i = 0; i < ngrid; ++i) {
+    rv[i] = std::exp(log_rmin + i*dlogr);
+    dv[i] = 4.0*M_PI*rv[i]*rv[i]*density(rv[i]);
+  }
+  
+  mv[0] = 0.0;
+  for (int i=1; i<ngrid; i++)
+    mv[i] = mv[i-1] + 0.5*(dv[i] + dv[i-1])*(rv[i] - rv[i-1]);
+}
+
+double ExpDeproj::density(double R)
+{
+  if (R < 0) {
+    throw std::invalid_argument("R must be non-negative");
+  }
+  return 0.5*std::cyl_bessel_k(0.0, R)/(M_PI*M_PI);
+}
+
+double ExpDeproj::mass(double R)
+{
+  // Initialize precomputed arrays?
+  if (mv.size() == 0) {
+    initialize();
+  }
+
+  if (R < 0) {
+    throw std::invalid_argument("R must be non-negative");
+  }
+
+  if (R < rmin) return 0.0;
+  if (R > rmax) return mv.back();
+
+  auto it = std::lower_bound(rv.begin(), rv.end(), R);
+  // If R is slightly larger than the largest grid point due to rounding,
+  // lower_bound may return rv.end(); in that case, use the last mass value.
+  if (it == rv.end()) {
+    return mv.back();
+  }
+
+  std::size_t idx = static_cast<std::size_t>(std::distance(rv.begin(), it));
+  if (idx >= rv.size()) {
+    // Defensive guard; should not happen after the it == rv.end() check.
+    return mv.back();
+  }
+  if (rv[idx] == R) {
+    return mv[idx];
+  } else {
+    // Ensure idx-1 is valid
+    if (idx == 0) idx++;
+    // Linear interpolation
+    double x0 = rv[idx-1];
+    double x1 = rv[idx  ];
+    double y0 = mv[idx-1];
+    double y1 = mv[idx  ];
+
+    return y0 + (y1 - y0)*(R - x0)/(x1 - x0);
+  }
+}
diff --git a/include/EmpCylSL.H b/include/EmpCylSL.H
index ebb2af947..f6f4d268f 100644
--- a/include/EmpCylSL.H
+++ b/include/EmpCylSL.H
@@ -17,6 +17,7 @@
 #include <dmalloc.h>
 #endif
 
+#include "ExpDeproj.H"
 #include "Particle.H"
 #include "SLGridMP2.H"
 #include "coef.H"
@@ -80,6 +81,9 @@ protected:
 
   int rank2, rank3;
 
+  //! Precomputed exponential deprojection profile
+  ExpDeproj expdeproj;
+
   //@{
   //! Storage buffers for MPI
   std::vector<double> MPIin, MPIout, MPIin2, MPIout2;
@@ -364,6 +368,7 @@ public:
   //! Type of density model to use
   enum EmpModel {
     Exponential,
+    ExpSphere,
     Gaussian, 
     Plummer,
     Power,
diff --git a/include/ExpDeproj.H b/include/ExpDeproj.H
new file mode 100644
index 000000000..1c9da8a02
--- /dev/null
+++ b/include/ExpDeproj.H
@@ -0,0 +1,28 @@
+#pragma once
+
+#include <vector>
+#include <cmath>
+
+class ExpDeproj
+{
+  const double rmin = 1.0e-4;
+  const double rmax = 30.0;
+
+  // Precomputed radial grid
+  int ngrid;
+  std::vector<double> rv, mv;
+  void initialize();
+
+public:
+  //! Constructor
+  ExpDeproj(int n=4000) : ngrid(n) {}
+
+  //! Destructor
+  virtual ~ExpDeproj() {}
+
+  //! Density profile at radius R
+  double density(double R);
+
+  //! Enclosed mass at radius R
+  double mass(double R);
+};
diff --git a/src/Cylinder.H b/src/Cylinder.H
index d7fcac153..a97f0da6e 100644
--- a/src/Cylinder.H
+++ b/src/Cylinder.H
@@ -34,6 +34,8 @@ class MixtureBasis;
 
     @param acyl is the scale length
 
+    @param bias is the O(1) bias correction for deprojected scale length
+
     @param hcyl is the scale height
 
     @param lmaxfid is the maximum spherical harmonic index for EOF construction
@@ -144,7 +146,7 @@ private:
   void initialize(void);
 
 				// Parameters
-  double rcylmin, rcylmax, zmax, acyl;
+  double rcylmin, rcylmax, zmax, acyl, bias;
   int nmaxfid, lmaxfid, mmax, mlim;
   int ncylnx, ncylny, ncylr;
   double hcyl, hexp, snr, rem;
@@ -334,6 +336,7 @@ public:
   //! 
   //! Parameters settable in the configuration line include:
   //! \param acyl is the disk scale length for the basis
+  //! \param bias is the disk scale length bias factor for the basis
   //! \param hcyl is the disk scale height for the basis
   //! \param nmaxfid is the radial order for the underlying spherical basis
   //! \param lmaxfid is the harmonic order for the underlying spherical basis
diff --git a/src/Cylinder.cc b/src/Cylinder.cc
index 448a5f769..1c01437c2 100644
--- a/src/Cylinder.cc
+++ b/src/Cylinder.cc
@@ -28,6 +28,7 @@ Cylinder::valid_keys = {
   "rcylmin",
   "rcylmax",
   "acyl",
+  "bias",
   "hcyl",
   "sech2",
   "hexp",
@@ -104,6 +105,7 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) :
   ncylr           = 2000;
 
   acyl            = 0.01;
+  bias            = 1.0;
   lmaxfid         = 128;
   nmaxfid         = 64;
   mmax            = 6;
@@ -176,7 +178,7 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) :
   // Make the empirical orthogonal basis instance
   //
   ortho = std::make_shared<CylEXP>
-    (nmaxfid, lmaxfid, mmax, nmax, acyl, hcyl, ncylodd, cachename);
+    (nmaxfid, lmaxfid, mmax, nmax, bias*acyl, hcyl, ncylodd, cachename);
   
   // Set azimuthal harmonic order restriction?
   //
@@ -332,6 +334,7 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) :
 		<< std::endl << sep << "rcylmin="     << rcylmin
 		<< std::endl << sep << "rcylmax="     << rcylmax
 		<< std::endl << sep << "acyl="        << acyl
+		<< std::endl << sep << "bias="        << bias
 		<< std::endl << sep << "hcyl="        << hcyl
 		<< std::endl << sep << "precond="     << precond
 		<< std::endl << sep << "pcavar="      << pcavar
@@ -361,6 +364,7 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) :
 	      << std::endl << sep << "rcylmin="     << rcylmin
 	      << std::endl << sep << "rcylmax="     << rcylmax
 	      << std::endl << sep << "acyl="        << acyl
+	      << std::endl << sep << "bias="        << bias
 	      << std::endl << sep << "hcyl="        << hcyl
 	      << std::endl << sep << "precond="     << precond
 	      << std::endl << sep << "pcavar="      << pcavar
@@ -417,6 +421,7 @@ void Cylinder::initialize()
     if (conf["rcylmax"   ])    rcylmax  = conf["rcylmax"   ].as<double>();
 
     if (conf["acyl"      ])       acyl  = conf["acyl"      ].as<double>();
+    if (conf["bias"      ])       bias  = conf["bias"      ].as<double>();
     if (conf["hcyl"      ])       hcyl  = conf["hcyl"      ].as<double>();
     if (conf["sech2"     ])      sech2  = conf["sech2"     ].as<bool>();
     if (conf["hexp"      ])       hexp  = conf["hexp"      ].as<double>();
diff --git a/utils/Test/CMakeLists.txt b/utils/Test/CMakeLists.txt
index 6d7ed3a45..5a8a6fc4b 100644
--- a/utils/Test/CMakeLists.txt
+++ b/utils/Test/CMakeLists.txt
@@ -1,5 +1,5 @@
 
-set(bin_PROGRAMS testBarrier expyaml orthoTest)
+set(bin_PROGRAMS testBarrier expyaml orthoTest testED)
 
 set(common_LINKLIB OpenMP::OpenMP_CXX MPI::MPI_CXX expui exputil
   yaml-cpp ${VTK_LIBRARIES})
@@ -32,6 +32,7 @@ endif()
 add_executable(testBarrier    test_barrier.cc)
 add_executable(expyaml        test_config.cc)
 add_executable(orthoTest      orthoTest.cc Biorth2Ortho.cc)
+add_executable(testED         testED.cc)
 
 foreach(program ${bin_PROGRAMS})
   target_link_libraries(${program} ${common_LINKLIB})
diff --git a/utils/Test/testED.cc b/utils/Test/testED.cc
new file mode 100644
index 000000000..9fb09d4dc
--- /dev/null
+++ b/utils/Test/testED.cc
@@ -0,0 +1,55 @@
+#include <iostream>
+#include <iomanip>
+#include "ExpDeproj.H"
+
+double projectedDensity(double R, double rmax, int nsteps, ExpDeproj & deproj)
+{
+  double dz = rmax / double(nsteps-1);
+  double sum = 0.0;
+  double curr = deproj.density(R), last = 0.0;
+  for (int i=1; i<nsteps; i++) {
+    last = curr;
+    double z = dz * i;
+    double r = sqrt(R*R + z*z);
+    curr = deproj.density(r);
+    sum += (curr + last) * dz;
+  }
+  return sum;
+}
+
+int main()
+{
+  ExpDeproj deproj;
+
+  double rmin = 1.0e-4;
+  double rmax = 30.0;
+  int  nsteps = 4000;
+
+  std::cout << "# Testing ExpDeproj class" << std::endl;
+  std::cout << "# The density and mass columns are the deprojected model" << std::endl;
+  std::cout << "# The projDensity column is the projected deprojected density calculated via numerical integration" << std::endl;
+  std::cout << "# The error column is the difference between the projected density and exponential surface density" << std::endl;
+  std::cout << "# " << std::endl;
+  std::cout << std::setw(18) << "# R"
+	    << std::setw(18) << "density"
+	    << std::setw(18) << "mass"
+	    << std::setw(18) << "projDensity"
+	    << std::setw(18) << "error"
+	    << std::endl;
+
+  for (int i=0; i<100; i++) {
+    double R = rmin + (0.5+i)*(rmax - rmin)/100.0;
+    double dens = deproj.density(R);
+    double mass = deproj.mass(R);
+    double proj = projectedDensity(R, rmax, nsteps, deproj);
+    double expt = 0.5 * exp(-R) / M_PI;
+    std::cout << std::setw(18) << R
+	      << std::setw(18) << dens
+	      << std::setw(18) << mass
+	      << std::setw(18) << proj
+	      << std::setw(18) << (proj - expt)/expt
+	      << std::endl;
+  }
+  
+  return 0;
+}