diff --git a/singularity-opac/constants/constants.hpp b/singularity-opac/constants/constants.hpp
index d08bb4a..f5bbede 100644
--- a/singularity-opac/constants/constants.hpp
+++ b/singularity-opac/constants/constants.hpp
@@ -137,10 +137,6 @@ struct PhysicalConstants {
                                                  (mass * time * time);
   static constexpr Real g_newt = gravitational_constant;
 
-  static constexpr Real acceleration_from_gravity =
-      BASE::acceleration_from_gravity * length / (time * time);
-  static constexpr Real g_accel = acceleration_from_gravity;
-
   static constexpr Real electron_mass = BASE::electron_mass * mass;
   static constexpr Real me = electron_mass;
 
@@ -163,19 +159,6 @@ struct PhysicalConstants {
   static constexpr Real radiation_constant = 4.0 * sb / c;
   static constexpr Real ar = radiation_constant;
 
-  // Faraday constant (derived)
-  static constexpr Real faraday_constant = na * qe;
-  static constexpr Real faraday = faraday_constant;
-
-  // Permeability of free space
-  static constexpr Real vacuum_permeability =
-      BASE::vacuum_permeability * force / (current * current);
-  static constexpr Real mu0 = vacuum_permeability;
-
-  // Permittivity of free space (derived)
-  static constexpr Real vacuum_permittivity = 1.0 / (mu0 * c * c);
-  static constexpr Real eps0 = vacuum_permittivity;
-
   // Electron volt (derived)
   static constexpr Real electron_volt = BASE::elementary_charge * energy;
   static constexpr Real eV = electron_volt;
@@ -195,14 +178,38 @@ struct PhysicalConstants {
 };
 
 struct RuntimePhysicalConstants {
+  template <typename T>
+  PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants(T pc)
+      : length(1.), time(1.), mass(1.), temp(1.), 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), me(pc.me), mp(pc.mp), mn(pc.mn), amu(pc.amu),
+        sb(pc.sb), ar(pc.ar), eV(pc.eV), Fc(pc.Fc), nu_sigma0(pc.nu_sigma0),
+        gA(pc.gA) {}
+
   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) {}
+  RuntimePhysicalConstants(T pc, const Real time_, const Real mass_,
+                           const Real length_, const Real temp_)
+      : time(time_), mass(mass_), length(length_), temp(temp_), na(pc.na),
+        alpha(pc.alpha), h(pc.h * time / mass * std::pow(length, -2)),
+        hbar(pc.hbar * time / mass * std::pow(length, -2)),
+        kb(pc.kb * std::pow(time, 2) / mass * std::pow(length, -2) * temp),
+        r_gas(pc.r_gas * std::pow(time, 2) / mass * std::pow(length, -2)),
+        qe(pc.qe), c(pc.c * time / length),
+        g_newt(pc.g_newt * std::pow(length, -3) / mass * std::pow(time, 2)),
+        me(pc.me / mass), mp(pc.mp / mass), mn(pc.mn / mass),
+        amu(pc.amu / mass),
+        sb(pc.sb * std::pow(time, 3) / mass * std::pow(temp, 4)),
+        ar(pc.ar * std::pow(time, 2) * length / mass * std::pow(temp, 4)),
+        eV(pc.eV * std::pow(time, 2) / mass * std::pow(length, -2)),
+        Fc(pc.Fc * std::pow(time, 2) * length / mass),
+        nu_sigma0(pc.nu_sigma0 / std::pow(length, 2)), gA(pc.gA) {}
+
+  // Conversion factors to the T pc unit system from code units
+  const Real time;
+  const Real mass;
+  const Real length;
+  const Real temp;
 
   const Real na;
   const Real alpha;
@@ -213,28 +220,18 @@ struct RuntimePhysicalConstants {
   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 9138084..6ac9cf1 100644
--- a/singularity-opac/neutrinos/brt_neutrinos.hpp
+++ b/singularity-opac/neutrinos/brt_neutrinos.hpp
@@ -173,6 +173,11 @@ class BRTOpacity {
     return dist_.NumberDensityFromTemperature(temp, type, lambda);
   }
 
+  PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
+  GetRuntimePhysicalConstants() const {
+    return RuntimePhysicalConstants(PC());
+  }
+
  private:
   PORTABLE_INLINE_FUNCTION
   Real GetSigmac(const RadiationType type, const Real nu) const {
diff --git a/singularity-opac/neutrinos/gray_opacity_neutrinos.hpp b/singularity-opac/neutrinos/gray_opacity_neutrinos.hpp
index 830615e..8dae38d 100644
--- a/singularity-opac/neutrinos/gray_opacity_neutrinos.hpp
+++ b/singularity-opac/neutrinos/gray_opacity_neutrinos.hpp
@@ -175,6 +175,11 @@ class GrayOpacity {
     return dist_.NumberDensityFromTemperature(temp, type, lambda);
   }
 
+  PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
+  GetRuntimePhysicalConstants() const {
+    return RuntimePhysicalConstants(PC());
+  }
+
  private:
   Real kappa_; // Opacity. Units of cm^2/g
   ThermalDistribution dist_;
diff --git a/singularity-opac/neutrinos/mean_neutrino_variant.hpp b/singularity-opac/neutrinos/mean_neutrino_variant.hpp
index 488085b..0675fed 100644
--- a/singularity-opac/neutrinos/mean_neutrino_variant.hpp
+++ b/singularity-opac/neutrinos/mean_neutrino_variant.hpp
@@ -72,11 +72,7 @@ class MeanVariant {
   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_);
+        [](auto &opac) { return opac.GetRuntimePhysicalConstants(); }, opac_);
   }
 
   PORTABLE_INLINE_FUNCTION Real PlanckMeanAbsorptionCoefficient(
diff --git a/singularity-opac/neutrinos/mean_opacity_neutrinos.hpp b/singularity-opac/neutrinos/mean_opacity_neutrinos.hpp
index 52e9485..9cd578e 100644
--- a/singularity-opac/neutrinos/mean_opacity_neutrinos.hpp
+++ b/singularity-opac/neutrinos/mean_opacity_neutrinos.hpp
@@ -37,9 +37,11 @@ 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:
+  using PC = pc;
+
   MeanOpacity() = default;
   template <typename Opacity>
   MeanOpacity(const Opacity &opac, const Real lRhoMin, const Real lRhoMax,
@@ -106,6 +108,11 @@ class MeanOpacity {
     lkappaRosseland_.finalize();
   }
 
+  PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
+  GetRuntimePhysicalConstants() const {
+    return RuntimePhysicalConstants(PC());
+  }
+
   PORTABLE_INLINE_FUNCTION
   Real PlanckMeanAbsorptionCoefficient(const Real rho, const Real temp,
                                        const Real Ye,
@@ -133,7 +140,9 @@ 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;
+    static_assert(std::is_same<typename Opacity::PC, PC>::value,
+                  "Error: MeanOpacity physical constants do not match those of "
+                  "Opacity used for construction.");
 
     lkappaPlanck_.resize(NRho, NT, NYe, NEUTRINO_NTYPES);
     // index 0 is the species and is not interpolatable
@@ -220,7 +229,7 @@ class MeanOpacity {
 
 } // namespace impl
 
-using MeanOpacityBase = impl::MeanOpacity;
+using MeanOpacityBase = impl::MeanOpacity<PhysicalConstantsCGS>;
 using MeanOpacity =
     impl::MeanVariant<MeanOpacityBase, MeanNonCGSUnits<MeanOpacityBase>>;
 
diff --git a/singularity-opac/neutrinos/neutrino_variant.hpp b/singularity-opac/neutrinos/neutrino_variant.hpp
index 5736077..b53551a 100644
--- a/singularity-opac/neutrinos/neutrino_variant.hpp
+++ b/singularity-opac/neutrinos/neutrino_variant.hpp
@@ -74,11 +74,7 @@ class Variant {
   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_);
+        [](auto &opac) { return opac.GetRuntimePhysicalConstants(); }, opac_);
   }
 
   // Directional absorption coefficient with units of 1/length
diff --git a/singularity-opac/neutrinos/non_cgs_neutrinos.hpp b/singularity-opac/neutrinos/non_cgs_neutrinos.hpp
index 400cb2b..16a7b7f 100644
--- a/singularity-opac/neutrinos/non_cgs_neutrinos.hpp
+++ b/singularity-opac/neutrinos/non_cgs_neutrinos.hpp
@@ -239,6 +239,12 @@ class NonCGSUnits {
     return NoH * mass_unit_ / rho_unit_;
   }
 
+  PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
+  GetRuntimePhysicalConstants() const {
+    return RuntimePhysicalConstants(PC(), time_unit_, mass_unit_, length_unit_,
+                                    temp_unit_);
+  }
+
  private:
   Opac opac_;
   Real time_unit_, mass_unit_, length_unit_, temp_unit_;
@@ -249,6 +255,8 @@ class NonCGSUnits {
 template <typename MeanOpac>
 class MeanNonCGSUnits {
  public:
+  using PC = typename MeanOpac::PC;
+
   MeanNonCGSUnits() = default;
   MeanNonCGSUnits(MeanOpac &&mean_opac, const Real time_unit,
                   const Real mass_unit, const Real length_unit,
@@ -298,6 +306,12 @@ class MeanNonCGSUnits {
     return alpha * length_unit_;
   }
 
+  PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
+  GetRuntimePhysicalConstants() const {
+    return RuntimePhysicalConstants(PhysicalConstantsCGS(), time_unit_,
+                                    mass_unit_, length_unit_, temp_unit_);
+  }
+
  private:
   MeanOpac mean_opac_;
   Real time_unit_, mass_unit_, length_unit_, temp_unit_;
diff --git a/singularity-opac/neutrinos/spiner_opac_neutrinos.hpp b/singularity-opac/neutrinos/spiner_opac_neutrinos.hpp
index 906bb20..99563c1 100644
--- a/singularity-opac/neutrinos/spiner_opac_neutrinos.hpp
+++ b/singularity-opac/neutrinos/spiner_opac_neutrinos.hpp
@@ -369,6 +369,11 @@ class SpinerOpacity {
     return dist_.NumberDensityFromTemperature(temp, type, lambda);
   }
 
+  PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
+  GetRuntimePhysicalConstants() const {
+    return RuntimePhysicalConstants(PC());
+  }
+
  private:
   // TODO(JMM): Offsets probably not necessary
   PORTABLE_INLINE_FUNCTION Real toLog_(const Real x, const Real offset) const {
diff --git a/singularity-opac/neutrinos/tophat_emissivity_neutrinos.hpp b/singularity-opac/neutrinos/tophat_emissivity_neutrinos.hpp
index bfb024b..a63e325 100644
--- a/singularity-opac/neutrinos/tophat_emissivity_neutrinos.hpp
+++ b/singularity-opac/neutrinos/tophat_emissivity_neutrinos.hpp
@@ -191,6 +191,11 @@ class TophatEmissivity {
     return dist_.NumberDensityFromTemperature(temp, type, lambda);
   }
 
+  PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
+  GetRuntimePhysicalConstants() const {
+    return RuntimePhysicalConstants(PC());
+  }
+
  private:
   PORTABLE_INLINE_FUNCTION
   Real GetYeF(RadiationType type, Real Ye) const {
diff --git a/singularity-opac/photons/epbremsstrahlung_opacity_photons.hpp b/singularity-opac/photons/epbremsstrahlung_opacity_photons.hpp
index 725681d..1c85da5 100644
--- a/singularity-opac/photons/epbremsstrahlung_opacity_photons.hpp
+++ b/singularity-opac/photons/epbremsstrahlung_opacity_photons.hpp
@@ -179,6 +179,11 @@ class EPBremsstrahlungOpacity {
     return dist_.NumberDensityFromTemperature(temp, lambda);
   }
 
+  PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
+  GetRuntimePhysicalConstants() const {
+    return RuntimePhysicalConstants(PC());
+  }
+
  private:
   Real mass_ion_;
   PlanckDistribution<pc> dist_;
diff --git a/singularity-opac/photons/gray_opacity_photons.hpp b/singularity-opac/photons/gray_opacity_photons.hpp
index dcca68c..a22a0f7 100644
--- a/singularity-opac/photons/gray_opacity_photons.hpp
+++ b/singularity-opac/photons/gray_opacity_photons.hpp
@@ -167,6 +167,11 @@ class GrayOpacity {
     return dist_.NumberDensityFromTemperature(temp, lambda);
   }
 
+  PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
+  GetRuntimePhysicalConstants() const {
+    return RuntimePhysicalConstants(PC());
+  }
+
  private:
   Real kappa_; // Opacity. Units of cm^2/g
   PlanckDistribution<pc> dist_;
diff --git a/singularity-opac/photons/mean_opacity_photons.hpp b/singularity-opac/photons/mean_opacity_photons.hpp
index d4ac76b..ee85645 100644
--- a/singularity-opac/photons/mean_opacity_photons.hpp
+++ b/singularity-opac/photons/mean_opacity_photons.hpp
@@ -37,9 +37,11 @@ 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 PC = pc;
+
   MeanOpacity() = default;
   template <typename Opacity>
   MeanOpacity(const Opacity &opac, const Real lRhoMin, const Real lRhoMax,
@@ -102,6 +104,11 @@ class MeanOpacity {
     lkappaRosseland_.finalize();
   }
 
+  PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
+  GetRuntimePhysicalConstants() const {
+    return RuntimePhysicalConstants(PC());
+  }
+
   PORTABLE_INLINE_FUNCTION
   Real PlanckMeanAbsorptionCoefficient(const Real rho, const Real temp) const {
     Real lRho = toLog_(rho);
@@ -123,7 +130,9 @@ 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) {
-    using PC = typename Opacity::PC;
+    static_assert(std::is_same<typename Opacity::PC, PC>::value,
+                  "Error: MeanOpacity physical constants do not match those of "
+                  "Opacity used for construction.");
 
     lkappaPlanck_.resize(NRho, NT);
     lkappaPlanck_.setRange(0, lTMin, lTMax, NT);
@@ -196,7 +205,7 @@ class MeanOpacity {
 
 } // namespace impl
 
-using MeanOpacityBase = impl::MeanOpacity;
+using MeanOpacityBase = impl::MeanOpacity<PhysicalConstantsCGS>;
 using MeanOpacity =
     impl::MeanVariant<MeanOpacityBase, MeanNonCGSUnits<MeanOpacityBase>>;
 
diff --git a/singularity-opac/photons/mean_photon_variant.hpp b/singularity-opac/photons/mean_photon_variant.hpp
index a259db4..ff860b0 100644
--- a/singularity-opac/photons/mean_photon_variant.hpp
+++ b/singularity-opac/photons/mean_photon_variant.hpp
@@ -72,11 +72,7 @@ class MeanVariant {
   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_);
+        [](auto &opac) { return opac.GetRuntimePhysicalConstants(); }, opac_);
   }
 
   PORTABLE_INLINE_FUNCTION Real
diff --git a/singularity-opac/photons/non_cgs_photons.hpp b/singularity-opac/photons/non_cgs_photons.hpp
index e128cf9..36857b0 100644
--- a/singularity-opac/photons/non_cgs_photons.hpp
+++ b/singularity-opac/photons/non_cgs_photons.hpp
@@ -222,9 +222,10 @@ class NonCGSUnits {
     return NoH * mass_unit_ / rho_unit_;
   }
 
-  template <typename T>
-  PORTABLE_INLINE_FUNCTION T GetPhysicalConstants() const {
-    return opac_.GetPhysicalConstants();
+  PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
+  GetRuntimePhysicalConstants() const {
+    return RuntimePhysicalConstants(PC(), time_unit_, mass_unit_, length_unit_,
+                                    temp_unit_);
   }
 
  private:
@@ -237,6 +238,8 @@ class NonCGSUnits {
 template <typename MeanOpac>
 class MeanNonCGSUnits {
  public:
+  using PC = typename MeanOpac::PC;
+
   MeanNonCGSUnits() = default;
   MeanNonCGSUnits(MeanOpac &&mean_opac, const Real time_unit,
                   const Real mass_unit, const Real length_unit,
@@ -283,6 +286,12 @@ class MeanNonCGSUnits {
     return alpha * length_unit_;
   }
 
+  PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
+  GetRuntimePhysicalConstants() const {
+    return RuntimePhysicalConstants(PhysicalConstantsCGS(), time_unit_,
+                                    mass_unit_, length_unit_, temp_unit_);
+  }
+
  private:
   MeanOpac mean_opac_;
   Real time_unit_, mass_unit_, length_unit_, temp_unit_;
diff --git a/singularity-opac/photons/photon_variant.hpp b/singularity-opac/photons/photon_variant.hpp
index 3386b8d..0b2a1e7 100644
--- a/singularity-opac/photons/photon_variant.hpp
+++ b/singularity-opac/photons/photon_variant.hpp
@@ -74,11 +74,7 @@ class Variant {
   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_);
+        [](auto &opac) { return opac.GetRuntimePhysicalConstants(); }, opac_);
   }
 
   // Directional absorption coefficient with units of 1/length
diff --git a/singularity-opac/photons/powerlaw_opacity_photons.hpp b/singularity-opac/photons/powerlaw_opacity_photons.hpp
index 33a4c3f..c971eea 100644
--- a/singularity-opac/photons/powerlaw_opacity_photons.hpp
+++ b/singularity-opac/photons/powerlaw_opacity_photons.hpp
@@ -174,6 +174,11 @@ class PowerLawOpacity {
     return dist_.NumberDensityFromTemperature(temp, lambda);
   }
 
+  PORTABLE_INLINE_FUNCTION RuntimePhysicalConstants
+  GetRuntimePhysicalConstants() const {
+    return RuntimePhysicalConstants(PC());
+  }
+
  private:
   Real kappa0_;   // Opacity scale. Units of cm^2/g
   Real rho_exp_;  // Power law index of density
diff --git a/test/test_gray_opacities.cpp b/test/test_gray_opacities.cpp
index aba0158..7c898a5 100644
--- a/test/test_gray_opacities.cpp
+++ b/test/test_gray_opacities.cpp
@@ -109,6 +109,38 @@ TEST_CASE("Gray neutrino opacities", "[GrayNeutrinos]") {
               neutrinos::Gray(1), time_unit, mass_unit, length_unit, temp_unit);
       auto funny_units = funny_units_host.GetOnDevice();
 
+      THEN("We can retrieve physical constants in code units") {
+        auto noncgs_rpc = funny_units.GetRuntimePhysicalConstants();
+        REQUIRE(FractionalDifference(noncgs_rpc.length, length_unit) <
+                EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.time, time_unit) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.mass, mass_unit) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.temp, temp_unit) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.na, 6.022141e+23) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.alpha, 7.297353e-03) <
+                EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.h, 2.871060e-33) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.hbar, 4.569434e-34) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.kb, 2.030877e-18) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.r_gas, 4.431243e+03) <
+                EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.qe, 4.803205e-10) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.c, 4.673571e+09) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.g_newt, 4.508065e-15) <
+                EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.me, 1.997672e-30) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.mp, 3.668030e-27) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.mn, 3.673086e-27) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.amu, 3.641533e-27) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.sb, 1.342760e+09) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.ar, 1.149237e+00) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.eV, 8.538896e-17) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.Fc, 1.189435e-09) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.nu_sigma0, 2.829094e-74) <
+                EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.gA, -1.23) < EPS_TEST);
+      }
+
       THEN("We can convert meaningfully into and out of funny units") {
         int n_wrong_h = 0;
 #ifdef PORTABILITY_STRATEGY_KOKKOS
@@ -269,6 +301,38 @@ TEST_CASE("Gray photon opacities", "[GrayPhotons]") {
           photons::Gray(1), time_unit, mass_unit, length_unit, temp_unit);
       auto funny_units = funny_units_host.GetOnDevice();
 
+      THEN("We can retrieve physical constants in code units") {
+        auto noncgs_rpc = funny_units.GetRuntimePhysicalConstants();
+        REQUIRE(FractionalDifference(noncgs_rpc.length, length_unit) <
+                EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.time, time_unit) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.mass, mass_unit) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.temp, temp_unit) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.na, 6.022141e+23) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.alpha, 7.297353e-03) <
+                EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.h, 2.871060e-33) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.hbar, 4.569434e-34) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.kb, 2.030877e-18) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.r_gas, 4.431243e+03) <
+                EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.qe, 4.803205e-10) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.c, 4.673571e+09) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.g_newt, 4.508065e-15) <
+                EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.me, 1.997672e-30) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.mp, 3.668030e-27) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.mn, 3.673086e-27) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.amu, 3.641533e-27) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.sb, 1.342760e+09) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.ar, 1.149237e+00) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.eV, 8.538896e-17) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.Fc, 1.189435e-09) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.nu_sigma0, 2.829094e-74) <
+                EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.gA, -1.23) < EPS_TEST);
+      }
+
       THEN("We can convert meaningfully into and out of funny units") {
         int n_wrong_h = 0;
 #ifdef PORTABILITY_STRATEGY_KOKKOS
diff --git a/test/test_mean_opacities.cpp b/test/test_mean_opacities.cpp
index 932d3d1..0d384f0 100644
--- a/test/test_mean_opacities.cpp
+++ b/test/test_mean_opacities.cpp
@@ -174,13 +174,48 @@ TEST_CASE("Mean neutrino opacities", "[MeanNeutrinos]") {
       constexpr Real rho_unit =
           mass_unit / (length_unit * length_unit * length_unit);
 
+      auto mean_opac_host_base =
+          neutrinos::MeanOpacityBase(opac_host, lRhoMin, lRhoMax, NRho, lTMin,
+                                     lTMax, NT, YeMin, YeMax, NYe);
       auto funny_units_host =
-          neutrinos::MeanNonCGSUnits<neutrinos::MeanOpacity>(
-              std::forward<neutrinos::MeanOpacity>(mean_opac_host), time_unit,
-              mass_unit, length_unit, temp_unit);
+          neutrinos::MeanNonCGSUnits<neutrinos::MeanOpacityBase>(
+              std::forward<neutrinos::MeanOpacityBase>(mean_opac_host_base),
+              time_unit, mass_unit, length_unit, temp_unit);
 
       auto funny_units = funny_units_host.GetOnDevice();
 
+      THEN("We can retrieve physical constants in code units") {
+        auto noncgs_rpc = funny_units.GetRuntimePhysicalConstants();
+        REQUIRE(FractionalDifference(noncgs_rpc.length, length_unit) <
+                EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.time, time_unit) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.mass, mass_unit) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.temp, temp_unit) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.na, 6.022141e+23) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.alpha, 7.297353e-03) <
+                EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.h, 2.871060e-33) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.hbar, 4.569434e-34) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.kb, 2.030877e-18) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.r_gas, 4.431243e+03) <
+                EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.qe, 4.803205e-10) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.c, 4.673571e+09) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.g_newt, 4.508065e-15) <
+                EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.me, 1.997672e-30) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.mp, 3.668030e-27) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.mn, 3.673086e-27) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.amu, 3.641533e-27) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.sb, 1.342760e+09) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.ar, 1.149237e+00) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.eV, 8.538896e-17) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.Fc, 1.189435e-09) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.nu_sigma0, 2.829094e-74) <
+                EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.gA, -1.23) < EPS_TEST);
+      }
+
       int n_wrong_h = 0;
 #ifdef PORTABILITY_STRATEGY_KOKKOS
       Kokkos::View<int, atomic_view> n_wrong_d("wrong");
@@ -495,12 +530,47 @@ TEST_CASE("Mean photon opacities", "[MeanPhotons]") {
       constexpr Real rho_unit =
           mass_unit / (length_unit * length_unit * length_unit);
 
-      auto funny_units_host = photons::MeanNonCGSUnits<photons::MeanOpacity>(
-          std::forward<photons::MeanOpacity>(mean_opac_host), time_unit,
-          mass_unit, length_unit, temp_unit);
+      auto mean_opac_host_base = photons::MeanOpacityBase(
+          opac_host, lRhoMin, lRhoMax, NRho, lTMin, lTMax, NT);
+      auto funny_units_host =
+          photons::MeanNonCGSUnits<photons::MeanOpacityBase>(
+              std::forward<photons::MeanOpacityBase>(mean_opac_host_base),
+              time_unit, mass_unit, length_unit, temp_unit);
 
       auto funny_units = funny_units_host.GetOnDevice();
 
+      THEN("We can retrieve physical constants in code units") {
+        auto noncgs_rpc = funny_units.GetRuntimePhysicalConstants();
+        REQUIRE(FractionalDifference(noncgs_rpc.length, length_unit) <
+                EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.time, time_unit) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.mass, mass_unit) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.temp, temp_unit) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.na, 6.022141e+23) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.alpha, 7.297353e-03) <
+                EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.h, 2.871060e-33) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.hbar, 4.569434e-34) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.kb, 2.030877e-18) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.r_gas, 4.431243e+03) <
+                EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.qe, 4.803205e-10) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.c, 4.673571e+09) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.g_newt, 4.508065e-15) <
+                EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.me, 1.997672e-30) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.mp, 3.668030e-27) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.mn, 3.673086e-27) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.amu, 3.641533e-27) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.sb, 1.342760e+09) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.ar, 1.149237e+00) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.eV, 8.538896e-17) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.Fc, 1.189435e-09) < EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.nu_sigma0, 2.829094e-74) <
+                EPS_TEST);
+        REQUIRE(FractionalDifference(noncgs_rpc.gA, -1.23) < EPS_TEST);
+      }
+
       int n_wrong_h = 0;
 #ifdef PORTABILITY_STRATEGY_KOKKOS
       Kokkos::View<int, atomic_view> n_wrong_d("wrong");