Merge pull request #103 from CarrieFilion/diskbulge

expui/BiorthBasis.H

@@ -8,6 +8,7 @@
 #include <unsupported/Eigen/CXX11/Tensor> // For 3d rectangular grids
 #include <yaml-cpp/yaml.h>
 
+#include <DiskDensityFunc.H>
 #include <ParticleReader.H>
 #include <Coefficients.H>
 #include <BiorthBess.H>
@@ -301,7 +302,7 @@ namespace BasisClasses
     {
       auto [ret, worst, lworst] = orthoCompute(orthoCheck(knots));
       // For the CTest log
-      std::cout << "Spherical::orthoTest: worst=" << worst << std::endl;
+      std::cout << "---- Spherical::orthoTest: worst=" << worst << std::endl;
       return ret;
     }
 
@@ -536,7 +537,7 @@ namespace BasisClasses
     {
       auto [ret, worst, lworst] = orthoCompute(orthoCheck());
       // For the CTest log
-      std::cout << "FlatDisk::orthoTest: worst=" << worst << std::endl;
+      std::cout << "---- FlatDisk::orthoTest: worst=" << worst << std::endl;
       return ret;
     }
 
@@ -741,18 +742,20 @@ namespace BasisClasses
 
     //@{
     //! Basis construction parameters
-    double aratio, hratio, dweight, rwidth, ashift, rfactor, rtrunc, ppow;
+    double aratio, hratio, dweight, rwidth, ashift, rfactor, rtrunc, ppow, Mfac, HERNA;
     bool Ignore, deproject;
+    std::string pyname;
 
     //! DiskType support
     //
-    enum DiskType { constant, gaussian, mn, exponential, doubleexpon };
+    enum DiskType { constant, gaussian, mn, exponential, doubleexpon, diskbulge, python };
     DiskType DTYPE;
     std::string mtype, dtype, dmodel;
 
     static const std::map<std::string, DiskType> dtlookup;
     double DiskDens(double R, double z, double phi);
     double dcond(double R, double z, double phi, int M);
+    std::shared_ptr<DiskDensityFunc> pyDens;
 
   protected:
 
@@ -837,7 +840,7 @@ namespace BasisClasses
     {
       auto [ret, worst, lworst] = orthoCompute(sl->orthoCheck());
       // For the CTest log
-      std::cout << "Cylindrical::orthoTest: worst=" << worst << std::endl;
+      std::cout << "---- Cylindrical::orthoTest: worst=" << worst << std::endl;
       return ret;
     }
   };
@@ -981,7 +984,7 @@ namespace BasisClasses
 	}
       }
 
-      std::cout << "Slab::orthoTest: worst=" << worst << std::endl;
+      std::cout << "---- Slab::orthoTest: worst=" << worst << std::endl;
       if (worst > __EXP__::orthoTol) return false;
       return true;
     }
@@ -1103,7 +1106,7 @@ namespace BasisClasses
 	}
       }
 
-      std::cout << "Cube::orthoTest: worst=" << worst << std::endl;
+      std::cout << "---- Cube::orthoTest: worst=" << worst << std::endl;
       if (worst > __EXP__::orthoTol) return false;
       return true;
     }

expui/BiorthBasis.cc

@@ -60,6 +60,7 @@ namespace BasisClasses
     "self_consistent",
     "cachename",
     "modelname",
+    "pyname",
     "rnum"
   };
 
@@ -860,7 +861,9 @@ namespace BasisClasses
       {"gaussian",    DiskType::gaussian},
       {"mn",          DiskType::mn},
       {"exponential", DiskType::exponential},
-      {"doubleexpon", DiskType::doubleexpon}
+      {"doubleexpon", DiskType::doubleexpon},
+      {"diskbulge",   DiskType::diskbulge},
+      {"python",      DiskType::python}
     };
 
   const std::set<std::string>
@@ -911,6 +914,8 @@ namespace BasisClasses
     "aratio",
     "hratio",
     "dweight",
+    "Mfac",
+    "HERNA",
     "rwidth",
     "rfactor",
     "rtrunc",
@@ -921,7 +926,8 @@ namespace BasisClasses
     "self_consistent",
     "playback",
     "coefCompute",
-    "coefMaster"
+    "coefMaster",
+    "pyname"
   };
 
   Cylindrical::Cylindrical(const YAML::Node& CONF) :
@@ -980,6 +986,22 @@ namespace BasisClasses
 	  w2*exp(-R/a2)/(4.0*M_PI*a2*a2*h2*f2*f2) ;
       }
       break;
+
+    case DiskType::diskbulge:
+      {
+	double f  = cosh(z/hcyl);
+	double rr = pow(pow(R, 2) + pow(z,2), 0.5);
+	double w1 = Mfac;
+	double w2 = 1.0 - Mfac;
+	double as = HERNA;
+
+	ans = w1*exp(-R/acyl)/(4.0*M_PI*acyl*acyl*hcyl*f*f) + 
+	  w2*pow(as, 4)/(2.0*M_PI*rr)*pow(rr+as,-3.0) ;
+      }
+      break;
+    case DiskType::python:
+      ans = (*pyDens)(R, z, phi);
+      break;
     case DiskType::exponential:
     default:
       {
@@ -1068,6 +1090,12 @@ namespace BasisClasses
     // Ratio of second disk relative to the first disk for disk basis construction with double-exponential
     dweight     = 1.0;              
 
+    // mass fraction for disk for diskbulge
+    Mfac      = 1.0; 
+
+    // Hernquist scale a disk basis construction with diskbulge
+    HERNA      = 0.10; 
+
     // Width for erf truncation for EOF conditioning density (ignored if zero)
     rwidth      = 0.0;             
 
@@ -1124,16 +1152,19 @@ namespace BasisClasses
       if (conf["dmodel"    ])      dmodel = conf["dmodel"    ].as<bool>();
 
       if (conf["aratio"    ])      aratio = conf["aratio"    ].as<double>();
-      if (conf["hratio"    ])      aratio = conf["hratio"    ].as<double>();
-      if (conf["dweight"   ])      aratio = conf["dweight"   ].as<double>();
-      if (conf["rwidth"    ])      aratio = conf["rwidth"    ].as<double>();
-      if (conf["ashift"    ])      aratio = conf["ashift"    ].as<double>();
+      if (conf["hratio"    ])      hratio = conf["hratio"    ].as<double>();
+      if (conf["dweight"   ])      dweight = conf["dweight"   ].as<double>();
+      if (conf["Mfac"   ])         Mfac = conf["Mfac"   ].as<double>();
+      if (conf["HERNA"   ])        HERNA = conf["HERNA"   ].as<double>();
+      if (conf["rwidth"    ])      rwidth = conf["rwidth"    ].as<double>();
+      if (conf["ashift"    ])      ashift = conf["ashift"    ].as<double>();
       if (conf["rfactor"   ])     rfactor = conf["rfactor"   ].as<double>();
       if (conf["rtrunc"    ])      rtrunc = conf["rtrunc"    ].as<double>();
       if (conf["pow"       ])        ppow = conf["ppow"      ].as<double>();
       if (conf["mtype"     ])       mtype = conf["mtype"     ].as<std::string>();
       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>();
 
       // Deprecation warning
       if (conf["density"   ]) {
@@ -1261,7 +1292,7 @@ namespace BasisClasses
 	DTYPE = dtlookup.at(dtype);	// key is not in the map.
 
 	if (myid==0)		// Report DiskType
-	  std::cout << "DiskType is <" << dtype << ">" << std::endl;
+	  std::cout << "---- DiskType is <" << dtype << ">" << std::endl;
       }
       catch (const std::out_of_range& err) {
 	if (myid==0) {
@@ -1273,6 +1304,12 @@ namespace BasisClasses
 	throw std::runtime_error("Cylindrical::initialize: invalid DiskType");
       }
 
+      // Check for and initialize the Python density type
+      //
+      if (DTYPE == DiskType::python) {
+	pyDens = std::make_shared<DiskDensityFunc>(pyname);
+      }
+
       // Use these user models to deproject for the EOF spherical basis
       //
       if (deproject) {

exputil/CMakeLists.txt

@@ -35,13 +35,14 @@ set(QPDISTF_SRC QPDistF.cc qld.c)
 set(SLEDGE_SRC sledge.f)
 set(PARTICLE_SRC Particle.cc ParticleReader.cc header.cc)
 set(CUDA_SRC cudaParticle.cu cudaSLGridMP2.cu)
+set(PYWRAP_SRC DiskDensityFunc.cc)
 
 set(exputil_SOURCES ${ODE_SRC} ${ROOT_SRC} ${QUAD_SRC}
   ${RANDOM_SRC} ${UTIL_SRC} ${SPECFUNC_SRC}
   ${PHASE_SRC}  ${SYMP_SRC} ${INTERP_SRC} ${MASSMODEL_SRC}
   ${ORBIT_SRC} ${BIORTH_SRC} ${POLY_SRC} ${GAUSS_SRC}
   ${QPDISTF_SRC} ${BESSEL_SRC} ${OPTIMIZATION_SRC}
-  ${SLEDGE_SRC} ${PARTICLE_SRC} ${CUDA_SRC})
+  ${SLEDGE_SRC} ${PARTICLE_SRC} ${CUDA_SRC} ${PYWRAP_SRC})
 
 set(common_INCLUDE_DIRS $<INSTALL_INTERFACE:include>
   $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/include/>
@@ -53,7 +54,7 @@ set(common_INCLUDE_DIRS $<INSTALL_INTERFACE:include>
 
 set(common_LINKLIB ${DEP_LIB} OpenMP::OpenMP_CXX MPI::MPI_CXX yaml-cpp
   ${VTK_LIBRARIES} ${HDF5_CXX_LIBRARIES} ${HDF5_LIBRARIES}
-  ${HDF5_HL_LIBRARIES} ${FFTW_DOUBLE_LIB})
+  ${HDF5_HL_LIBRARIES} ${FFTW_DOUBLE_LIB} pybind11::embed)
 
 if(ENABLE_CUDA)
   list(APPEND common_LINKLIB CUDA::toolkit CUDA::cudart)

exputil/DiskDensityFunc.cc

@@ -0,0 +1,31 @@
+#include <DiskDensityFunc.H>
+
+namespace py = pybind11;
+
+DiskDensityFunc::DiskDensityFunc(const std::string& modulename,
+				 const std::string& funcname)
+  : funcname(funcname)
+{
+  // Check if a Python interpreter exists
+  if (Py_IsInitialized() == 0) {
+    py::initialize_interpreter();
+    // Mark the interpreter as started by this instance
+    started = true;
+  }
+
+  // Bind the disk_density function from Python
+  disk_density =
+    py::reinterpret_borrow<py::function>
+    ( py::module::import(modulename.c_str()).attr(funcname.c_str()) );
+}
+
+DiskDensityFunc::~DiskDensityFunc()
+{
+  // Only end the interpreter if it was started by this instance
+  if (started) py::finalize_interpreter();
+}
+
+double DiskDensityFunc::operator() (double R, double z, double phi)
+{
+  return disk_density(R, z, phi).cast<double>();
+}

include/DiskDensityFunc.H

@@ -0,0 +1,61 @@
+#ifndef _DiskDensityFunc_H_
+#define _DiskDensityFunc_H_
+
+#include <functional>
+
+#include <pybind11/pybind11.h>
+#include <pybind11/embed.h>
+
+/** Python disk-density function wrapper
+
+    The wrapper class will take a user supplied Python module that
+    must supply the 'disk-density' function
+
+    An example disk-density Python module:
+    --------------------------cut here--------------------------------
+    import math
+
+    def disk_density(R, z, phi):
+        """Computes the disk density at a given radius, height, and
+        azimuth.  This would be the user's target density for basis
+        creation.
+
+        """
+        a = 1.0                         # Scale radius
+        h = 0.2                         # Scale height
+        f = math.exp(-math.fabs(z)/h)   # Prevent overflows
+        sech = 2.0*f / (1.0 + f*f)      # 
+        return math.exp(-R/a)*sech*sech/(4*math.pi*h*a*a)
+	
+    --------------------------cut here--------------------------------
+
+ */
+class __attribute__((visibility("default")))
+DiskDensityFunc
+{
+private:
+
+  //! The disk-density function
+  pybind11::function disk_density;
+
+  //! Interpreter started?
+  bool started = false;
+
+  //! Name of density function
+  std::string funcname;
+
+public:
+
+  //! Constructor
+  DiskDensityFunc(const std::string& modulename,
+		  const std::string& funcname="disk_density");
+
+  //! Destructor
+  ~DiskDensityFunc();
+
+  //! Evaluate the density
+  double operator() (double R, double z, double phi);
+
+};
+
+#endif

include/DiskModels.H

@@ -2,7 +2,7 @@
 #define _DISK_MODELS_H
 
 #include <EmpCylSL.H>
-
+#include <DiskDensityFunc.H>
 
 //! A Ferrers Ellipsoid + Evacuated Exponential Disc (semi-realistic
 //! bar+disk model)
@@ -223,5 +223,58 @@ public:
 
 };
 
-#endif
 
+//! The usual exponential disk + a Hernquist bulge
+class diskbulge : public EmpCylSL::AxiDisk
+{
+private:
+  double a, h, as, Mfac, d1, d2;
+
+public:
+
+  diskbulge(double a, double h, double as, double Mfac, double M=1) :
+    a(a), h(h), as(as), Mfac(Mfac), EmpCylSL::AxiDisk(M, "diskbulge")
+  {
+    params.push_back(a);
+    params.push_back(h);
+    params.push_back(as);
+    params.push_back(Mfac);
+  }
+
+  double operator()(double R, double z, double phi=0.)
+  {
+    double sh = 1.0/cosh(z/h);
+    d1 = 0.25*M*Mfac/(M_PI*a*a*h) * exp(-R/a) * sh * sh;
+
+    double rr = pow(pow(R, 2) + pow(z,2), 0.5);
+    d2 = M*(1-Mfac)*pow(as, 4)/(2.0*M_PI*rr)*pow(rr+as,-3.0);
+
+    return d1 + d2;
+
+  }
+
+};
+
+
+//! A user-defined Python model
+class PyModel : public EmpCylSL::AxiDisk
+{
+private:
+
+  std::shared_ptr<DiskDensityFunc> pyDens;
+
+public:
+
+  PyModel(std::string& pyname)
+  {
+    pyDens = std::make_shared<DiskDensityFunc>(pyname);
+  }
+
+  double operator()(double R, double z, double phi=0.)
+  {
+    return (*pyDens)(R, z, phi);
+  }
+
+};
+
+#endif