IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 5c9290c7 authored by ralfulrich's avatar ralfulrich Committed by Ralf Ulrich
Browse files

Max comments

parent c3f780c8
No related branches found
No related tags found
1 merge request!265Add medium type, as new property (air, water, rock...)
......@@ -74,37 +74,26 @@ int main(int argc, char** argv) {
EnvType env;
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
Point const center{rootCS, 0_m, 0_m, 0_m};
environment::LayeredSphericalAtmosphereBuilder<setup::EnvironmentInterface> builder{
center};
environment::LayeredSphericalAtmosphereBuilder<setup::EnvironmentInterface, MEnv>
builder{center};
builder.setNuclearComposition(
{{particles::Code::Nitrogen, particles::Code::Oxygen},
{0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
builder.addExponentialLayer<
environment::MediumPropertyModel<environment::UniformMagneticField<
SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km, environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.addExponentialLayer<
environment::MediumPropertyModel<environment::UniformMagneticField<
SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km, environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.addExponentialLayer<
environment::MediumPropertyModel<environment::UniformMagneticField<
SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km, environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.addExponentialLayer<
environment::MediumPropertyModel<environment::UniformMagneticField<
SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km, environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.addLinearLayer<environment::MediumPropertyModel<
environment::UniformMagneticField<HomogeneousMedium<setup::EnvironmentInterface>>>>(
1e9_cm, 112.8_km, environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km,
environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.addExponentialLayer(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km,
environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km,
environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km,
environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.addLinearLayer(1e9_cm, 112.8_km, environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.assemble(env);
// setup particle stack, and add primary particle
......
......@@ -73,7 +73,7 @@ void registerRandomStreams(const int seed) {
}
template <typename T>
using MEnv = environment::UniformMediumType<environment::UniformMagneticField<T>>;
using MEnv = environment::MediumPropertyModel<environment::UniformMagneticField<T>>;
int main(int argc, char** argv) {
......@@ -98,37 +98,26 @@ int main(int argc, char** argv) {
EnvType env;
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
Point const center{rootCS, 0_m, 0_m, 0_m};
environment::LayeredSphericalAtmosphereBuilder<setup::EnvironmentInterface> builder{
center};
environment::LayeredSphericalAtmosphereBuilder<setup::EnvironmentInterface, MEnv>
builder{center};
builder.setNuclearComposition(
{{particles::Code::Nitrogen, particles::Code::Oxygen},
{0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
builder.addExponentialLayer<
environment::MediumPropertyModel<environment::UniformMagneticField<
SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km, environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.addExponentialLayer<
environment::MediumPropertyModel<environment::UniformMagneticField<
SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km, environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.addExponentialLayer<
environment::MediumPropertyModel<environment::UniformMagneticField<
SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km, environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.addExponentialLayer<
environment::MediumPropertyModel<environment::UniformMagneticField<
SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km, environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.addLinearLayer<environment::MediumPropertyModel<
environment::UniformMagneticField<HomogeneousMedium<setup::EnvironmentInterface>>>>(
1e9_cm, 112.8_km, environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km,
environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.addExponentialLayer(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km,
environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km,
environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km,
environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.addLinearLayer(1e9_cm, 112.8_km, environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.assemble(env);
// setup particle stack, and add primary particle
......
......@@ -76,7 +76,7 @@ void registerRandomStreams(const int seed) {
}
template <typename T>
using MEnv = environment::UniformMediumType<environment::UniformMagneticField<T>>;
using MEnv = environment::MediumPropertyModel<environment::UniformMagneticField<T>>;
int main(int argc, char** argv) {
......@@ -101,37 +101,26 @@ int main(int argc, char** argv) {
EnvType env;
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
Point const center{rootCS, 0_m, 0_m, 0_m};
environment::LayeredSphericalAtmosphereBuilder<setup::EnvironmentInterface> builder{
center};
environment::LayeredSphericalAtmosphereBuilder<setup::EnvironmentInterface, MEnv>
builder{center};
builder.setNuclearComposition(
{{particles::Code::Nitrogen, particles::Code::Oxygen},
{0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
builder.addExponentialLayer<
environment::MediumPropertyModel<environment::UniformMagneticField<
SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km, environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.addExponentialLayer<
environment::MediumPropertyModel<environment::UniformMagneticField<
SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km, environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.addExponentialLayer<
environment::MediumPropertyModel<environment::UniformMagneticField<
SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km, environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.addExponentialLayer<
environment::MediumPropertyModel<environment::UniformMagneticField<
SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km, environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.addLinearLayer<environment::MediumPropertyModel<
environment::UniformMagneticField<HomogeneousMedium<setup::EnvironmentInterface>>>>(
1e9_cm, 112.8_km, environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km,
environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.addExponentialLayer(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km,
environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km,
environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km,
environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.addLinearLayer(1e9_cm, 112.8_km, environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.assemble(env);
// setup particle stack, and add primary particle
......
......@@ -20,6 +20,7 @@
#include <memory>
#include <stack>
#include <type_traits>
namespace corsika::environment {
......@@ -35,7 +36,23 @@ namespace corsika::environment {
*
*/
template <typename TMediumInterface = environment::IMediumModel>
namespace detail {
struct NoExtraModelInner {};
template <typename M>
struct NoExtraModel {};
template <template <typename> typename M>
struct has_extra_models : std::true_type {};
template <>
struct has_extra_models<NoExtraModel> : std::false_type {};
} // namespace detail
template <typename TMediumInterface = environment::IMediumModel,
template <typename> typename TMediumModelExtra = detail::NoExtraModel>
class LayeredSphericalAtmosphereBuilder {
std::unique_ptr<NuclearComposition> composition_;
geometry::Point center_;
......@@ -51,6 +68,12 @@ namespace corsika::environment {
}
}
LayeredSphericalAtmosphereBuilder() = delete;
LayeredSphericalAtmosphereBuilder(const LayeredSphericalAtmosphereBuilder&) = delete;
LayeredSphericalAtmosphereBuilder(const LayeredSphericalAtmosphereBuilder&&) = delete;
LayeredSphericalAtmosphereBuilder& operator=(
const LayeredSphericalAtmosphereBuilder&) = delete;
public:
LayeredSphericalAtmosphereBuilder(
corsika::geometry::Point center,
......@@ -62,9 +85,7 @@ namespace corsika::environment {
composition_ = std::make_unique<NuclearComposition>(composition);
}
template <
typename TMediumModel = environment::SlidingPlanarExponential<TMediumInterface>,
typename... TArgs>
template <typename... TArgs>
void addExponentialLayer(units::si::GrammageType b, units::si::LengthType c,
units::si::LengthType upperBoundary, TArgs&&... args) {
using namespace units::si;
......@@ -79,14 +100,19 @@ namespace corsika::environment {
auto const rho0 = b / c;
std::cout << "rho0 = " << rho0 << ", c = " << c << std::endl;
node->template SetModelProperties<TMediumModel>(args..., center_, rho0, -c,
*composition_, earthRadius_);
if constexpr (detail::has_extra_models<TMediumModelExtra>::value)
node->template SetModelProperties<
TMediumModelExtra<environment::SlidingPlanarExponential<TMediumInterface>>>(
args..., center_, rho0, -c, *composition_, earthRadius_);
else
node->template SetModelProperties<
environment::SlidingPlanarExponential<TMediumInterface>>(
center_, rho0, -c, *composition_, earthRadius_);
layers_.push(std::move(node));
}
template <typename TMediumModel = environment::HomogeneousMedium<TMediumInterface>,
typename... TArgs>
template <typename... TArgs>
void addLinearLayer(units::si::LengthType c, units::si::LengthType upperBoundary,
TArgs&&... args) {
using namespace units::si;
......@@ -103,7 +129,14 @@ namespace corsika::environment {
auto node = std::make_unique<VolumeTreeNode<TMediumInterface>>(
std::make_unique<geometry::Sphere>(center_, radius));
node->template SetModelProperties<TMediumModel>(args..., rho0, *composition_);
if constexpr (detail::has_extra_models<TMediumModelExtra>::value)
node->template SetModelProperties<
TMediumModelExtra<environment::HomogeneousMedium<TMediumInterface>>>(
args..., rho0, *composition_);
else
node->template SetModelProperties<
environment::HomogeneousMedium<TMediumInterface>>(args..., rho0,
*composition_);
layers_.push(std::move(node));
}
......
......@@ -415,6 +415,7 @@ def inc_start():
string = """
// generated by readProperties.py
// MANUAL EDITS ON OWN RISK. THEY WILL BE OVERWRITTEN.
// since this is automatic code, we don't attempt to generate automatic unit testing, too: LCOV_EXCL_START
namespace corsika::environment {
"""
......@@ -440,7 +441,10 @@ def detail_end():
#
#
def inc_end():
string = "\n} // end namespace corsika::environment"
string = """
\n} // end namespace corsika::environment
// since this was automatic code, we didn't attempt to generate automatic unit testing, too: LCOV_EXCL_STOP
"""
return string
......
......@@ -202,7 +202,7 @@ TEST_CASE("InhomogeneousMedium") {
}
TEST_CASE("LayeredSphericalAtmosphereBuilder") {
LayeredSphericalAtmosphereBuilder builder(gOrigin);
LayeredSphericalAtmosphereBuilder builder(gOrigin, units::constants::EarthRadius::Mean);
builder.setNuclearComposition(
{{{particles::Code::Nitrogen, particles::Code::Oxygen}}, {{.6, .4}}});
......@@ -306,17 +306,13 @@ TEST_CASE("LayeredSphericalAtmosphereBuilder w/ magnetic field") {
// create magnetic field vectors
Vector B0(gCS, 0_T, 0_T, 1_T);
Vector B1(gCS, 1_T, 1_T, 0_T);
LayeredSphericalAtmosphereBuilder<ModelInterface> builder(gOrigin);
LayeredSphericalAtmosphereBuilder<ModelInterface, UniformMagneticField> builder{
gOrigin};
builder.setNuclearComposition(
{{{particles::Code::Nitrogen, particles::Code::Oxygen}}, {{.6, .4}}});
builder.addLinearLayer<UniformMagneticField<HomogeneousMedium<ModelInterface>>>(
1_km, 10_km, B0);
builder.addExponentialLayer<
UniformMagneticField<SlidingPlanarExponential<ModelInterface>>>(
1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 20_km, B1);
builder.addLinearLayer(1_km, 10_km, B0);
builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 20_km, B0);
CHECK(builder.size() == 2);
......@@ -334,7 +330,7 @@ TEST_CASE("LayeredSphericalAtmosphereBuilder w/ magnetic field") {
.GetMagneticField(pTest)
.GetComponents(gCS));
const Point pTest2(gCS, 10_m, -4_m, R + 15_km);
CHECK(B1.GetComponents(gCS) == univ->GetContainingNode(pTest2)
CHECK(B0.GetComponents(gCS) == univ->GetContainingNode(pTest2)
->GetModelProperties()
.GetMagneticField(pTest2)
.GetComponents(gCS));
......@@ -409,6 +405,7 @@ TEST_CASE("MediumProperties") {
CHECK(air.Cbar() == 10.5961);
CHECK(air.x0() == 1.7418);
CHECK(air.x1() == 4.2759);
CHECK(air.aa() == 0.10914);
CHECK(air.sk() == 3.3994);
CHECK(air.dlt0() == 0.0);
}
......
......@@ -35,8 +35,8 @@ using namespace corsika::environment;
using namespace corsika::geometry;
using namespace corsika::units::si;
// template <typename T>
// using MEnv = environment::UniformMediumType<environment::UniformMagneticField<T>>;
template <typename T>
using MEnv = environment::MediumPropertyModel<environment::UniformMagneticField<T>>;
TEST_CASE("CONEXSourceCut") {
random::RNGManager::GetInstance().RegisterRandomStream("cascade");
......@@ -49,37 +49,27 @@ TEST_CASE("CONEXSourceCut") {
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
Point const center{rootCS, 0_m, 0_m, 0_m};
environment::LayeredSphericalAtmosphereBuilder<setup::EnvironmentInterface> builder{
center, conex::earthRadius};
environment::LayeredSphericalAtmosphereBuilder<setup::EnvironmentInterface, MEnv>
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
builder.addExponentialLayer<
environment::MediumPropertyModel<environment::UniformMagneticField<
SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km, environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.addExponentialLayer<
environment::MediumPropertyModel<environment::UniformMagneticField<
SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km, environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.addExponentialLayer<
environment::MediumPropertyModel<environment::UniformMagneticField<
SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km, environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.addExponentialLayer<
environment::MediumPropertyModel<environment::UniformMagneticField<
SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km, environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.addLinearLayer<environment::MediumPropertyModel<
environment::UniformMagneticField<HomogeneousMedium<setup::EnvironmentInterface>>>>(
1e9_cm, 112.8_km, environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km,
environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.addExponentialLayer(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km,
environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km,
environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km,
environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.addLinearLayer(1e9_cm, 112.8_km, environment::Medium::AirDry1Atm,
geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.assemble(env);
......
......@@ -49,8 +49,9 @@ namespace corsika::process::proposal {
}
media[comp.hash()] =
PROPOSAL::Medium(medium.name(), medium.Ieff(), -medium.Cbar(), medium.aa(), medium.sk(),
medium.x0(), medium.x1(), medium.dlt0(), medium.corrected_density(), comp_vec);
PROPOSAL::Medium(medium.name(), medium.Ieff(), -medium.Cbar(), medium.aa(),
medium.sk(), medium.x0(), medium.x1(), medium.dlt0(),
medium.corrected_density(), comp_vec);
}
});
......
......@@ -81,12 +81,6 @@ TEST_CASE("Pythia", "[processes]") {
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/environment/UniformMagneticField.h>
#include <corsika/environment/UniformMediumType.h>
using namespace corsika;
using namespace corsika::units::si;
......
......@@ -114,12 +114,6 @@ TEST_CASE("QgsjetII", "[processes]") {
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/environment/UniformMagneticField.h>
#include <corsika/environment/UniformMediumType.h>
using namespace corsika::units::si;
using namespace corsika::units;
using namespace corsika;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment