Missing instantiation in 'initial.cc'

EXP-code · Dec 26, 2024 · b578949 · b578949
1 parent 86c8c4e
commit b578949
Show file tree

Hide file tree

Showing 2 changed files with 21 additions and 4 deletions.
diff --git a/src/Cylinder.cc b/src/Cylinder.cc
@@ -265,13 +265,18 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) :
       std::shared_ptr<DiskDensityFunc> dfunc;
       if (pyname.size()) dfunc = std::make_shared<DiskDensityFunc>(pyname);
 
+      // Use either the user-supplied Python target or the default
+      // exponential disk model
       auto DiskDens = [&](double R, double z, double phi)
       {
 	if (dfunc) return (*dfunc)(R, z, phi);
-	double f = cosh(z/hcyl);
-	return exp(-R/acyl)/(4.0*M_PI*acyl*acyl*hcyl*f*f);
+	double f = exp(-fabs(z)/hcyl); // Overflow prevention
+	double s = 2.0*f/(1.0 + f*f);  // in sech computation
+	return exp(-R/acyl)*s*s/(4.0*M_PI*acyl*acyl*hcyl);
       };
 
+      // The conditioning function for the EOF with an optional shift
+      // for M>0
       auto dcond = [&](double R, double z, double phi, int M)
       {
 	//

diff --git a/utils/ICs/initial.cc b/utils/ICs/initial.cc
@@ -92,6 +92,7 @@
 #include <numerical.H>
 #include <gaussQ.H>
 #include <isothermal.H>
+#include <DiskDensityFunc.H>
 #include <hernquist_model.H>
 #include <model3d.H>
 #include <biorth.H>
@@ -239,18 +240,20 @@ const double f_H = 0.76;
 
 // Global variables
 //
-enum DiskType { constant, gaussian, mn, exponential, doubleexpon, diskbulge };
+enum DiskType { constant, gaussian, mn, exponential, doubleexpon, diskbulge, python };
 
 std::map<std::string, DiskType> dtlookup =
   { {"constant",    DiskType::constant},
     {"gaussian",    DiskType::gaussian},
     {"mn",          DiskType::mn},
     {"exponential", DiskType::exponential},
     {"doubleexpon", DiskType::doubleexpon},
-    {"diskbulge",   DiskType::diskbulge}
+    {"diskbulge",   DiskType::diskbulge},
+    {"python",      DiskType::python}
   };
 
 DiskType     DTYPE;
+std::string  pyname;
 double       ASCALE;
 double       ASHIFT;
 double       HSCALE;
@@ -267,6 +270,7 @@ double       HERNA = 0.10;
 double DiskDens(double R, double z, double phi)
 {
   double ans = 0.0;
+  static std::shared_ptr<DiskDensityFunc> dfunc;
 
   switch (DTYPE) {
 
@@ -322,6 +326,12 @@ double DiskDens(double R, double z, double phi)
             w2*pow(as, 4)/(2.0*M_PI*rr)*pow(rr+as,-3.0);
     }
     break;
+  case DiskType::python:
+    {
+      if (!dfunc) dfunc = std::make_shared<DiskDensityFunc>(pyname);
+      ans = (*dfunc)(R, z, phi);
+    }
+    break;
   case DiskType::exponential:
   default:
     {
@@ -598,6 +608,8 @@ main(int ac, char **av)
     ("allow", "Allow multimass algorithm to generature negative masses for testing")
     ("nomono", "Allow non-monotonic mass interpolation")
     ("diskmodel", "Table describing the model for the disk plane")
+    ("pyname", "Name of module with the user-specified target disk density",
+     cxxopts::value<std::string>(pyname))
     ;
 
   cxxopts::ParseResult vm;