IAP GITLAB

Skip to content
Snippets Groups Projects
Commit a3ed8076 authored by Ralf M Ulrich's avatar Ralf M Ulrich Committed by Ralf Ulrich
Browse files

added media type property

parent 03aa719a
No related branches found
No related tags found
1 merge request!265Add medium type, as new property (air, water, rock...)
Hide whitespace changes
Inline Side-by-side
Showing
with 437 additions and 114 deletions
Documentation/Examples/cascade_example.cc
+ 14
− 12
View file @ a3ed8076
......@@ -45,7 +45,6 @@ using namespace corsika::process;
using namespace corsika::units;
using namespace corsika::particles;
using namespace corsika::random;
using namespace corsika::setup;
using namespace corsika::geometry;
using namespace corsika::environment;
......@@ -68,27 +67,30 @@ int main() {
random::RNGManager::GetInstance().RegisterRandomStream("cascade");
// setup environment, geometry
using EnvType = environment::Environment<setup::IEnvironmentModel>;
EnvType env;
setup::Environment env;
auto& universe = *(env.GetUniverse());
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
auto outerMedium = EnvType::CreateNode<Sphere>(
auto outerMedium = setup::EnvironmentType::CreateNode<Sphere>(
Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
using EnvironmentModel = environment::UniformMediumType<environment::UniformMagneticField<environment::HomogeneousMedium<setup::IEnvironment>>>;
// fraction of oxygen
const float fox = 0.20946;
auto const props =
outerMedium
->SetModelProperties<environment::HomogeneousMedium<setup::IEnvironmentModel>>(
1_kg / (1_m * 1_m * 1_m),
environment::NuclearComposition(
std::vector<particles::Code>{particles::Code::Nitrogen,
particles::Code::Oxygen},
std::vector<float>{1.f - fox, fox}));
auto innerMedium = EnvType::CreateNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 5000_m);
->SetModelProperties<EnvironmentModel>(environment::EMediumType::eAir,
Vector(rootCS, 0_T, 0_T, 0_T),
1_kg / (1_m * 1_m * 1_m),
environment::NuclearComposition(
std::vector<particles::Code>{particles::Code::Nitrogen,
particles::Code::Oxygen},
std::vector<float>{1.f - fox, fox}));
auto innerMedium = setup::Environment::CreateNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 5000_m);
innerMedium->SetModelProperties(props);
......
Environment/CMakeLists.txt
+ 3
− 0
View file @ a3ed8076
......@@ -25,6 +25,9 @@ set (
UniformMagneticField.h
IRefractiveIndexModel.h
UniformRefractiveIndex.h
MediumTypes.h
IMediumTypeModel.h
UniformMediumType.h
)
set (
......
Environment/Environment.h
+ 0
− 5
View file @ a3ed8076
......@@ -38,8 +38,6 @@ namespace corsika::environment {
, fUniverse(std::make_unique<BaseNodeType>(
std::make_unique<Universe>(fCoordinateSystem))) {}
// using IEnvironmentModel = corsika::setup::IEnvironmentModel;
auto& GetUniverse() { return fUniverse; }
auto const& GetUniverse() const { return fUniverse; }
......@@ -61,7 +59,4 @@ namespace corsika::environment {
typename BaseNodeType::VTNUPtr fUniverse;
};
// using SetupBaseNodeType = VolumeTreeNode<corsika::setup::IEnvironmentModel>;
// using SetupEnvironment = Environment<corsika::setup::IEnvironmentModel>;
} // namespace corsika::environment
Environment/IMediumType.h 0 → 100644
+ 41
− 0
View file @ a3ed8076
/*
* (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.
*/
#pragma once
#include <corsika/geometry/Point.h>
#include <corsika/units/PhysicalUnits.h>
namespace corsika::environment {
/**
* An interface for type of media, needed e.g. to determine energy losses
*
* This is the base interface for media types.
*
*/
template <typename Model>
class IMediumType : public Model {
public:
/**
* Evaluate the medium type at a given location.
*
* @param point The location to evaluate at.
* @returns The media type
*/
virtual double medium_type(corsika::geometry::Point const&) const = 0;
/**
* A virtual default destructor.
*/
virtual ~IMediumType() = default;
}; // END: class IMediumType
} // namespace corsika::environment
Environment/IMediumTypeModel.h 0 → 100644
+ 43
− 0
View file @ a3ed8076
/*
* (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.
*/
#pragma once
#include <corsika/environment/MediumTypes.h>
#include <corsika/geometry/Point.h>
#include <corsika/units/PhysicalUnits.h>
namespace corsika::environment {
/**
* An interface for type of media, needed e.g. to determine energy losses
*
* This is the base interface for media types.
*
*/
template <typename Model>
class IMediumTypeModel : public Model {
public:
/**
* Evaluate the medium type at a given location.
*
* @param point The location to evaluate at.
* @returns The media type
*/
virtual EMediumType medium_type(corsika::geometry::Point const&) const = 0;
/**
* A virtual default destructor.
*/
virtual ~IMediumTypeModel() = default;
}; // END: class IMediumTypeModel
} // namespace corsika::environment
Environment/MediumTypes.h 0 → 100644
+ 15
− 0
View file @ a3ed8076
/*
* (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.
*/
#pragma once
namespace corsika::environment {
enum class EMediumType { eUnknown, eAir, eWater, eIce, eRock };
}
Environment/UniformMediaType.h 0 → 100644
+ 61
− 0
View file @ a3ed8076
/*
* (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.
*/
#pragma once
#include <corsika/environment/IMediaTypeModel.h>
namespace corsika::environment {
/**
* A uniform refractive index.
*
* This class returns the same refractive index
* for all evaluated locations.
*
*/
template <typename T>
class UniformMediaType : public T {
EMediaType double n_; ///< The constant refractive index that we use.
public:
/**
* Construct a UniformMediaType.
*
* This is initialized with a fixed refractive index
* and returns this refractive index at all locations.
*
* @param field The refractive index to return everywhere.
*/
template <typename... Args>
UniformMediaType(double const n, Args&&... args)
: T(std::forward<Args>(args)...)
, n_(n) {}
/**
* Evaluate the refractive index at a given location.
*
* @param point The location to evaluate at.
* @returns The refractive index at this point.
*/
double GetMediaType(corsika::geometry::Point const&) const final override {
return n_;
}
/**
* Set the refractive index returned by this instance.
*
* @param point The location to evaluate at.
* @returns The refractive index at this location.
*/
void SetMediaType(double const& n) { n_ = n; }
}; // END: class MediaType
} // namespace corsika::environment
Environment/UniformMediumType.h 0 → 100644
+ 62
− 0
View file @ a3ed8076
/*
* (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.
*/
#pragma once
#include <corsika/environment/IMediumTypeModel.h>
#include <corsika/environment/MediumTypes.h>
namespace corsika::environment {
/**
* A uniform refractive index.
*
* This class returns the same refractive index
* for all evaluated locations.
*
*/
template <typename T>
class UniformMediumType : public T {
EMediumType type_; ///< The constant medium type
public:
/**
* Construct a UniformMediumType.
*
* This is initialized with a fixed medium type
* and returns this at all locations.
*
* @param field The medium type to return everywhere.
*/
template <typename... Args>
UniformMediumType(EMediumType const type, Args&&... args)
: T(std::forward<Args>(args)...)
, type_(type) {}
/**
* Evaluate the refractive index at a given location.
*
* @param point The location to evaluate at.
* @returns The refractive index at this point.
*/
EMediumType medium_type(corsika::geometry::Point const&) const final override {
return type_;
}
/**
* Set the refractive index returned by this instance.
*
* @param point The location to evaluate at.
* @returns The refractive index at this location.
*/
void set_medium_type(EMediumType const& type) { type_ = type; }
}; // END: class MediumType
} // namespace corsika::environment
Environment/testEnvironment.cc
+ 123
− 62
View file @ a3ed8076
......@@ -11,6 +11,7 @@
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/IMagneticFieldModel.h>
#include <corsika/environment/IMediumModel.h>
#include <corsika/environment/IMediumTypeModel.h>
#include <corsika/environment/IRefractiveIndexModel.h>
#include <corsika/environment/InhomogeneousMedium.h>
#include <corsika/environment/LayeredSphericalAtmosphereBuilder.h>
......@@ -18,6 +19,7 @@
#include <corsika/environment/NuclearComposition.h>
#include <corsika/environment/SlidingPlanarExponential.h>
#include <corsika/environment/UniformMagneticField.h>
#include <corsika/environment/UniformMediumType.h>
#include <corsika/environment/UniformRefractiveIndex.h>
#include <corsika/environment/VolumeTreeNode.h>
#include <corsika/geometry/Line.h>
......@@ -61,8 +63,8 @@ TEST_CASE("FlatExponential") {
gCS, {20_cm / second, 0_m / second, 0_m / second}));
Trajectory<Line> const trajectory(line, tEnd);
REQUIRE((medium.IntegratedGrammage(trajectory, 2_m) / (rho0 * 2_m)) == Approx(1));
REQUIRE((medium.ArclengthFromGrammage(trajectory, rho0 * 5_m) / 5_m) == Approx(1));
CHECK((medium.IntegratedGrammage(trajectory, 2_m) / (rho0 * 2_m)) == Approx(1));
CHECK((medium.ArclengthFromGrammage(trajectory, rho0 * 5_m) / 5_m) == Approx(1));
}
SECTION("vertical") {
......@@ -72,8 +74,8 @@ TEST_CASE("FlatExponential") {
LengthType const length = 2 * lambda;
GrammageType const exact = rho0 * lambda * (exp(length / lambda) - 1);
REQUIRE((medium.IntegratedGrammage(trajectory, length) / exact) == Approx(1));
REQUIRE((medium.ArclengthFromGrammage(trajectory, exact) / length) == Approx(1));
CHECK((medium.IntegratedGrammage(trajectory, length) / exact) == Approx(1));
CHECK((medium.ArclengthFromGrammage(trajectory, exact) / length) == Approx(1));
}
SECTION("escape grammage") {
......@@ -83,9 +85,9 @@ TEST_CASE("FlatExponential") {
GrammageType const escapeGrammage = rho0 * lambda;
REQUIRE(trajectory.NormalizedDirection().dot(axis).magnitude() < 0);
REQUIRE(medium.ArclengthFromGrammage(trajectory, 1.2 * escapeGrammage) ==
std::numeric_limits<typename GrammageType::value_type>::infinity() * 1_m);
CHECK(trajectory.NormalizedDirection().dot(axis).magnitude() < 0);
CHECK(medium.ArclengthFromGrammage(trajectory, 1.2 * escapeGrammage) ==
std::numeric_limits<typename GrammageType::value_type>::infinity() * 1_m);
}
SECTION("inclined") {
......@@ -96,8 +98,8 @@ TEST_CASE("FlatExponential") {
LengthType const length = 2 * lambda;
GrammageType const exact =
rho0 * lambda * (exp(cosTheta * length / lambda) - 1) / cosTheta;
REQUIRE((medium.IntegratedGrammage(trajectory, length) / exact) == Approx(1));
REQUIRE((medium.ArclengthFromGrammage(trajectory, exact) / length) == Approx(1));
CHECK((medium.IntegratedGrammage(trajectory, length) / exact) == Approx(1));
CHECK((medium.ArclengthFromGrammage(trajectory, exact) / length) == Approx(1));
}
}
......@@ -170,9 +172,9 @@ TEST_CASE("InhomogeneousMedium") {
DensityFunction<decltype(e), LinearApproximationIntegrator> const rho(e);
SECTION("DensityFunction") {
REQUIRE(e.Derivative<1>(gOrigin, direction) / (1_kg / 1_m / 1_m / 1_m / 1_m) ==
Approx(1));
REQUIRE(rho.EvaluateAt(gOrigin) == e(gOrigin));
CHECK(e.Derivative<1>(gOrigin, direction) / (1_kg / 1_m / 1_m / 1_m / 1_m) ==
Approx(1));
CHECK(rho.EvaluateAt(gOrigin) == e(gOrigin));
}
auto const exactGrammage = [](auto l) { return 1_m * rho0 * (exp(l / 1_m) - 1); };
......@@ -184,18 +186,18 @@ TEST_CASE("InhomogeneousMedium") {
InhomogeneousMedium<IMediumModel, decltype(rho)> const inhMedium(composition, rho);
SECTION("Integration") {
REQUIRE(rho.IntegrateGrammage(trajectory, l) / exactGrammage(l) ==
Approx(1).epsilon(1e-2));
REQUIRE(rho.ArclengthFromGrammage(trajectory, exactGrammage(l)) /
exactLength(exactGrammage(l)) ==
Approx(1).epsilon(1e-2));
REQUIRE(rho.MaximumLength(trajectory, 1e-2) >
l); // todo: write reasonable test when implementation is working
REQUIRE(rho.IntegrateGrammage(trajectory, l) ==
inhMedium.IntegratedGrammage(trajectory, l));
REQUIRE(rho.ArclengthFromGrammage(trajectory, 20_g / (1_cm * 1_cm)) ==
inhMedium.ArclengthFromGrammage(trajectory, 20_g / (1_cm * 1_cm)));
CHECK(rho.IntegrateGrammage(trajectory, l) / exactGrammage(l) ==
Approx(1).epsilon(1e-2));
CHECK(rho.ArclengthFromGrammage(trajectory, exactGrammage(l)) /
exactLength(exactGrammage(l)) ==
Approx(1).epsilon(1e-2));
CHECK(rho.MaximumLength(trajectory, 1e-2) >
l); // todo: write reasonable test when implementation is working
CHECK(rho.IntegrateGrammage(trajectory, l) ==
inhMedium.IntegratedGrammage(trajectory, l));
CHECK(rho.ArclengthFromGrammage(trajectory, 20_g / (1_cm * 1_cm)) ==
inhMedium.ArclengthFromGrammage(trajectory, 20_g / (1_cm * 1_cm)));
}
}
......@@ -208,27 +210,27 @@ TEST_CASE("LayeredSphericalAtmosphereBuilder") {
builder.addLinearLayer(2_km, 20_km);
builder.addLinearLayer(3_km, 30_km);
REQUIRE(builder.size() == 3);
CHECK(builder.size() == 3);
auto const builtEnv = builder.assemble();
auto const& univ = builtEnv.GetUniverse();
REQUIRE(builder.size() == 0);
CHECK(builder.size() == 0);
auto const R = builder.getEarthRadius();
REQUIRE(univ->GetChildNodes().size() == 1);
REQUIRE(univ->GetContainingNode(Point(gCS, 0_m, 0_m, R + 35_km)) == univ.get());
REQUIRE(dynamic_cast<Sphere const&>(
univ->GetContainingNode(Point(gCS, 0_m, 0_m, R + 8_km))->GetVolume())
.GetRadius() == R + 10_km);
REQUIRE(dynamic_cast<Sphere const&>(
univ->GetContainingNode(Point(gCS, 0_m, 0_m, R + 12_km))->GetVolume())
.GetRadius() == R + 20_km);
REQUIRE(dynamic_cast<Sphere const&>(
univ->GetContainingNode(Point(gCS, 0_m, 0_m, R + 24_km))->GetVolume())
.GetRadius() == R + 30_km);
CHECK(univ->GetChildNodes().size() == 1);
CHECK(univ->GetContainingNode(Point(gCS, 0_m, 0_m, R + 35_km)) == univ.get());
CHECK(dynamic_cast<Sphere const&>(
univ->GetContainingNode(Point(gCS, 0_m, 0_m, R + 8_km))->GetVolume())
.GetRadius() == R + 10_km);
CHECK(dynamic_cast<Sphere const&>(
univ->GetContainingNode(Point(gCS, 0_m, 0_m, R + 12_km))->GetVolume())
.GetRadius() == R + 20_km);
CHECK(dynamic_cast<Sphere const&>(
univ->GetContainingNode(Point(gCS, 0_m, 0_m, R + 24_km))->GetVolume())
.GetRadius() == R + 30_km);
}
TEST_CASE("UniformMagneticField w/ Homogeneous Medium") {
......@@ -251,13 +253,13 @@ TEST_CASE("UniformMagneticField w/ Homogeneous Medium") {
AtmModel medium(B0, density, protonComposition);
// and test at several locations
REQUIRE(B0.GetComponents(gCS) ==
medium.GetMagneticField(Point(gCS, -10_m, 4_m, 35_km)).GetComponents(gCS));
REQUIRE(
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));
REQUIRE(B0.GetComponents(gCS) ==
medium.GetMagneticField(Point(gCS, 0_m, 0_m, 0_m)).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);
......@@ -266,16 +268,16 @@ TEST_CASE("UniformMagneticField w/ Homogeneous Medium") {
medium.SetMagneticField(B1);
// and test at several locations
REQUIRE(B1.GetComponents(gCS) ==
medium.GetMagneticField(Point(gCS, -10_m, 4_m, 35_km)).GetComponents(gCS));
REQUIRE(
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)).GetComponents(gCS));
REQUIRE(B1.GetComponents(gCS) ==
medium.GetMagneticField(Point(gCS, 0_m, 0_m, 0_m)).GetComponents(gCS));
CHECK(B1.GetComponents(gCS) ==
medium.GetMagneticField(Point(gCS, 0_m, 0_m, 0_m)).GetComponents(gCS));
// check the density and nuclear composition
REQUIRE(density == medium.GetMassDensity(Point(gCS, 0_m, 0_m, 0_m)));
CHECK(density == medium.GetMassDensity(Point(gCS, 0_m, 0_m, 0_m)));
medium.GetNuclearComposition();
// create a line of length 1 m
......@@ -289,8 +291,8 @@ TEST_CASE("UniformMagneticField w/ Homogeneous Medium") {
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));
CHECK((medium.IntegratedGrammage(trajectory, 3_m) / (density * 3_m)) == Approx(1));
CHECK((medium.ArclengthFromGrammage(trajectory, density * 5_m) / 5_m) == Approx(1));
}
TEST_CASE("UniformRefractiveIndex w/ Homogeneous") {
......@@ -313,10 +315,10 @@ TEST_CASE("UniformRefractiveIndex w/ Homogeneous") {
AtmModel medium(n, density, protonComposition);
// and require that it is constant
REQUIRE(n == medium.GetRefractiveIndex(Point(gCS, -10_m, 4_m, 35_km)));
REQUIRE(n == medium.GetRefractiveIndex(Point(gCS, +210_m, 0_m, 7_km)));
REQUIRE(n == medium.GetRefractiveIndex(Point(gCS, 0_m, 0_m, 0_km)));
REQUIRE(n == medium.GetRefractiveIndex(Point(gCS, 100_km, 400_km, 350_km)));
CHECK(n == medium.GetRefractiveIndex(Point(gCS, -10_m, 4_m, 35_km)));
CHECK(n == medium.GetRefractiveIndex(Point(gCS, +210_m, 0_m, 7_km)));
CHECK(n == medium.GetRefractiveIndex(Point(gCS, 0_m, 0_m, 0_km)));
CHECK(n == medium.GetRefractiveIndex(Point(gCS, 100_km, 400_km, 350_km)));
// a new refractive index
const double n2{2.3472123};
......@@ -325,16 +327,75 @@ TEST_CASE("UniformRefractiveIndex w/ Homogeneous") {
medium.SetRefractiveIndex(n2);
// check that the returned refractive index is correct
REQUIRE(n2 == medium.GetRefractiveIndex(Point(gCS, -10_m, 4_m, 35_km)));
REQUIRE(n2 == medium.GetRefractiveIndex(Point(gCS, +210_m, 0_m, 7_km)));
REQUIRE(n2 == medium.GetRefractiveIndex(Point(gCS, 0_m, 0_m, 0_km)));
REQUIRE(n2 == medium.GetRefractiveIndex(Point(gCS, 100_km, 400_km, 350_km)));
CHECK(n2 == medium.GetRefractiveIndex(Point(gCS, -10_m, 4_m, 35_km)));
CHECK(n2 == medium.GetRefractiveIndex(Point(gCS, +210_m, 0_m, 7_km)));
CHECK(n2 == medium.GetRefractiveIndex(Point(gCS, 0_m, 0_m, 0_km)));
CHECK(n2 == medium.GetRefractiveIndex(Point(gCS, 100_km, 400_km, 350_km)));
// define our axis vector
Vector const axis(gCS, QuantityVector<dimensionless_d>(0, 0, 1));
// check the density and nuclear composition
CHECK(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
CHECK((medium.IntegratedGrammage(trajectory, 3_m) / (density * 3_m)) == Approx(1));
CHECK((medium.ArclengthFromGrammage(trajectory, density * 5_m) / 5_m) == Approx(1));
}
TEST_CASE("UniformMediumType w/ Homogeneous") {
// setup our interface types
using IModelInterface = IMediumTypeModel<IMediumModel>;
using AtmModel = UniformMediumType<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
const EMediumType type = EMediumType::eAir;
// create the atmospheric model
AtmModel medium(type, density, protonComposition);
// and require that it is constant
CHECK(type == medium.medium_type(Point(gCS, -10_m, 4_m, 35_km)));
CHECK(type == medium.medium_type(Point(gCS, +210_m, 0_m, 7_km)));
CHECK(type == medium.medium_type(Point(gCS, 0_m, 0_m, 0_km)));
CHECK(type == medium.medium_type(Point(gCS, 100_km, 400_km, 350_km)));
// a new refractive index
const EMediumType type2 = EMediumType::eRock;
// update the refractive index of this atmospheric model
medium.set_medium_type(type2);
// check that the returned refractive index is correct
CHECK(type2 == medium.medium_type(Point(gCS, -10_m, 4_m, 35_km)));
CHECK(type2 == medium.medium_type(Point(gCS, +210_m, 0_m, 7_km)));
CHECK(type2 == medium.medium_type(Point(gCS, 0_m, 0_m, 0_km)));
CHECK(type2 == medium.medium_type(Point(gCS, 100_km, 400_km, 350_km)));
// define our axis vector
Vector const axis(gCS, QuantityVector<dimensionless_d>(0, 0, 1));
// check the density and nuclear composition
REQUIRE(density == medium.GetMassDensity(Point(gCS, 0_m, 0_m, 0_m)));
CHECK(density == medium.GetMassDensity(Point(gCS, 0_m, 0_m, 0_m)));
medium.GetNuclearComposition();
// create a line of length 1 m
......@@ -348,6 +409,6 @@ TEST_CASE("UniformRefractiveIndex w/ Homogeneous") {
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));
CHECK((medium.IntegratedGrammage(trajectory, 3_m) / (density * 3_m)) == Approx(1));
CHECK((medium.ArclengthFromGrammage(trajectory, density * 5_m) / 5_m) == Approx(1));
}
Processes/CONEXSourceCut/testCONEXSourceCut.cc
+ 15
− 2
View file @ a3ed8076
......@@ -6,18 +6,28 @@
* the license.
*/
#include <corsika/setup/SetupEnvironment.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/LayeredSphericalAtmosphereBuilder.h>
#include <corsika/environment/UniformMediumType.h>
#include <corsika/environment/UniformMagneticField.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/geometry/Vector.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/conex_source_cut/CONEXSourceCut.h>
#include <corsika/process/sibyll/Interaction.h>
#include <corsika/process/sibyll/NuclearInteraction.h>
#include <corsika/random/RNGManager.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/utl/CorsikaFenv.h>
#include <catch2/catch.hpp>
using namespace corsika;
......@@ -32,11 +42,14 @@ TEST_CASE("CONEXSourceCut") {
feenableexcept(FE_INVALID);
// setup environment, geometry
using EnvType = Environment<setup::IEnvironmentModel>;
EnvType env;
setup::Environment env;
auto& universe = *(env.GetUniverse());
using EnvironmentModel = environment::UniformMediumType<environment::UniformMagneticField<environment::InhomogeneousMedium<setup::IEnvironment>>>;
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
Point const center{rootCS, 0_m, 0_m, 0_m};
environment::LayeredSphericalAtmosphereBuilder builder{center, conex::earthRadius};
builder.setNuclearComposition(
{{particles::Code::Nitrogen, particles::Code::Oxygen},
{0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
......
Processes/Proposal/CMakeLists_PROPOSAL.txt
+ 2
− 2
View file @ a3ed8076
......@@ -55,8 +55,8 @@ else (IN_CORSIKA8)
endif (IN_CORSIKA8)
add_library(PROPOSAL::PROPOSAL ALIAS PROPOSAL)
target_compile_features(PROPOSAL PUBLIC cxx_std_11)
set_target_properties(PROPOSAL PROPERTIES CXX_EXTENSIONS OFF)
#target_compile_features(PROPOSAL PUBLIC cxx_std_11)
#set_target_properties(PROPOSAL PROPERTIES CXX_EXTENSIONS OFF)
target_include_directories(
PROPOSAL PUBLIC
$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>
......
Processes/Pythia/testPythia8.cc
+ 19
− 10
View file @ a3ed8076
......@@ -15,9 +15,11 @@
#include <corsika/particles/ParticleProperties.h>
#include <corsika/geometry/Point.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/utl/CorsikaFenv.h>
#include <catch2/catch.hpp>
TEST_CASE("Pythia", "[processes]") {
......@@ -77,10 +79,13 @@ TEST_CASE("Pythia", "[processes]") {
#include <corsika/particles/ParticleProperties.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/setup/SetupEnvironment.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/environment/UniformMediumType.h>
#include <corsika/environment/UniformMagneticField.h>
using namespace corsika;
using namespace corsika::units::si;
......@@ -97,24 +102,28 @@ auto sumMomentum(TStackView const& view, geometry::CoordinateSystem const& vCS)
TEST_CASE("pythia process") {
// setup environment, geometry
environment::Environment<environment::IMediumModel> env;
geometry::CoordinateSystem const& cs = env.GetCoordinateSystem();
setup::Environment env;
auto& universe = *(env.GetUniverse());
using EnvironmentModel = environment::UniformMediumType<environment::UniformMagneticField<environment::HomogeneousMedium<setup::IEnvironment>>>;
auto theMedium =
environment::Environment<environment::IMediumModel>::CreateNode<geometry::Sphere>(
geometry::Point{cs, 0_m, 0_m, 0_m},
setup::Environment::CreateNode<geometry::Sphere>(
geometry::Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
1_km * std::numeric_limits<double>::infinity());
using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
theMedium->SetModelProperties<MyHomogeneousModel>(
theMedium->SetModelProperties<EnvironmentModel>(
environment::EMediumType::eAir,
geometry::Vector(env.GetCoordinateSystem(), 0_T, 0_T, 0_T),
1_kg / (1_m * 1_m * 1_m),
environment::NuclearComposition(
std::vector<particles::Code>{particles::Code::Hydrogen},
std::vector<float>{1.}));
std::vector<particles::Code>{particles::Code::Oxygen}, std::vector<float>{1.}));
auto const* nodePtr = theMedium.get();
universe.AddChild(std::move(theMedium));
auto const* nodePtr = theMedium.get(); // save the medium for later use before moving it
const geometry::CoordinateSystem& cs = env.GetCoordinateSystem();
SECTION("pythia decay") {
feenableexcept(FE_INVALID);
setup::Stack stack;
......
Processes/QGSJetII/testQGSJetII.cc
+ 9
− 5
View file @ a3ed8076
......@@ -111,12 +111,14 @@ TEST_CASE("QgsjetII", "[processes]") {
#include <corsika/particles/ParticleProperties.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupEnvironment.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/process/qgsjetII/qgsjet-II-04.h>
#include <corsika/environment/UniformMediumType.h>
#include <corsika/environment/UniformMagneticField.h>
using namespace corsika::units::si;
using namespace corsika::units;
......@@ -124,16 +126,18 @@ using namespace corsika::units;
TEST_CASE("QgsjetIIInterface", "[processes]") {
// setup environment, geometry
environment::Environment<environment::IMediumModel> env;
setup::Environment env;
auto& universe = *(env.GetUniverse());
using EnvironmentModel = environment::UniformMediumType<environment::UniformMagneticField<environment::HomogeneousMedium<setup::IEnvironment>>>;
auto theMedium =
environment::Environment<environment::IMediumModel>::CreateNode<geometry::Sphere>(
setup::Environment::CreateNode<geometry::Sphere>(
geometry::Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
1_km * std::numeric_limits<double>::infinity());
using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
theMedium->SetModelProperties<MyHomogeneousModel>(
theMedium->SetModelProperties<EnvironmentModel>(
environment::EMediumType::eAir,
geometry::Vector(env.GetCoordinateSystem(), 0_T, 0_T, 0_T),
1_kg / (1_m * 1_m * 1_m),
environment::NuclearComposition(
std::vector<particles::Code>{particles::Code::Oxygen}, std::vector<float>{1.}));
......
Processes/Sibyll/NuclearInteraction.cc
+ 7
− 7
View file @ a3ed8076
......@@ -43,7 +43,7 @@ namespace corsika::process::sibyll {
}
template <>
void NuclearInteraction<SetupEnvironment>::PrintCrossSectionTable(
void NuclearInteraction<setup::Environment>::PrintCrossSectionTable(
corsika::particles::Code pCode) {
using namespace corsika::particles;
const int k = targetComponentsIndex_.at(pCode);
......@@ -68,7 +68,7 @@ namespace corsika::process::sibyll {
}
template <>
void NuclearInteraction<SetupEnvironment>::InitializeNuclearCrossSections() {
void NuclearInteraction<setup::Environment>::InitializeNuclearCrossSections() {
using namespace corsika::particles;
using namespace units::si;
......@@ -134,7 +134,7 @@ namespace corsika::process::sibyll {
}
template <>
units::si::CrossSectionType NuclearInteraction<SetupEnvironment>::ReadCrossSectionTable(
units::si::CrossSectionType NuclearInteraction<setup::Environment>::ReadCrossSectionTable(
const int ia, particles::Code pTarget, units::si::HEPEnergyType elabnuc) {
using namespace corsika::particles;
using namespace units::si;
......@@ -154,7 +154,7 @@ namespace corsika::process::sibyll {
template <>
template <>
tuple<units::si::CrossSectionType, units::si::CrossSectionType>
NuclearInteraction<SetupEnvironment>::GetCrossSection(Particle const& vP,
NuclearInteraction<setup::Environment>::GetCrossSection(Particle const& vP,
const particles::Code TargetId) {
using namespace units::si;
if (vP.GetPID() != particles::Code::Nucleus)
......@@ -195,7 +195,7 @@ namespace corsika::process::sibyll {
template <>
template <>
units::si::GrammageType NuclearInteraction<SetupEnvironment>::GetInteractionLength(
units::si::GrammageType NuclearInteraction<setup::Environment>::GetInteractionLength(
Particle const& vP) {
using namespace units;
......@@ -615,8 +615,8 @@ namespace corsika::process::sibyll {
}
template <>
NuclearInteraction<SetupEnvironment>::NuclearInteraction(
process::sibyll::Interaction& hadint, SetupEnvironment const& env)
NuclearInteraction<setup::Environment>::NuclearInteraction(
process::sibyll::Interaction& hadint, setup::Environment const& env)
: environment_(env)
, hadronicInteraction_(hadint) {
......
Processes/Sibyll/testSibyll.cc
+ 10
− 6
View file @ a3ed8076
......@@ -75,12 +75,14 @@ TEST_CASE("Sibyll", "[processes]") {
#include <corsika/particles/ParticleProperties.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupEnvironment.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/process/sibyll/sibyll2.3d.h>
#include <corsika/environment/UniformMediumType.h>
#include <corsika/environment/UniformMagneticField.h>
using namespace corsika::units::si;
using namespace corsika::units;
......@@ -95,17 +97,19 @@ auto sumMomentum(TStackView const& view, geometry::CoordinateSystem const& vCS)
TEST_CASE("SibyllInterface", "[processes]") {
// setup environment, geometry
environment::Environment<environment::IMediumModel> env;
setup::Environment env;
auto& universe = *(env.GetUniverse());
using EnvironmentModel = environment::UniformMediumType<environment::UniformMagneticField<environment::HomogeneousMedium<setup::IEnvironment>>>;
auto theMedium =
environment::Environment<environment::IMediumModel>::CreateNode<geometry::Sphere>(
setup::Environment::CreateNode<geometry::Sphere>(
geometry::Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
1_km * std::numeric_limits<double>::infinity());
using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
theMedium->SetModelProperties<MyHomogeneousModel>(
1_kg / (1_m * 1_m * 1_m),
theMedium->SetModelProperties<EnvironmentModel>(
environment::EMediumType::eAir,
geometry::Vector(env.GetCoordinateSystem(), 0_T, 0_T, 0_T),
1_kg / (1_m * 1_m * 1_m),
environment::NuclearComposition(
std::vector<particles::Code>{particles::Code::Oxygen}, std::vector<float>{1.}));
......
Setup/SetupEnvironment.h
+ 13
− 3
View file @ a3ed8076
......@@ -9,10 +9,20 @@
#pragma once
#include <corsika/environment/Environment.h>
#include <corsika/environment/IMagneticFieldModel.h>
#include <corsika/environment/IMediumModel.h>
#include <corsika/environment/NameModel.h>
#include <corsika/environment/IMediumTypeModel.h>
#include <corsika/environment/IRefractiveIndexModel.h>
#include <corsika/environment/InhomogeneousMedium.h>
namespace corsika::setup {
using IEnvironmentModel = environment::IMediumModel;
using SetupEnvironment = environment::Environment<IEnvironmentModel>;
/**
Definition of the default environemnt model interface. Each model
interface provides properties of the environment in a position
bdependent way.
*/
using IEnvironment = environment::IMediumTypeModel<environment::IMagneticFieldModel<environment::IMediumModel>>;
using Environment = environment::Environment<IEnvironment>;
} // namespace corsika::setup
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

data privacy - impressum