Adds an spherical basis that is an exact deprojection of an exponential surface density

.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

_codeql_detected_source_root

@@ -0,0 +1 @@
+.

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;

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);
   }
 

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)

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;

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);
+  }
+}

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,

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);
+};