Added a spherical histogram1d with logarithmic scaling to pyEXP; sinc…

…e this is a new function, there is no API change
EXP-code · Dec 17, 2024 · febd3b7 · febd3b7
1 parent dba7ade
commit febd3b7
Show file tree

Hide file tree

Showing 3 changed files with 83 additions and 0 deletions.
diff --git a/expui/FieldGenerator.H b/expui/FieldGenerator.H
@@ -115,6 +115,11 @@ namespace Field
     histogram1d(PR::PRptr reader, double rmax, int nbins, std::string proj,
 		std::vector<double> center={0.0, 0.0, 0.0});
 
+    //! Compute spherical log histogram from particles
+    std::tuple<Eigen::VectorXf, Eigen::VectorXf>
+    histo1dlog(PR::PRptr reader, double rmin, double rmax, int nbins,
+		std::vector<double> center={0.0, 0.0, 0.0});
+
     //! Write field slices to files.  This will be VTK your build is
     //! compiled with VTK and ascii tables otherwise.
     void file_slices(BasisClasses::BasisPtr basis, CoefClasses::CoefsPtr coefs,

diff --git a/expui/FieldGenerator.cc b/expui/FieldGenerator.cc
@@ -894,6 +894,58 @@ namespace Field
   }
   // END histogram1d
 
+  std::tuple<Eigen::VectorXf, Eigen::VectorXf>
+  FieldGenerator::histo1dlog(PR::PRptr reader, double rmin, double rmax,
+			     int nbins, std::vector<double> ctr)
+  {
+    if (rmin <= 0.0)  throw std::runtime_error("FieldGenerator::histo1dlog: rmin must be > 0.0");
+    if (rmax <= rmin) throw std::runtime_error("FieldGenerator::histo1dlog: rmax must be > rmin");
+
+    const double pi = 3.14159265358979323846;
+
+    Eigen::VectorXf rad = Eigen::VectorXf::Zero(nbins);
+    Eigen::VectorXf ret = Eigen::VectorXf::Zero(nbins);
+
+    double lrmin = log(rmin), lrmax = log(rmax);
+    double del = (lrmax - lrmin)/nbins;
+
+    // Make the histogram
+    //
+    for (auto p=reader->firstParticle(); p!=0; p=reader->nextParticle()) {
+      double rad = 0.0;
+      for (int k=0; k<3; k++) {
+	double pp = p->pos[k] - ctr[k];
+	rad += pp*pp;
+      }
+
+      int indx = floor((log(sqrt(rad)) - lrmin)/del);
+      if (indx>=0 and indx<nbins) ret[indx] += p->mass;
+    }
+
+    // Accumulate between MPI nodes; return value to root node
+    //
+    if (use_mpi) {
+      if (myid==0) 
+	MPI_Reduce(MPI_IN_PLACE, ret.data(), ret.size(),
+		   MPI_FLOAT, MPI_SUM, 0, MPI_COMM_WORLD);
+      else
+	MPI_Reduce(ret.data(), NULL, ret.size(),
+		   MPI_FLOAT, MPI_SUM, 0, MPI_COMM_WORLD);
+    }
+
+    // Compute density
+    //
+    double d3 = exp(3.0*del);
+    double rf = 4.0*pi/3.0*(d3 - 1.0);
+    for (int i=0; i<nbins; i++) {
+      rad[i] = exp(lrmin + del*(0.5 + i));
+      ret[i] /= exp(3.0*(lrmin + del*i)) * rf;
+    }
+
+    return {rad, ret};
+  }
+  // END histo1dlog
+
   std::map<double, std::map<std::string, Eigen::VectorXf>>
   FieldGenerator::points(BasisClasses::BasisPtr basis,
 			 CoefClasses::CoefsPtr coefs)

diff --git a/pyEXP/FieldWrappers.cc b/pyEXP/FieldWrappers.cc
@@ -304,6 +304,32 @@ void FieldGeneratorClasses(py::module &m) {
 	py::arg("projection"),
 	py::arg("center") = std::vector<double>(3, 0.0));
 
+  f.def("histo1dlog", &Field::FieldGenerator::histo1dlog,
+	R"(
+        Make a 1d density histogram (array) in spherical projection with
+        logarithmic scaling
+
+        Parameters
+        ----------
+        reader : ParticleReader
+            particle reader instance
+        rmin : float
+            minimum radius of the histogram (>0)
+        rmax : float
+            linear extent of the histogram
+        nbins : int
+            number of bins
+        center : list(float, float, float), default=[0, 0, 0]
+            origin for computing the histogram
+
+        Returns
+        -------
+        tuple(numpy.ndarray, numpy.ndarray)
+            the evaluation grid and computed 1d histogram
+        )",
+	py::arg("reader"), py::arg("rmin"), py::arg("rmax"), py::arg("nbins"),
+	py::arg("center") = std::vector<double>(3, 0.0));
+
   f.def("file_lines", &Field::FieldGenerator::file_lines,
 	R"(
         Write field arrays to files using the supplied string prefix.