BiorthBasis::Cylinder deprojection improvements

expui/BiorthBasis.H

@@ -996,11 +996,12 @@ namespace BasisClasses
     //! Basis construction parameters
     double aratio, hratio, dweight, rwidth, ashift, rfactor, rtrunc, ppow, Mfac, HERNA;
     bool Ignore, deproject;
-    std::string pyname;
+    std::string pyname, pyproj;
 
     //! DiskType support
-    //
-    enum DiskType { constant, gaussian, mn, exponential, doubleexpon, diskbulge, python };
+    //!
+    enum class DiskType
+      { constant, gaussian, mn, exponential, doubleexpon, diskbulge, python };
     DiskType DTYPE;
     std::string mtype, dtype, dmodel;
 
@@ -1009,6 +1010,20 @@ namespace BasisClasses
     double dcond(double R, double z, double phi, int M);
     std::shared_ptr<DiskDensityFunc> pyDens;
 
+    //@{
+    //! DeprojType support
+
+    //! Disk models used for deprojection
+    enum class DeprojType
+      { mn, toomre, python, exponential};
+
+    //! Current model
+    DeprojType PTYPE;
+
+    //! Look up by string
+    static const std::map<std::string, DeprojType> dplookup;
+    //@}
+
   protected:
 
     //! Evaluate basis in spherical coordinates

expui/BiorthBasis.cc

@@ -1198,6 +1198,14 @@ namespace BasisClasses
       {"python",      DiskType::python}
     };
 
+  // Deprojection model for cylindrical basis construction
+  const std::map<std::string, Cylindrical::DeprojType> Cylindrical::dplookup =
+    { {"mn",          DeprojType::mn},
+      {"exponential", DeprojType::exponential},
+      {"toomre",      DeprojType::toomre},
+      {"python",      DeprojType::python}
+    };
+
   const std::set<std::string>
   Cylindrical::valid_keys = {
     "tk_type",
@@ -1233,6 +1241,7 @@ namespace BasisClasses
     "expcond",
     "ignore",
     "deproject",
+    "dmodel",
     "logr",
     "pcavar",
     "pcaeof",
@@ -1261,7 +1270,8 @@ namespace BasisClasses
     "playback",
     "coefCompute",
     "coefMaster",
-    "pyname"
+    "pyname",
+    "pyproj"
   };
 
   Cylindrical::Cylindrical(const YAML::Node& CONF) :
@@ -1494,7 +1504,7 @@ namespace BasisClasses
       if (conf["cmapz"     ])      cmapZ  = conf["cmapz"     ].as<int>();
       if (conf["ignore"    ])      Ignore = conf["ignore"    ].as<bool>();
       if (conf["deproject" ])   deproject = conf["deproject" ].as<bool>();
-      if (conf["dmodel"    ])      dmodel = conf["dmodel"    ].as<bool>();
+      if (conf["dmodel"    ])      dmodel = conf["dmodel"    ].as<std::string>();
 
       if (conf["aratio"    ])      aratio = conf["aratio"    ].as<double>();
       if (conf["hratio"    ])      hratio = conf["hratio"    ].as<double>();
@@ -1510,6 +1520,7 @@ namespace BasisClasses
       if (conf["dtype"     ])       dtype = conf["dtype"     ].as<std::string>();
       if (conf["vflag"     ])       vflag = conf["vflag"     ].as<int>();
       if (conf["pyname"    ])      pyname = conf["pyname"    ].as<std::string>();
+      if (conf["pyproj"    ])      pyproj = conf["pyproj"    ].as<std::string>();
       if (conf["pcavar"]    )      pcavar = conf["pcavar"    ].as<bool>();
       if (conf["subsamp"]   )      sampT  = conf["subsamp"   ].as<int>();
 
@@ -1697,16 +1708,50 @@ namespace BasisClasses
 	//
 	EmpCylSL::AxiDiskPtr model;
 
-	if (dmodel.compare("MN")==0) // Miyamoto-Nagai
+	// Convert dmodel string to lower case
+	//
+	std::transform(dmodel.begin(), dmodel.end(), dmodel.begin(),
+		       [](unsigned char c){ return std::tolower(c); });
+
+	// Map legacy/short model names to canonical keys expected by dplookup
+	//
+	if (dmodel == "exp") {
+	  dmodel = "exponential";
+	}
+
+	// Check for map entry
+	//
+	try {
+	  PTYPE = dplookup.at(dmodel);
+
+	  // Report DeprojType
+	  if (myid==0) {
+	    std::cout << "---- Deprojection type is <" << dmodel
+		      << ">" << std::endl;
+	  }
+	}
+	catch (const std::out_of_range& err) {
+	  if (myid==0) {
+	    std::cout << "DeprojType error in configuration file" << std::endl;
+	    std::cout << "Valid options are: ";
+	    for (auto v : dplookup) std::cout << v.first << " ";
+	    std::cout << std::endl;
+	  }
+	  throw std::runtime_error("Cylindrical::initialize: invalid DiskModel");
+	}
+
+	if (PTYPE == DeprojType::mn) // Miyamoto-Nagai
 	  model = std::make_shared<MNdisk>(1.0, H);
-	else if (DTYPE == DiskType::python) {
-	  model = std::make_shared<AxiSymPyModel>(pyname, acyl);
+	else if (PTYPE == DeprojType::toomre) {
+	  model = std::make_shared<Toomre>(1.0, H, 5.0);
+	} else if (PTYPE == DeprojType::python) {
+	  model = std::make_shared<AxiSymPyModel>(pyproj, acyl);
 	  std::cout << "Using AxiSymPyModel for deprojection from Python function <"
 		    << pyname << ">" << std::endl;
-	}
-	else			// Default to exponential
+	} else {		// Default to exponential
 	  model = std::make_shared<Exponential>(1.0, H);
-
+	}
+
 	if (rwidth>0.0) {
 	  model = std::make_shared<Truncated>(rtrunc/acyl,
 					      rwidth/acyl,

exputil/EmpCylSL.cc

@@ -487,28 +487,45 @@ void EmpCylSL::create_deprojection(double H, double Rf, int NUMR, int NINT,
   // Now, compute Abel inverion integral
   //
   for (int i=0; i<NUMR; i++) {
+
+    double surfS = abel_type==AbelType::Subtraction ? surf.eval(rl[i]) : 0.0;
     double r = rr[i];
 
     // Interval by Legendre
     //
     rhoI[i] = 0.0;
+
     for (int n=0; n<NINT; n++) {
+
       double x   = lq.knot(n);
       double x12 = 1.0 - x*x;
-      double z   = x/sqrt(x12);
-      double R   = sqrt(z*z + r*r);
+      double R   = r/sqrt(x12);
       double lR  = log(R);
 
-      rhoI[i]   += lq.weight(n)*2.0*pow(x12, -1.5)*surf.eval(lR);
+      switch (abel_type) {
+      case AbelType::Derivative:
+	rhoI[i]   += lq.weight(n)/x12 * surf.deriv(lR)/R;
+	break;
+      case AbelType::Subtraction:
+	rhoI[i]   += lq.weight(n)/(x*x*sqrt(x12)*r) * ( surf.eval(lR) - surfS );
+	break;
+      case AbelType::IBP:
+	rhoI[i]   += lq.weight(n)/x12 * surf.eval(lR);
+	break;
+      }
     }
   }
 
-  std::vector<double> rho(NUMR), mass(NUMR);
-
-  Linear1d intgr(rl, rhoI);
+  Linear1d intgr;
+  if (abel_type == AbelType::IBP) intgr = Linear1d(rl, rhoI);
 
-  for (int i=0; i<NUMR; i++)
-    rho[i] = -intgr.deriv(rl[i])/(2.0*M_PI*rr[i]*rr[i]);
+  std::vector<double> rho(NUMR), mass(NUMR);
+  for (int i=0; i<NUMR; i++) {
+    if (abel_type == AbelType::IBP)
+      rho[i] = -intgr.deriv(rl[i])/rr[i]/M_PI;
+    else
+      rho[i] = -rhoI[i]/M_PI;
+  }
 
   mass[0] = 0.0;
   for (int i=1; i<NUMR; i++) {

include/DiskModels.H

@@ -3,6 +3,7 @@
 
 #include "EmpCylSL.H"
 #include "DiskDensityFunc.H"
+#include <stdexcept>
 
 //! A Ferrers Ellipsoid + Evacuated Exponential Disc (semi-realistic
 //! bar+disk model)
@@ -196,6 +197,34 @@ public:
   }
 };
 
+//! Toomre disk
+class Toomre : public EmpCylSL::AxiDisk
+{
+private:
+  double a, h, n;
+  double norm;
+
+public:
+
+  Toomre(double a, double h, double n=3, double M=1) :
+    a(a), h(h), n(n), EmpCylSL::AxiDisk(M, "Toomre")
+  {
+    params.push_back(a);
+    params.push_back(h);
+    params.push_back(n);
+    if (n <= 2) throw std::runtime_error("Toomre index must be > 2");
+    norm = M*(n - 2.0)/(2.0*M_PI*a*a*4.0*h);
+  }
+
+  double operator()(double R, double z, double phi=0.)
+  {
+    double sigma = norm * std::pow(1.0 + (R/a)*(R/a), -0.5*n);
+    double Q = std::exp(-0.5*std::fabs(z)/h);
+    double sech = 2.0*Q/(1.0 + Q*Q); // Prevent overflow
+    return sigma * sech*sech;
+  }
+};
+
 //! Truncate a AxiDisk
 class Truncated : public EmpCylSL::AxiDisk
 {

include/EmpCylSL.H

@@ -351,6 +351,12 @@ protected:
   //! with the new Eigen3 API
   bool allow_old_cache = false;
 
+  //! Abel type for deprojection
+  enum class AbelType { Derivative, Subtraction, IBP };
+
+  //! Deprojection method (default: derivative)
+  AbelType abel_type = AbelType::Derivative;
+
 public:
 
   /*! Enum listing the possible selection algorithms for coefficient

utils/Test/CMakeLists.txt

@@ -1,5 +1,6 @@
 
-set(bin_PROGRAMS testBarrier expyaml orthoTest)
+set(bin_PROGRAMS testBarrier expyaml orthoTest testDeproj
+		 testDeprojPlum testEmpDeproj testEmp)
 
 set(common_LINKLIB OpenMP::OpenMP_CXX MPI::MPI_CXX expui exputil
   yaml-cpp ${VTK_LIBRARIES})
@@ -32,6 +33,12 @@ endif()
 add_executable(testBarrier    test_barrier.cc)
 add_executable(expyaml        test_config.cc)
 add_executable(orthoTest      orthoTest.cc Biorth2Ortho.cc)
+add_executable(testDeproj     testDeproject.cc CubicSpline.cc Deprojector.cc)
+add_executable(testDeprojPlum testDeprojPlummer.cc CubicSpline.cc
+			      Deprojector.cc)
+add_executable(testEmpDeproj  testEmpDeproj.cc CubicSpline.cc
+			      Deprojector.cc EmpDeproj.cc)
+add_executable(testEmp        testEmp.cc EmpDeproj.cc)
 
 foreach(program ${bin_PROGRAMS})
   target_link_libraries(${program} ${common_LINKLIB})
@@ -40,4 +47,7 @@ foreach(program ${bin_PROGRAMS})
   # install(TARGETS ${program} DESTINATION bin)
 endforeach()
 
+# For Python interpreter...
+target_link_libraries(testEmp ${common_LINKLIB} pybind11::embed)
+
 install(TARGETS expyaml DESTINATION bin)

utils/Test/CubicSpline.H

@@ -0,0 +1,31 @@
+#ifndef CUBIC_SPLINE_H
+#define CUBIC_SPLINE_H
+
+#include <vector>
+
+namespace Deproject
+{
+
+  class CubicSpline {
+  public:
+    CubicSpline() = default;
+    CubicSpline(const std::vector<double>& x_in, const std::vector<double>& y_in);
+
+    // set data (x must be strictly increasing)
+    void set_data(const std::vector<double>& x_in, const std::vector<double>& y_in);
+
+    // evaluate spline and its derivative (xx should lie within [xmin(), xmax()], but endpoints are clamped)
+    double eval(double xx) const;
+    double deriv(double xx) const;
+
+    double xmin() const;
+    double xmax() const;
+
+  private:
+    std::vector<double> x_, y_, y2_;
+    int locate(double xx) const;
+  };
+
+}
+
+#endif // CUBIC_SPLINE_H

utils/Test/CubicSpline.cc

@@ -0,0 +1,73 @@
+#include <stdexcept>
+
+#include "CubicSpline.H"
+
+namespace Deproject
+{
+
+  CubicSpline::CubicSpline(const std::vector<double>& x_in, const std::vector<double>& y_in) {
+    set_data(x_in, y_in);
+  }
+
+  void CubicSpline::set_data(const std::vector<double>& x_in, const std::vector<double>& y_in) {
+    x_ = x_in;
+    y_ = y_in;
+    int n = (int)x_.size();
+    if (n < 2 || (int)y_.size() != n) throw std::runtime_error("CubicSpline: need at least two points and equal-sized x,y.");
+
+    y2_.assign(n, 0.0);
+    std::vector<double> u(n - 1, 0.0);
+
+    // natural spline boundary conditions (second derivatives at endpoints = 0)
+    for (int i = 1; i < n - 1; ++i) {
+      double sig = (x_[i] - x_[i-1]) / (x_[i+1] - x_[i-1]);
+      double p = sig * y2_[i-1] + 2.0;
+      y2_[i] = (sig - 1.0) / p;
+      double dY1 = (y_[i+1] - y_[i]) / (x_[i+1] - x_[i]);
+      double dY0 = (y_[i]   - y_[i-1]) / (x_[i]   - x_[i-1]);
+      u[i] = (6.0 * (dY1 - dY0) / (x_[i+1] - x_[i-1]) - sig * u[i-1]) / p;
+    }
+
+    for (int k = n - 2; k >= 0; --k) y2_[k] = y2_[k] * y2_[k+1] + u[k];
+  }
+
+  int CubicSpline::locate(double xx) const {
+    int n = (int)x_.size();
+    if (xx <= x_.front()) return 0;
+    if (xx >= x_.back()) return n - 2;
+    int lo = 0, hi = n - 1;
+    while (hi - lo > 1) {
+      int mid = (lo + hi) >> 1;
+      if (x_[mid] > xx) hi = mid; else lo = mid;
+    }
+    return lo;
+  }
+
+  double CubicSpline::eval(double xx) const {
+    int klo = locate(xx);
+    int khi = klo + 1;
+    double h = x_[khi] - x_[klo];
+    if (h <= 0.0) throw std::runtime_error("CubicSpline::eval: non-increasing x.");
+    double A = (x_[khi] - xx) / h;
+    double B = (xx - x_[klo]) / h;
+    double val = A * y_[klo] + B * y_[khi]
+      + ((A*A*A - A) * y2_[klo] + (B*B*B - B) * y2_[khi]) * (h*h) / 6.0;
+    return val;
+  }
+
+  double CubicSpline::deriv(double xx) const {
+    int klo = locate(xx);
+    int khi = klo + 1;
+    double h = x_[khi] - x_[klo];
+    if (h <= 0.0) throw std::runtime_error("CubicSpline::deriv: non-increasing x.");
+    double A = (x_[khi] - xx) / h;
+    double B = (xx - x_[klo]) / h;
+    double dy = (y_[khi] - y_[klo]) / h
+      + ( - (3.0*A*A - 1.0) * y2_[klo] + (3.0*B*B - 1.0) * y2_[khi] ) * (h / 6.0);
+    return dy;
+  }
+
+  double CubicSpline::xmin() const { return x_.front(); }
+  double CubicSpline::xmax() const { return x_.back(); }
+
+}