diff --git a/corsika/detail/framework/geometry/Path.inl b/corsika/detail/framework/geometry/Path.inl new file mode 100644 index 0000000000000000000000000000000000000000..c367531883c0bce37a488fb2901c6cc13c0bad09 --- /dev/null +++ b/corsika/detail/framework/geometry/Path.inl @@ -0,0 +1,66 @@ +/* + * (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 diff --git a/corsika/detail/framework/geometry/Point.inl b/corsika/detail/framework/geometry/Point.inl index c5c0a98a64491d3c6c752b647d60cdb1c2090e73..fe6e8b162383352523f7ba2399ef8685b0b5a2cb 100644 --- a/corsika/detail/framework/geometry/Point.inl +++ b/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 diff --git a/corsika/detail/media/ExponentialRefractiveIndex.inl b/corsika/detail/media/ExponentialRefractiveIndex.inl new file mode 100644 index 0000000000000000000000000000000000000000..760abb9e2410b3e9f1a1fced6560b3a75998779b --- /dev/null +++ b/corsika/detail/media/ExponentialRefractiveIndex.inl @@ -0,0 +1,29 @@ +/* + * (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 diff --git a/corsika/framework/geometry/Path.hpp b/corsika/framework/geometry/Path.hpp new file mode 100644 index 0000000000000000000000000000000000000000..620cbeeab2210141e8a51f423e9422688313e4d2 --- /dev/null +++ b/corsika/framework/geometry/Path.hpp @@ -0,0 +1,85 @@ +/* + * (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 diff --git a/corsika/framework/geometry/Point.hpp b/corsika/framework/geometry/Point.hpp index 2c684438438455985d0bebfbb95ed6bdafa4c86b..6eb8786c9df6accf7e654585a5bd4a8d8a957313 100644 --- a/corsika/framework/geometry/Point.hpp +++ b/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 diff --git a/corsika/media/ExponentialRefractiveIndex.hpp b/corsika/media/ExponentialRefractiveIndex.hpp new file mode 100644 index 0000000000000000000000000000000000000000..cc3ec838e9531a57d2ebd45d3434e00538481a06 --- /dev/null +++ b/corsika/media/ExponentialRefractiveIndex.hpp @@ -0,0 +1,55 @@ +/* + * (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> diff --git a/tests/framework/testGeometry.cpp b/tests/framework/testGeometry.cpp index e8847ee5e824242344147f46afeb37dfe0eb3708..809ae5e1391071b948a6307600cabb4bfc201136 100644 --- a/tests/framework/testGeometry.cpp +++ b/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 diff --git a/tests/media/testRefractiveIndex.cpp b/tests/media/testRefractiveIndex.cpp index fb1e71fec7b914fcb449f8a2bd2d06dce952222e..1aee5973e8dca12f3f45b743e47658ef324e5600 100644 --- a/tests/media/testRefractiveIndex.cpp +++ b/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)); +}