testMagneticField.cpp

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

#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
#include <corsika/media/HomogeneousMedium.hpp>
#include <corsika/media/IMediumModel.hpp>
#include <corsika/media/UniformMagneticField.hpp>
#include <corsika/media/IMagneticFieldModel.hpp>
#include <corsika/media/VolumeTreeNode.hpp>
#include <corsika/media/GeomagneticModel.hpp>

#include <SetupTestTrajectory.hpp>
#include <corsika/setup/SetupTrajectory.hpp>

#include <catch2/catch.hpp>

using namespace corsika;

TEST_CASE("UniformMagneticField w/ Homogeneous Medium") {

  logging::set_level(logging::level::trace);

  CoordinateSystemPtr const& gCS = get_root_CoordinateSystem();
  Point const gOrigin(gCS, {0_m, 0_m, 0_m});

  SECTION("UniformMagneticField interface") {

    // setup our interface types
    using IModelInterface = IMagneticFieldModel<IMediumModel>;
    using AtmModel = UniformMagneticField<HomogeneousMedium<IModelInterface>>;

    // the composition we use for the homogenous medium
    NuclearComposition const protonComposition({Code::Proton}, {1.});

    // create a magnetic field vector
    Vector B0(gCS, 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);

    // and test at several locations
    CHECK(B0.getComponents(gCS) ==
          medium.getMagneticField(Point(gCS, -10_m, 4_m, 35_km)).getComponents(gCS));
    CHECK(B0.getComponents(gCS) ==
          medium.getMagneticField(Point(gCS, 1000_km, -1000_km, 1000_km))
              .getComponents(gCS));
    CHECK(B0.getComponents(gCS) ==
          medium.getMagneticField(Point(gCS, 0_m, 0_m, 0_m)).getComponents(gCS));

    // create a new magnetic field vector
    Vector B1(gCS, 23_T, 57_T, -4_T);

    // and update this atmospheric model
    medium.setMagneticField(B1);

    // and test at several locations
    CHECK(B1.getComponents(gCS) ==
          medium.getMagneticField(Point(gCS, -10_m, 4_m, 35_km)).getComponents(gCS));
    CHECK(B1.getComponents(gCS) ==
          medium.getMagneticField(Point(gCS, 1000_km, -1000_km, 1000_km))