diff --git a/Environment/CMakeLists.txt b/Environment/CMakeLists.txt index d5e69be6dd383bd5c1ce126331a59c66bbdb3be6..e90160a67a2c74c55c4a40fbad479cfe86d05416 100644 --- a/Environment/CMakeLists.txt +++ b/Environment/CMakeLists.txt @@ -23,6 +23,8 @@ set ( ShowerAxis.h IMagneticFieldModel.h UniformMagneticField.h + IRefractiveIndexModel.h + UniformRefractiveIndex.h ) set ( diff --git a/Environment/IRefractiveIndexModel.h b/Environment/IRefractiveIndexModel.h new file mode 100644 index 0000000000000000000000000000000000000000..465dd43e16aead86c4f9a0de8d2a5e10a2eaf492 --- /dev/null +++ b/Environment/IRefractiveIndexModel.h @@ -0,0 +1,42 @@ +/* + * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu + * + * See file AUTHORS for a list of contributors. + * + * 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 <corsika/geometry/Point.h> +#include <corsika/units/PhysicalUnits.h> + +namespace corsika::environment { + + /** + * An interface for refractive index models. + * + * This is the base interface for refractive index mixins. + * + */ + template <typename Model> + class IRefractiveIndexModel : public Model { + + public: + /** + * Evaluate the refractive index at a given location. + * + * @param point The location to evaluate at. + * @returns The refractive index at this point. + */ + virtual double GetRefractiveIndex(corsika::geometry::Point const&) const = 0; + + /** + * A virtual default destructor. + */ + virtual ~IRefractiveIndexModel() = default; + + }; // END: class IRefractiveIndexModel + +} // namespace corsika::environment diff --git a/Environment/UniformRefractiveIndex.h b/Environment/UniformRefractiveIndex.h new file mode 100644 index 0000000000000000000000000000000000000000..030ef593798437aba2e9fb0047ac2d5ba4229386 --- /dev/null +++ b/Environment/UniformRefractiveIndex.h @@ -0,0 +1,62 @@ +/* + * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu + * + * See file AUTHORS for a list of contributors. + * + * 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 <corsika/environment/IRefractiveIndexModel.h> + +namespace corsika::environment { + + /** + * A uniform refractive index. + * + * This class returns the same refractive index + * for all evaluated locations. + * + */ + template <typename T> + class UniformRefractiveIndex : public T { + + double n_; ///< The constant refractive index that we use. + + public: + /** + * Construct a UniformRefractiveIndex. + * + * This is initialized with a fixed refractive index + * and returns this refractive index at all locations. + * + * @param field The refractive index to return everywhere. + */ + template <typename... Args> + UniformRefractiveIndex(double const n, Args&&... args) + : T(std::forward<Args>(args)...) + , n_(n) {} + + /** + * Evaluate the refractive index at a given location. + * + * @param point The location to evaluate at. + * @returns The refractive index at this point. + */ + double GetRefractiveIndex(corsika::geometry::Point const&) const final override { + return n_; + } + + /** + * Set the refractive index returned by this instance. + * + * @param point The location to evaluate at. + * @returns The refractive index at this location. + */ + void SetRefractiveIndex(double const& n) { n_ = n; } + + }; // END: class RefractiveIndex + +} // namespace corsika::environment diff --git a/Environment/testEnvironment.cc b/Environment/testEnvironment.cc index 55ef5fba4f2473e6e59141953c1b191d2ceadaec..31955e15ac8618e431315fb56685429f6a07f046 100644 --- a/Environment/testEnvironment.cc +++ b/Environment/testEnvironment.cc @@ -13,12 +13,14 @@ #include <corsika/environment/HomogeneousMedium.h> #include <corsika/environment/IMagneticFieldModel.h> #include <corsika/environment/IMediumModel.h> +#include <corsika/environment/IRefractiveIndexModel.h> #include <corsika/environment/InhomogeneousMedium.h> #include <corsika/environment/LayeredSphericalAtmosphereBuilder.h> #include <corsika/environment/LinearApproximationIntegrator.h> #include <corsika/environment/NuclearComposition.h> #include <corsika/environment/SlidingPlanarExponential.h> #include <corsika/environment/UniformMagneticField.h> +#include <corsika/environment/UniformRefractiveIndex.h> #include <corsika/environment/VolumeTreeNode.h> #include <corsika/geometry/Line.h> #include <corsika/geometry/RootCoordinateSystem.h> @@ -251,6 +253,11 @@ TEST_CASE("UniformMagneticField w/ Homogeneous Medium") { // create our atmospheric model AtmModel medium(B0, density, protonComposition); + // and test at several locations + REQUIRE(B0 == medium.GetMagneticField(Point(gCS, -10_m, 4_m, 35_km))); + REQUIRE(B0 == medium.GetMagneticField(Point(gCS, 1000_km, -1000_km, 1000_km))); + REQUIRE(B0 == medium.GetMagneticField(Point(gCS, 0_m, 0_m, 0_m))); + // create a new magnetic field vector QuantityVector B1(23_T, 57_T, -4_T); @@ -280,3 +287,62 @@ TEST_CASE("UniformMagneticField w/ Homogeneous Medium") { REQUIRE((medium.IntegratedGrammage(trajectory, 3_m) / (density * 3_m)) == Approx(1)); REQUIRE((medium.ArclengthFromGrammage(trajectory, density * 5_m) / 5_m) == Approx(1)); } + +TEST_CASE("UniformRefractiveIndex w/ Homogeneous") { + + // setup our interface types + using IModelInterface = IRefractiveIndexModel<IMediumModel>; + using AtmModel = UniformRefractiveIndex<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}); + + // the refrative index that we use + const double n{1.000327}; + + // create the atmospheric model + AtmModel medium(n, density, protonComposition); + + // and require that it is constant + REQUIRE(n == medium.GetRefractiveIndex(Point(gCS, -10_m, 4_m, 35_km))); + REQUIRE(n == medium.GetRefractiveIndex(Point(gCS, +210_m, 0_m, 7_km))); + REQUIRE(n == medium.GetRefractiveIndex(Point(gCS, 0_m, 0_m, 0_km))); + REQUIRE(n == medium.GetRefractiveIndex(Point(gCS, 100_km, 400_km, 350_km))); + + // a new refractive index + const double n2{2.3472123}; + + // update the refractive index of this atmospheric model + medium.SetRefractiveIndex(n2); + + // check that the returned refractive index is correct + REQUIRE(n2 == medium.GetRefractiveIndex(Point(gCS, -10_m, 4_m, 35_km))); + REQUIRE(n2 == medium.GetRefractiveIndex(Point(gCS, +210_m, 0_m, 7_km))); + REQUIRE(n2 == medium.GetRefractiveIndex(Point(gCS, 0_m, 0_m, 0_km))); + REQUIRE(n2 == medium.GetRefractiveIndex(Point(gCS, 100_km, 400_km, 350_km))); + + // 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(); + + // 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 + Trajectory<Line> const trajectory(line, tEnd); + + // and check the integrated grammage + REQUIRE((medium.IntegratedGrammage(trajectory, 3_m) / (density * 3_m)) == Approx(1)); + REQUIRE((medium.ArclengthFromGrammage(trajectory, density * 5_m) / 5_m) == Approx(1)); +}