uniform real distribution with units

Framework/Random/CMakeLists.txt

@@ -7,6 +7,7 @@ set (
 set (
   CORSIKArandom_HEADERS
   RNGManager.h
+  UniformRealDistribution.h
   )
 
 set (

Framework/Random/UniformRealDistribution.h

+#ifndef _include_random_distributions_h
+#define _include_random_distributions_h
+
+#include <corsika/units/PhysicalUnits.h>
+#include <random>
+
+namespace corsika::random {
+
+  template <class TQuantity>
+  class UniformRealDistribution {
+    using RealType = typename TQuantity::value_type;
+    std::uniform_real_distribution<RealType> dist{RealType(0.), RealType(1.)};
+
+    TQuantity const a, b;
+
+  public:
+    UniformRealDistribution(TQuantity b)
+        : a{TQuantity(phys::units::detail::magnitude_tag, 0)}
+        , b(b) {}
+    UniformRealDistribution(TQuantity a, TQuantity b)
+        : a(a)
+        , b(b) {}
+
+    template <class Generator>
+    TQuantity operator()(Generator& g) {
+      return a + dist(g) * (b - a);
+    }
+  };
+
+} // namespace corsika::random
+
+#endif

Framework/Random/testRandom.cc

-
 /**
  * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
  *
@@ -14,7 +13,11 @@
 #include <catch2/catch.hpp>
 
 #include <corsika/random/RNGManager.h>
+#include <corsika/random/UniformRealDistribution.h>
+#include <corsika/units/PhysicalUnits.h>
 #include <iostream>
+#include <limits>
+#include <random>
 
 using namespace corsika::random;
 
@@ -37,3 +40,46 @@ SCENARIO("random-number streams can be registered and retrieved") {
     }
   }
 }
+
+TEST_CASE("UniformRealDistribution") {
+  using namespace corsika::units::si;
+  std::mt19937 rng;
+
+  corsika::random::UniformRealDistribution<LengthType> dist(1_m, 2_m);
+
+  SECTION("range") {
+    corsika::random::UniformRealDistribution<LengthType> dist(1_m, 2_m);
+
+    LengthType min =
+        +1_m * std::numeric_limits<typename LengthType::value_type>::infinity();
+    LengthType max =
+        -1_m * std::numeric_limits<typename LengthType::value_type>::infinity();
+
+    for (int i{0}; i < 1'000'000; ++i) {
+      LengthType x = dist(rng);
+      min = std::min(min, x);
+      max = std::max(max, x);
+    }
+
+    CHECK(min / 1_m == Approx(1.));
+    CHECK(max / 2_m == Approx(1.));
+  }
+
+  SECTION("range") {
+    corsika::random::UniformRealDistribution<LengthType> dist(18_cm);
+
+    LengthType min =
+        +1_m * std::numeric_limits<typename LengthType::value_type>::infinity();
+    LengthType max =
+        -1_m * std::numeric_limits<typename LengthType::value_type>::infinity();
+
+    for (int i{0}; i < 1'000'000; ++i) {
+      LengthType x = dist(rng);
+      min = std::min(min, x);
+      max = std::max(max, x);
+    }
+
+    CHECK(min / 1_m == Approx(0.).margin(1e-3));
+    CHECK(max / 18_cm == Approx(1.));
+  }
+}