conversion function HEP to SI based on dim. analysis

Framework/Units/PhysicalConstants.h

@@ -37,7 +37,7 @@ namespace corsika::units::constants {
       phys::units::square(phys::units::second)};
 
   // Avogadro constant
-  constexpr quantity<dimensions<0, 0, 0, 0, 0, -1> > N_sub_A{Rep(6.02214199e+23L) / mole};
+  constexpr quantity<dimensions<0, 0, 0, 0, 0, -1>> N_sub_A{Rep(6.02214199e+23L) / mole};
 
   // elementary charge
   constexpr quantity<electric_charge_d> e{Rep(1.6021766208e-19L) * coulomb};
@@ -46,12 +46,16 @@ namespace corsika::units::constants {
   // constexpr quantity<hepenergy_d> eV{e / coulomb * joule};
 
   // Planck constant
-  constexpr quantity<dimensions<2, 1, -1> > h{Rep(6.62606876e-34L) * joule * second};
+  constexpr quantity<dimensions<2, 1, -1>> h{Rep(6.62606876e-34L) * joule * second};
+  constexpr quantity<dimensions<2, 1, -1>> hBar{h / (2 * M_PI)};
 
   // speed of light in a vacuum
   constexpr quantity<speed_d> c{Rep(299792458L) * meter / second};
   constexpr auto cSquared = c * c;
 
+  constexpr quantity<dimensions<1, 0, 0, 0, 0, 0, 0, 1>> hBarC{
+      Rep(1.973'269'78e-7L) * electronvolt * meter}; // from RPP 2018
+
   // unified atomic mass unit
   constexpr quantity<mass_d> u{Rep(1.6605402e-27L) * kilogram};
 

Framework/Units/PhysicalUnits.h

@@ -74,6 +74,51 @@ namespace corsika::units::si {
       phys::units::quantity<phys::units::dimensions<0, 0, -1>, double>;
   using InverseGrammageType =
       phys::units::quantity<phys::units::dimensions<2, -1, 0>, double>;
+
+  namespace detail {
+    template <int N, typename T>
+    auto constexpr static_pow([[maybe_unused]] T x) {
+      if constexpr (N == 0) {
+        return 1;
+      } else if constexpr (N > 0) {
+        return x * static_pow<N - 1, T>(x);
+      } else {
+        return 1 / static_pow<-N, T>(x);
+      }
+    }
+  } // namespace detail
+
+  template <typename DimFrom, typename DimTo>
+  auto constexpr ConversionFactorHEPToSI() {
+    static_assert(DimFrom::dim1 == 0 && DimFrom::dim2 == 0 && DimFrom::dim3 == 0 &&
+                      DimFrom::dim4 == 0 && DimFrom::dim5 == 0 && DimFrom::dim6 == 0 &&
+                      DimFrom::dim7 == 0,
+                  "must be a pure HEP type");
+
+    static_assert(
+        DimTo::dim4 == 0 && DimTo::dim5 == 0 && DimTo::dim6 == 0 && DimTo::dim7 == 0,
+        "conversion possible only into L, M, T dimensions");
+
+    int constexpr e = DimFrom::dim8; // HEP dim.
+
+    int constexpr l = DimTo::dim1; // SI length dim.
+    int constexpr m = DimTo::dim2; // SI mass dim.
+    int constexpr t = DimTo::dim3; // SI time dim.
+
+    int constexpr p = m;
+    int constexpr q = -m - t;
+    static_assert(q == l + e - 2 * m, "HEP/SI dimension mismatch!");
+
+    using namespace detail;
+    return static_pow<-e>(corsika::units::constants::hBarC) *
+           static_pow<p>(corsika::units::constants::hBar) *
+           static_pow<q>(corsika::units::constants::c);
+  }
+
+  template <typename DimFrom, typename DimTo>
+  auto constexpr ConvertHEPToSI(quantity<DimFrom> q) {
+    return ConversionFactorHEPToSI<DimFrom, DimTo>() * q;
+  }
 } // end namespace corsika::units::si
 
 /**

Framework/Units/testUnits.cc

@@ -115,4 +115,35 @@ TEST_CASE("PhysicalUnits", "[Units]") {
     REQUIRE(farAway > 100000_m);
     REQUIRE_FALSE(farAway < 1e19 * meter);
   }
+
+  SECTION("static_pow") {
+    using namespace corsika::units::si::detail;
+    double x = 235.7913;
+    REQUIRE(1 == static_pow<0, double>(x));
+    REQUIRE(x == static_pow<1, double>(x));
+    REQUIRE(x * x == static_pow<2, double>(x));
+    REQUIRE(1 / x == static_pow<-1, double>(x));
+    REQUIRE(1 / x / x == static_pow<-2, double>(x));
+  }
+
+  SECTION("HEP/SI conversion") {
+    auto const invEnergy = 1 / 197.326978_MeV; // should be convertible to length or time
+
+    LengthType const length =
+        ConvertHEPToSI<decltype(invEnergy)::dimension_type, LengthType::dimension_type>(
+            invEnergy);
+    REQUIRE((length / 1_fm) == Approx(1));
+
+    TimeType const time =
+        ConvertHEPToSI<decltype(invEnergy)::dimension_type, TimeType::dimension_type>(
+            invEnergy);
+    REQUIRE((time / (1_fm / corsika::units::constants::c)) == Approx(1));
+
+    auto const protonMass = 938.272'081'3_MeV; // convertible to mass or SI energy
+    MassType protonMassSI =
+        ConvertHEPToSI<decltype(protonMass)::dimension_type, MassType::dimension_type>(
+            protonMass);
+    REQUIRE((protonMassSI / 1.672'621'898e-27_kg) == Approx(1));
+    REQUIRE((protonMassSI / (1.007'276 * corsika::units::constants::u)) == Approx(1));
+  }
 }