IAP GITLAB

testMedium.cpp 3.76 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
/*
 * (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/PhysicalUnits.hpp>
#include <corsika/framework/geometry/Line.hpp>
#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
#include <corsika/framework/geometry/Vector.hpp>
#include <corsika/media/FlatExponential.hpp>
#include <corsika/media/HomogeneousMedium.hpp>
#include <corsika/media/IMediumModel.hpp>
#include <corsika/media/InhomogeneousMedium.hpp>
#include <corsika/media/NuclearComposition.hpp>
#include <corsika/media/MediumPropertyModel.hpp>
19
#include <corsika/media/MediumProperties.hpp>
20 21 22

#include <catch2/catch.hpp>

Ralf Ulrich's avatar
Ralf Ulrich committed
23 24
#include <SetupTestTrajectory.hpp>

25
using namespace corsika;
26 27 28

TEST_CASE("MediumProperties") {

Ralf Ulrich's avatar
Ralf Ulrich committed
29 30 31
  logging::set_level(logging::level::info);
  corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");

32 33 34 35 36 37 38 39 40 41 42 43 44
  // test access of medium properties via enum and class types

  const Medium type = Medium::AirDry1Atm;
  const MediumData& air = mediumData(type);
  CHECK(air.getIeff() == 85.7);
  CHECK(air.getCbar() == 10.5961);
  CHECK(air.getX0() == 1.7418);
  CHECK(air.getX1() == 4.2759);
  CHECK(air.getAA() == 0.10914);
  CHECK(air.getSK() == 3.3994);
  CHECK(air.getDlt0() == 0.0);
}

45 46
TEST_CASE("MediumPropertyModel w/ Homogeneous") {

Ralf Ulrich's avatar
Ralf Ulrich committed
47 48 49
  logging::set_level(logging::level::info);
  corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");

Ralf Ulrich's avatar
Ralf Ulrich committed
50
  CoordinateSystemPtr gCS = get_root_CoordinateSystem();
51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

  Point const gOrigin(gCS, {0_m, 0_m, 0_m});

  // setup our interface types
  using IModelInterface = IMediumPropertyModel<IMediumModel>;
  using AtmModel = MediumPropertyModel<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
66
  const Medium type = corsika::Medium::AirDry1Atm;
67 68 69 70 71

  // create the atmospheric model
  AtmModel medium(type, density, protonComposition);

  // and require that it is constant
72 73 74 75
  CHECK(type == medium.getMedium(Point(gCS, -10_m, 4_m, 35_km)));
  CHECK(type == medium.getMedium(Point(gCS, +210_m, 0_m, 7_km)));
  CHECK(type == medium.getMedium(Point(gCS, 0_m, 0_m, 0_km)));
  CHECK(type == medium.getMedium(Point(gCS, 100_km, 400_km, 350_km)));
76 77

  // a new refractive index
78
  const Medium type2 = corsika::Medium::StandardRock;
79 80

  // update the refractive index of this atmospheric model
81
  medium.setMedium(type2);
82 83

  // check that the returned refractive index is correct
84 85 86 87
  CHECK(type2 == medium.getMedium(Point(gCS, -10_m, 4_m, 35_km)));
  CHECK(type2 == medium.getMedium(Point(gCS, +210_m, 0_m, 7_km)));
  CHECK(type2 == medium.getMedium(Point(gCS, 0_m, 0_m, 0_km)));
  CHECK(type2 == medium.getMedium(Point(gCS, 100_km, 400_km, 350_km)));
88 89 90 91 92

  // define our axis vector
  Vector const axis(gCS, QuantityVector<dimensionless_d>(0, 0, 1));

  // check the density and nuclear composition
93 94
  CHECK(density == medium.getMassDensity(Point(gCS, 0_m, 0_m, 0_m)));
  medium.getNuclearComposition();
95 96

  // create a line of length 1 m
Ralf Ulrich's avatar
Ralf Ulrich committed
97 98
  Line const line(gOrigin,
                  VelocityVector(gCS, {1_m / second, 0_m / second, 0_m / second}));
99 100 101 102 103

  // the end time of our line
  auto const tEnd = 1_s;

  // and the associated trajectory
Ralf Ulrich's avatar
Ralf Ulrich committed
104
  setup::Trajectory const trajectory(line, tEnd);
105 106

  // and check the integrated grammage
Ralf Ulrich's avatar
Ralf Ulrich committed
107 108
  CHECK((medium.getIntegratedGrammage(trajectory, 3_m) / (density * 3_m)) == Approx(1));
  CHECK((medium.getArclengthFromGrammage(trajectory, density * 5_m) / 5_m) == Approx(1));
109
}