diff --git a/CMakeLists.txt b/CMakeLists.txt
index e34cb89..b34b014 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -53,8 +53,8 @@ target_include_directories(singularity-opac::flags
                            INTERFACE
                            $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}>)
 
-include (SetupDeps)
 include (SetupOptions)
+include (SetupDeps)
 include (SetupCompilers)
 include (SetupFlags)
 
diff --git a/README.md b/README.md
index 6a34156..7f24b2a 100644
--- a/README.md
+++ b/README.md
@@ -57,6 +57,21 @@ A number of options are avaialable for compiling:
 | --------------------------------- | ------- | ------------------------------------------------------------------------------------ |
 | SINGULARITY_BUILD_TESTS           | OFF     | Build test infrastructure.                                                           |
 | SINGULARITY_USE_HDF5              | ON      | Enables HDF5. Required for Spiner opacities.                                         |
+| SINGULARITY_KOKKOS_IN_TREE        | OFF     | Force cmake to use Kokkos source included in tree.                                   |
+
+### Loading ASCII Data
+
+Currently, the MeanOpacity class, defined in singularity-opac/photons/mean_opacity_photons.hpp, supports
+loading grey Rosseland and Planck opacity data in an ASCII format.  An example of this format is
+provided by singularity-opac/photons/example_ascii/kap_plaw.txt.  The 1st row of the header has the
+number of density points, NRho, then the number of temperature points, NT.  The 2nd (3rd) row of the header
+has min and max density (temperature) bounds.  These bounds are inclusive, so the opacity data in the file
+should have evaluations at these min and max values.  The rest of the ASCII file is a two-column table, where
+the 1st (2nd) column is Rosseland (Planck) opacity.  The number of rows in each column is NRhoxNT, where
+density is the slow index and temperature is the fast index (thus the row index = temperature index
++ NT x (density index), indexing from 0).  Each opacity is assumed to be evaluated on log-spaced
+density and temperature grids, where these grids are defined by NRho, NT, and the (again inclusive) min and
+max bounds the header.
 
 ## Copyright
 
diff --git a/cmake/SetupDeps.cmake b/cmake/SetupDeps.cmake
index 16455f5..99e1ba0 100644
--- a/cmake/SetupDeps.cmake
+++ b/cmake/SetupDeps.cmake
@@ -10,6 +10,25 @@ else()
   message(status "CUDA::toolkit provided by parent package")
 endif()
 
+#=======================================
+# Setup Kokkos
+# - provides Kokkos::kokkos
+#=======================================
+if (SINGULARITY_USE_KOKKOS)
+  if (NOT TARGET Kokkos::kokkos)
+    message(status "Kokkos::kokkos must be found")
+    if (SINGULARITY_KOKKOS_IN_TREE)
+      message(status "Using in-tree Kokkos")
+      add_subdirectory(utils/kokkos)
+    else()
+      message(status "Using system Kokkos if available")
+      find_package(Kokkos REQUIRED)
+    endif()
+  else()
+    message(status "Kokkos::kokkos provided by parent package")
+  endif()
+endif()
+
 #=======================================
 # Setup ports of call
 # - provides PortsofCall::PortsofCall
@@ -17,16 +36,6 @@ endif()
 find_package(PortsofCall REQUIRED)
 target_link_libraries(singularity-opac::flags INTERFACE PortsofCall::PortsofCall)
 
-#=======================================
-# Setup Kokkos
-# - provides Kokkos::kokkos
-#=======================================
-if (NOT TARGET Kokkos::kokkos)
-  find_package(Kokkos QUIET)
-else()
-  message(status "Kokkos::kokkos provided by parent package")
-endif()
-
 #=======================================
 # Find HDF5
 # - cmake@3.20+ provides HDF5::HDF5, but
diff --git a/cmake/SetupOptions.cmake b/cmake/SetupOptions.cmake
index 21a857c..9c9e53d 100644
--- a/cmake/SetupOptions.cmake
+++ b/cmake/SetupOptions.cmake
@@ -17,11 +17,18 @@ option (SINGULARITY_USE_FMATH "Enable fast-math logarithms" ON)
 # Dependency options
 #=======================================
 # check for using in-tree third-party dependencies
-option (SINULARITY_KOKKOS_IN_TREE "Use in-tree dependencies" OFF)
-
 option (SINGULARITY_USE_KOKKOS "Use Kokkos for portability" OFF)
 option (SINGULARITY_USE_CUDA "Enable cuda support" OFF)
 option (SINGULARITY_USE_HDF5 "Pull in hdf5" OFF)
+cmake_dependent_option(SINULARITY_KOKKOS_IN_TREE
+  "Use in-tree dependencies" OFF
+  ${SINGULARITY_USE_KOKKOS} OFF)
+
+# Kill cmake's package registry because it can interfere
+if (SINGULARITY_KOKKOS_IN_TREE)
+  set(CMAKE_FIND_PACKAGE_NO_PACKAGE_REGISTRY ON CACHE BOOL "" FORCE)
+  set(CMAKE_FIND_PACKAGE_NO_SYSTEM_PACKAGE_REGISTRY ON CACHE BOOL "" FORCE)
+endif()
 
 # If the conditional is TRUE, it's the first default, else it's the
 # second.
diff --git a/singularity-opac/constants/constants.hpp b/singularity-opac/constants/constants.hpp
index 580037a..d08bb4a 100644
--- a/singularity-opac/constants/constants.hpp
+++ b/singularity-opac/constants/constants.hpp
@@ -1,5 +1,5 @@
 // ======================================================================
-// © 2021. Triad National Security, LLC. All rights reserved.  This
+// © 2021-2024. Triad National Security, LLC. All rights reserved.  This
 // program was produced under U.S. Government contract
 // 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
 // is operated by Triad National Security, LLC for the U.S.
@@ -194,6 +194,47 @@ struct PhysicalConstants {
   static constexpr Real gA = axial_vector_coupling_constant;
 };
 
+struct RuntimePhysicalConstants {
+  template <typename T>
+  PORTABLE_INLINE_FUNCTION
+  RuntimePhysicalConstants(T pc)
+      : na(pc.na), alpha(pc.alpha), h(pc.h), hbar(pc.hbar), kb(pc.kb),
+        r_gas(pc.r_gas), qe(pc.qe), c(pc.c), g_newt(pc.g_newt),
+        g_accel(pc.g_accel), me(pc.me), mp(pc.mp), mn(pc.mn), amu(pc.amu),
+        sb(pc.sb), ar(pc.ar), faraday(pc.faraday), mu0(pc.mu0), eps0(pc.eps0),
+        eV(pc.eV), Fc(pc.Fc), nu_sigma0(pc.nu_sigma0), gA(pc.gA) {}
+
+  const Real na;
+  const Real alpha;
+  const Real h;
+  const Real hbar;
+  const Real kb;
+  const Real r_gas;
+  const Real qe;
+  const Real c;
+  const Real g_newt;
+  const Real g_accel;
+  const Real me;
+  const Real mp;
+  const Real mn;
+  const Real amu;
+  const Real sb;
+  const Real ar;
+  const Real faraday;
+  const Real mu0;
+  const Real eps0;
+  const Real eV;
+  const Real Fc;
+  const Real nu_sigma0;
+  const Real gA;
+};
+
+template <typename T>
+PORTABLE_INLINE_FUNCTION
+RuntimePhysicalConstants GetRuntimePhysicalConstants(T phys_constants) {
+  return RuntimePhysicalConstants(phys_constants);
+}
+
 using PhysicalConstantsUnity =
     PhysicalConstants<BaseUnity, UnitConversionDefault>;
 using PhysicalConstantsSI =
diff --git a/singularity-opac/neutrinos/brt_neutrinos.hpp b/singularity-opac/neutrinos/brt_neutrinos.hpp
index 5721021..9138084 100644
--- a/singularity-opac/neutrinos/brt_neutrinos.hpp
+++ b/singularity-opac/neutrinos/brt_neutrinos.hpp
@@ -1,5 +1,5 @@
 // ======================================================================
-// © 2021. Triad National Security, LLC. All rights reserved.  This
+// © 2021-2024. Triad National Security, LLC. All rights reserved.  This
 // program was produced under U.S. Government contract
 // 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
 // is operated by Triad National Security, LLC for the U.S.
@@ -32,6 +32,8 @@ namespace neutrinos {
 template <typename ThermalDistribution, typename pc = PhysicalConstantsCGS>
 class BRTOpacity {
  public:
+  using PC = pc;
+
   BRTOpacity() = default;
   BRTOpacity(const ThermalDistribution &dist) : dist_(dist) {}
   BRTOpacity GetOnDevice() { return *this; }
diff --git a/singularity-opac/neutrinos/gray_opacity_neutrinos.hpp b/singularity-opac/neutrinos/gray_opacity_neutrinos.hpp
index fc9a57e..830615e 100644
--- a/singularity-opac/neutrinos/gray_opacity_neutrinos.hpp
+++ b/singularity-opac/neutrinos/gray_opacity_neutrinos.hpp
@@ -1,5 +1,5 @@
 // ======================================================================
-// © 2021. Triad National Security, LLC. All rights reserved.  This
+// © 2021-2024. Triad National Security, LLC. All rights reserved.  This
 // program was produced under U.S. Government contract
 // 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
 // is operated by Triad National Security, LLC for the U.S.
@@ -30,6 +30,8 @@ namespace neutrinos {
 template <typename ThermalDistribution, typename pc = PhysicalConstantsCGS>
 class GrayOpacity {
  public:
+  using PC = pc;
+
   GrayOpacity() = default;
   GrayOpacity(const Real kappa) : kappa_(kappa) {}
   GrayOpacity(const ThermalDistribution &dist, const Real kappa)
diff --git a/singularity-opac/neutrinos/mean_neutrino_variant.hpp b/singularity-opac/neutrinos/mean_neutrino_variant.hpp
index fd50081..488085b 100644
--- a/singularity-opac/neutrinos/mean_neutrino_variant.hpp
+++ b/singularity-opac/neutrinos/mean_neutrino_variant.hpp
@@ -1,5 +1,5 @@
 // ======================================================================
-// © 2021. Triad National Security, LLC. All rights reserved.  This
+// © 2021-2024. Triad National Security, LLC. All rights reserved.  This
 // program was produced under U.S. Government contract
 // 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
 // is operated by Triad National Security, LLC for the U.S.
@@ -69,6 +69,16 @@ class MeanVariant {
         opac_);
   }
 
+  PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
+  GetRuntimePhysicalConstants() const {
+    return mpark::visit(
+        [=](const auto &opac) {
+          using PC = typename std::decay_t<decltype(opac)>::PC;
+          return singularity::GetRuntimePhysicalConstants(PC());
+        },
+        opac_);
+  }
+
   PORTABLE_INLINE_FUNCTION Real PlanckMeanAbsorptionCoefficient(
       const Real rho, const Real temp, const Real Ye,
       const RadiationType type) const {
@@ -91,6 +101,12 @@ class MeanVariant {
   inline void Finalize() noexcept {
     return mpark::visit([](auto &opac) { return opac.Finalize(); }, opac_);
   }
+
+#ifdef SPINER_USE_HDF
+  void Save(const std::string &filename) const {
+    return mpark::visit([=](auto &opac) { return opac.Save(filename); }, opac_);
+  }
+#endif
 };
 
 } // namespace impl
diff --git a/singularity-opac/neutrinos/mean_opacity_neutrinos.hpp b/singularity-opac/neutrinos/mean_opacity_neutrinos.hpp
index 61a6396..52e9485 100644
--- a/singularity-opac/neutrinos/mean_opacity_neutrinos.hpp
+++ b/singularity-opac/neutrinos/mean_opacity_neutrinos.hpp
@@ -1,5 +1,5 @@
 // ======================================================================
-// © 2022. Triad National Security, LLC. All rights reserved.  This
+// © 2022-2024. Triad National Security, LLC. All rights reserved.  This
 // program was produced under U.S. Government contract
 // 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
 // is operated by Triad National Security, LLC for the U.S.
@@ -37,7 +37,6 @@ namespace impl {
 // TODO(BRR) Note: It is assumed that lambda is constant for all densities,
 // temperatures, and Ye
 
-template <typename pc = PhysicalConstantsCGS>
 class MeanOpacity {
 
  public:
@@ -134,6 +133,8 @@ class MeanOpacity {
                         const Real lTMax, const int NT, const Real YeMin,
                         const Real YeMax, const int NYe, Real lNuMin,
                         Real lNuMax, const int NNu, Real *lambda = nullptr) {
+    using PC = typename Opacity::PC;
+
     lkappaPlanck_.resize(NRho, NT, NYe, NEUTRINO_NTYPES);
     // index 0 is the species and is not interpolatable
     lkappaPlanck_.setRange(1, YeMin, YeMax, NYe);
@@ -159,8 +160,8 @@ class MeanOpacity {
             // Choose default temperature-specific frequency grid if frequency
             // grid not specified
             if (AUTOFREQ) {
-              lNuMin = toLog_(1.e-3 * pc::kb * fromLog_(lTMin) / pc::h);
-              lNuMax = toLog_(1.e3 * pc::kb * fromLog_(lTMax) / pc::h);
+              lNuMin = toLog_(1.e-3 * PC::kb * fromLog_(lTMin) / PC::h);
+              lNuMax = toLog_(1.e3 * PC::kb * fromLog_(lTMax) / PC::h);
             }
             const Real dlnu = (lNuMax - lNuMin) / (NNu - 1);
             // Integrate over frequency
@@ -219,10 +220,9 @@ class MeanOpacity {
 
 } // namespace impl
 
-using MeanOpacityScaleFree = impl::MeanOpacity<PhysicalConstantsUnity>;
-using MeanOpacityCGS = impl::MeanOpacity<PhysicalConstantsCGS>;
-using MeanOpacity = impl::MeanVariant<MeanOpacityScaleFree, MeanOpacityCGS,
-                                      MeanNonCGSUnits<MeanOpacityCGS>>;
+using MeanOpacityBase = impl::MeanOpacity;
+using MeanOpacity =
+    impl::MeanVariant<MeanOpacityBase, MeanNonCGSUnits<MeanOpacityBase>>;
 
 } // namespace neutrinos
 } // namespace singularity
diff --git a/singularity-opac/neutrinos/neutrino_variant.hpp b/singularity-opac/neutrinos/neutrino_variant.hpp
index ae14525..5736077 100644
--- a/singularity-opac/neutrinos/neutrino_variant.hpp
+++ b/singularity-opac/neutrinos/neutrino_variant.hpp
@@ -1,5 +1,5 @@
 // ======================================================================
-// © 2021. Triad National Security, LLC. All rights reserved.  This
+// © 2021-2024. Triad National Security, LLC. All rights reserved.  This
 // program was produced under U.S. Government contract
 // 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
 // is operated by Triad National Security, LLC for the U.S.
@@ -71,6 +71,16 @@ class Variant {
         opac_);
   }
 
+  PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
+  GetRuntimePhysicalConstants() const {
+    return mpark::visit(
+        [=](const auto &opac) {
+          using PC = typename std::decay_t<decltype(opac)>::PC;
+          return singularity::GetRuntimePhysicalConstants(PC());
+        },
+        opac_);
+  }
+
   // Directional absorption coefficient with units of 1/length
   // Signature should be at least
   // rho, temp, Ye, type, nu, lambda
diff --git a/singularity-opac/neutrinos/non_cgs_neutrinos.hpp b/singularity-opac/neutrinos/non_cgs_neutrinos.hpp
index 7a760c7..400cb2b 100644
--- a/singularity-opac/neutrinos/non_cgs_neutrinos.hpp
+++ b/singularity-opac/neutrinos/non_cgs_neutrinos.hpp
@@ -1,5 +1,5 @@
 // ======================================================================
-// © 2021. Triad National Security, LLC. All rights reserved.  This
+// © 2021-2024. Triad National Security, LLC. All rights reserved.  This
 // program was produced under U.S. Government contract
 // 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
 // is operated by Triad National Security, LLC for the U.S.
@@ -29,6 +29,8 @@ namespace neutrinos {
 template <typename Opac>
 class NonCGSUnits {
  public:
+  using PC = typename Opac::PC;
+
   NonCGSUnits() = default;
   NonCGSUnits(Opac &&opac, const Real time_unit, const Real mass_unit,
               const Real length_unit, const Real temp_unit)
@@ -47,6 +49,7 @@ class NonCGSUnits {
     return NonCGSUnits<Opac>(opac_.GetOnDevice(), time_unit_, mass_unit_,
                              length_unit_, temp_unit_);
   }
+
   inline void Finalize() noexcept { opac_.Finalize(); }
 
   PORTABLE_INLINE_FUNCTION
@@ -66,10 +69,11 @@ class NonCGSUnits {
   }
 
   template <typename FrequencyIndexer, typename DataIndexer>
-  PORTABLE_INLINE_FUNCTION void AbsorptionCoefficient(
-      const Real rho, const Real temp, const Real Ye, RadiationType type,
-      FrequencyIndexer &nu_bins, DataIndexer &coeffs, const int nbins,
-      Real *lambda = nullptr) const {
+  PORTABLE_INLINE_FUNCTION void
+  AbsorptionCoefficient(const Real rho, const Real temp, const Real Ye,
+                        RadiationType type, FrequencyIndexer &nu_bins,
+                        DataIndexer &coeffs, const int nbins,
+                        Real *lambda = nullptr) const {
     for (int i = 0; i < nbins; ++i) {
       nu_bins[i] *= freq_unit_;
     }
@@ -262,6 +266,12 @@ class MeanNonCGSUnits {
   PORTABLE_INLINE_FUNCTION
   int nlambda() const noexcept { return mean_opac_.nlambda(); }
 
+#ifdef SPINER_USE_HDF
+  void Save(const std::string &filename) const {
+    return mean_opac_.Save(filename);
+  }
+#endif
+
   PORTABLE_INLINE_FUNCTION
   Real PlanckMeanAbsorptionCoefficient(const Real rho, const Real temp,
                                        const Real Ye,
diff --git a/singularity-opac/neutrinos/spiner_opac_neutrinos.hpp b/singularity-opac/neutrinos/spiner_opac_neutrinos.hpp
index 6ee5bb0..039125f 100644
--- a/singularity-opac/neutrinos/spiner_opac_neutrinos.hpp
+++ b/singularity-opac/neutrinos/spiner_opac_neutrinos.hpp
@@ -1,5 +1,5 @@
 // ======================================================================
-// © 2021. Triad National Security, LLC. All rights reserved.  This
+// © 2021-2024. Triad National Security, LLC. All rights reserved.  This
 // program was produced under U.S. Government contract
 // 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
 // is operated by Triad National Security, LLC for the U.S.
@@ -63,7 +63,9 @@ enum class DataStatus { Deallocated, OnDevice, OnHost };
 template <typename ThermalDistribution, typename pc = PhysicalConstantsCGS>
 class SpinerOpacity {
  public:
+  using PC = pc;
   using DataBox = Spiner::DataBox<Real>;
+
   static constexpr Real MEV = 1e6 * pc::eV;
   static constexpr Real EPS = 10.0 * std::numeric_limits<Real>::min();
   static constexpr Real Hz2MeV = pc::h / (1e6 * pc::eV);
diff --git a/singularity-opac/neutrinos/tophat_emissivity_neutrinos.hpp b/singularity-opac/neutrinos/tophat_emissivity_neutrinos.hpp
index ee7b9b1..bfb024b 100644
--- a/singularity-opac/neutrinos/tophat_emissivity_neutrinos.hpp
+++ b/singularity-opac/neutrinos/tophat_emissivity_neutrinos.hpp
@@ -1,5 +1,5 @@
 // ======================================================================
-// © 2021. Triad National Security, LLC. All rights reserved.  This
+// © 2021-2024. Triad National Security, LLC. All rights reserved.  This
 // program was produced under U.S. Government contract
 // 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
 // is operated by Triad National Security, LLC for the U.S.
@@ -32,6 +32,8 @@ namespace neutrinos {
 template <typename ThermalDistribution, typename pc = PhysicalConstantsCGS>
 class TophatEmissivity {
  public:
+  using PC = pc;
+
   TophatEmissivity(const Real C, const Real numin, const Real numax)
       : C_(C), numin_(numin), numax_(numax) {}
   TophatEmissivity(const ThermalDistribution &dist, const Real C,
@@ -46,6 +48,7 @@ class TophatEmissivity {
     printf("Tophat emissivity. C, numin, numax = %g, %g, %g\n", C_, numin_,
            numax_);
   }
+
   inline void Finalize() noexcept {}
 
   PORTABLE_INLINE_FUNCTION
diff --git a/singularity-opac/photons/epbremsstrahlung_opacity_photons.hpp b/singularity-opac/photons/epbremsstrahlung_opacity_photons.hpp
index e4e9ccc..725681d 100644
--- a/singularity-opac/photons/epbremsstrahlung_opacity_photons.hpp
+++ b/singularity-opac/photons/epbremsstrahlung_opacity_photons.hpp
@@ -1,5 +1,5 @@
 // ======================================================================
-// © 2022. Triad National Security, LLC. All rights reserved.  This
+// © 2022-2024. Triad National Security, LLC. All rights reserved.  This
 // program was produced under U.S. Government contract
 // 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
 // is operated by Triad National Security, LLC for the U.S.
@@ -32,6 +32,8 @@ namespace photons {
 template <typename pc = PhysicalConstantsCGS>
 class EPBremsstrahlungOpacity {
  public:
+  using PC = pc;
+
   EPBremsstrahlungOpacity() = default;
   EPBremsstrahlungOpacity(const PlanckDistribution<pc> &dist) : dist_(dist) {}
 
@@ -42,6 +44,10 @@ class EPBremsstrahlungOpacity {
   void PrintParams() const noexcept {
     printf("Electron-proton bremsstrahlung opacity.\n");
   }
+
+  PORTABLE_INLINE_FUNCTION
+  pc GetPhysicalConstants() const { return pc(); }
+
   inline void Finalize() noexcept {}
 
   PORTABLE_INLINE_FUNCTION
diff --git a/singularity-opac/photons/example_ascii/kap_plaw.txt b/singularity-opac/photons/example_ascii/kap_plaw.txt
new file mode 100644
index 0000000..59568c1
--- /dev/null
+++ b/singularity-opac/photons/example_ascii/kap_plaw.txt
@@ -0,0 +1,19 @@
+ NRho 4 NT 4
+ rho_min 1.00000000e-14 rho_max 7.94328235e-01
+ T_min 1.00000000e+00 T_max 7.94328235e+06
+1.00000000e-03 1.00000000e-01
+1.00000000e-03 1.00000000e-01
+1.00000000e-03 1.00000000e-01
+1.00000000e-03 1.00000000e-01
+1.00000000e-03 1.00000000e-01
+1.00000000e-03 1.00000000e-01
+1.00000000e-03 1.00000000e-01
+1.00000000e-03 1.00000000e-01
+1.00000000e-03 1.00000000e-01
+1.00000000e-03 1.00000000e-01
+1.00000000e-03 1.00000000e-01
+1.00000000e-03 1.00000000e-01
+1.00000000e-03 1.00000000e-01
+1.00000000e-03 1.00000000e-01
+1.00000000e-03 1.00000000e-01
+1.00000000e-03 1.00000000e-01
diff --git a/singularity-opac/photons/example_ascii/preproc_ascii_opac.py b/singularity-opac/photons/example_ascii/preproc_ascii_opac.py
new file mode 100644
index 0000000..295c0ad
--- /dev/null
+++ b/singularity-opac/photons/example_ascii/preproc_ascii_opac.py
@@ -0,0 +1,103 @@
+# ======================================================================
+# © 2024. Triad National Security, LLC. All rights reserved.  This
+# program was produced under U.S. Government contract
+# 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
+# is operated by Triad National Security, LLC for the U.S.
+# Department of Energy/National Nuclear Security Administration. All
+# rights in the program are reserved by Triad National Security, LLC,
+# and the U.S. Department of Energy/National Nuclear Security
+# Administration. The Government is granted for itself and others
+# acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+# license in this material to reproduce, prepare derivative works,
+# distribute copies to the public, perform publicly and display
+# publicly, and to permit others to do so.
+# ======================================================================
+
+#--------------------------------------------------------------------------------------------------#
+# Regrid (interpolate) Rosseland and Planck opacity data to a log-log rho-T grid.
+#--------------------------------------------------------------------------------------------------#
+import numpy as np
+from scipy.interpolate import interp2d
+import argparse
+
+#-- parse command line arguments
+parser = argparse.ArgumentParser(description="Re-grid opacity data to be log-log in density-temperature.")
+parser.add_argument("fname", type=str, help="Opacity file to re-grid.")
+parser.add_argument("--is_loglog", action="store_true", help="Avoid regridding if already log-log.")
+parser.add_argument("--fname_new", type=str, default="kappa_new.txt", help="File name for new file created.")
+parser.add_argument("--plot", action="store_true", help="If not log-log, plot colormap of original and regridded data.")
+args = parser.parse_args()
+
+NRho = -1
+NT = -1
+with open(args.fname, "r") as ff:
+    svec = ff.readline().split(" ")
+    NRho = int(svec[2])
+    NT = int(svec[4])
+
+print("NRho = ", NRho)
+print("NT = ", NT)
+
+opac = np.loadtxt(args.fname, skiprows=1)
+
+print("opac.shape = ", opac.shape)
+assert np.size(opac, 0) == NT * NRho, "np.size(opac, 0) != NT * NRho"
+
+#-- density, temperature grids
+Rho = np.unique(opac[:,0])
+T = np.unique(opac[:,1])
+
+assert np.size(Rho) == NRho, "np.size(Rho) != Rho"
+assert np.size(T) == NT, "np.size(T) != NT"
+
+r = np.logspace(np.log10(Rho[0]), np.log10(Rho[NRho - 1]), NRho)
+tt = np.logspace(np.log10(T[0]), np.log10(T[NT - 1]), NT)
+
+if (args.is_loglog):
+    assert np.max(abs(r - Rho) / Rho) < 1e-6, "np.max(abs(r - Rho) / Rho) >= 1e-6"
+    assert np.max(abs(tt - T) / T) < 1e-6, "np.max(abs(tt - T) / T) >= 1e-6"
+else:
+    print()
+    print("Interpolating data to log-log rho-T grid...")
+    #-- reshape to rho-T grid
+    ross_old_opac = np.reshape(opac[:,2], (NRho, NT))
+    plnk_old_opac = np.reshape(opac[:,3], (NRho, NT))
+    #-- linearly interpolate in log-log rho-T
+    ross_f2d = interp2d(np.log10(T), np.log10(Rho), ross_old_opac, kind='linear')
+    plnk_f2d = interp2d(np.log10(T), np.log10(Rho), plnk_old_opac, kind='linear')
+    ross_new_opac = ross_f2d(np.log10(tt), np.log10(r))
+    plnk_new_opac = plnk_f2d(np.log10(tt), np.log10(r))
+    #-- plot 4-panel of data
+    if (args.plot):
+        import matplotlib.pyplot as plt
+        XOLD, YOLD = np.meshgrid(np.log10(T), np.log10(Rho))
+        XNEW, YNEW = np.meshgrid(np.log10(tt), np.log10(r))
+        fig, axes = plt.subplots(2, 2)
+        im1 = axes[0, 0].pcolormesh(XOLD, YOLD, ross_old_opac)
+        axes[0, 0].set_title('Old Rosseland')
+        fig.colorbar(im1, ax=axes[0, 0])
+        im2 = axes[0, 1].pcolormesh(XOLD, YOLD, plnk_old_opac)
+        axes[0, 1].set_title('Old Planck')
+        fig.colorbar(im2, ax=axes[0, 1])
+        im3 = axes[1, 0].pcolormesh(XNEW, YNEW, ross_new_opac)
+        axes[1, 0].set_title('New Rosseland')
+        fig.colorbar(im3, ax=axes[1, 0])
+        im4 = axes[1, 1].pcolormesh(XNEW, YNEW, plnk_new_opac)
+        axes[1, 1].set_title('New Planck')
+        fig.colorbar(im4, ax=axes[1, 1])
+        plt.tight_layout()
+        plt.show()
+    #-- reset opacity array
+    for i in range(NRho):
+        for j in range(NT):
+            k = j + NT * i
+            opac[k,0] = r[i]
+            opac[k,1] = tt[j]
+            opac[k,2] = ross_new_opac[i,j]
+            opac[k,3] = plnk_new_opac[i,j]
+
+#-- save with new header containing min/max rho, T
+hdr = "NRho "+str(NRho)+" NT "+str(NT) + "\n" \
+    "rho_min {:.8e}".format(r[0])+" rho_max {:.8e}".format(r[NRho-1]) + "\n" \
+    "T_min {:.8e}".format(tt[0]) + " T_max {:.8e}".format(tt[NT-1])
+np.savetxt(args.fname_new, opac[:,2:], delimiter=" ", header=hdr, comments=" ", fmt="%.8e")
diff --git a/singularity-opac/photons/gray_opacity_photons.hpp b/singularity-opac/photons/gray_opacity_photons.hpp
index 9bc0358..dcca68c 100644
--- a/singularity-opac/photons/gray_opacity_photons.hpp
+++ b/singularity-opac/photons/gray_opacity_photons.hpp
@@ -1,5 +1,5 @@
 // ======================================================================
-// © 2021. Triad National Security, LLC. All rights reserved.  This
+// © 2021-2024. Triad National Security, LLC. All rights reserved.  This
 // program was produced under U.S. Government contract
 // 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
 // is operated by Triad National Security, LLC for the U.S.
@@ -30,6 +30,8 @@ namespace photons {
 template <typename pc = PhysicalConstantsCGS>
 class GrayOpacity {
  public:
+  using PC = pc;
+
   GrayOpacity() = default;
   GrayOpacity(const Real kappa) : kappa_(kappa) {}
   GrayOpacity(const PlanckDistribution<pc> &dist, const Real kappa)
diff --git a/singularity-opac/photons/mean_opacity_photons.hpp b/singularity-opac/photons/mean_opacity_photons.hpp
index 8a58859..a6a7511 100644
--- a/singularity-opac/photons/mean_opacity_photons.hpp
+++ b/singularity-opac/photons/mean_opacity_photons.hpp
@@ -1,5 +1,5 @@
 // ======================================================================
-// © 2022. Triad National Security, LLC. All rights reserved.  This
+// © 2022-2024. Triad National Security, LLC. All rights reserved.  This
 // program was produced under U.S. Government contract
 // 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
 // is operated by Triad National Security, LLC for the U.S.
@@ -17,6 +17,8 @@
 #define SINGULARITY_OPAC_PHOTONS_MEAN_OPACITY_PHOTONS_
 
 #include <cmath>
+#include <filesystem>
+#include <fstream>
 
 #include <ports-of-call/portability.hpp>
 #include <singularity-opac/base/radiation_types.hpp>
@@ -37,11 +39,9 @@ namespace impl {
 // TODO(BRR) Note: It is assumed that lambda is constant for all densities and
 // temperatures
 
-template <typename pc = PhysicalConstantsCGS>
 class MeanOpacity {
 
  public:
-  using DataBox = Spiner::DataBox<Real>;
   MeanOpacity() = default;
   template <typename Opacity>
   MeanOpacity(const Opacity &opac, const Real lRhoMin, const Real lRhoMax,
@@ -59,27 +59,46 @@ class MeanOpacity {
                                      NT, lNuMin, lNuMax, NNu, lambda);
   }
 
+  // construct Planck/Rosseland DataBox from ascii file
+  MeanOpacity(const std::string &filename)  : filename_(filename.c_str()) {
+
+    // get number of density and temperature points
+    std::ifstream ff(filename.c_str());
+    const bool fexists = ff.good();
+
+    if (fexists) {
+
+      std::filesystem::path filePath(filename);
+      std::string extension = filePath.extension().string();
+
+      if (extension == ".txt") {
+        loadASCII(ff);
 #ifdef SPINER_USE_HDF
-  MeanOpacity(const std::string &filename) : filename_(filename.c_str()) {
-    herr_t status = H5_SUCCESS;
-    hid_t file = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
-    status += lkappaPlanck_.loadHDF(file, SP5::MeanOpac::PlanckMeanOpacity);
-    status +=
-        lkappaRosseland_.loadHDF(file, SP5::MeanOpac::RosselandMeanOpacity);
-    status += H5Fclose(file);
+      } else if (extension == ".hdf5" || extension == ".h5" || extension == ".sp5") {
+        herr_t status = H5_SUCCESS;
+        hid_t file = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
+        status += lkappa_.loadHDF(file, "Rosseland & Planck Mean Opacities");
+        status += H5Fclose(file);
 
-    if (status != H5_SUCCESS) {
-      OPAC_ERROR("photons::MeanOpacity: HDF5 error\n");
+        if (status != H5_SUCCESS) {
+          OPAC_ERROR("photons::MeanOpacity: HDF5 error\n");
+        }
+#endif
+      } else {
+        OPAC_ERROR("photons::MeanOpacity: unrecognized file extension");
+      }
+
+    } else {
+      OPAC_ERROR("photons::MeanOpacity: file does not exist");
     }
   }
 
+#ifdef SPINER_USE_HDF
   void Save(const std::string &filename) const {
     herr_t status = H5_SUCCESS;
     hid_t file =
         H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
-    status += lkappaPlanck_.saveHDF(file, SP5::MeanOpac::PlanckMeanOpacity);
-    status +=
-        lkappaRosseland_.saveHDF(file, SP5::MeanOpac::RosselandMeanOpacity);
+    status += lkappa_.saveHDF(file, "Rosseland & Planck Mean Opacities");
     status += H5Fclose(file);
 
     if (status != H5_SUCCESS) {
@@ -94,21 +113,19 @@ class MeanOpacity {
 
   MeanOpacity GetOnDevice() {
     MeanOpacity other;
-    other.lkappaPlanck_ = Spiner::getOnDeviceDataBox(lkappaPlanck_);
-    other.lkappaRosseland_ = Spiner::getOnDeviceDataBox(lkappaRosseland_);
+    other.lkappa_ = Spiner::getOnDeviceDataBox(lkappa_);
     return other;
   }
 
   void Finalize() {
-    lkappaPlanck_.finalize();
-    lkappaRosseland_.finalize();
+    lkappa_.finalize();
   }
 
   PORTABLE_INLINE_FUNCTION
   Real PlanckMeanAbsorptionCoefficient(const Real rho, const Real temp) const {
     Real lRho = toLog_(rho);
     Real lT = toLog_(temp);
-    return rho * fromLog_(lkappaPlanck_.interpToReal(lRho, lT));
+    return rho * fromLog_(lkappa_.interpToReal(lRho, lT, Planck));
   }
 
   PORTABLE_INLINE_FUNCTION
@@ -116,7 +133,7 @@ class MeanOpacity {
                                           const Real temp) const {
     Real lRho = toLog_(rho);
     Real lT = toLog_(temp);
-    return rho * fromLog_(lkappaRosseland_.interpToReal(lRho, lT));
+    return rho * fromLog_(lkappa_.interpToReal(lRho, lT, Rosseland));
   }
 
  private:
@@ -125,25 +142,26 @@ class MeanOpacity {
                         const Real lRhoMax, const int NRho, const Real lTMin,
                         const Real lTMax, const int NT, Real lNuMin,
                         Real lNuMax, const int NNu, Real *lambda = nullptr) {
-    lkappaPlanck_.resize(NRho, NT);
-    lkappaPlanck_.setRange(0, lTMin, lTMax, NT);
-    lkappaPlanck_.setRange(1, lRhoMin, lRhoMax, NRho);
-    lkappaRosseland_.copyMetadata(lkappaPlanck_);
+    using PC = typename Opacity::PC;
+
+    lkappa_.resize(NRho, NT, 2);
+    lkappa_.setRange(1, lTMin, lTMax, NT);
+    lkappa_.setRange(2, lRhoMin, lRhoMax, NRho);
 
     // Fill tables
     for (int iRho = 0; iRho < NRho; ++iRho) {
-      Real lRho = lkappaPlanck_.range(1).x(iRho);
+      Real lRho = lkappa_.range(2).x(iRho);
       Real rho = fromLog_(lRho);
       for (int iT = 0; iT < NT; ++iT) {
-        Real lT = lkappaPlanck_.range(0).x(iT);
+        Real lT = lkappa_.range(1).x(iT);
         Real T = fromLog_(lT);
         Real kappaPlanckNum = 0.;
         Real kappaPlanckDenom = 0.;
         Real kappaRosselandNum = 0.;
         Real kappaRosselandDenom = 0.;
         if (AUTOFREQ) {
-          lNuMin = toLog_(1.e-3 * pc::kb * fromLog_(lTMin) / pc::h);
-          lNuMax = toLog_(1.e3 * pc::kb * fromLog_(lTMax) / pc::h);
+          lNuMin = toLog_(1.e-3 * PC::kb * fromLog_(lTMin) / PC::h);
+          lNuMax = toLog_(1.e3 * PC::kb * fromLog_(lTMax) / PC::h);
         }
         const Real dlnu = (lNuMax - lNuMin) / (NNu - 1);
         // Integrate over frequency
@@ -172,34 +190,121 @@ class MeanOpacity {
                                   ? singularity_opac::robust::ratio(
                                         kappaRosselandDenom, kappaRosselandNum)
                                   : 0.;
-        lkappaPlanck_(iRho, iT) = toLog_(kappaPlanck);
-        lkappaRosseland_(iRho, iT) = toLog_(kappaRosseland);
-        if (std::isnan(lkappaPlanck_(iRho, iT)) ||
-            std::isnan(lkappaRosseland_(iRho, iT))) {
+
+        lkappa_(iRho, iT, Rosseland) = toLog_(kappaRosseland);
+        lkappa_(iRho, iT, Planck) = toLog_(kappaPlanck);
+        if (std::isnan(lkappa_(iRho, iT, Rosseland)) ||
+            std::isnan(lkappa_(iRho, iT, Planck))) {
           OPAC_ERROR("photons::MeanOpacity: NAN in opacity evaluations");
         }
       }
     }
   }
+
+  // ASCII opacity file reader
+  void loadASCII(std::ifstream &ff) {
+
+    int NRho = -1;
+    int NT = -1;
+
+    // line read from file
+    std::string fline;
+
+    // read 1st line of header to get sizes
+    std::getline(ff, fline);
+
+    // tokenize fline
+    char* cfline = const_cast<char*>(fline.c_str());
+    char* fl_tok = std::strtok(cfline, " ");
+
+    // move to next token to get number of density points
+    fl_tok = std::strtok(nullptr, " ");
+    NRho = std::stoi(fl_tok);
+
+    // move to next token to get number of temperature points
+    fl_tok = std::strtok(nullptr, " ");
+    fl_tok = std::strtok(nullptr, " ");
+    NT = std::stoi(fl_tok);
+
+    // read 2nd line of header to get min/max density
+    std::getline(ff, fline);
+    // tokenize fline
+    cfline = const_cast<char*>(fline.c_str());
+    fl_tok = std::strtok(cfline, " ");
+    fl_tok = std::strtok(nullptr, " ");
+    const Real RhoMin = std::stod(fl_tok);
+    fl_tok = std::strtok(nullptr, " ");
+    fl_tok = std::strtok(nullptr, " ");
+    const Real RhoMax = std::stod(fl_tok);
+
+    // read 3nd line of header to get min/max temperature
+    std::getline(ff, fline);
+    // tokenize fline
+    cfline = const_cast<char*>(fline.c_str());
+    fl_tok = std::strtok(cfline, " ");
+    fl_tok = std::strtok(nullptr, " ");
+    const Real TMin = std::stod(fl_tok);
+    fl_tok = std::strtok(nullptr, " ");
+    fl_tok = std::strtok(nullptr, " ");
+    const Real TMax = std::stod(fl_tok);
+
+    // reseize the Planck and Rosseland databox
+    lkappa_.resize(NRho, NT, 2);
+
+    // set rho-T ranges
+    const Real lTMin = toLog_(TMin);
+    const Real lTMax = toLog_(TMax);
+    const Real lRhoMin = toLog_(RhoMin);
+    const Real lRhoMax = toLog_(RhoMax);
+    lkappa_.setRange(1, lTMin, lTMax, NT);
+    lkappa_.setRange(2, lRhoMin, lRhoMax, NRho);
+
+    // fill tables
+    for (int iRho = 0; iRho < NRho; ++iRho) {
+      const Real lRho_i = lkappa_.range(2).x(iRho);
+      for (int iT = 0; iT < NT; ++iT) {
+
+        // get new file-line
+        std::getline(ff, fline);
+        cfline = const_cast<char*>(fline.c_str());
+        fl_tok = std::strtok(cfline, " ");
+
+        // populate Rosseland opacity [cm^2/g]
+        lkappa_(iRho, iT, Rosseland) = toLog_(std::stod(fl_tok));
+
+        // populate Planck opacity [cm^2/g]
+        fl_tok = std::strtok(nullptr, " ");
+        lkappa_(iRho, iT, Planck) = toLog_(std::stod(fl_tok));
+
+        if (std::isnan(lkappa_(iRho, iT, Rosseland)) ||
+            std::isnan(lkappa_(iRho, iT, Planck))) {
+          OPAC_ERROR("photons::MeanOpacity: NAN in parsed ASCII opacity");
+        }
+      }
+    }
+  }
+
   PORTABLE_INLINE_FUNCTION Real toLog_(const Real x) const {
     return std::log10(std::abs(x) + EPS);
   }
   PORTABLE_INLINE_FUNCTION Real fromLog_(const Real lx) const {
     return std::pow(10., lx);
   }
-  DataBox lkappaPlanck_;
-  DataBox lkappaRosseland_;
+  Spiner::DataBox<Real> lkappa_;
   const char *filename_;
+  enum OpacityAveraging {
+    Rosseland = 0,
+    Planck = 1
+  };
 };
 
 #undef EPS
 
 } // namespace impl
 
-using MeanOpacityScaleFree = impl::MeanOpacity<PhysicalConstantsUnity>;
-using MeanOpacityCGS = impl::MeanOpacity<PhysicalConstantsCGS>;
-using MeanOpacity = impl::MeanVariant<MeanOpacityScaleFree, MeanOpacityCGS,
-                                      MeanNonCGSUnits<MeanOpacityCGS>>;
+using MeanOpacityBase = impl::MeanOpacity;
+using MeanOpacity =
+    impl::MeanVariant<MeanOpacityBase, MeanNonCGSUnits<MeanOpacityBase>>;
 
 } // namespace photons
 } // namespace singularity
diff --git a/singularity-opac/photons/mean_photon_variant.hpp b/singularity-opac/photons/mean_photon_variant.hpp
index af292ea..a259db4 100644
--- a/singularity-opac/photons/mean_photon_variant.hpp
+++ b/singularity-opac/photons/mean_photon_variant.hpp
@@ -1,5 +1,5 @@
 // ======================================================================
-// © 2022. Triad National Security, LLC. All rights reserved.  This
+// © 2022-2024. Triad National Security, LLC. All rights reserved.  This
 // program was produced under U.S. Government contract
 // 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
 // is operated by Triad National Security, LLC for the U.S.
@@ -69,6 +69,16 @@ class MeanVariant {
         opac_);
   }
 
+  PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
+  GetRuntimePhysicalConstants() const {
+    return mpark::visit(
+        [=](const auto &opac) {
+          using PC = typename std::decay_t<decltype(opac)>::PC;
+          return singularity::GetRuntimePhysicalConstants(PC());
+        },
+        opac_);
+  }
+
   PORTABLE_INLINE_FUNCTION Real
   PlanckMeanAbsorptionCoefficient(const Real rho, const Real temp) const {
     return mpark::visit(
@@ -89,6 +99,12 @@ class MeanVariant {
   inline void Finalize() noexcept {
     return mpark::visit([](auto &opac) { return opac.Finalize(); }, opac_);
   }
+
+#ifdef SPINER_USE_HDF
+  void Save(const std::string &filename) const {
+    return mpark::visit([=](auto &opac) { return opac.Save(filename); }, opac_);
+  }
+#endif
 };
 
 } // namespace impl
diff --git a/singularity-opac/photons/non_cgs_photons.hpp b/singularity-opac/photons/non_cgs_photons.hpp
index a63df14..e128cf9 100644
--- a/singularity-opac/photons/non_cgs_photons.hpp
+++ b/singularity-opac/photons/non_cgs_photons.hpp
@@ -29,6 +29,8 @@ namespace photons {
 template <typename Opac>
 class NonCGSUnits {
  public:
+  using PC = typename Opac::PC;
+
   NonCGSUnits() = default;
   NonCGSUnits(Opac &&opac, const Real time_unit, const Real mass_unit,
               const Real length_unit, const Real temp_unit)
@@ -65,9 +67,10 @@ class NonCGSUnits {
   }
 
   template <typename FrequencyIndexer, typename DataIndexer>
-  PORTABLE_INLINE_FUNCTION void AbsorptionCoefficient(
-      const Real rho, const Real temp, FrequencyIndexer &nu_bins,
-      DataIndexer &coeffs, const int nbins, Real *lambda = nullptr) const {
+  PORTABLE_INLINE_FUNCTION void
+  AbsorptionCoefficient(const Real rho, const Real temp,
+                        FrequencyIndexer &nu_bins, DataIndexer &coeffs,
+                        const int nbins, Real *lambda = nullptr) const {
     for (int i = 0; i < nbins; ++i) {
       nu_bins[i] *= freq_unit_;
     }
@@ -219,6 +222,11 @@ class NonCGSUnits {
     return NoH * mass_unit_ / rho_unit_;
   }
 
+  template <typename T>
+  PORTABLE_INLINE_FUNCTION T GetPhysicalConstants() const {
+    return opac_.GetPhysicalConstants();
+  }
+
  private:
   Opac opac_;
   Real time_unit_, mass_unit_, length_unit_, temp_unit_;
@@ -246,6 +254,12 @@ class MeanNonCGSUnits {
   PORTABLE_INLINE_FUNCTION
   int nlambda() const noexcept { return mean_opac_.nlambda(); }
 
+#ifdef SPINER_USE_HDF
+  void Save(const std::string &filename) const {
+    return mean_opac_.Save(filename);
+  }
+#endif
+
   PORTABLE_INLINE_FUNCTION
   Real PlanckMeanAbsorptionCoefficient(const Real rho, const Real temp) const {
     const Real alpha = mean_opac_.PlanckMeanAbsorptionCoefficient(
diff --git a/singularity-opac/photons/opac_photons.hpp b/singularity-opac/photons/opac_photons.hpp
index 2128a2c..50918a5 100644
--- a/singularity-opac/photons/opac_photons.hpp
+++ b/singularity-opac/photons/opac_photons.hpp
@@ -20,6 +20,7 @@
 #include <singularity-opac/photons/gray_opacity_photons.hpp>
 #include <singularity-opac/photons/non_cgs_photons.hpp>
 #include <singularity-opac/photons/photon_variant.hpp>
+#include <singularity-opac/photons/powerlaw_opacity_photons.hpp>
 #include <singularity-opac/photons/thermal_distributions_photons.hpp>
 
 #include <singularity-opac/photons/mean_opacity_photons.hpp>
@@ -29,10 +30,13 @@ namespace photons {
 
 using ScaleFree = GrayOpacity<PhysicalConstantsUnity>;
 using Gray = GrayOpacity<PhysicalConstantsCGS>;
+using PowerLawScaleFree = PowerLawOpacity<PhysicalConstantsUnity>;
+using PowerLaw = PowerLawOpacity<PhysicalConstantsCGS>;
 using EPBremss = EPBremsstrahlungOpacity<PhysicalConstantsCGS>;
 
-using Opacity = impl::Variant<ScaleFree, Gray, EPBremss, NonCGSUnits<Gray>,
-                              NonCGSUnits<EPBremss>>;
+using Opacity = impl::Variant<ScaleFree, Gray, PowerLawScaleFree, PowerLaw,
+                              EPBremss, NonCGSUnits<Gray>,
+                              NonCGSUnits<PowerLaw>, NonCGSUnits<EPBremss>>;
 
 } // namespace photons
 } // namespace singularity
diff --git a/singularity-opac/photons/photon_variant.hpp b/singularity-opac/photons/photon_variant.hpp
index 440e5e8..3386b8d 100644
--- a/singularity-opac/photons/photon_variant.hpp
+++ b/singularity-opac/photons/photon_variant.hpp
@@ -1,5 +1,5 @@
 // ======================================================================
-// © 2021. Triad National Security, LLC. All rights reserved.  This
+// © 2021-2024. Triad National Security, LLC. All rights reserved.  This
 // program was produced under U.S. Government contract
 // 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
 // is operated by Triad National Security, LLC for the U.S.
@@ -71,6 +71,16 @@ class Variant {
         opac_);
   }
 
+  PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
+  GetRuntimePhysicalConstants() const {
+    return mpark::visit(
+        [=](const auto &opac) {
+          using PC = typename std::decay_t<decltype(opac)>::PC;
+          return singularity::GetRuntimePhysicalConstants(PC());
+        },
+        opac_);
+  }
+
   // Directional absorption coefficient with units of 1/length
   // Signature should be at least
   // rho, temp, nu, lambda
diff --git a/singularity-opac/photons/powerlaw_opacity_photons.hpp b/singularity-opac/photons/powerlaw_opacity_photons.hpp
new file mode 100644
index 0000000..33a4c3f
--- /dev/null
+++ b/singularity-opac/photons/powerlaw_opacity_photons.hpp
@@ -0,0 +1,187 @@
+// ======================================================================
+// © 2024. Triad National Security, LLC. All rights reserved.  This
+// program was produced under U.S. Government contract
+// 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
+// is operated by Triad National Security, LLC for the U.S.
+// Department of Energy/National Nuclear Security Administration. All
+// rights in the program are reserved by Triad National Security, LLC,
+// and the U.S. Department of Energy/National Nuclear Security
+// Administration. The Government is granted for itself and others
+// acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+// license in this material to reproduce, prepare derivative works,
+// distribute copies to the public, perform publicly and display
+// publicly, and to permit others to do so.
+// ======================================================================
+
+#ifndef SINGULARITY_OPAC_PHOTONS_POWERLAW_OPACITY_PHOTONS_
+#define SINGULARITY_OPAC_PHOTONS_POWERLAW_OPACITY_PHOTONS_
+
+#include <cassert>
+#include <cmath>
+#include <cstdio>
+
+#include <ports-of-call/portability.hpp>
+#include <singularity-opac/base/opac_error.hpp>
+#include <singularity-opac/photons/thermal_distributions_photons.hpp>
+
+namespace singularity {
+namespace photons {
+
+template <typename pc = PhysicalConstantsCGS>
+class PowerLawOpacity {
+ public:
+  using PC = pc;
+
+  PowerLawOpacity() = default;
+  PowerLawOpacity(const Real kappa0, const Real rho_exp, const Real temp_exp)
+      : kappa0_(kappa0), rho_exp_(rho_exp), temp_exp_(temp_exp) {}
+  PowerLawOpacity(const PlanckDistribution<pc> &dist, const Real kappa0,
+                  const Real rho_exp, const Real temp_exp)
+      : dist_(dist), kappa0_(kappa0), rho_exp_(rho_exp), temp_exp_(temp_exp) {}
+
+  PowerLawOpacity GetOnDevice() { return *this; }
+  PORTABLE_INLINE_FUNCTION
+  int nlambda() const noexcept { return 0; }
+  PORTABLE_INLINE_FUNCTION
+  void PrintParams() const noexcept {
+    printf("Power law opacity. kappa0 = %g rho_exp = %g temp_exp = %g\n",
+           kappa0_, rho_exp_, temp_exp_);
+  }
+  inline void Finalize() noexcept {}
+
+  PORTABLE_INLINE_FUNCTION
+  Real AbsorptionCoefficient(const Real rho, const Real temp, const Real nu,
+                             Real *lambda = nullptr) const {
+    return dist_.AbsorptionCoefficientFromKirkhoff(*this, rho, temp, nu,
+                                                   lambda);
+  }
+
+  template <typename FrequencyIndexer, typename DataIndexer>
+  PORTABLE_INLINE_FUNCTION void
+  AbsorptionCoefficient(const Real rho, const Real temp,
+                        FrequencyIndexer &nu_bins, DataIndexer &coeffs,
+                        const int nbins, Real *lambda = nullptr) const {
+    for (int i = 0; i < nbins; ++i) {
+      coeffs[i] = AbsorptionCoefficient(rho, temp, nu_bins[i], lambda);
+    }
+  }
+
+  PORTABLE_INLINE_FUNCTION
+  Real AngleAveragedAbsorptionCoefficient(const Real rho, const Real temp,
+                                          const Real nu,
+                                          Real *lambda = nullptr) const {
+    return dist_.AngleAveragedAbsorptionCoefficientFromKirkhoff(
+        *this, rho, temp, nu, lambda);
+  }
+
+  template <typename FrequencyIndexer, typename DataIndexer>
+  PORTABLE_INLINE_FUNCTION void AngleAveragedAbsorptionCoefficient(
+      const Real rho, const Real temp, FrequencyIndexer &nu_bins,
+      DataIndexer &coeffs, const int nbins, Real *lambda = nullptr) const {
+    for (int i = 0; i < nbins; ++i) {
+      coeffs[i] =
+          AngleAveragedAbsorptionCoefficient(rho, temp, nu_bins[i], lambda);
+    }
+  }
+
+  PORTABLE_INLINE_FUNCTION
+  Real EmissivityPerNuOmega(const Real rho, const Real temp, const Real nu,
+                            Real *lambda = nullptr) const {
+    Real Bnu = dist_.ThermalDistributionOfTNu(temp, nu, lambda);
+    return rho *
+           (kappa0_ * std::pow(rho, rho_exp_) * std::pow(temp, temp_exp_)) *
+           Bnu;
+  }
+
+  template <typename FrequencyIndexer, typename DataIndexer>
+  PORTABLE_INLINE_FUNCTION void
+  EmissivityPerNuOmega(const Real rho, const Real temp,
+                       FrequencyIndexer &nu_bins, DataIndexer &coeffs,
+                       const int nbins, Real *lambda = nullptr) const {
+    for (int i = 0; i < nbins; ++i) {
+      coeffs[i] = EmissivityPerNuOmega(rho, temp, nu_bins[i], lambda);
+    }
+  }
+
+  PORTABLE_INLINE_FUNCTION
+  Real EmissivityPerNu(const Real rho, const Real temp, const Real nu,
+                       Real *lambda = nullptr) const {
+    return 4 * M_PI * EmissivityPerNuOmega(rho, temp, nu, lambda);
+  }
+
+  template <typename FrequencyIndexer, typename DataIndexer>
+  PORTABLE_INLINE_FUNCTION void
+  EmissivityPerNu(const Real rho, const Real temp, FrequencyIndexer &nu_bins,
+                  DataIndexer &coeffs, const int nbins,
+                  Real *lambda = nullptr) const {
+    for (int i = 0; i < nbins; ++i) {
+      coeffs[i] = EmissivityPerNu(rho, temp, nu_bins[i], lambda);
+    }
+  }
+
+  PORTABLE_INLINE_FUNCTION
+  Real Emissivity(const Real rho, const Real temp,
+                  Real *lambda = nullptr) const {
+    Real B = dist_.ThermalDistributionOfT(temp, lambda);
+    return rho *
+           (kappa0_ * std::pow(rho, rho_exp_) * std::pow(temp, temp_exp_)) * B;
+  }
+
+  PORTABLE_INLINE_FUNCTION
+  Real NumberEmissivity(const Real rho, const Real temp,
+                        Real *lambda = nullptr) const {
+    return (kappa0_ * std::pow(rho, rho_exp_) * std::pow(temp, temp_exp_)) *
+           dist_.ThermalNumberDistributionOfT(temp, lambda);
+  }
+
+  PORTABLE_INLINE_FUNCTION
+  Real ThermalDistributionOfTNu(const Real temp, const Real nu,
+                                Real *lambda = nullptr) const {
+    return dist_.ThermalDistributionOfTNu(temp, nu, lambda);
+  }
+
+  PORTABLE_INLINE_FUNCTION
+  Real DThermalDistributionOfTNuDT(const Real temp, const Real nu,
+                                   Real *lambda = nullptr) const {
+    return dist_.DThermalDistributionOfTNuDT(temp, nu, lambda);
+  }
+
+  PORTABLE_INLINE_FUNCTION
+  Real ThermalDistributionOfT(const Real temp, Real *lambda = nullptr) const {
+    return dist_.ThermalDistributionOfT(temp, lambda);
+  }
+
+  PORTABLE_INLINE_FUNCTION Real
+  ThermalNumberDistributionOfT(const Real temp, Real *lambda = nullptr) const {
+    return dist_.ThermalNumberDistributionOfT(temp, lambda);
+  }
+
+  PORTABLE_INLINE_FUNCTION
+  Real EnergyDensityFromTemperature(const Real temp,
+                                    Real *lambda = nullptr) const {
+    return dist_.EnergyDensityFromTemperature(temp, lambda);
+  }
+
+  PORTABLE_INLINE_FUNCTION
+  Real TemperatureFromEnergyDensity(const Real er,
+                                    Real *lambda = nullptr) const {
+    return dist_.TemperatureFromEnergyDensity(er, lambda);
+  }
+
+  PORTABLE_INLINE_FUNCTION
+  Real NumberDensityFromTemperature(const Real temp,
+                                    Real *lambda = nullptr) const {
+    return dist_.NumberDensityFromTemperature(temp, lambda);
+  }
+
+ private:
+  Real kappa0_;   // Opacity scale. Units of cm^2/g
+  Real rho_exp_;  // Power law index of density
+  Real temp_exp_; // Power law index of temperature
+  PlanckDistribution<pc> dist_;
+};
+
+} // namespace photons
+} // namespace singularity
+
+#endif // SINGULARITY_OPAC_PHOTONS_POWERLAW_OPACITY_PHOTONS_
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 06f6a86..22d3860 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -39,9 +39,10 @@ PRIVATE
   test_thermal_dist.cpp
   test_scalefree_opacities.cpp
   test_gray_opacities.cpp
-  test_epbremsstrahlung_opacities.cpp
   test_gray_s_opacities.cpp
+  test_epbremsstrahlung_opacities.cpp
   test_brt_opacities.cpp
+  test_powerlaw_opacities.cpp
   test_thomson_s_opacities.cpp
   test_chebyshev.cpp
   test_spiner_opac_neutrinos.cpp
@@ -49,6 +50,11 @@ PRIVATE
   test_variant.cpp
 )
 
+# Link data file directory for Zhu table test
+file(CREATE_LINK
+    ${CMAKE_SOURCE_DIR}/singularity-opac/photons/example_ascii/kap_plaw.txt
+    ${CMAKE_CURRENT_BINARY_DIR}/kap_plaw.txt SYMBOLIC)
+
 target_link_libraries(${PROJECT_NAME}_unit_tests
 PRIVATE
   catch2_define
@@ -58,7 +64,7 @@ PRIVATE
 # Ensure code works with C++11 and earlier
 # TODO(MM): Remove this later when it's not needed.
 set_target_properties(${PROJECT_NAME}_unit_tests
-  PROPERTIES CXX_STANDARD 14
+  PROPERTIES CXX_STANDARD 17
   CXX_STANDARD_REQUIRED YES
   CXX_EXTENSIONS NO)
 
diff --git a/test/test_gray_opacities.cpp b/test/test_gray_opacities.cpp
index 464dbd1..aba0158 100644
--- a/test/test_gray_opacities.cpp
+++ b/test/test_gray_opacities.cpp
@@ -1,5 +1,5 @@
 // ======================================================================
-// © 2021. Triad National Security, LLC. All rights reserved.  This
+// © 2021-2024. Triad National Security, LLC. All rights reserved.  This
 // program was produced under U.S. Government contract
 // 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
 // is operated by Triad National Security, LLC for the U.S.
@@ -54,9 +54,25 @@ TEST_CASE("Gray neutrino opacities", "[GrayNeutrinos]") {
     constexpr RadiationType type = RadiationType::NU_ELECTRON;
     constexpr Real nu = 1.25 * MeV2Hz; // 1 MeV
 
-    neutrinos::Gray opac_host(1);
+    neutrinos::Opacity opac_host = neutrinos::Gray(1.);
     neutrinos::Opacity opac = opac_host.GetOnDevice();
 
+    // Check constants from opacity
+    THEN("Check constants from mean opacity for consistency") {
+      auto constants = opac_host.GetRuntimePhysicalConstants();
+      int n_wrong = 0;
+      if (FractionalDifference(pc::eV, constants.eV) > EPS_TEST) {
+        n_wrong += 1;
+      }
+      if (FractionalDifference(pc::kb, constants.kb) > EPS_TEST) {
+        n_wrong += 1;
+      }
+      if (FractionalDifference(pc::h, constants.h) > EPS_TEST) {
+        n_wrong += 1;
+      }
+      REQUIRE(n_wrong == 0);
+    }
+
     THEN("The emissivity per nu omega is consistent with the emissity per nu") {
       int n_wrong_h = 0;
 #ifdef PORTABILITY_STRATEGY_KOKKOS
@@ -203,6 +219,23 @@ TEST_CASE("Gray photon opacities", "[GrayPhotons]") {
 
     photons::Opacity opac_host = photons::Gray(1);
     photons::Opacity opac = opac_host.GetOnDevice();
+
+    // Check constants from mean opacity
+    THEN("Check constants from mean opacity for consistency") {
+      auto constants = opac_host.GetRuntimePhysicalConstants();
+      int n_wrong = 0;
+      if (FractionalDifference(pc::eV, constants.eV) > EPS_TEST) {
+        n_wrong += 1;
+      }
+      if (FractionalDifference(pc::kb, constants.kb) > EPS_TEST) {
+        n_wrong += 1;
+      }
+      if (FractionalDifference(pc::h, constants.h) > EPS_TEST) {
+        n_wrong += 1;
+      }
+      REQUIRE(n_wrong == 0);
+    }
+
     THEN("The emissivity per nu omega is consistent with the emissity per nu") {
       int n_wrong_h = 0;
 #ifdef PORTABILITY_STRATEGY_KOKKOS
@@ -283,8 +316,7 @@ TEST_CASE("Gray photon opacities", "[GrayPhotons]") {
 
       Real *nu_bins = (Real *)PORTABLE_MALLOC(nbins * sizeof(Real));
       Real *temp_bins = (Real *)PORTABLE_MALLOC(ntemps * sizeof(Real));
-      DataBox loglin_bins(Spiner::AllocationTarget::Device, ntemps,
-                                  nbins);
+      DataBox loglin_bins(Spiner::AllocationTarget::Device, ntemps, nbins);
 
       portableFor(
           "set nu bins", 0, nbins, PORTABLE_LAMBDA(const int &i) {
diff --git a/test/test_mean_opacities.cpp b/test/test_mean_opacities.cpp
index b6ad91f..af2d291 100644
--- a/test/test_mean_opacities.cpp
+++ b/test/test_mean_opacities.cpp
@@ -1,5 +1,5 @@
 // ======================================================================
-// © 2022. Triad National Security, LLC. All rights reserved.  This
+// © 2022-2024. Triad National Security, LLC. All rights reserved.  This
 // program was produced under U.S. Government contract
 // 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
 // is operated by Triad National Security, LLC for the U.S.
@@ -46,7 +46,7 @@ using atomic_view = Kokkos::MemoryTraits<Kokkos::Atomic>;
 
 template <typename T>
 PORTABLE_INLINE_FUNCTION T FractionalDifference(const T &a, const T &b) {
-  return 2 * std::abs(b - a) / (std::abs(a) + std::abs(b) + 1e-20);
+  return 2 * std::abs(b - a) / (std::abs(a) + std::abs(b) + 1e-100);
 }
 constexpr Real EPS_TEST = 1e-3;
 template <typename T>
@@ -87,9 +87,8 @@ TEST_CASE("Mean neutrino opacities", "[MeanNeutrinos]") {
     neutrinos::Gray opac_host(kappa);
     neutrinos::Opacity opac = opac_host.GetOnDevice();
 
-    neutrinos::MeanOpacityCGS mean_opac_host(
+    neutrinos::MeanOpacity mean_opac_host = neutrinos::MeanOpacityBase(
         opac_host, lRhoMin, lRhoMax, NRho, lTMin, lTMax, NT, YeMin, YeMax, NYe);
-    // neutrinos::MeanOpacity mean_opac = mean_opac_host.GetOnDevice();
     auto mean_opac = mean_opac_host.GetOnDevice();
 
     THEN("The emissivity per nu omega is consistent with the emissity per nu") {
@@ -123,7 +122,7 @@ TEST_CASE("Mean neutrino opacities", "[MeanNeutrinos]") {
 #ifdef SPINER_USE_HDF
     THEN("We can save to disk and reload") {
       mean_opac.Save(grayname);
-      neutrinos::MeanOpacityCGS mean_opac_host_load(grayname);
+      neutrinos::MeanOpacity mean_opac_host_load(grayname);
       AND_THEN("The reloaded table matches the gray opacities") {
 
         auto mean_opac_load = mean_opac_host_load.GetOnDevice();
@@ -253,7 +252,6 @@ TEST_CASE("Mean neutrino scattering opacities", "[MeanNeutrinosS]") {
 
     neutrinos::MeanSOpacityCGS mean_opac_host(
         opac_host, lRhoMin, lRhoMax, NRho, lTMin, lTMax, NT, YeMin, YeMax, NYe);
-    // neutrinos::MeanOpacity mean_opac = mean_opac_host.GetOnDevice();
     auto mean_opac = mean_opac_host.GetOnDevice();
 
     THEN("The emissivity per nu omega is consistent with the emissity per nu") {
@@ -414,8 +412,8 @@ TEST_CASE("Mean photon opacities", "[MeanPhotons]") {
     photons::Gray opac_host(kappa);
     photons::Opacity opac = opac_host.GetOnDevice();
 
-    photons::MeanOpacityCGS mean_opac_host(opac_host, lRhoMin, lRhoMax, NRho,
-                                           lTMin, lTMax, NT);
+    photons::MeanOpacity mean_opac_host = photons::MeanOpacityBase(
+        opac_host, lRhoMin, lRhoMax, NRho, lTMin, lTMax, NT);
     auto mean_opac = mean_opac_host.GetOnDevice();
 
     THEN("The emissivity per nu omega is consistent with the emissity per nu") {
@@ -449,7 +447,8 @@ TEST_CASE("Mean photon opacities", "[MeanPhotons]") {
 #ifdef SPINER_USE_HDF
     THEN("We can save to disk and reload") {
       mean_opac.Save(grayname);
-      photons::MeanOpacityCGS mean_opac_host_load(grayname);
+      photons::MeanOpacity mean_opac_host_load =
+          photons::MeanOpacityBase(grayname);
       AND_THEN("The reloaded table matches the gray opacities") {
 
         auto mean_opac_load = mean_opac_host_load.GetOnDevice();
@@ -611,7 +610,8 @@ TEST_CASE("Mean photon scattering opacities", "[MeanPhotonS]") {
         int n_wrong = 0;
         portableReduce(
             "rebuilt table vs gray", 0, NRho, 0, NT, 0, 0,
-            PORTABLE_LAMBDA(const int iRho, const int iT, const int igarbage, int &accumulate) {
+            PORTABLE_LAMBDA(const int iRho, const int iT, const int igarbage,
+                            int &accumulate) {
               const Real lRho =
                   lRhoMin + (lRhoMax - lRhoMin) / (NRho - 1) * iRho;
               const Real rho = std::pow(10, lRho);
@@ -695,3 +695,62 @@ TEST_CASE("Mean photon scattering opacities", "[MeanPhotonS]") {
     opac.Finalize();
   }
 }
+
+TEST_CASE("ASCII-parsed Mean photon opacities", "[MeanPhotons]") {
+  const std::string grayname = "kap_plaw.txt";
+
+  WHEN("We initialize a mean photon opacity from an ASCII table") {
+
+    constexpr Real rho_min = 1e-14;  // g/cc.
+    constexpr Real temp_min = 1.0; // Kelvin.
+    constexpr Real ross_at_min = 0.001; // cm^2/g
+    constexpr Real plnk_at_min = 0.1; // cm^2/g
+    constexpr Real rho_max = 0.7943282347241912;  // g/cc.
+    constexpr Real temp_max = 7943282.347242886; // Kelvin.
+    constexpr Real ross_at_max = 0.001; // cm^2/g
+    constexpr Real plnk_at_max = 0.1; // cm^2/g
+
+    photons::MeanOpacity mean_opac_host = photons::MeanOpacity(grayname);
+    auto mean_opac = mean_opac_host.GetOnDevice();
+
+    THEN("The emissivity per nu omega is consistent with the emissity per nu") {
+      int n_wrong_h = 0;
+#ifdef PORTABILITY_STRATEGY_KOKKOS
+      Kokkos::View<int, atomic_view> n_wrong_d("wrong");
+#else
+      PortableMDArray<int> n_wrong_d(&n_wrong_h, 1);
+#endif
+
+      // get a test value at min rho-T point
+      Real mross = mean_opac.RosselandMeanAbsorptionCoefficient(rho_min, temp_min);
+      Real mplnk = mean_opac.PlanckMeanAbsorptionCoefficient(rho_min, temp_min);
+
+      // check min rho-T point
+      if (FractionalDifference(rho_min * ross_at_min, mross) > EPS_TEST) {
+        n_wrong_d() += 1;
+      }
+      if (FractionalDifference(rho_min * plnk_at_min, mplnk) > EPS_TEST) {
+        n_wrong_d() += 1;
+      }
+
+      // get a test value at max rho-T point
+      mross = mean_opac.RosselandMeanAbsorptionCoefficient(rho_max, temp_max);
+      mplnk = mean_opac.PlanckMeanAbsorptionCoefficient(rho_max, temp_max);
+
+      // check min rho-T point
+      if (FractionalDifference(rho_max * ross_at_max, mross) > EPS_TEST) {
+        n_wrong_d() += 1;
+      }
+      if (FractionalDifference(rho_max * plnk_at_max, mplnk) > EPS_TEST) {
+        n_wrong_d() += 1;
+      }
+
+#ifdef PORTABILITY_STRATEGY_KOKKOS
+      Kokkos::deep_copy(n_wrong_h, n_wrong_d);
+#endif
+      REQUIRE(n_wrong_h == 0);
+    }
+
+    mean_opac.Finalize();
+  }
+}
diff --git a/test/test_powerlaw_opacities.cpp b/test/test_powerlaw_opacities.cpp
new file mode 100644
index 0000000..73e39b1
--- /dev/null
+++ b/test/test_powerlaw_opacities.cpp
@@ -0,0 +1,198 @@
+// ======================================================================
+// © 2024. Triad National Security, LLC. All rights reserved.  This
+// program was produced under U.S. Government contract
+// 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
+// is operated by Triad National Security, LLC for the U.S.
+// Department of Energy/National Nuclear Security Administration. All
+// rights in the program are reserved by Triad National Security, LLC,
+// and the U.S. Department of Energy/National Nuclear Security
+// Administration. The Government is granted for itself and others
+// acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+// license in this material to reproduce, prepare derivative works,
+// distribute copies to the public, perform publicly and display
+// publicly, and to permit others to do so.
+// ======================================================================
+
+#include <cmath>
+#include <iostream>
+
+#include <catch2/catch.hpp>
+
+#include <ports-of-call/portability.hpp>
+#include <ports-of-call/portable_arrays.hpp>
+#include <spiner/databox.hpp>
+
+#include <singularity-opac/photons/opac_photons.hpp>
+
+using namespace singularity;
+
+using DataBox = Spiner::DataBox<Real>;
+
+#ifdef PORTABILITY_STRATEGY_KOKKOS
+using atomic_view = Kokkos::MemoryTraits<Kokkos::Atomic>;
+#endif
+
+template <typename T>
+PORTABLE_INLINE_FUNCTION T FractionalDifference(const T &a, const T &b) {
+  return 2 * std::abs(b - a) / (std::abs(a) + std::abs(b) + 1e-20);
+}
+PORTABLE_INLINE_FUNCTION Real CalcFrequency(const int n, const Real nu_min,
+                                            const Real nu_max, const int n_nu) {
+  const Real lnu_min = std::log(nu_min);
+  const Real lnu_max = std::log(nu_max);
+  const Real dlnu = (lnu_max - lnu_min) / static_cast<Real>(n_nu);
+  return std::exp(lnu_min + (n + 0.5) * dlnu);
+}
+
+constexpr Real EPS_TEST = 1e-3;
+
+constexpr Real rho_exp = 2.;
+constexpr Real temp_exp = 2.5;
+
+TEST_CASE("Scale free power law photon opacities",
+          "[PowerLawScaleFreePhotonOpacities]") {
+  WHEN("We initialize a scale-free power law photon opacity") {
+    constexpr Real rho = 1.5e0;
+    constexpr Real temp = 3.2e0;
+    constexpr Real nu_min = 1.e-1;
+    constexpr Real nu_max = 1.e1;
+    constexpr int n_nu = 100;
+    constexpr Real kappa0 = 1.5;
+
+    photons::PowerLawScaleFree opac_host(kappa0, rho_exp, temp_exp);
+    photons::Opacity opac = opac_host.GetOnDevice();
+
+    THEN("The emissivity per nu omega is consistent with the emissity per nu") {
+      int n_wrong_h = 0;
+#ifdef PORTABILITY_STRATEGY_KOKKOS
+      Kokkos::View<int, atomic_view> n_wrong_d("wrong");
+#else
+      PortableMDArray<int> n_wrong_d(&n_wrong_h, 1);
+#endif
+
+      portableFor(
+          "calc emissivities", 0, n_nu, PORTABLE_LAMBDA(const int &i) {
+            const Real nu = CalcFrequency(i, nu_min, nu_max, n_nu);
+            const Real jnu = opac.EmissivityPerNuOmega(rho, temp, nu);
+            const Real Jnu = opac.EmissivityPerNu(rho, temp, nu);
+            const Real kappa =
+                kappa0 * std::pow(rho, rho_exp) * std::pow(temp, temp_exp);
+            if (FractionalDifference(jnu, rho * kappa *
+                                              opac.ThermalDistributionOfTNu(
+                                                  temp, nu)) > EPS_TEST) {
+              n_wrong_d() += 1;
+            }
+            if (FractionalDifference(Jnu, 4 * M_PI * jnu) > EPS_TEST) {
+              n_wrong_d() += 1;
+            }
+          });
+
+#ifdef PORTABILITY_STRATEGY_KOKKOS
+      Kokkos::deep_copy(n_wrong_h, n_wrong_d);
+#endif
+      REQUIRE(n_wrong_h == 0);
+    }
+
+    opac.Finalize();
+  }
+}
+
+TEST_CASE("CGS power law photon opacities", "[PowerLawCGSPhotonOpacities]") {
+  constexpr Real rho = 1.5e0;    // g/cc
+  constexpr Real temp = 1.e3;    // K
+  constexpr Real nu_min = 1.e10; // Hz
+  constexpr Real nu_max = 1.e14; // Hz
+  constexpr int n_nu = 100;
+  constexpr Real kappa0 = 1.5; // cm^2 / g
+
+  WHEN("We initialize a CGS power law photon opacity") {
+    photons::PowerLaw opac_host(kappa0, rho_exp, temp_exp);
+    photons::Opacity opac = opac_host.GetOnDevice();
+
+    THEN("The emissivity per nu omega is consistent with the emissity per nu") {
+      int n_wrong_h = 0;
+#ifdef PORTABILITY_STRATEGY_KOKKOS
+      Kokkos::View<int, atomic_view> n_wrong_d("wrong");
+#else
+      PortableMDArray<int> n_wrong_d(&n_wrong_h, 1);
+#endif
+
+      portableFor(
+          "calc emissivities", 0, n_nu, PORTABLE_LAMBDA(const int &i) {
+            const Real nu = CalcFrequency(i, nu_min, nu_max, n_nu);
+            const Real jnu = opac.EmissivityPerNuOmega(rho, temp, nu);
+            const Real Jnu = opac.EmissivityPerNu(rho, temp, nu);
+            const Real kappa =
+                kappa0 * std::pow(rho, rho_exp) * std::pow(temp, temp_exp);
+            if (FractionalDifference(jnu, rho * kappa *
+                                              opac.ThermalDistributionOfTNu(
+                                                  temp, nu)) > EPS_TEST) {
+              n_wrong_d() += 1;
+            }
+            if (FractionalDifference(Jnu, 4 * M_PI * jnu) > EPS_TEST) {
+              n_wrong_d() += 1;
+            }
+          });
+
+#ifdef PORTABILITY_STRATEGY_KOKKOS
+      Kokkos::deep_copy(n_wrong_h, n_wrong_d);
+#endif
+      REQUIRE(n_wrong_h == 0);
+    }
+
+    opac.Finalize();
+  }
+
+  WHEN("We initialize a CGS power law photon opacity with non-CGS units") {
+    constexpr Real time_unit = 123.;
+    constexpr Real mass_unit = 456.;
+    constexpr Real length_unit = 789.;
+    constexpr Real temp_unit = 276.;
+    constexpr Real rho_unit =
+        mass_unit / (length_unit * length_unit * length_unit);
+    constexpr Real j_unit = mass_unit / (length_unit * time_unit * time_unit);
+
+    photons::NonCGSUnits<photons::PowerLaw> opac_host(
+        photons::PowerLaw(kappa0, rho_exp, temp_exp), time_unit, mass_unit,
+        length_unit, temp_unit);
+    photons::Opacity opac = opac_host.GetOnDevice();
+    photons::PowerLaw opac_cgs_host(kappa0, rho_exp, temp_exp);
+    photons::Opacity opac_cgs = opac_cgs_host.GetOnDevice();
+
+    THEN("The emissivity per nu omega is consistent with the emissity per nu") {
+      int n_wrong_h = 0;
+#ifdef PORTABILITY_STRATEGY_KOKKOS
+      Kokkos::View<int, atomic_view> n_wrong_d("wrong");
+#else
+      PortableMDArray<int> n_wrong_d(&n_wrong_h, 1);
+#endif
+
+      portableFor(
+          "calc emissivities", 0, n_nu, PORTABLE_LAMBDA(const int &i) {
+            const Real nu = CalcFrequency(i, nu_min, nu_max, n_nu);
+            const Real jnu = opac.EmissivityPerNuOmega(
+                rho / rho_unit, temp / temp_unit, nu * time_unit);
+            const Real Jnu = opac.EmissivityPerNu(
+                rho / rho_unit, temp / temp_unit, nu * time_unit);
+            const Real kappa =
+                kappa0 * std::pow(rho, rho_exp) * std::pow(temp, temp_exp);
+            if (FractionalDifference(
+                    jnu * j_unit,
+                    opac_cgs.EmissivityPerNuOmega(rho, temp, nu)) > EPS_TEST) {
+              n_wrong_d() += 1;
+            }
+            if (FractionalDifference(Jnu, 4 * M_PI * jnu) > EPS_TEST) {
+              n_wrong_d() += 1;
+            }
+          });
+
+#ifdef PORTABILITY_STRATEGY_KOKKOS
+      Kokkos::deep_copy(n_wrong_h, n_wrong_d);
+#endif
+      REQUIRE(n_wrong_h == 0);
+    }
+
+    opac.Finalize();
+  }
+}
+
diff --git a/utils/kokkos b/utils/kokkos
index 1fb0c28..08ceff9 160000
--- a/utils/kokkos
+++ b/utils/kokkos
@@ -1 +1 @@
-Subproject commit 1fb0c284d458c75370094921d9f202c287502325
+Subproject commit 08ceff92bcf3a828844480bc1e6137eb74028517
diff --git a/utils/ports-of-call b/utils/ports-of-call
index 993d820..58ce118 160000
--- a/utils/ports-of-call
+++ b/utils/ports-of-call
@@ -1 +1 @@
-Subproject commit 993d8209de98da186c5e0b445a4d6c359afa8c81
+Subproject commit 58ce1181b2d835bd32673ad70550c9130381f91b
diff --git a/utils/spiner b/utils/spiner
index 33862f8..bd57161 160000
--- a/utils/spiner
+++ b/utils/spiner
@@ -1 +1 @@
-Subproject commit 33862f831be65e2dd2284b42f05c477b6fe7367a
+Subproject commit bd57161576a62a13341dada183f2b336a3e99b08
diff --git a/utils/variant b/utils/variant
index 3c7fc82..23cb94f 160000
--- a/utils/variant
+++ b/utils/variant
@@ -1 +1 @@
-Subproject commit 3c7fc8266bb46046b42c2dc2663f9f505f0cec28
+Subproject commit 23cb94f027d4ef33bf48133acc2695c7e5c6f1e7