Merge branch 'Toomre' of github.com:EXP-code/EXP into Toomre

Author: Martin D. Weinberg

expui/BiorthBasis.cc

@@ -1198,7 +1198,7 @@ namespace BasisClasses
       {"python",      DiskType::python}
     };
 
-  // Dprojection model for cylindrical basis construction
+  // Deprojection model for cylindrical basis construction
   const std::map<std::string, Cylindrical::DeprojType> Cylindrical::dplookup =
     { {"mn",          DeprojType::mn},
       {"exponential", DeprojType::exponential},
@@ -1503,7 +1503,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<string>();
+      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>();
@@ -1711,6 +1711,10 @@ namespace BasisClasses
 	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 {
@@ -1724,7 +1728,7 @@ namespace BasisClasses
 	}
 	catch (const std::out_of_range& err) {
 	  if (myid==0) {
-	    std::cout << "DeprojType error in configuraton file" << std::endl;
+	    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;
@@ -1735,7 +1739,7 @@ namespace BasisClasses
 	if (PTYPE == DeprojType::mn) // Miyamoto-Nagai
 	  model = std::make_shared<MNdisk>(1.0, H);
 	else if (PTYPE == DeprojType::toomre) {
-	  model = std::make_shared<Toomre>(1.0, H);
+	  model = std::make_shared<Toomre>(1.0, H, 5.0);
 	} else if (PTYPE == DeprojType::python and 
 		   DTYPE == DiskType::python) {
 	  model = std::make_shared<AxiSymPyModel>(pyname, acyl);

include/DiskModels.H

@@ -219,9 +219,8 @@ public:
   double operator()(double R, double z, double phi=0.)
   {
     double sigma = norm * std::pow(1.0 + (R/a)*(R/a), -0.5*n);
-    double Z = std::fabs(z);
-    double Q = std::exp(-0.5*Z/h);
-    double sech = 2.0*Q/(1.0 + Q*Q);
+    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;
   }
 };

utils/Test/Deprojector.cc

@@ -46,10 +46,14 @@ namespace Deproject
     for (size_t i=0;i<pairs.size();++i) { R[i] = pairs[i].first; S[i] = pairs[i].second; }
 
     // ensure positive radii and strictly increasing
+    double eps = 1e-12;
     if (R.front() <= 0.0) {
-      double eps = 1e-12;
-      if (R.front() <= 0.0) R[0] = eps;
-      for (size_t i=1;i<R.size();++i) if (R[i] <= R[i-1]) R[i] = R[i-1] + eps;
+      R[0] = eps;
+    }
+    for (size_t i=1;i<R.size();++i) {
+      if (R[i] <= R[i-1]) {
+        R[i] = R[i-1] + eps;
+      }
     }
 
     spline_.set_data(R, S);

utils/Test/EmpDeproj.cc

@@ -104,18 +104,20 @@ EmpDeproj::EmpDeproj(double H, double RMIN, double RMAX, int NUMR, int NINT,
 
   // Debug
   //
-  if (true) {
+#ifdef EMPDEPROJ_DEBUG
+  {
     std::string fname("deproject_sl.txt");
     std::ofstream out(fname);
     if (out) {
       for (int i=0; i<NUMR; i++)
-	out << std::setw(18) << rl[i]
-	    << std::setw(18) << rr[i]
-	    << std::setw(18) << rho[i]
-	    << std::setw(18) << mass[i]
-	    << std::endl;
+        out << std::setw(18) << rl[i]
+            << std::setw(18) << rr[i]
+            << std::setw(18) << rho[i]
+            << std::setw(18) << mass[i]
+            << std::endl;
     }
   }
+#endif
 
   // Finalize
   //

utils/Test/testDeprojPlummer.cc

@@ -64,12 +64,13 @@ int main()
   auto rho = D.rho(r_eval);
 
   std::ofstream ofs("rho_test.txt");
-  for (size_t i = 0; i < r_eval.size(); ++i)
+  for (size_t i = 0; i < r_eval.size(); ++i) {
     ofs << std::setw(16) << r_eval[i]
 	<< std::setw(16) << rho[i]
 	<< std::setw(16) << RhoFunc(r_eval[i])
 	<< std::endl;
-      ofs.close();
+  }
+  ofs.close();
   std::cout << "Wrote rho_test.txt\n";
 
   return 0;

utils/Test/testDeproject.cc

@@ -18,7 +18,7 @@ int main()
       double t = (double)i / (Ndata - 1);
       double r = Rmin + t * (Rmax - Rmin);
       Rdata.push_back(r);
-      Sigma.push_back(1.0 / std::pow(1.0 + r*r, -1.5));
+      Sigma.push_back(std::pow(1.0 + r*r, -1.5));
     }
 
     Deprojector D(Rdata, Sigma, /*R_max_extend=*/50.0, /*tail_power=*/-4.0, /*Ngrid=*/6000);
@@ -39,7 +39,7 @@ int main()
 
   // Example B: construct from analytic functor + analytic derivative
   {
-    auto SigmaFunc = [](double R)->double { return 1.0 / std::pow(1.0 + R*R, -1.5); };
+    auto SigmaFunc = [](double R)->double { return std::pow(1.0 + R*R, -1.5); };
     auto dSigmaFunc = [](double R)->double { return -3.0 * R / std::pow(1.0 + R*R, -2.5); }; // analytic derivative
 
     Deprojector D(SigmaFunc, dSigmaFunc, /*R_data_min=*/0.01, /*R_data_max=*/10.0,

utils/Test/testEmp.cc

@@ -3,6 +3,7 @@
 #include <fstream>
 #include <vector>
 #include <cmath>
+#include <stdexcept>
 
 #include "Deprojector.H"
 #include "EmpDeproj.H"