From e515f7c3dc848a88c6053aa2d97e0c2bf9e87b21 Mon Sep 17 00:00:00 2001
From: Philipp Grete <pgrete@hs.uni-hamburg.de>
Date: Tue, 17 Oct 2023 15:24:26 +0200
Subject: [PATCH] Add support for calculating lin and log bin edges

---
 src/outputs/histogram.cpp                     | 31 +++++++++++++++++--
 .../output_hdf5/parthinput.advection          |  6 ++--
 2 files changed, 32 insertions(+), 5 deletions(-)

diff --git a/src/outputs/histogram.cpp b/src/outputs/histogram.cpp
index 94158b55635e..046dcaec3fbc 100644
--- a/src/outputs/histogram.cpp
+++ b/src/outputs/histogram.cpp
@@ -64,13 +64,38 @@ namespace HistUtil {
 
 ParArray1D<Real> GetEdges(ParameterInput *pin, const std::string &block_name,
                           const std::string &prefix) {
-
   std::vector<Real> edges_in;
 
   const auto edge_type_str = pin->GetString(block_name, prefix + "type");
-  if (edge_type_str == "lin") {
+  if (edge_type_str == "lin" || edge_type_str == "log") {
+    const auto edge_min = pin->GetReal(block_name, prefix + "min");
+    const auto edge_max = pin->GetReal(block_name, prefix + "max");
+    PARTHENON_REQUIRE_THROWS(edge_max > edge_min,
+                             "Histogram max needs to be larger than min.")
+
+    const auto edge_num_bins = pin->GetReal(block_name, prefix + "num_bins");
+    PARTHENON_REQUIRE_THROWS(edge_num_bins >= 1, "Need at least one bin for histogram.");
+
+    if (edge_type_str == "lin") {
+      auto dbin = (edge_max - edge_min) / (edge_num_bins);
+      for (int i = 0; i < edge_num_bins; i++) {
+        edges_in.emplace_back(edge_min + i * dbin);
+      }
+      edges_in.emplace_back(edge_max);
+    } else if (edge_type_str == "log") {
+      PARTHENON_REQUIRE_THROWS(
+          edge_min > 0.0 && edge_max > 0.0,
+          "Log binning for negative values not implemented. However, you can specify "
+          "arbitrary bin edges through the 'list' edge type.")
 
-  } else if (edge_type_str == "log") {
+      auto dbin = (std::log10(edge_max) - std::log10(edge_min)) / (edge_num_bins);
+      for (int i = 0; i < edge_num_bins; i++) {
+        edges_in.emplace_back(std::pow(10., std::log10(edge_min) + i * dbin));
+      }
+      edges_in.emplace_back(edge_max);
+    } else {
+      PARTHENON_FAIL("Not sure how I got here...")
+    }
 
   } else if (edge_type_str == "list") {
     edges_in = pin->GetVector<Real>(block_name, prefix + "list");
diff --git a/tst/regression/test_suites/output_hdf5/parthinput.advection b/tst/regression/test_suites/output_hdf5/parthinput.advection
index 677a2ebf2e18..347ff6d75644 100644
--- a/tst/regression/test_suites/output_hdf5/parthinput.advection
+++ b/tst/regression/test_suites/output_hdf5/parthinput.advection
@@ -88,8 +88,10 @@ hist0_binned_variable_component = 0
 # 2D histogram of a variable, binned by a coordinate and a different variable
 hist1_ndim = 2
 hist1_x_variable = HIST_COORD_X1
-hist1_x_edges_type = list
-hist1_x_edges_list = -0.5, -0.25, 0.0, 0.25, 0.5
+hist1_x_edges_type = lin
+hist1_x_edges_num_bins = 4
+hist1_x_edges_min = -0.5
+hist1_x_edges_max = 0.5
 hist1_y_variable = one_minus_advected_sq
 hist1_y_variable_component = 0
 hist1_y_edges_type = list