diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index ddc484ae59537ff64a9bdf058156b18bd0179467..bdad6b8e933657e926a21d1028f809827a3f792a 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -64,3 +64,5 @@ CORSIKA_REGISTER_EXAMPLE (corsika RUN_OPTIONS --energy=100 --zenith=10 --nevent= add_executable (mars mars.cpp) target_link_libraries (mars CORSIKA8::CORSIKA8) +add_executable (environment environment.cpp) +target_link_libraries (environment CORSIKA8::CORSIKA8) diff --git a/examples/environment.cpp b/examples/environment.cpp new file mode 100644 index 0000000000000000000000000000000000000000..406b475ae92426cdfc6140fa358b46b600863c94 --- /dev/null +++ b/examples/environment.cpp @@ -0,0 +1,77 @@ +/* + * (c) Copyright 2022 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/media/Environment.hpp> +#include <corsika/media/HomogeneousMedium.hpp> +#include <corsika/media/IMediumModel.hpp> +#include <corsika/media/MediumProperties.hpp> +#include <corsika/media/MediumPropertyModel.hpp> + +using namespace corsika; + +// Example to show how to construct an environment in CORSIKA + +// MyMediumInterface is constructed via "mixin inheritance" +// Here, we construct our base class with IMediumModel and IMediumPropertyModel: +// the former allow us define a density profile and nuclear composite, +// the later is used to determine energy losses parameters. +using MyMediumInterface = IMediumPropertyModel<IMediumModel>; + +using MyEnvType = Environment<MyMediumInterface>; + +int main() { + // create environment and universe, empty so far + MyEnvType env; + + // create a geometry object: a sphere with radius of 1000Â m + CoordinateSystemPtr const& rootCS = env.getCoordinateSystem(); + Point const center{rootCS, 0_m, 0_m, 0_m}; + LengthType const radius = 1000_m; + auto sphere = std::make_unique<Sphere>(center, radius); + + // create node from geometry object + auto node = std::make_unique<VolumeTreeNode<MyMediumInterface>>(std::move(sphere)); + + // set medium properties to our node, say it is water + NuclearComposition const nucl_comp{{Code::Hydrogen, Code::Oxygen}, {0.11, 0.89}}; + MassDensityType const density = 1_g / (1_cm * 1_cm * 1_cm); + + // create concrete implementation of MyMediumInterface by combining parts that + // implement the corresponding interfaces. HomogeneousMedium implements IMediumModel, + // MediumPropertyModel implements IMediumPropertyModel. + auto const medium = + std::make_shared<MediumPropertyModel<HomogeneousMedium<MyMediumInterface>>>( + Medium::WaterLiquid, density, nucl_comp); + node->setModelProperties(medium); + + // put our node into universe + // node: this has to be down after setting node model properties, since + // std::move will make our previous defined node, which is a unique pointer + // un-referencable in the context + VolumeTreeNode<MyMediumInterface>* const universe = env.getUniverse().get(); + universe->addChild(std::move(node)); + + // example to explore the media properties of the node + for (auto h = 0_m; h < radius; h += 100_m) { + Point const ptest{rootCS, 0_m, 0_m, h}; + MyMediumInterface const& media_prop = + env.getUniverse()->getContainingNode(ptest)->getModelProperties(); + MassDensityType const rho = media_prop.getMassDensity(ptest); + NuclearComposition const& nuc_comp = media_prop.getNuclearComposition(); + std::vector<Code> const& nuclei = nuc_comp.getComponents(); + std::vector<double> const& fractions = nuc_comp.getFractions(); + CORSIKA_LOG_INFO( + "radius: {:.2f} m, density: {:.2f} g/cm^3, nuclei: ({:d}, {:d}), fractions: " + "({:.2f}, " + "{:.2f})", + h / 1_m, rho / (1_g / static_pow<3>(1_cm)), get_PDG(nuclei.at(0)), + get_PDG(nuclei.at(1)), fractions.at(0), fractions.at(1)); + } + + return 0; +}