Merge branch 'rprechelt-magnetic-field-mixins' into 'master'

Implementation of a templated magnetic field interface and UniformMagneticField type.

See merge request AirShowerPhysics/corsika!211

Environment/CMakeLists.txt

@@ -21,6 +21,8 @@ set (
   SlidingPlanarExponential.h
   LayeredSphericalAtmosphereBuilder.h
   ShowerAxis.h
+  IMagneticFieldModel.h
+  UniformMagneticField.h
   )
 
 set (

Environment/IMagneticFieldModel.h

+/*
+ * (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

Environment/UniformMagneticField.h

+/*
+ * (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

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));
+}