From f57e4951dcabc4d5d485976ad674bd25968017ac Mon Sep 17 00:00:00 2001
From: Remy Prechelt <prechelt@hawaii.edu>
Date: Sat, 27 Jun 2020 14:05:20 -1000
Subject: [PATCH] Implementation of magnetic field interface.
---
Environment/CMakeLists.txt | 2 +
Environment/IMagneticFieldModel.h | 47 +++++++++++++++++++++
Environment/UniformMagneticField.h | 67 ++++++++++++++++++++++++++++++
Environment/testEnvironment.cc | 46 ++++++++++++++++++++
4 files changed, 162 insertions(+)
create mode 100644 Environment/IMagneticFieldModel.h
create mode 100644 Environment/UniformMagneticField.h
diff --git a/Environment/CMakeLists.txt b/Environment/CMakeLists.txt
index e01188125..d5e69be6d 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 000000000..5941f3ba4
--- /dev/null
+++ b/Environment/IMagneticFieldModel.h
@@ -0,0 +1,47 @@
+/*
+ * (c) Copyright 2018 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 000000000..4e3ee298a
--- /dev/null
+++ b/Environment/UniformMagneticField.h
@@ -0,0 +1,67 @@
+/*
+ * (c) Copyright 2018 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 574907e59..d2434aac2 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,47 @@ 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);
+
+ // create our atmospheric model
+ AtmModel const medium(B0, 19.2_g / cube(1_cm), protonComposition);
+}
+
+TEST_CASE("UniformMagneticField w/ FlatExponential") {
+
+ // setup our interface types
+ using IModelInterface = IMagneticFieldModel<IMediumModel>;
+ using AtmModel = UniformMagneticField<FlatExponential<IModelInterface>>;
+
+ // the composition we use for the homogenous medium
+ NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
+ std::vector<float>{1.f});
+
+ // define our quantity vector
+ Vector const axis(gCS, QuantityVector<dimensionless_d>(0, 0, 1));
+
+ // the parameters of our exponential model
+ LengthType const lambda = 3_m;
+ auto const rho0 = 1_g / units::si::detail::static_pow<3>(1_cm);
+
+ // create a magnetic field vector
+ QuantityVector B0(23_T, 57_T, -4_T);
+
+ // create our atmospheric model
+ AtmModel const medium(B0, gOrigin, axis, rho0, lambda, protonComposition);
+
+ // check that the returned magnetic field is correct
+ REQUIRE(B0 == medium.GetMagneticField(Point(gCS, -10_m, 4_m, 35_km)));
+}
--
GitLab