Merge branch...

Merge branch '367-geometry-and-environment-feature-updates-merge-to-refactored-version' into 'master'

Resolve "Geometry and environment feature updates - merge to refactored version"

Closes #367

See merge request !313

corsika/detail/framework/geometry/Path.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <deque>
+#include <corsika/framework/geometry/Point.hpp>
+
+namespace corsika {
+
+  Path::Path(Point const& point) { points_.push_front(point); }
+
+  Path::Path(std::deque<Point> const& points)
+      : points_(points) {
+    int dequesize_ = points.size();
+    if (dequesize_ == 0 || dequesize_ == 1) {
+      length_ = LengthType::zero();
+    } else if (dequesize_ == 2) {
+      length_ = (points.back() - points.front()).getNorm();
+    } else {
+      for (auto point = points.begin(); point != points.end() - 1; ++point) {
+        auto point_next = *(point + 1);
+        auto point_now = *(point);
+        length_ += (point_next - point_now).getNorm();
+      }
+    }
+  }
+
+  inline void Path::addToEnd(Point const& point) {
+    length_ += (point - points_.back()).getNorm();
+    points_.push_back(point);
+  }
+
+  inline void Path::removeFromEnd() {
+    auto lastpoint_ = points_.back();
+    points_.pop_back();
+    int dequesize_ = points_.size();
+    if (dequesize_ == 0 || dequesize_ == 1) {
+      length_ = LengthType::zero();
+    } else if (dequesize_ == 2) {
+      length_ = (points_.back() - points_.front()).getNorm();
+    } else {
+      length_ -= (lastpoint_ - points_.back()).getNorm();
+    }
+  }
+
+  inline LengthType Path::getLength() const { return length_; }
+
+  inline Point Path::getStart() const { return points_.front(); }
+
+  inline Point Path::getEnd() const { return points_.back(); }
+
+  inline Point Path::getPoint(std::size_t const index) const { return points_.at(index); }
+
+  inline auto Path::begin() { return points_.begin(); }
+
+  inline auto Path::end() { return points_.end(); }
+
+  inline int Path::getNSegments() const { return points_.size() - 1; }
+
+} // namespace corsika
\ No newline at end of file

corsika/detail/framework/geometry/Point.inl

@@ -102,4 +102,8 @@ namespace corsika {
     return os;
   }
 
-} // namespace corsika
+  inline LengthType distance(Point const& p1, Point const& p2) {
+    return (p1 - p2).getNorm();
+  }
+
+} // namespace corsika
\ No newline at end of file

corsika/detail/media/ExponentialRefractiveIndex.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <bits/stdc++.h>
+#include <corsika/media/IRefractiveIndexModel.hpp>
+
+namespace corsika {
+
+  template <typename T>
+  template <typename... Args>
+  ExponentialRefractiveIndex<T>::ExponentialRefractiveIndex(
+      double const n0, InverseLengthType const lambda, Args&&... args)
+      : T(std::forward<Args>(args)...)
+      , n_0(n0)
+      , lambda_(lambda) {}
+
+  template <typename T>
+  double ExponentialRefractiveIndex<T>::getRefractiveIndex(Point const& point) const {
+    return n_0 * exp((-lambda_) * point.getCoordinates().getZ());
+  }
+
+} // namespace corsika

corsika/framework/geometry/Path.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <deque>
+#include <corsika/framework/geometry/Point.hpp>
+
+namespace corsika {
+
+  /**
+   * This class represents a (potentially) curved path between two
+   * points using N >= 1 straight-line segments.
+   */
+  class Path {
+    std::deque<Point> points_;               ///< The points that make up this path.
+    LengthType length_ = LengthType::zero(); ///< The length of the path.
+  public:
+    /**
+     * Create a Path with a given starting Point.
+     */
+    Path(Point const& point);
+
+    /**
+     * Initialize a Path from an existing collection of Points.
+     */
+    Path(std::deque<Point> const& points);
+
+    /**
+     * Add a new Point to the end of the path.
+     */
+    inline void addToEnd(Point const& point);
+
+    /**
+     * Remove a point from the end of the path.
+     */
+    inline void removeFromEnd();
+
+    /**
+     * Get the total length of the path.
+     */
+    inline LengthType getLength() const;
+
+    /**
+     * Get the starting point of the path.
+     */
+    inline Point getStart() const;
+
+    /**
+     * Get the end point of the path.
+     */
+    inline Point getEnd() const;
+
+    /**
+     * Get a specific point of the path.
+     */
+    inline Point getPoint(std::size_t const index) const;
+
+    /**
+     * Return an iterator to the start of the Path.
+     */
+    inline auto begin();
+
+    /**
+     * Return an iterator to the end of the Path.
+     */
+    inline auto end();
+
+    /**
+     * Get the number of steps in the path.
+     * This is one less than the number of points that
+     * defines the path.
+     */
+    inline int getNSegments() const;
+
+  }; // class Path
+
+} // namespace corsika
+
+#include <corsika/detail/framework/geometry/Path.inl>
\ No newline at end of file

corsika/framework/geometry/Point.hpp

@@ -31,8 +31,8 @@ namespace corsika {
     /** \todo TODO: this should be private or protected, we don NOT want to expose numbers
      * without reference to outside:
      */
-    QuantityVector<length_d> const& getCoordinates() const;
-    QuantityVector<length_d>& getCoordinates();
+    inline QuantityVector<length_d> const& getCoordinates() const;
+    inline QuantityVector<length_d>& getCoordinates();
 
     /**
        this always returns a QuantityVector as triple
@@ -40,7 +40,7 @@ namespace corsika {
        \returns A value type QuantityVector, since it may have to create a temporary
        object to transform to pCS.
     **/
-    QuantityVector<length_d> getCoordinates(CoordinateSystemPtr const& pCS) const;
+    inline QuantityVector<length_d> getCoordinates(CoordinateSystemPtr const& pCS) const;
 
     /**
      * this always returns a QuantityVector as triple
@@ -49,7 +49,7 @@ namespace corsika {
      *   is actually transformed to pCS, if needed. Thus, there may be an implicit call to
      *   \ref rebase.
      **/
-    QuantityVector<length_d>& getCoordinates(CoordinateSystemPtr const& pCS);
+    inline QuantityVector<length_d>& getCoordinates(CoordinateSystemPtr const& pCS);
 
     /**
      * \name access coordinate components
@@ -60,25 +60,30 @@ namespace corsika {
      * created and destroyed each call. This can be avoided by using
      * \ref rebase first.
      **/
-    LengthType getX(CoordinateSystemPtr const& pCS) const;
-    LengthType getY(CoordinateSystemPtr const& pCS) const;
-    LengthType getZ(CoordinateSystemPtr const& pCS) const;
+    inline LengthType getX(CoordinateSystemPtr const& pCS) const;
+    inline LengthType getY(CoordinateSystemPtr const& pCS) const;
+    inline LengthType getZ(CoordinateSystemPtr const& pCS) const;
     /** \} **/
 
     /*!
      * transforms the Point into another CoordinateSystem by changing its
      * coordinates interally
      */
-    void rebase(CoordinateSystemPtr const& pCS);
+    inline void rebase(CoordinateSystemPtr const& pCS);
 
-    Point operator+(Vector<length_d> const& pVec) const;
+    inline Point operator+(Vector<length_d> const& pVec) const;
 
     /*!
      * returns the distance Vector between two points
      */
-    Vector<length_d> operator-(Point const& pB) const;
+    inline Vector<length_d> operator-(Point const& pB) const;
   };
 
+  /*
+   * calculates the distance between two points
+   */
+  inline LengthType distance(Point const& p1, Point const& p2);
+
 } // namespace corsika
 
-#include <corsika/detail/framework/geometry/Point.inl>
+#include <corsika/detail/framework/geometry/Point.inl>
\ No newline at end of file

corsika/media/ExponentialRefractiveIndex.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <bits/stdc++.h>
+#include <corsika/media/IRefractiveIndexModel.hpp>
+
+namespace corsika {
+
+  /**
+   * An exponential refractive index.
+   *
+   * This class returns the value of an exponential refractive index
+   * for all evaluated locations.
+   *
+   */
+  template <typename T>
+  class ExponentialRefractiveIndex : public T {
+
+    double n_0;                ///< n0 constant.
+    InverseLengthType lambda_; ///< lambda parameter.
+
+  public:
+    /**
+     * Construct an ExponentialRefractiveIndex.
+     *
+     * This is initialized with two parameters n_0 and lambda
+     * and returns the value of the exponential refractive index
+     * at a given point location.
+     *
+     * @param field    The refractive index to return to a given point.
+     */
+    template <typename... Args>
+    ExponentialRefractiveIndex(double const n0, InverseLengthType const lambda,
+                               Args&&... args);
+
+    /**
+     * Evaluate the refractive index at a given location using its z-coordinate.
+     *
+     * @param  point    The location to evaluate at.
+     * @returns    The refractive index at this point.
+     */
+    double getRefractiveIndex(Point const& point) const final override;
+
+  }; // END: class ExponentialRefractiveIndex
+
+} // namespace corsika
+
+#include <corsika/detail/media/ExponentialRefractiveIndex.inl>

tests/framework/testGeometry.cpp

@@ -14,6 +14,7 @@
 #include <corsika/framework/geometry/Line.hpp>
 #include <corsika/framework/geometry/Helix.hpp>
 #include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/Path.hpp>
 #include <corsika/framework/geometry/RootCoordinateSystem.hpp>
 #include <corsika/framework/geometry/Sphere.hpp>
 #include <corsika/framework/geometry/StraightTrajectory.hpp>
@@ -312,3 +313,78 @@ TEST_CASE("Geometry Trajectories") {
               .magnitude() == Approx(0).margin(absMargin));
   }
 }
+
+TEST_CASE("Distance between points") {
+  // define a known CS
+  CoordinateSystemPtr root = get_root_CoordinateSystem();
+
+  // define known points
+  Point p1(root, {0_m, 0_m, 0_m});
+  Point p2(root, {0_m, 0_m, 5_m});
+  Point p3(root, {1_m, 0_m, 0_m});
+  Point p4(root, {5_m, 0_m, 0_m});
+  Point p5(root, {0_m, 4_m, 0_m});
+  Point p6(root, {0_m, 5_m, 0_m});
+
+  // check distance() function
+  CHECK(distance(p1, p2) / 1_m == Approx(5));
+  CHECK(distance(p3, p4) / 1_m == Approx(4));
+  CHECK(distance(p5, p6) / 1_m == Approx(1));
+}
+
+TEST_CASE("Path") {
+  // define a known CS
+  CoordinateSystemPtr root = get_root_CoordinateSystem();
+
+  // define known points
+  Point p1(root, {0_m, 0_m, 0_m});
+  Point p2(root, {0_m, 0_m, 1_m});
+  Point p3(root, {0_m, 0_m, 2_m});
+  Point p4(root, {0_m, 0_m, 3_m});
+  Point p5(root, {0_m, 0_m, 4_m});
+  // define paths
+  Path P1(p1);
+  Path P2({p1, p2});
+  Path P3({p1, p2, p3});
+  // define deque that include point(s)
+  std::deque<Point> l1 = {p1};
+  std::deque<Point> l2 = {p1, p2};
+  std::deque<Point> l3 = {p1, p2, p3};
+
+  // test the various path constructors
+  SECTION("Test Constructors") {
+    // check constructor for one point
+    CHECK(std::equal(P1.begin(), P1.end(), l1.begin(),
+                     [](Point a, Point b) { return (a - b).getNorm() / 1_m < 1e-5; }));
+    // check constructor for collection of points
+    CHECK(std::equal(P3.begin(), P3.end(), l3.begin(),
+                     [](Point a, Point b) { return (a - b).getNorm() / 1_m < 1e-5; }));
+  }
+
+  // test the length and access methods
+  SECTION("Test getLength() and modifications to Path") {
+    P1.addToEnd(p2);
+    P2.removeFromEnd();
+    // Check modifications to path
+    CHECK(std::equal(P1.begin(), P1.end(), l2.begin(),
+                     [](Point a, Point b) { return (a - b).getNorm() / 1_m < 1e-5; }));
+    CHECK(std::equal(P2.begin(), P2.end(), l1.begin(),
+                     [](Point a, Point b) { return (a - b).getNorm() / 1_m < 1e-5; }));
+    // Check GetStart(), GetEnd(), GetPoint()
+    CHECK((P3.getEnd() - P3.getStart()).getNorm() / 1_m == Approx(2));
+    CHECK((P1.getPoint(1) - p2).getNorm() / 1_m == Approx(0));
+    // Check GetLength()
+    CHECK(P1.getLength() / 1_m == Approx(1));
+    CHECK(P2.getLength() / 1_m == Approx(0));
+    CHECK(P3.getLength() / 1_m == Approx(2));
+    P2.removeFromEnd();
+    CHECK(P2.getLength() / 1_m == Approx(0)); // Check the length of an empty path
+    P3.addToEnd(p4);
+    P3.addToEnd(p5);
+    CHECK(P3.getLength() / 1_m == Approx(4));
+    P3.removeFromEnd();
+    CHECK(P3.getLength() / 1_m == Approx(3)); // Check RemoveFromEnd() else case
+    // Check GetNSegments()
+    CHECK(P3.getNSegments() - 3 == Approx(0));
+  }
+}
\ No newline at end of file

tests/media/testRefractiveIndex.cpp

@@ -16,6 +16,7 @@
 #include <corsika/media/InhomogeneousMedium.hpp>
 #include <corsika/media/NuclearComposition.hpp>
 #include <corsika/media/UniformRefractiveIndex.hpp>
+#include <corsika/media/ExponentialRefractiveIndex.hpp>
 
 #include <SetupTestTrajectory.hpp>
 
@@ -89,3 +90,81 @@ TEST_CASE("UniformRefractiveIndex w/ Homogeneous") {
   CHECK((medium.getIntegratedGrammage(track, 3_m) / (density * 3_m)) == Approx(1));
   CHECK((medium.getArclengthFromGrammage(track, density * 5_m) / 5_m) == Approx(1));
 }
+
+TEST_CASE("ExponentialRefractiveIndex w/ Homogeneous medium") {
+
+  logging::set_level(logging::level::info);
+  corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
+
+  // get a CS and a point
+  CoordinateSystemPtr const& gCS = get_root_CoordinateSystem();
+
+  Point const gOrigin(gCS, {0_m, 0_m, 0_m});
+
+  // setup our interface types
+  using IModelInterface = IRefractiveIndexModel<IMediumModel>;
+  using AtmModel = ExponentialRefractiveIndex<HomogeneousMedium<IModelInterface>>;
+
+  // the constant density
+  const auto density{19.2_g / cube(1_cm)};
+
+  // the composition we use for the homogenous medium
+  NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
+                                             std::vector<float>{1.f});
+
+  // a new refractive index
+  const double n0{2};
+  const InverseLengthType lambda{56 / 1_m};
+
+  // create the atmospheric model and check refractive index
+  AtmModel medium(n0, lambda, density, protonComposition);
+  CHECK(n0 == medium.getRefractiveIndex(Point(gCS, 10_m, 14_m, 0_m)));
+
+  // another refractive index
+  const double n0_{1};
+  const InverseLengthType lambda_{1 / 1_km};
+
+  // create the atmospheric model and check refractive index
+  AtmModel medium_(n0_, lambda_, density, protonComposition);
+  CHECK(medium_.getRefractiveIndex(Point(gCS, 12_m, 56_m, 1_km)) == Approx(0.3678794));
+
+  // another refractive index
+  const double n0__{4};
+  const InverseLengthType lambda__{2 / 1_m};
+
+  // create the atmospheric model and check refractive index
+  AtmModel medium__(n0__, lambda__, density, protonComposition);
+  CHECK(medium__.getRefractiveIndex(Point(gCS, 0_m, 0_m, 35_km)) == Approx(0));
+
+  // define our axis vector
+  Vector const axis(gCS, QuantityVector<dimensionless_d>(0, 0, 1));
+
+  // check the density and nuclear composition
+  REQUIRE(density == medium.getMassDensity(Point(gCS, 0_m, 0_m, 0_m)));
+  medium.getNuclearComposition();
+  REQUIRE(density == medium_.getMassDensity(Point(gCS, 0_m, 0_m, 0_m)));
+  medium_.getNuclearComposition();
+  REQUIRE(density == medium__.getMassDensity(Point(gCS, 0_m, 0_m, 0_m)));
+  medium__.getNuclearComposition();
+
+  // create a line of length 1 m
+  Line const line(gOrigin, Vector<SpeedType::dimension_type>(
+                               gCS, {1_m / second, 0_m / second, 0_m / second}));
+
+  // the end time of our line
+  auto const tEnd = 1_s;
+
+  // and the associated trajectory
+  setup::Trajectory const track =
+      setup::testing::make_track<setup::Trajectory>(line, tEnd);
+  //  // and the associated trajectory
+  //  Trajectory<Line> const trajectory(line, tEnd);
+
+  // and check the integrated grammage
+  REQUIRE((medium.getIntegratedGrammage(track, 3_m) / (density * 3_m)) == Approx(1));
+  REQUIRE((medium.getArclengthFromGrammage(track, density * 5_m) / 5_m) == Approx(1));
+  REQUIRE((medium_.getIntegratedGrammage(track, 3_m) / (density * 3_m)) == Approx(1));
+  REQUIRE((medium_.getArclengthFromGrammage(track, density * 5_m) / 5_m) == Approx(1));
+  REQUIRE((medium__.getIntegratedGrammage(track, 3_m) / (density * 3_m)) == Approx(1));
+  REQUIRE((medium__.getArclengthFromGrammage(track, density * 5_m) / 5_m) == Approx(1));
+}