diff --git a/Environment/CMakeLists.txt b/Environment/CMakeLists.txt index e01188125c647e1f802edf844ce69b7014e70a49..d5e69be6dd383bd5c1ce126331a59c66bbdb3be6 100644 --- a/Environment/CMakeLists.txt +++ b/Environment/CMakeLists.txt @@ -21,6 +21,8 @@ set ( SlidingPlanarExponential.h LayeredSphericalAtmosphereBuilder.h ShowerAxis.h + IMagneticFieldModel.h + UniformMagneticField.h ) set ( diff --git a/Environment/IMagneticFieldModel.h b/Environment/IMagneticFieldModel.h new file mode 100644 index 0000000000000000000000000000000000000000..a2666ba267f530e7a56346f1ad7faf141ca9d7c1 --- /dev/null +++ b/Environment/IMagneticFieldModel.h @@ -0,0 +1,47 @@ +/* + * (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 magnetic field models. + * + * This is the base interface for magnetic field model mixins. + * + */ + template <typename Model> + class IMagneticFieldModel : public Model { + + // a type-alias for a magnetic field vector + using MagneticFieldVector = + corsika::geometry::QuantityVector<corsika::units::si::magnetic_flux_density_d>; + + public: + /** + * Evaluate the magnetic field at a given location. + * + * @param point The location to evaluate the field at. + * @returns The magnetic field vector at that point. + */ + virtual auto GetMagneticField(corsika::geometry::Point const&) const + -> MagneticFieldVector = 0; + + /** + * A virtual default destructor. + */ + virtual ~IMagneticFieldModel() = default; + + }; // END: class MagneticField + +} // namespace corsika::environment diff --git a/Environment/UniformMagneticField.h b/Environment/UniformMagneticField.h new file mode 100644 index 0000000000000000000000000000000000000000..31f2a1b94b78cb6ef5e7e318019c48401f87beb9 --- /dev/null +++ b/Environment/UniformMagneticField.h @@ -0,0 +1,67 @@ +/* + * (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/IMagneticFieldModel.h> + +namespace corsika::environment { + + /** + * A uniform (constant) magnetic field. + * + * This class returns the same magnetic field vector + * for all evaluated locations. + * + */ + template <typename T> + class UniformMagneticField : public T { + + // a type-alias for a magnetic field vector + using MagneticFieldVector = + corsika::geometry::QuantityVector<corsika::units::si::magnetic_flux_density_d>; + + MagneticFieldVector B_; ///< The constant magnetic field we use. + + public: + /** + * Construct a UniformMagneticField. + * + * This is initialized with a fixed magnetic field + * and returns this magnetic field at all locations. + * + * @param field The fixed magnetic field to return. + */ + template <typename... Args> + UniformMagneticField(MagneticFieldVector const& B, Args&&... args) + : T(std::forward<Args>(args)...) + , B_(B) {} + + /** + * Evaluate the magnetic field at a given location. + * + * @param point The location to evaluate the field at. + * @returns The magnetic field vector. + */ + MagneticFieldVector GetMagneticField( + corsika::geometry::Point const&) const final override { + return B_; + } + + /** + * Set the magnetic field returned by this instance. + * + * @param point The location to evaluate the field at. + * @returns The magnetic field vector. + */ + auto SetMagneticField(MagneticFieldVector const& Bfield) -> void { B_ = Bfield; } + + }; // END: class MagneticField + +} // namespace corsika::environment diff --git a/Environment/testEnvironment.cc b/Environment/testEnvironment.cc index 574907e5914ec7e97628023483f4dd7572dde836..55ef5fba4f2473e6e59141953c1b191d2ceadaec 100644 --- a/Environment/testEnvironment.cc +++ b/Environment/testEnvironment.cc @@ -11,12 +11,14 @@ #include <corsika/environment/DensityFunction.h> #include <corsika/environment/FlatExponential.h> #include <corsika/environment/HomogeneousMedium.h> +#include <corsika/environment/IMagneticFieldModel.h> #include <corsika/environment/IMediumModel.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/VolumeTreeNode.h> #include <corsika/geometry/Line.h> #include <corsika/geometry/RootCoordinateSystem.h> @@ -229,3 +231,52 @@ TEST_CASE("LayeredSphericalAtmosphereBuilder") { univ->GetContainingNode(Point(gCS, 0_m, 0_m, R + 24_km))->GetVolume()) .GetRadius() == R + 30_km); } + +TEST_CASE("UniformMagneticField w/ Homogeneous Medium") { + + // setup our interface types + using IModelInterface = IMagneticFieldModel<IMediumModel>; + using AtmModel = UniformMagneticField<HomogeneousMedium<IModelInterface>>; + + // the composition we use for the homogenous medium + NuclearComposition const protonComposition(std::vector<Code>{Code::Proton}, + std::vector<float>{1.f}); + + // create a magnetic field vector + QuantityVector B0(0_T, 0_T, 0_T); + + // the constant density + const auto density{19.2_g / cube(1_cm)}; + + // create our atmospheric model + AtmModel medium(B0, density, protonComposition); + + // create a new magnetic field vector + QuantityVector B1(23_T, 57_T, -4_T); + + // and update this atmospheric model + medium.SetMagneticField(B1); + + // and test at several locations + REQUIRE(B1 == medium.GetMagneticField(Point(gCS, -10_m, 4_m, 35_km))); + REQUIRE(B1 == medium.GetMagneticField(Point(gCS, 1000_km, -1000_km, 1000_km))); + REQUIRE(B1 == medium.GetMagneticField(Point(gCS, 0_m, 0_m, 0_m))); + + // 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)); +}