From b5789498b77ff998821e5954c165b069a55e806e Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Thu, 26 Dec 2024 10:53:54 -0500
Subject: [PATCH] Missing instantiation in 'initial.cc'

---
 src/Cylinder.cc      |  9 +++++++--
 utils/ICs/initial.cc | 16 ++++++++++++++--
 2 files changed, 21 insertions(+), 4 deletions(-)

diff --git a/src/Cylinder.cc b/src/Cylinder.cc
index ad369eb6..98773a41 100644
--- 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
index 7cbc9f5b..f0ad4b6d 100644
--- 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,7 +240,7 @@ 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},
@@ -247,10 +248,12 @@ std::map<std::string, DiskType> dtlookup =
     {"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;