Environment/FlatAtmosphere/FlatAtmosphere.h

-
-/**
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * 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.
- */

Environment/HomogeneousMedium.h

-/**
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * 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.
- */
-
-#ifndef _include_HomogeneousMedium_h_
-#define _include_HomogeneousMedium_h_
-
-#include <corsika/environment/NuclearComposition.h>
-#include <corsika/geometry/Line.h>
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Trajectory.h>
-#include <corsika/particles/ParticleProperties.h>
-#include <corsika/units/PhysicalUnits.h>
-
-/**
- * a homogeneous medium
- */
-
-namespace corsika::environment {
-
-  template <class T>
-  class HomogeneousMedium : public T {
-    corsika::units::si::MassDensityType const fDensity;
-    NuclearComposition const fNuclComp;
-
-  public:
-    HomogeneousMedium(corsika::units::si::MassDensityType pDensity,
-                      NuclearComposition pNuclComp)
-        : fDensity(pDensity)
-        , fNuclComp(pNuclComp){};
-
-    corsika::units::si::MassDensityType GetMassDensity(
-        corsika::geometry::Point const&) const override {
-      return fDensity;
-    }
-    NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
-
-    corsika::units::si::GrammageType IntegratedGrammage(
-        corsika::geometry::Trajectory<corsika::geometry::Line> const& pTraj,
-        corsika::units::si::TimeType pTo) const override {
-      using namespace corsika::units::si;
-      return pTraj.ArcLength(0_s, pTo) * fDensity;
-    }
-
-    corsika::units::si::TimeType FromGrammage(
-        corsika::geometry::Trajectory<corsika::geometry::Line> const& pTraj,
-        corsika::units::si::GrammageType pGrammage) const override {
-      return pTraj.TimeFromArclength(pGrammage / fDensity);
-    }
-  };
-
-} // namespace corsika::environment
-#endif

Environment/IMediumModel.h

-#ifndef _include_IMediumModel_h
-#define _include_IMediumModel_h
-
-#include <corsika/environment/NuclearComposition.h>
-#include <corsika/geometry/Line.h>
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Trajectory.h>
-#include <corsika/units/PhysicalUnits.h>
-
-namespace corsika::environment {
-
-  class IMediumModel {
-  public:
-    virtual ~IMediumModel() = default;
-
-    virtual corsika::units::si::MassDensityType GetMassDensity(
-        corsika::geometry::Point const&) const = 0;
-
-    // todo: think about the mixin inheritance of the trajectory vs the BaseTrajectory
-    // approach for now, only lines are supported
-    virtual corsika::units::si::GrammageType IntegratedGrammage(
-        corsika::geometry::Trajectory<corsika::geometry::Line> const&,
-        corsika::units::si::TimeType) const = 0;
-
-    virtual corsika::units::si::TimeType FromGrammage(
-        corsika::geometry::Trajectory<corsika::geometry::Line> const&,
-        corsika::units::si::GrammageType) const = 0;
-
-    virtual NuclearComposition const& GetNuclearComposition() const = 0;
-  };
-
-} // namespace corsika::environment
-
-#endif

Environment/NuclearComposition.h

-#ifndef _include_NuclearComposition_h
-#define _include_NuclearComposition_h
-
-#include <corsika/particles/ParticleProperties.h>
-#include <numeric>
-#include <stdexcept>
-#include <vector>
-
-namespace corsika::environment {
-  class NuclearComposition {
-    std::vector<float> const fNumberFractions; //!< relative fractions of number density
-    std::vector<corsika::particles::Code> const
-        fComponents; //!< particle codes of consitutents
-
-  public:
-    NuclearComposition(std::vector<corsika::particles::Code> pComponents,
-                       std::vector<float> pFractions)
-        : fNumberFractions(pFractions)
-        , fComponents(pComponents) {
-      auto const sumFractions =
-          std::accumulate(pFractions.cbegin(), pFractions.cend(), 0.f);
-
-      if (!(0.999f < sumFractions && sumFractions < 1.001f)) {
-        throw std::runtime_error("element fractions do not add up to 1");
-      }
-    }
-
-    auto size() const { return fNumberFractions.size(); }
-
-    auto const& GetFractions() const { return fNumberFractions; }
-    auto const& GetComponents() const { return fComponents; }
-  };
-
-} // namespace corsika::environment
-
-#endif

Environment/VolumeTreeNode.h

-/**
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * 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.
- */
-
-#ifndef _include_VolumeTreeNode_H
-#define _include_VolumeTreeNode_H
-
-#include <corsika/geometry/Volume.h>
-#include <memory>
-#include <vector>
-
-namespace corsika::environment {
-
-  class Empty {}; //<! intended for usage as default template argument
-
-  template <typename IModelProperties = Empty>
-  class VolumeTreeNode {
-  public:
-    using VTNUPtr = std::unique_ptr<VolumeTreeNode<IModelProperties>>;
-    using IMPSharedPtr = std::shared_ptr<IModelProperties>;
-    using VolUPtr = std::unique_ptr<corsika::geometry::Volume>;
-
-    VolumeTreeNode(VolUPtr pVolume = nullptr)
-        : fGeoVolume(std::move(pVolume)) {}
-
-    bool Contains(corsika::geometry::Point const& p) const {
-      return fGeoVolume->Contains(p);
-    }
-
-    VolumeTreeNode<IModelProperties> const* Excludes(
-        corsika::geometry::Point const& p) const {
-      auto exclContainsIter =
-          std::find_if(fExcludedNodes.cbegin(), fExcludedNodes.cend(),
-                       [&](auto const& s) { return bool(s->Contains(p)); });
-
-      return exclContainsIter != fExcludedNodes.cend() ? *exclContainsIter : nullptr;
-    }
-
-    /** returns a pointer to the sub-VolumeTreeNode which is "responsible" for the given
-     * \class Point \arg p, or nullptr iff \arg p is not contained in this volume.
-     */
-    VolumeTreeNode<IModelProperties> const* GetContainingNode(
-        corsika::geometry::Point const& p) const {
-      if (!Contains(p)) { return nullptr; }
-
-      if (auto const childContainsIter =
-              std::find_if(fChildNodes.cbegin(), fChildNodes.cend(),
-                           [&](auto const& s) { return bool(s->Contains(p)); });
-          childContainsIter == fChildNodes.cend()) // not contained in any of the children
-      {
-        if (auto const exclContainsIter = Excludes(p)) // contained in any excluded nodes
-        {
-          return exclContainsIter->GetContainingNode(p);
-        } else {
-          return this;
-        }
-      } else {
-        return (*childContainsIter)->GetContainingNode(p);
-      }
-    }
-
-    void AddChild(VTNUPtr pChild) {
-      pChild->fParentNode = this;
-      fChildNodes.push_back(std::move(pChild));
-      // It is a bad idea to return an iterator to the inserted element
-      // because it might get invalidated when the vector needs to grow
-      // later and the caller won't notice.
-    }
-
-    void ExcludeOverlapWith(VTNUPtr const& pNode) {
-      fExcludedNodes.push_back(pNode.get());
-    }
-
-    auto* GetParent() const { return fParentNode; };
-
-    auto const& GetChildNodes() const { return fChildNodes; }
-
-    auto const& GetExcludedNodes() const { return fExcludedNodes; }
-
-    auto const& GetVolume() const { return *fGeoVolume; }
-
-    auto const& GetModelProperties() const { return *fModelProperties; }
-
-    template <typename TModelProperties, typename... Args>
-    auto SetModelProperties(Args&&... args) {
-      static_assert(std::is_base_of_v<IModelProperties, TModelProperties>,
-                    "unusable type provided");
-
-      fModelProperties = std::make_shared<TModelProperties>(std::forward<Args>(args)...);
-      return fModelProperties;
-    }
-
-    void SetModelProperties(IMPSharedPtr ptr) { fModelProperties = ptr; }
-
-    template <class MediumType, typename... Args>
-    static auto CreateMedium(Args&&... args) {
-      static_assert(std::is_base_of_v<IMediumModel, MediumType>,
-                    "unusable type provided, needs to be derived from \"IMediumModel\"");
-
-      return std::make_shared<MediumType>(std::forward<Args>(args)...);
-    }
-
-  private:
-    std::vector<VTNUPtr> fChildNodes;
-    std::vector<VolumeTreeNode<IModelProperties> const*> fExcludedNodes;
-    VolumeTreeNode<IModelProperties> const* fParentNode = nullptr;
-    VolUPtr fGeoVolume;
-    IMPSharedPtr fModelProperties;
-  };
-
-} // namespace corsika::environment
-
-#endif

Environment/testEnvironment.cc

-/**
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * 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.
- */
-
-#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
-                          // cpp file
-
-#include <corsika/environment/HomogeneousMedium.h>
-#include <corsika/environment/IMediumModel.h>
-#include <corsika/environment/NuclearComposition.h>
-#include <corsika/environment/VolumeTreeNode.h>
-#include <corsika/particles/ParticleProperties.h>
-#include <catch2/catch.hpp>
-
-using namespace corsika::geometry;
-using namespace corsika::environment;
-using namespace corsika::units::si;
-
-TEST_CASE("HomogeneousMedium") {
-  NuclearComposition const protonComposition(
-      std::vector<corsika::particles::Code>{corsika::particles::Code::Proton},
-      std::vector<float>{1.f});
-  HomogeneousMedium<IMediumModel> const medium(19.2_g / cube(1_cm), protonComposition);
-}

Framework/CMakeLists.txt

-
-add_subdirectory (Utilities)
-add_subdirectory (Units)
-add_subdirectory (Geometry)
-add_subdirectory (Particles)
-add_subdirectory (Logging)
-add_subdirectory (StackInterface)
-add_subdirectory (ProcessSequence)
-add_subdirectory (Cascade)
-add_subdirectory (Random)

Framework/Cascade/CMakeLists.txt

-
-# namespace of library -> location of header files
-set (
-  CORSIKAcascade_NAMESPACE
-  corsika/cascade
-  )
-
-# header files of this library
-set (
-  CORSIKAcascade_HEADERS
-  Cascade.h
-  )
-
-add_library (CORSIKAcascade INTERFACE)
-
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAcascade ${CORSIKAcascade_NAMESPACE} ${CORSIKAcascade_HEADERS})
-
-#target_link_libraries (
-#  CORSIKAcascade
-#  CORSIKAparticles
-#  CORSIKAunits
-#  CORSIKAthirdparty # for catch2
-#  )
-
-# include directive for upstream code
-target_include_directories (
-  CORSIKAcascade
-  INTERFACE
-  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-  $<INSTALL_INTERFACE:include/>
-  )
-
-# install library
-install (
-  FILES ${CORSIKAcascade_HEADERS}
-  DESTINATION include/${CORSIKAcascade_NAMESPACE}
-  )
-
-# ----------------
-# code unit testing
-add_executable (
-  testCascade
-  testCascade.cc
-  )
-
-target_link_libraries (
-  testCascade
-  #  CORSIKAutls
-  CORSIKArandom
-  ProcessSibyll
-  CORSIKAcascade
-  ProcessStackInspector
-  CORSIKAstackinterface
-  CORSIKAprocesses
-  CORSIKAparticles
-  CORSIKAgeometry
-  CORSIKAenvironment
-  CORSIKAprocesssequence
-  CORSIKAunits
-  CORSIKAthirdparty # for catch2
-  )
-
-add_test (
-  NAME testCascade
-  COMMAND testCascade
-  )

Framework/Cascade/Cascade.h

-
-/**
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * 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.
- */
-
-#ifndef _include_corsika_cascade_Cascade_h_
-#define _include_corsika_cascade_Cascade_h_
-
-#include <corsika/process/ProcessReturn.h>
-#include <corsika/random/RNGManager.h>
-#include <corsika/setup/SetupTrajectory.h>
-#include <corsika/units/PhysicalUnits.h>
-
-#include <type_traits>
-
-using namespace corsika;
-using namespace corsika::units::si;
-
-namespace corsika::cascade {
-
-  template <typename Tracking, typename ProcessList, typename Stack>
-  class Cascade {
-
-    typedef typename Stack::ParticleType Particle;
-
-    Cascade() = delete;
-
-  public:
-    Cascade(Tracking& tr, ProcessList& pl, Stack& stack)
-        : fTracking(tr)
-        , fProcesseList(pl)
-        , fStack(stack) {}
-
-    void Init() {
-      fTracking.Init();
-      fProcesseList.Init();
-      fStack.Init();
-    }
-
-    void Run() {
-      while (!fStack.IsEmpty()) {
-        while (!fStack.IsEmpty()) {
-          Particle& pNext = *fStack.GetNextParticle();
-          Step(pNext);
-        }
-        // do cascade equations, which can put new particles on Stack,
-        // thus, the double loop
-        // DoCascadeEquations();
-      }
-    }
-
-    void Step(Particle& particle) {
-
-      // get access to random number generator
-      static corsika::random::RNG& rmng =
-          corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
-
-      // determine geometric tracking
-      corsika::setup::Trajectory step = fTracking.GetTrack(particle);
-
-      // determine combined total interaction length (inverse)
-      const double total_inv_lambda =
-          fProcesseList.GetTotalInverseInteractionLength(particle, step);
-
-      // sample random exponential step length
-      const double sample_step = rmng() / (double)rmng.max();
-      const double next_interact = -log(sample_step) / total_inv_lambda;
-      std::cout << "total_inv_lambda=" << total_inv_lambda
-                << ", next_interact=" << next_interact << std::endl;
-
-      // determine the maximum geometric step length
-      const double distance_max = fProcesseList.MaxStepLength(particle, step);
-      std::cout << "distance_max=" << distance_max << std::endl;
-
-      // determine combined total inverse decay time
-      const double total_inv_lifetime = fProcesseList.GetTotalInverseLifetime(particle);
-
-      // sample random exponential decay time
-      const double sample_decay = rmng() / (double)rmng.max();
-      const double next_decay = -log(sample_decay) / total_inv_lifetime;
-      std::cout << "total_inv_lifetime=" << total_inv_lifetime
-                << ", next_decay=" << next_decay << std::endl;
-
-      // convert next_step from grammage to length [m]
-      // Environment::GetDistance(step, next_step);
-      const double distance_interact = next_interact;
-      // ....
-
-      // convert next_decay from time to length [m]
-      // Environment::GetDistance(step, next_decay);
-      const double distance_decay = next_decay;
-      // ....
-
-      // take minimum of geometry, interaction, decay for next step
-      const double distance_decay_interact = std::min(next_decay, next_interact);
-      const double distance_next = std::min(distance_decay_interact, distance_max);
-
-      /// here the particle is actually moved along the trajectory to new position:
-      // std::visit(corsika::setup::ParticleUpdate<Particle>{particle}, step);
-      particle.SetPosition(step.GetPosition(1));
-      // .... also update time, momentum, direction, ...
-
-      // apply all continuous processes on particle + track
-      corsika::process::EProcessReturn status =
-          fProcesseList.DoContinuous(particle, step, fStack);
-
-      if (status == corsika::process::EProcessReturn::eParticleAbsorbed) {
-        // fStack.Delete(particle); // TODO: check if this is really needed
-      } else {
-
-        std::cout << "select process " << (distance_decay_interact < distance_max)
-                  << std::endl;
-        // check if geometry limits step, then quit this step
-        if (distance_decay_interact < distance_max) {
-          // check weather decay or interaction limits this step
-          if (distance_decay > distance_interact) {
-            std::cout << "collide" << std::endl;
-            const double actual_inv_length =
-                fProcesseList.GetTotalInverseInteractionLength(particle, step);
-            const double sample_process = rmng() / (double)rmng.max();
-            double inv_lambda_count = 0;
-            fProcesseList.SelectInteraction(particle, fStack, actual_inv_length,
-                                            sample_process, inv_lambda_count);
-          } else {
-            std::cout << "decay" << std::endl;
-            const double actual_decay_time =
-                fProcesseList.GetTotalInverseLifetime(particle);
-            const double sample_process = rmng() / (double)rmng.max();
-            double inv_decay_count = 0;
-            fProcesseList.SelectDecay(particle, fStack, actual_decay_time, sample_process,
-                                      inv_decay_count);
-          }
-        }
-      }
-    }
-
-  private:
-    Tracking& fTracking;
-    ProcessList& fProcesseList;
-    Stack& fStack;
-  };
-
-} // namespace corsika::cascade
-
-#endif

Framework/Cascade/Step.cc

-
-/**
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * 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.
- */
-
-namespace cascade;
-
-void Cascade::Step(auto& sequence, Particle& particle) {
-  double nextStep = sequence.MinStepLength(particle);
-  Trajectory trajectory = sequence.Transport(particle, nextStep);
-  sequence.DoContinuous(particle, trajectory);
-  sequence.DoDiscrete(particle);
-}

Framework/Cascade/testCascade.cc

-/**
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * 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 <limits>
-
-#include <corsika/environment/Environment.h>
-
-#include <corsika/cascade/Cascade.h>
-
-#include <corsika/process/ProcessSequence.h>
-#include <corsika/process/stack_inspector/StackInspector.h>
-#include <corsika/process/tracking_line/TrackingLine.h>
-
-#include <corsika/stack/super_stupid/SuperStupidStack.h>
-
-#include <corsika/particles/ParticleProperties.h>
-
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/RootCoordinateSystem.h>
-#include <corsika/geometry/Vector.h>
-
-#include <corsika/setup/SetupStack.h>
-#include <corsika/setup/SetupTrajectory.h>
-using corsika::setup::Trajectory;
-
-#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
-                          // cpp file
-#include <catch2/catch.hpp>
-
-using namespace corsika;
-using namespace corsika::process;
-using namespace corsika::units;
-using namespace corsika::geometry;
-
-#include <iostream>
-using namespace std;
-using namespace corsika::units::si;
-
-static int fCount = 0;
-
-class ProcessSplit : public corsika::process::ContinuousProcess<ProcessSplit> {
-public:
-  ProcessSplit() {}
-
-  template <typename Particle, typename T>
-  double MaxStepLength(Particle&, T&) const {
-    return 1;
-  }
-
-  template <typename Particle, typename T, typename Stack>
-  void DoContinuous(Particle& p, T&, Stack& s) const {
-    EnergyType E = p.GetEnergy();
-    if (E < 85_MeV) {
-      p.Delete();
-      fCount++;
-    } else {
-      p.SetEnergy(E / 2);
-      auto pnew = s.NewParticle();
-      // s.Copy(p, pnew); fix that .... todo
-      pnew.SetPID(p.GetPID());
-      pnew.SetTime(p.GetTime());
-      pnew.SetEnergy(E / 2);
-      pnew.SetPosition(p.GetPosition());
-      pnew.SetMomentum(p.GetMomentum());
-    }
-  }
-
-  void Init() { fCount = 0; }
-
-  int GetCount() const { return fCount; }
-
-private:
-};
-
-TEST_CASE("Cascade", "[Cascade]") {
-  corsika::random::RNGManager& rmng = corsika::random::RNGManager::GetInstance();
-  rmng.RegisterRandomStream("s_rndm");
-
-  corsika::environment::Environment env; // dummy environment
-  auto& universe = *(env.GetUniverse());
-  auto const radius = 1_m * std::numeric_limits<double>::infinity();
-  ;
-  auto theMedium = corsika::environment::Environment::CreateNode<Sphere>(
-      Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, radius);
-  universe.AddChild(std::move(theMedium));
-
-  tracking_line::TrackingLine<setup::Stack> tracking(env);
-
-  stack_inspector::StackInspector<setup::Stack> p0(true);
-  ProcessSplit p1;
-  const auto sequence = p0 + p1;
-  setup::Stack stack;
-
-  corsika::cascade::Cascade EAS(tracking, sequence, stack);
-  CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
-
-  stack.Clear();
-  auto particle = stack.NewParticle();
-  EnergyType E0 = 100_GeV;
-  particle.SetPID(particles::Code::Electron);
-  particle.SetEnergy(E0);
-  particle.SetPosition(Point(rootCS, {0_m, 0_m, 10_km}));
-  particle.SetMomentum(corsika::stack::super_stupid::MomentumVector(
-      rootCS, {0 * newton * second, 0 * newton * second, -1 * newton * second}));
-  particle.SetTime(0_ns);
-  EAS.Init();
-  EAS.Run();
-
-  /*
-  SECTION("sectionTwo") {
-    for (int i = 0; i < 0; ++i) {
-      stack.Clear();
-      auto particle = stack.NewParticle();
-      EnergyType E0 = 100_GeV * pow(10, i);
-      particle.SetEnergy(E0);
-      EAS.Init();
-      EAS.Run();
-
-      // cout << "Result: E0=" << E0 / 1_GeV << "GeV, count=" << p1.GetCount() << endl;
-    }
-  }
-  */
-}

Framework/Geometry/BaseTrajectory.h

-
-/**
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * 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.
- */
-
-#ifndef _include_BASETRAJECTORY_H
-#define _include_BASETRAJECTORY_H
-
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Vector.h>
-#include <corsika/units/PhysicalUnits.h>
-#include <string>
-
-namespace corsika::geometry {
-
-  /*!
-   * Interface / base class for trajectories.
-   */
-  class BaseTrajectory {
-
-    BaseTrajectory() = delete;
-
-  public:
-    BaseTrajectory(corsika::units::si::TimeType start, corsika::units::si::TimeType end)
-        : fTStart(start)
-        , fTEnd(end) {}
-
-    //!< for \f$ t = 0 \f$, the starting Point shall be returned.
-    virtual Point GetPosition(corsika::units::si::TimeType) const = 0;
-
-    //!< the Point is return from u=0 (start) to u=1 (end)
-    virtual Point GetPosition(double u) const = 0;
-
-    /*!
-     * returns the length between two points of the trajectory
-     * parameterized by \arg t1 and \arg t2. Requires \arg t2 > \arg t1.
-     */
-
-    virtual corsika::units::si::TimeType TimeFromArclength(
-        corsika::units::si::LengthType) const = 0;
-
-    virtual LengthType ArcLength(corsika::units::si::TimeType t1,
-                                 corsika::units::si::TimeType t2) const = 0;
-
-    virtual corsika::units::si::TimeType GetDuration(
-        corsika::units::si::TimeType t1, corsika::units::si::TimeType t2) const {
-      return t2 - t1;
-    }
-
-    virtual Point GetEndpoint() const { return GetPosition(fTEnd); }
-    virtual Point GetStartpoint() const { return GetPosition(fTStart); }
-
-  protected:
-    corsika::units::si::TimeType const fTStart, fTEnd;
-  };
-
-} // namespace corsika::geometry
-
-#endif

Framework/Geometry/BaseVector.h

-
-/**
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * 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.
- */
-
-#ifndef _include_BASEVECTOR_H_
-#define _include_BASEVECTOR_H_
-
-#include <corsika/geometry/CoordinateSystem.h>
-#include <corsika/geometry/QuantityVector.h>
-
-namespace corsika::geometry {
-
-  /*!
-   * Common base class for Vector and Point. Currently it does basically nothing.
-   */
-
-  template <typename dim>
-  class BaseVector {
-  protected:
-    QuantityVector<dim> qVector;
-    CoordinateSystem const* cs;
-
-  public:
-    BaseVector(CoordinateSystem const& pCS, QuantityVector<dim> pQVector)
-        : qVector(pQVector)
-        , cs(&pCS) {}
-  };
-
-} // namespace corsika::geometry
-
-#endif

Framework/Geometry/CMakeLists.txt

-
-set (
-  GEOMETRY_SOURCES
-  CoordinateSystem.cc
-  )
-
-set (
-  GEOMETRY_HEADERS
-  Vector.h
-  Point.h
-  Line.h
-  Sphere.h
-  Volume.h
-  CoordinateSystem.h
-  RootCoordinateSystem.h
-  Helix.h
-  BaseVector.h
-  QuantityVector.h
-  Trajectory.h
-  # BaseTrajectory.h
-  )
-
-set (
-  GEOMETRY_NAMESPACE
-  corsika/geometry
-  )
-
-add_library (CORSIKAgeometry STATIC ${GEOMETRY_SOURCES})
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAgeometry ${GEOMETRY_NAMESPACE} ${GEOMETRY_HEADERS})
-
-set_target_properties (
-  CORSIKAgeometry
-  PROPERTIES
-  VERSION ${PROJECT_VERSION}
-  SOVERSION 1
-  PUBLIC_HEADER "${GEOMETRY_HEADERS}"
-  )
-
-# target dependencies on other libraries (also the header onlys)
-target_link_libraries (
-  CORSIKAgeometry
-  CORSIKAunits
-  )
-
-target_include_directories (
-  CORSIKAgeometry
-  PUBLIC    ${EIGEN3_INCLUDE_DIR}
-  INTERFACE ${EIGEN3_INCLUDE_DIR}
-  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-  $<INSTALL_INTERFACE:include/include>
-  )
-
-install (
-  TARGETS CORSIKAgeometry
-  LIBRARY DESTINATION lib
-  ARCHIVE DESTINATION lib
-  PUBLIC_HEADER DESTINATION include/${GEOMETRY_NAMESPACE}
-  )
-
-
-# --------------------
-# code unit testing
-add_executable (testGeometry testGeometry.cc)
-
-target_link_libraries (
-  testGeometry
-  CORSIKAgeometry
-  CORSIKAunits
-  CORSIKAthirdparty # for catch2
-  )
-
-add_test (NAME testGeometry COMMAND testGeometry)
-

Framework/Geometry/CoordinateSystem.cc

-
-/**
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * 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/geometry/CoordinateSystem.h>
-#include <stdexcept>
-
-using namespace corsika::geometry;
-
-/**
- * returns the transformation matrix necessary to transform primitives with coordinates
- * in \a pFrom to \a pTo, e.g.
- * \f$ \vec{v}^{\text{(to)}} = \mathcal{M} \vec{v}^{\text{(from)}} \f$
- * (\f$ \vec{v}^{(.)} \f$ denotes the coordinates/components of the component in
- * the indicated CoordinateSystem).
- */
-EigenTransform CoordinateSystem::GetTransformation(CoordinateSystem const& pFrom,
-                                                   CoordinateSystem const& pTo) {
-  CoordinateSystem const* a{&pFrom};
-  CoordinateSystem const* b{&pTo};
-  CoordinateSystem const* commonBase{nullptr};
-
-  while (a != b && b != nullptr) {
-    a = &pFrom;
-
-    while (a != b && a != nullptr) { a = a->GetReference(); }
-
-    if (a == b) break;
-
-    b = b->GetReference();
-  }
-
-  if (a == b && a != nullptr) {
-    commonBase = a;
-
-  } else {
-    throw std::runtime_error("no connection between coordinate systems found!");
-  }
-
-  EigenTransform t = EigenTransform::Identity();
-  auto* p = &pFrom;
-
-  while (p != commonBase) {
-    t = p->GetTransform() * t;
-    p = p->GetReference();
-  }
-
-  p = &pTo;
-
-  while (p != commonBase) {
-    t = t * p->GetTransform().inverse(Eigen::TransformTraits::Isometry);
-    p = p->GetReference();
-  }
-
-  return t;
-}

Framework/Geometry/CoordinateSystem.h

-
-/**
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * 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.
- */
-
-#ifndef _include_COORDINATESYSTEM_H_
-#define _include_COORDINATESYSTEM_H_
-
-#include <corsika/geometry/QuantityVector.h>
-#include <corsika/units/PhysicalUnits.h>
-#include <Eigen/Dense>
-#include <stdexcept>
-
-typedef Eigen::Transform<double, 3, Eigen::Affine> EigenTransform;
-typedef Eigen::Translation<double, 3> EigenTranslation;
-
-namespace corsika::geometry {
-
-  class RootCoordinateSystem;
-
-  using corsika::units::si::length_d;
-
-  class CoordinateSystem {
-    CoordinateSystem const* reference = nullptr;
-    EigenTransform transf;
-
-    CoordinateSystem(CoordinateSystem const& reference, EigenTransform const& transf)
-        : reference(&reference)
-        , transf(transf) {}
-
-    CoordinateSystem()
-        : // for creating the root CS
-        transf(EigenTransform::Identity()) {}
-
-  protected:
-    static auto CreateCS() { return CoordinateSystem(); }
-    friend corsika::geometry::RootCoordinateSystem; /// this is the only class that can
-                                                    /// creat ONE unique root CS
-
-  public:
-    static EigenTransform GetTransformation(CoordinateSystem const& c1,
-                                            CoordinateSystem const& c2);
-
-    auto& operator=(const CoordinateSystem& pCS) {
-      reference = pCS.reference;
-      transf = pCS.transf;
-      return *this;
-    }
-
-    auto translate(QuantityVector<length_d> vector) const {
-      EigenTransform const translation{EigenTranslation(vector.eVector)};
-
-      return CoordinateSystem(*this, translation);
-    }
-
-    auto rotate(QuantityVector<phys::units::length_d> axis, double angle) const {
-      if (axis.eVector.isZero()) {
-        throw std::runtime_error("null-vector given as axis parameter");
-      }
-
-      EigenTransform const rotation{Eigen::AngleAxisd(angle, axis.eVector.normalized())};
-
-      return CoordinateSystem(*this, rotation);
-    }
-
-    auto translateAndRotate(QuantityVector<phys::units::length_d> translation,
-                            QuantityVector<phys::units::length_d> axis, double angle) {
-      if (axis.eVector.isZero()) {
-        throw std::runtime_error("null-vector given as axis parameter");
-      }
-
-      EigenTransform const transf{Eigen::AngleAxisd(angle, axis.eVector.normalized()) *
-                                  EigenTranslation(translation.eVector)};
-
-      return CoordinateSystem(*this, transf);
-    }
-
-    auto const* GetReference() const { return reference; }
-
-    auto const& GetTransform() const { return transf; }
-  };
-
-} // namespace corsika::geometry
-
-#endif

Framework/Geometry/Helix.h

-
-/**
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * 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.
- */
-
-#ifndef _include_HELIX_H_
-#define _include_HELIX_H_
-
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Vector.h>
-#include <corsika/units/PhysicalUnits.h>
-#include <cmath>
-
-namespace corsika::geometry {
-  /*!
-   * A Helix is defined by the cyclotron frequency \f$ \omega_c \f$, the initial
-   * Point r0 and
-   * the velocity vectors \f$ \vec{v}_{\parallel} \f$ and \f$ \vec{v}_{\perp} \f$
-   * denoting the projections of the initial velocity \f$ \vec{v}_0 \f$ parallel
-   * and perpendicular to the axis \f$ \vec{B} \f$, respectively, i.e.
-   * \f{align*}{
-        \vec{v}_{\parallel} &= \frac{\vec{v}_0 \cdot \vec{B}}{\vec{B}^2} \vec{B} \\
-        \vec{v}_{\perp} &= \vec{v}_0 - \vec{v}_{\parallel}
-     \f}
-   */
-
-  class Helix {
-
-    using VelocityVec = Vector<corsika::units::si::SpeedType::dimension_type>;
-
-    Point const r0;
-    corsika::units::si::FrequencyType const omegaC;
-    VelocityVec const vPar;
-    VelocityVec const vPerp, uPerp;
-
-    corsika::units::si::LengthType const radius;
-
-  public:
-    Helix(Point const& pR0, corsika::units::si::FrequencyType pOmegaC,
-          VelocityVec const& pvPar, VelocityVec const& pvPerp)
-        : r0(pR0)
-        , omegaC(pOmegaC)
-        , vPar(pvPar)
-        , vPerp(pvPerp)
-        , uPerp(vPerp.cross(vPar.normalized()))
-        , radius(pvPar.norm() / abs(pOmegaC)) {}
-
-    Point GetPosition(corsika::units::si::TimeType t) const {
-      return r0 + vPar * t +
-             (vPerp * (cos(omegaC * t) - 1) + uPerp * sin(omegaC * t)) / omegaC;
-    }
-
-    auto GetRadius() const { return radius; }
-
-    corsika::units::si::LengthType ArcLength(corsika::units::si::TimeType t1,
-                                             corsika::units::si::TimeType t2) const {
-      return (vPar + vPerp).norm() * (t2 - t1);
-    }
-
-    corsika::units::si::TimeType TimeFromArclength(
-        corsika::units::si::LengthType t) const {
-      return t / (vPar + vPerp).norm();
-    }
-  };
-
-} // namespace corsika::geometry
-
-#endif

Framework/Geometry/Line.h

-
-/**
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * 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.
- */
-
-#ifndef _include_LINETRAJECTORY_H
-#define _include_LINETRAJECTORY_H
-
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Vector.h>
-#include <corsika/units/PhysicalUnits.h>
-
-namespace corsika::geometry {
-
-  class Line {
-
-    using VelocityVec = Vector<corsika::units::si::SpeedType::dimension_type>;
-
-    Point const r0;
-    VelocityVec const v0;
-
-  public:
-    Line(Point const& pR0, VelocityVec const& pV0)
-        : r0(pR0)
-        , v0(pV0) {}
-
-    Point GetPosition(corsika::units::si::TimeType t) const { return r0 + v0 * t; }
-
-    LengthType ArcLength(corsika::units::si::TimeType t1,
-                         corsika::units::si::TimeType t2) const {
-      return v0.norm() * (t2 - t1);
-    }
-
-    corsika::units::si::TimeType TimeFromArclength(
-        corsika::units::si::LengthType t) const {
-      return t / v0.norm();
-    }
-
-    auto GetR0() const { return r0; }
-    auto GetV0() const { return v0; }
-  };
-
-} // namespace corsika::geometry
-
-#endif

Framework/Geometry/Point.h

-
-/**
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * 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.
- */
-
-#ifndef _include_POINT_H_
-#define _include_POINT_H_
-
-#include <corsika/geometry/BaseVector.h>
-#include <corsika/geometry/QuantityVector.h>
-#include <corsika/geometry/Vector.h>
-#include <corsika/units/PhysicalUnits.h>
-
-namespace corsika::geometry {
-
-  using corsika::units::si::length_d;
-  using corsika::units::si::LengthType;
-
-  /*!
-   * A Point represents a point in position space. It is defined by its
-   * coordinates with respect to some CoordinateSystem.
-   */
-  class Point : public BaseVector<length_d> {
-  public:
-    Point(CoordinateSystem const& pCS, QuantityVector<length_d> pQVector)
-        : BaseVector<length_d>(pCS, pQVector) {}
-
-    Point(CoordinateSystem const& cs, LengthType x, LengthType y, LengthType z)
-        : BaseVector<length_d>(cs, {x, y, z}) {}
-
-    // TODO: this should be private or protected, we don NOT want to expose numbers
-    // without reference to outside:
-    auto GetCoordinates() const { return BaseVector<length_d>::qVector; }
-
-    /// this always returns a QuantityVector as triple
-    auto GetCoordinates(CoordinateSystem const& pCS) const {
-      if (&pCS == BaseVector<length_d>::cs) {
-        return BaseVector<length_d>::qVector;
-      } else {
-        return QuantityVector<length_d>(
-            CoordinateSystem::GetTransformation(*BaseVector<length_d>::cs, pCS) *
-            BaseVector<length_d>::qVector.eVector);
-      }
-    }
-
-    /*!
-     * transforms the Point into another CoordinateSystem by changing its
-     * coordinates interally
-     */
-    void rebase(CoordinateSystem const& pCS) {
-      BaseVector<length_d>::qVector = GetCoordinates(pCS);
-      BaseVector<length_d>::cs = &pCS;
-    }
-
-    Point operator+(Vector<length_d> const& pVec) const {
-      return Point(*BaseVector<length_d>::cs,
-                   GetCoordinates() + pVec.GetComponents(*BaseVector<length_d>::cs));
-    }
-
-    /*!
-     * returns the distance Vector between two points
-     */
-    Vector<length_d> operator-(Point const& pB) const {
-      auto& cs = *BaseVector<length_d>::cs;
-      return Vector<length_d>(cs, GetCoordinates() - pB.GetCoordinates(cs));
-    }
-  };
-
-} // namespace corsika::geometry
-
-#endif

Framework/Geometry/QuantityVector.h

-
-/**
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * 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.
- */
-
-#ifndef _include_QUANTITYVECTOR_H_
-#define _include_QUANTITYVECTOR_H_
-
-#include <corsika/units/PhysicalUnits.h>
-
-#include <Eigen/Dense>
-
-#include <iostream>
-#include <utility>
-
-namespace corsika::geometry {
-
-  /*!
-   * A QuantityVector is a three-component container based on Eigen::Vector3d
-   * with a phys::units::si::dimension. Arithmethic operators are defined that
-   * propagate the dimensions by dimensional analysis.
-   */
-
-  template <typename dim>
-  class QuantityVector {
-  protected:
-    // todo: check if we need to move "quantity" into namespace corsika::units
-    using Quantity = phys::units::quantity<dim, double>; //< the phys::units::quantity
-                                                         // corresponding to the dimension
-
-  public:
-    Eigen::Vector3d eVector; //!< the actual container where the raw numbers are stored
-
-    typedef dim dimension; //!< should be a phys::units::dimension
-
-    QuantityVector(Quantity a, Quantity b, Quantity c)
-        : eVector{a.magnitude(), b.magnitude(), c.magnitude()} {}
-
-    QuantityVector(Eigen::Vector3d pBareVector)
-        : eVector(pBareVector) {}
-
-    auto operator[](size_t index) const {
-      return Quantity(phys::units::detail::magnitude_tag, eVector[index]);
-    }
-
-    auto norm() const {
-      return Quantity(phys::units::detail::magnitude_tag, eVector.norm());
-    }
-
-    auto squaredNorm() const {
-      using QuantitySquared =
-          decltype(std::declval<Quantity>() * std::declval<Quantity>());
-      return QuantitySquared(phys::units::detail::magnitude_tag, eVector.squaredNorm());
-    }
-
-    auto operator+(QuantityVector<dim> const& pQVec) const {
-      return QuantityVector<dim>(eVector + pQVec.eVector);
-    }
-
-    auto operator-(QuantityVector<dim> const& pQVec) const {
-      return QuantityVector<dim>(eVector - pQVec.eVector);
-    }
-
-    template <typename ScalarDim>
-    auto operator*(phys::units::quantity<ScalarDim, double> const p) const {
-      using ResQuantity = phys::units::detail::Product<ScalarDim, dim, double, double>;
-
-      if constexpr (std::is_same<ResQuantity, double>::value) // result dimensionless, not
-                                                              // a "Quantity" anymore
-      {
-        return QuantityVector<phys::units::dimensionless_d>(eVector * p.magnitude());
-      } else {
-        return QuantityVector<typename ResQuantity::dimension_type>(eVector *
-                                                                    p.magnitude());
-      }
-    }
-
-    template <typename ScalarDim>
-    auto operator/(phys::units::quantity<ScalarDim, double> const p) const {
-      return (*this) * (1 / p);
-    }
-
-    auto operator*(double const p) const { return QuantityVector<dim>(eVector * p); }
-
-    auto operator/(double const p) const { return QuantityVector<dim>(eVector / p); }
-
-    auto& operator/=(double const p) {
-      eVector /= p;
-      return *this;
-    }
-
-    auto& operator*=(double const p) {
-      eVector *= p;
-      return *this;
-    }
-
-    auto& operator+=(QuantityVector<dim> const& pQVec) {
-      eVector += pQVec.eVector;
-      return *this;
-    }
-
-    auto& operator-=(QuantityVector<dim> const& pQVec) {
-      eVector -= pQVec.eVector;
-      return *this;
-    }
-
-    auto& operator-() const { return QuantityVector<dim>(-eVector); }
-
-    auto normalized() const { return (*this) * (1 / norm()); }
-
-    auto operator==(QuantityVector<dim> const& p) const { return eVector == p.eVector; }
-  };
-
-} // namespace corsika::geometry
-
-template <typename dim>
-auto& operator<<(std::ostream& os, corsika::geometry::QuantityVector<dim> qv) {
-  using Quantity = phys::units::quantity<dim, double>;
-
-  os << '(' << qv.eVector(0) << ' ' << qv.eVector(1) << ' ' << qv.eVector(2) << ") "
-     << phys::units::to_unit_symbol<dim, double>(
-            Quantity(phys::units::detail::magnitude_tag, 1));
-  return os;
-}
-
-#endif