diff --git a/src/analysis/history.cpp b/src/analysis/history.cpp
index 5b044939..d88e3d65 100644
--- a/src/analysis/history.cpp
+++ b/src/analysis/history.cpp
@@ -12,10 +12,12 @@
 // publicly, and to permit others to do so.
 
 #include "history.hpp"
+#include "analysis/analysis.hpp"
 #include "geometry/geometry.hpp"
 #include "geometry/geometry_utils.hpp"
 #include "history_utils.hpp"
 #include "phoebus_utils/relativity_utils.hpp"
+#include <interface/sparse_pack.hpp>
 
 namespace History {
 
@@ -235,4 +237,46 @@ Real ReduceMagneticFluxPhi(MeshData<Real> *md) {
   return 0.5 * result; // 0.5 \int detg B^r dx2 dx3
 } // Phi
 
+// SN analysis
+// ReduceLocalizationFunction is not used currently. However this function returns
+// normalization function that is used for localizing quantities at the center, or at some
+// particular case. For example SN diagnostics oftec computes quantities at 400 km.
+void ReduceLocalizationFunction(MeshData<Real> *md) {
+  const auto ib = md->GetBoundsI(IndexDomain::interior);
+  const auto jb = md->GetBoundsJ(IndexDomain::interior);
+  const auto kb = md->GetBoundsK(IndexDomain::interior);
+  namespace diag = diagnostic_variables;
+  using parthenon::MakePackDescriptor;
+  auto *pmb = md->GetParentPointer();
+  Mesh *pmesh = md->GetMeshPointer();
+  auto &resolved_pkgs = pmesh->resolved_packages;
+  const int ndim = pmesh->ndim;
+
+  static auto desc = MakePackDescriptor<diag::localization_function>(resolved_pkgs.get());
+  auto v = desc.GetPack(md);
+  const int nblocks = v.GetNBlocks();
+  auto geom = Geometry::GetCoordinateSystem(md);
+
+  parthenon::par_for(
+      parthenon::LoopPatternMDRange(), "Central Density for SN", DevExecSpace(), 0,
+      nblocks - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
+      KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
+        auto &coords = v.GetCoordinates(b);
+        auto analysis = pmb->packages.Get("analysis").get();
+        const Real x[3] = {coords.Xc<1>(k, j, i), coords.Xc<2>(k, j, i),
+                           coords.Xc<3>(k, j, i)};
+        const Real sigma = analysis->Param<Real>("sigma");
+        Real gam[3][3];
+        Real r2 = 0;
+        geom.Metric(CellLocation::Cent, 0, k, j, i, gam);
+        for (int n = 0; n < 3; ++n) {
+          for (int m = 0; m < 3; ++m) {
+            r2 += gam[n][m] * x[n] * x[m];
+          }
+        }
+        v(b, diag::localization_function(), k, j, i) = std::exp(-r2 / sigma / sigma);
+      });
+
+} // exp
+
 } // namespace History
diff --git a/src/analysis/history.hpp b/src/analysis/history.hpp
index 04a61d60..8dbfb164 100644
--- a/src/analysis/history.hpp
+++ b/src/analysis/history.hpp
@@ -17,12 +17,14 @@
 #include <string>
 #include <vector>
 
+#include "geometry/geometry.hpp"
+#include "geometry/geometry_utils.hpp"
+#include "phoebus_utils/variables.hpp"
 #include <kokkos_abstraction.hpp>
+#include <parthenon/driver.hpp>
 #include <parthenon/package.hpp>
 #include <utils/error_checking.hpp>
 
-#include "phoebus_utils/variables.hpp"
-
 using namespace parthenon::package::prelude;
 
 /*
@@ -41,6 +43,7 @@ Real ReduceMassAccretionRate(MeshData<Real> *md);
 Real ReduceJetEnergyFlux(MeshData<Real> *md);
 Real ReduceJetMomentumFlux(MeshData<Real> *md);
 Real ReduceMagneticFluxPhi(MeshData<Real> *md);
+void ReduceLocalizationFunction(MeshData<Real> *md);
 
 template <typename Reducer_t>
 Real ReduceOneVar(MeshData<Real> *md, const std::string &varname, int idx = 0) {
@@ -58,11 +61,14 @@ Real ReduceOneVar(MeshData<Real> *md, const std::string &varname, int idx = 0) {
   // We choose to apply volume weighting when using the sum reduction.
   // This assumes that for "Sum" variables, the variable is densitized, meaning
   // it already contains a factor of the measure: sqrt(det(gamma))
-  const bool volume_weighting =
-      std::is_same<Reducer_t, Kokkos::Sum<Real, HostExecSpace>>::value;
-
   Real result = 0.0;
   Reducer_t reducer(result);
+
+  const bool volume_weighting =
+      std::is_same<Reducer_t, Kokkos::Sum<Real, HostExecSpace>>::value ||
+      std::is_same<Reducer_t, Kokkos::Sum<Real, DevExecSpace>>::value ||
+      std::is_same<Reducer_t, Kokkos::Sum<Real, Kokkos::DefaultExecutionSpace>>::value;
+
   parthenon::par_reduce(
       parthenon::LoopPatternMDRange(), "Phoebus History for " + varname, DevExecSpace(),
       0, pack.GetDim(5) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
@@ -78,6 +84,44 @@ Real ReduceOneVar(MeshData<Real> *md, const std::string &varname, int idx = 0) {
   return result;
 }
 
+template <typename Varname>
+Real ReduceInGain(MeshData<Real> *md, int idx = 0) {
+  const auto ib = md->GetBoundsI(IndexDomain::interior);
+  const auto jb = md->GetBoundsJ(IndexDomain::interior);
+  const auto kb = md->GetBoundsK(IndexDomain::interior);
+
+  namespace iv = internal_variables;
+  using parthenon::MakePackDescriptor;
+  Mesh *pmesh = md->GetMeshPointer();
+  auto &resolved_pkgs = pmesh->resolved_packages;
+  const int ndim = pmesh->ndim;
+  static auto desc =
+      MakePackDescriptor<Varname, iv::GcovHeat, iv::GcovCool>(resolved_pkgs.get());
+  auto v = desc.GetPack(md);
+  PARTHENON_REQUIRE_THROWS(v.ContainsHost(0, iv::GcovHeat(), iv::GcovCool()),
+                           "Must be doing SN simulation");
+  const int nblocks = v.GetNBlocks();
+
+  Real result = 0.0;
+
+  auto geom = Geometry::GetCoordinateSystem(md);
+
+  parthenon::par_reduce(
+      parthenon::LoopPatternMDRange(),
+      "Phoebus History for integrals in gain region (SN)", DevExecSpace(), 0, nblocks - 1,
+      kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
+      KOKKOS_LAMBDA(const int b, const int k, const int j, const int i, Real &lresult) {
+        Real gdet = geom.DetGamma(CellLocation::Cent, 0, k, j, i);
+        bool is_netheat =
+            (v(b, iv::GcovHeat(), k, j, i) - v(b, iv::GcovCool(), k, j, i) > 0);
+        const auto &coords = v.GetCoordinates(b);
+        const Real vol = coords.CellVolume(k, j, i);
+        lresult += gdet * vol * is_netheat * v(b, Varname(idx), k, j, i);
+      },
+      Kokkos::Sum<Real>(result));
+  return result;
+}
+
 } // namespace History
 
 #endif // ANALYSIS_HISTORY_HPP_
diff --git a/src/fluid/fluid.cpp b/src/fluid/fluid.cpp
index d8264a13..cddff16f 100644
--- a/src/fluid/fluid.cpp
+++ b/src/fluid/fluid.cpp
@@ -13,6 +13,7 @@
 
 #include "fluid.hpp"
 
+#include "analysis/analysis.hpp"
 #include "analysis/history.hpp"
 #include "con2prim.hpp"
 #include "con2prim_robust.hpp"
@@ -207,6 +208,8 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   physics->AddField(p::entropy::name(), mprim_scalar);
   physics->AddField(p::cs::name(), mprim_scalar);
   physics->AddField(diag::ratio_divv_cs::name(), mprim_scalar);
+  // physics->AddField(diag::localization_function::name(), mprim_scalar);
+  // (MG) currently not in use, turn on when needed.
   physics->AddField(diag::entropy_z_0::name(), mprim_scalar);
   physics->AddField(p::gamma1::name(), mprim_scalar);
   if (ye) {
@@ -309,6 +312,8 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   // Reductions
   // By default compute integrated value of scalar conserved vars
   auto HstSum = parthenon::UserHistoryOperation::sum;
+  auto HstMax = parthenon::UserHistoryOperation::max;
+  using History::ReduceInGain;
   using History::ReduceOneVar;
   using parthenon::HistoryOutputVar;
   parthenon::HstVar_list hst_vars = {};
@@ -319,6 +324,11 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   auto ReduceEn = [](MeshData<Real> *md) {
     return ReduceOneVar<Kokkos::Sum<Real>>(md, fluid_cons::energy::name(), 0);
   };
+  auto MaxDensity = [](MeshData<Real> *md) {
+    return ReduceOneVar<Kokkos::Max<Real>>(md, fluid_prim::density::name(), 0);
+  };
+
+  hst_vars.emplace_back(HistoryOutputVar(HstMax, MaxDensity, "maximum density"));
   hst_vars.emplace_back(HistoryOutputVar(HstSum, ReduceMass, "total baryon number"));
   hst_vars.emplace_back(HistoryOutputVar(HstSum, ReduceEn, "total conserved energy tau"));
 
diff --git a/src/pgen/progenitor.cpp b/src/pgen/progenitor.cpp
index 093343d5..48f076fe 100644
--- a/src/pgen/progenitor.cpp
+++ b/src/pgen/progenitor.cpp
@@ -16,18 +16,20 @@
 #include <string>
 #include <vector>
 
+#include "analysis/history.hpp"
 #include "geometry/geometry.hpp"
 #include "microphysics/eos_phoebus/eos_phoebus.hpp"
 #include "monopole_gr/monopole_gr.hpp"
 #include "pgen/pgen.hpp"
 #include "progenitor/progenitordata.hpp"
+#include <parthenon/package.hpp>
+#include <utils/error_checking.hpp>
 
 using DataBox = Spiner::DataBox<Real>;
 
 // Homologously collapsing star.
 namespace homologous {
 void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
-  std::cout << "Entered into ProbGen" << std::endl;
   const bool is_monopole_cart =
       (typeid(PHOEBUS_GEOMETRY) == typeid(Geometry::MonopoleCart));
   const bool is_monopole_sph =
diff --git a/src/phoebus_utils/variables.hpp b/src/phoebus_utils/variables.hpp
index 5f675402..2c127ea2 100644
--- a/src/phoebus_utils/variables.hpp
+++ b/src/phoebus_utils/variables.hpp
@@ -122,6 +122,7 @@ namespace diagnostic_variables {
 VARIABLE_NONS(divb);
 VARIABLE_NONS(ratio_divv_cs);
 VARIABLE_NONS(entropy_z_0);
+VARIABLE_NONS(localization_function);
 VARIABLE_CUSTOM(divf, flux_divergence);
 VARIABLE_NONS(src_terms);
 VARIABLE_CUSTOM(r_divf, r.flux_divergence);
diff --git a/src/progenitor/progenitordata.cpp b/src/progenitor/progenitordata.cpp
index 3fdac85f..92043be8 100644
--- a/src/progenitor/progenitordata.cpp
+++ b/src/progenitor/progenitordata.cpp
@@ -1,6 +1,7 @@
 #include <memory>
 #include <vector>
 
+#include "analysis/history.hpp"
 #include "ascii_reader.hpp"
 #include "geometry/geometry.hpp"
 #include "microphysics/eos_phoebus/eos_phoebus.hpp"
@@ -107,6 +108,24 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   params.Add("S_adm_dev", S_adm_dev);
   params.Add("Srr_adm_dev", Srr_adm_dev);
 
+  // Reductions
+  auto HstSum = parthenon::UserHistoryOperation::sum;
+  auto HstMax = parthenon::UserHistoryOperation::max;
+  using History::ReduceInGain;
+  using History::ReduceOneVar;
+  using parthenon::HistoryOutputVar;
+  parthenon::HstVar_list hst_vars = {};
+  auto Mgain = [](MeshData<Real> *md) {
+    return ReduceInGain<class fluid_prim::density>(md, 0);
+  };
+  auto Qgain = [](MeshData<Real> *md) {
+    return ReduceInGain<class internal_variables::GcovHeat>(md, 0) -
+           ReduceInGain<class internal_variables::GcovCool>(md, 0);
+  };
+  hst_vars.emplace_back(HistoryOutputVar(HstSum, Mgain, "Mgain"));
+  hst_vars.emplace_back(HistoryOutputVar(HstSum, Qgain, "total net heat"));
+  params.Add(parthenon::hist_param_key, hst_vars);
+
   return progenitor_pkg;
 } // Initialize
 
diff --git a/src/radiation/cooling_function.cpp b/src/radiation/cooling_function.cpp
index b6171ec6..06b39a26 100644
--- a/src/radiation/cooling_function.cpp
+++ b/src/radiation/cooling_function.cpp
@@ -279,10 +279,8 @@ TaskStatus CoolingFunctionCalculateFourForce(MeshData<Real> *rc, const double dt
             // detg included above
             Kokkos::atomic_add(&(v(b, iv::Gcov(mu), k, j, i)), -Gcov_coord[mu]);
           }
-          v(b, iv::GcovHeat(), k, j, i) =
-              v(b, p::density(), k, j, i) * density_conversion_factor * heat;
-          v(b, iv::GcovCool(), k, j, i) =
-              v(b, p::density(), k, j, i) * density_conversion_factor * cool;
+          v(b, iv::GcovHeat(), k, j, i) = cdensity * H;
+          v(b, iv::GcovCool(), k, j, i) = cdensity * C;
           Kokkos::atomic_add(&(v(b, iv::Gye(), k, j, i)), Jye);
         });
 #else