Add AbsorptionCoefficient and Emissivity to MeanOpacity variants.

+ Use gmode input to toggle Rosseland (default) to Planck.
+ Add PlanckDistribution function as member to MeanOpacity.
+ Add the new functions to generic and non-CGS MeanOpacity variants.

Note: the motivation is to use one of the frequency-integrated
opacity values as a grey opacity in transport, which requires using
it in both absorption and emissivity, for any choice.

singularity-opac/photons/mean_opacity_photons.hpp

@@ -29,6 +29,7 @@
 
 #include <singularity-opac/photons/mean_photon_variant.hpp>
 #include <singularity-opac/photons/non_cgs_photons.hpp>
+#include <singularity-opac/photons/thermal_distributions_photons.hpp>
 
 namespace singularity {
 namespace photons {
@@ -143,6 +144,25 @@ class MeanOpacity {
     return rho * fromLog_(lkappa_.interpToReal(lRho, lT, Rosseland));
   }
 
+  // standard grey opacity functions
+  PORTABLE_INLINE_FUNCTION
+  Real AbsorptionCoefficient(const Real rho, const Real temp,
+                            const int gmode = Rosseland) const {
+    Real lRho = toLog_(rho);
+    Real lT = toLog_(temp);
+    return rho * fromLog_(lkappa_.interpToReal(lRho, lT, gmode));
+  }
+
+  PORTABLE_INLINE_FUNCTION
+  Real Emissivity(const Real rho, const Real temp,
+                  const int gmode = Rosseland,
+                  Real *lambda = nullptr) const {
+    Real B = dist_.ThermalDistributionOfT(temp, lambda);
+    Real lRho = toLog_(rho);
+    Real lT = toLog_(temp);
+    return rho * fromLog_(lkappa_.interpToReal(lRho, lT, gmode)) * B;
+  }
+
  private:
   template <typename Opacity, bool AUTOFREQ>
   void MeanOpacityImpl_(const Opacity &opac, const Real lRhoMin,
@@ -305,6 +325,7 @@ class MeanOpacity {
     Rosseland = 0,
     Planck = 1
   };
+  PlanckDistribution<PC> dist_;
 };
 
 #undef EPS

singularity-opac/photons/mean_photon_variant.hpp

@@ -92,6 +92,27 @@ class MeanVariant {
         opac_);
   }
 
+  PORTABLE_INLINE_FUNCTION
+  Real AbsorptionCoefficient(const Real rho, const Real temp,
+                            const int gmode = 0) const {
+    return mpark::visit(
+        [=](const auto &opac) {
+          return opac.AbsorptionCoefficient(rho, temp, gmode);
+        },
+        opac_);
+  }
+
+  PORTABLE_INLINE_FUNCTION
+  Real Emissivity(const Real rho, const Real temp,
+                  const int gmode = 0,
+                  Real *lambda = nullptr) const {
+    return mpark::visit(
+        [=](const auto &opac) {
+          return opac.Emissivity(rho, temp, gmode, lambda);
+        },
+        opac_);
+  }
+
   inline void Finalize() noexcept {
     return mpark::visit([](auto &opac) { return opac.Finalize(); }, opac_);
   }

singularity-opac/photons/non_cgs_photons.hpp

@@ -246,7 +246,8 @@ class MeanNonCGSUnits {
                   const Real temp_unit)
       : mean_opac_(std::forward<MeanOpac>(mean_opac)), time_unit_(time_unit),
         mass_unit_(mass_unit), length_unit_(length_unit), temp_unit_(temp_unit),
-        rho_unit_(mass_unit_ / (length_unit_ * length_unit_ * length_unit_)) {}
+        rho_unit_(mass_unit_ / (length_unit_ * length_unit_ * length_unit_)),
+        inv_emiss_unit_(length_unit_ * time_unit_ * time_unit_ / mass_unit_) {}
 
   auto GetOnDevice() {
     return MeanNonCGSUnits<MeanOpac>(mean_opac_.GetOnDevice(), time_unit_,
@@ -286,6 +287,21 @@ class MeanNonCGSUnits {
     return alpha * length_unit_;
   }
 
+  PORTABLE_INLINE_FUNCTION
+  Real AbsorptionCoefficient(const Real rho, const Real temp,
+                            const int gmode = 0) const {
+    const Real alpha = mean_opac_.AbsorptionCoefficient(rho, temp, gmode);
+    return alpha * length_unit_;
+  }
+
+  PORTABLE_INLINE_FUNCTION
+  Real Emissivity(const Real rho, const Real temp,
+                  const int gmode = 0,
+                  Real *lambda = nullptr) const {
+    const Real J = mean_opac_.Emissivity(rho, temp, gmode);
+    return J * inv_emiss_unit_ * time_unit_;
+  }
+
   PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
   GetRuntimePhysicalConstants() const {
     return RuntimePhysicalConstants(PhysicalConstantsCGS(), time_unit_,
@@ -295,7 +311,7 @@ class MeanNonCGSUnits {
  private:
   MeanOpac mean_opac_;
   Real time_unit_, mass_unit_, length_unit_, temp_unit_;
-  Real rho_unit_;
+  Real rho_unit_, inv_emiss_unit_;
 };
 
 } // namespace photons

test/test_mean_opacities.cpp

@@ -815,6 +815,33 @@ TEST_CASE("ASCII-parsed Mean photon opacities", "[MeanPhotons]") {
         n_wrong_d() += 1;
       }
 
+      // test absorption and emission in Rossland and Planck modes
+      Real mabs_ross = mean_opac.AbsorptionCoefficient(rho_max, temp_max, 0);
+      Real mabs_plnk = mean_opac.AbsorptionCoefficient(rho_max, temp_max, 1);
+
+      // compare to Rossland and Planck
+      if (FractionalDifference(mross, mabs_ross) > EPS_TEST) {
+        n_wrong_d() += 1;
+      }
+      if (FractionalDifference(mplnk, mabs_plnk) > EPS_TEST) {
+        n_wrong_d() += 1;
+      }
+
+      Real *lambda = nullptr;
+      singularity::photons::PlanckDistribution<PhysicalConstantsCGS> dist;
+      Real B = dist.ThermalDistributionOfT(temp_max, lambda);
+
+      Real memiss_ross = mean_opac.Emissivity(rho_max, temp_max, 0);
+      Real memiss_plnk = mean_opac.Emissivity(rho_max, temp_max, 0);
+
+      // compare to Rossland and Planck
+      if (FractionalDifference(mross * B, memiss_ross) > EPS_TEST) {
+        n_wrong_d() += 1;
+      }
+      if (FractionalDifference(mplnk * B, memiss_plnk) > EPS_TEST) {
+        n_wrong_d() += 1;
+      }
+
 #ifdef PORTABILITY_STRATEGY_KOKKOS
       Kokkos::deep_copy(n_wrong_h, n_wrong_d);
 #endif