diff --git a/CMakeLists.txt b/CMakeLists.txt
index 11d78541572968bcd3c2c7d338b568473c80addc..62d488fa2713a44964460d939aa6a5742b1840d8 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -65,6 +65,7 @@ find_package (Eigen3 REQUIRED)
 add_subdirectory (ThirdParty)
 #add_subdirectory (Utilities)
 add_subdirectory (Framework)
+add_subdirectory (Environment)
 add_subdirectory (Stack)
 add_subdirectory (Setup)
 add_subdirectory (Processes)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 34b55476d9a4f869bd9a29c51688d0f49e680c37..c3a8c95191effe7928ac470281617c4f34d486dd 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -1,4 +1,3 @@
-
 /**
  * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
  *
@@ -17,22 +16,33 @@
 #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/random/RNGManager.h>
 
 #include <corsika/cascade/SibStack.h>
 #include <corsika/cascade/sibyll2.3c.h>
+#include <corsika/geometry/Sphere.h>
 #include <corsika/process/sibyll/ParticleConversion.h>
 
 #include <corsika/units/PhysicalUnits.h>
+
+#include <iostream>
+#include <limits>
+#include <typeinfo>
+
 using namespace corsika;
 using namespace corsika::process;
 using namespace corsika::units;
 using namespace corsika::particles;
 using namespace corsika::random;
+using namespace corsika::geometry;
+using namespace corsika::environment;
 
-#include <iostream>
-#include <typeinfo>
 using namespace std;
+using namespace corsika::units::si;
 
 static int fCount = 0;
 
@@ -40,9 +50,8 @@ class ProcessSplit : public corsika::process::BaseProcess<ProcessSplit> {
 public:
   ProcessSplit() {}
 
-  template <typename Particle>
-  double MinStepLength(Particle& p, setup::Trajectory&) const {
-
+  template <typename Particle, typename Track>
+  double MinStepLength(Particle& p, Track&) const {
     // beam particles for sibyll : 1, 2, 3 for p, pi, k
     // read from cross section code table
     int kBeam = 1;
@@ -92,9 +101,8 @@ public:
     return next_step;
   }
 
-  template <typename Particle, typename Stack>
-  EProcessReturn DoContinuous(Particle&, Trajectory&, Stack&) const {
-    // corsika::utls::ignore(p);
+  template <typename Particle, typename Track, typename Stack>
+  EProcessReturn DoContinuous(Particle&, Track&, Stack&) const {
     return EProcessReturn::eOk;
   }
 
@@ -186,9 +194,8 @@ public:
     //    rnd_ini_();
 
     corsika::random::RNGManager& rmng = corsika::random::RNGManager::GetInstance();
-    ;
-    const std::string str_name = "s_rndm";
-    rmng.RegisterRandomStream(str_name);
+
+    rmng.RegisterRandomStream("s_rndm");
 
     // //    corsika::random::RNG srng;
     // auto srng = rmng.GetRandomStream("s_rndm");
@@ -236,9 +243,28 @@ double s_rndm_(int&) {
   return rmng() / (double)rmng.max();
 }
 
+class A {};
+class B : public A {};
+
 int main() {
+  corsika::environment::Environment env; // dummy environment
+  auto& universe = *(env.GetUniverse());
+
+  auto theMedium = corsika::environment::Environment::CreateNode<Sphere>(
+      Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
+      1_km * std::numeric_limits<double>::infinity());
+
+  using MyHomogeneousModel =
+      corsika::environment::HomogeneousMedium<corsika::environment::IMediumModel>;
+  theMedium->SetModelProperties<MyHomogeneousModel>(
+      1_g / (1_m * 1_m * 1_m),
+      corsika::environment::NuclearComposition(
+          std::vector<corsika::particles::Code>{corsika::particles::Code::Proton},
+          std::vector<float>{1.}));
+
+  universe.AddChild(std::move(theMedium));
 
-  tracking_line::TrackingLine<setup::Stack> tracking;
+  tracking_line::TrackingLine<setup::Stack> tracking(env);
   stack_inspector::StackInspector<setup::Stack> p0(true);
   ProcessSplit p1;
   const auto sequence = p0 + p1;
@@ -248,7 +274,7 @@ int main() {
 
   stack.Clear();
   auto particle = stack.NewParticle();
-  EnergyType E0 = 100_GeV;
+  EnergyType E0 = 100_TeV;
   particle.SetEnergy(E0);
   particle.SetPID(Code::Proton);
   EAS.Init();
diff --git a/Documentation/Examples/geometry_example.cc b/Documentation/Examples/geometry_example.cc
index abba9d9fc0350752fe211bd092490b17e9da923e..67f43344c55d99fadfca0483e824f862c533db4d 100644
--- a/Documentation/Examples/geometry_example.cc
+++ b/Documentation/Examples/geometry_example.cc
@@ -50,12 +50,12 @@ int main() {
   std::cout << "p2-p1 components in cs3: " << diff.GetComponents(cs3)
             << std::endl; // but not under rotations
   std::cout << "p2-p1 norm^2: " << norm << std::endl;
-  assert(norm == 1 * meter*meter);
-  
+  assert(norm == 1 * meter * meter);
+
   Sphere s(p1, 10_m); // define a sphere around a point with a radius
   std::cout << "p1 inside s:  " << s.Contains(p2) << std::endl;
   assert(s.Contains(p2) == 1);
-  
+
   Sphere s2(p1, 3_um); // another sphere
   std::cout << "p1 inside s2: " << s2.Contains(p2) << std::endl;
   assert(s2.Contains(p2) == 0);
diff --git a/Documentation/Examples/stack_example.cc b/Documentation/Examples/stack_example.cc
index 4e2aae20b9fcd34ee0d1d7acc75c157558d4e027..fbe4cd9da69440c48d345c17a1affc52d6c0f150 100644
--- a/Documentation/Examples/stack_example.cc
+++ b/Documentation/Examples/stack_example.cc
@@ -11,9 +11,9 @@
 
 #include <corsika/particles/ParticleProperties.h>
 #include <corsika/stack/super_stupid/SuperStupidStack.h>
+#include <cassert>
 #include <iomanip>
 #include <iostream>
-#include <cassert>
 
 using namespace corsika::units::si;
 using namespace corsika::stack;
@@ -28,7 +28,7 @@ void fill(corsika::stack::super_stupid::SuperStupidStack& s) {
 }
 
 void read(corsika::stack::super_stupid::SuperStupidStack& s) {
-  assert(s.GetSize() == 11);  // stack has 11 particles
+  assert(s.GetSize() == 11); // stack has 11 particles
 
   EnergyType total_energy;
   int i = 0;
@@ -38,7 +38,7 @@ void read(corsika::stack::super_stupid::SuperStupidStack& s) {
     assert(p.GetPID() == corsika::particles::Code::Electron);
     assert(p.GetEnergy() == 1.5_GeV * (i++));
   }
-  //assert(total_energy == 82.5_GeV);
+  // assert(total_energy == 82.5_GeV);
 }
 
 int main() {
diff --git a/Documentation/Examples/staticsequence_example.cc b/Documentation/Examples/staticsequence_example.cc
index 53ebc73e44704f5d38ef01cd33ca67a9c1fac3d7..522c17114d05182c51237340f1afa63becc9e8d6 100644
--- a/Documentation/Examples/staticsequence_example.cc
+++ b/Documentation/Examples/staticsequence_example.cc
@@ -15,20 +15,23 @@
 
 #include <corsika/process/ProcessSequence.h>
 
-#include <corsika/setup/SetupTrajectory.h> // TODO: try to break this dependence later
-using corsika::setup::Trajectory;
-#include <corsika/units/PhysicalUnits.h> // dito
-using namespace corsika::units::si;
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
+#include <corsika/geometry/Vector.h>
 
-using namespace std;
+using namespace corsika;
+using namespace corsika::units::si;
 using namespace corsika::process;
+using namespace std;
+
+const int nData = 10;
 
 class Process1 : public BaseProcess<Process1> {
 public:
   Process1() {}
   template <typename D, typename T, typename S>
   EProcessReturn DoContinuous(D& d, T&, S&) const {
-    for (int i = 0; i < 10; ++i) d.p[i] += 1;
+    for (int i = 0; i < nData; ++i) d.p[i] += 1;
     return EProcessReturn::eOk;
   }
 };
@@ -39,66 +42,64 @@ public:
 
   template <typename D, typename T, typename S>
   inline EProcessReturn DoContinuous(D& d, T&, S&) const {
-    for (int i=0; i<10; ++i) d.p[i] -= 0.1*i;
+    for (int i = 0; i < nData; ++i) d.p[i] -= 0.1 * i;
     return EProcessReturn::eOk;
   }
 };
 
 class Process3 : public BaseProcess<Process3> {
 public:
-  // Process3(const int v) :fV(v) {}
   Process3() {}
 
   template <typename D, typename T, typename S>
-  inline EProcessReturn DoContinuous(D& /*d*/, T& /*t*/, S& /*s*/) const {
-    // for (int i=0; i<10; ++i) d.p[i] += fV;
+  inline EProcessReturn DoContinuous(D&, T&, S&) const {
     return EProcessReturn::eOk;
   }
-
-private:
-  // int fV;
 };
 
 class Process4 : public BaseProcess<Process4> {
 public:
-  // Process4(const int v) : fV(v) {}
-  Process4() {}
+  Process4(const double v)
+      : fV(v) {}
   template <typename D, typename T, typename S>
-  inline EProcessReturn DoContinuous(D& /*d*/, T& /*t*/, S& /*s*/) const {
-    // for (int i=0; i<10; ++i) d.p[i] /= fV;
+  inline EProcessReturn DoContinuous(D& d, T&, S&) const {
+    for (int i = 0; i < nData; ++i) d.p[i] *= fV;
     return EProcessReturn::eOk;
   }
 
 private:
-  // int fV;
+  double fV;
 };
 
 struct DummyData {
-  double p[10] = {0,0,0,0,0,0,0,0,0,0};
+  double p[nData] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+};
+struct DummyStack {
+  void clear() {}
 };
-struct DummyStack {};
+struct DummyTrajectory {};
 
 void modular() {
 
   Process1 m1;
   Process2 m2;
   Process3 m3;
-  Process4 m4;
+  Process4 m4(0.9);
 
   const auto sequence = m1 + m2 + m3 + m4;
-  
+
   DummyData p;
   DummyStack s;
-  Trajectory t;
+  DummyTrajectory t;
 
   const int n = 1000;
   for (int i = 0; i < n; ++i) { sequence.DoContinuous(p, t, s); }
-  
-  for(int i=0; i<10; ++i) {
-    //cout << p.p[i] << endl;
-    //assert(p.p[i] == n-i*100);
+
+  for (int i = 0; i < nData; ++i) {
+    // cout << p.p[i] << endl;
+    // assert(p.p[i] == n-i*100);
   }
-  
+
   cout << " done (nothing...) " << endl;
 }
 
diff --git a/Environment/CMakeLists.txt b/Environment/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..420fc20f087d0cb8ba3ab884f4ae0db7efc8c469
--- /dev/null
+++ b/Environment/CMakeLists.txt
@@ -0,0 +1,52 @@
+set (
+  ENVIRONMENT_HEADERS
+  VolumeTreeNode.h
+  IMediumModel.h
+  NuclearComposition.h
+  HomogeneousMedium.h
+  Environment.h
+  )
+
+set (
+  ENVIRONMENT_NAMESPACE
+  corsika/environment
+  )
+
+add_library (CORSIKAenvironment INTERFACE)
+CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAenvironment ${ENVIRONMENT_NAMESPACE} ${ENVIRONMENT_HEADERS})
+
+# target dependencies on other libraries (also the header onlys)
+target_link_libraries (
+  CORSIKAenvironment
+  INTERFACE
+  CORSIKAgeometry
+  CORSIKAparticles
+  CORSIKAunits
+  )
+
+target_include_directories (
+  CORSIKAenvironment
+  INTERFACE
+  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
+  $<INSTALL_INTERFACE:include>
+  )
+
+install (
+  TARGETS CORSIKAenvironment
+  LIBRARY DESTINATION lib
+  ARCHIVE DESTINATION lib
+  PUBLIC_HEADER DESTINATION include/${ENVIRONMENT_NAMESPACE}
+  )
+
+# --------------------
+# code unit testing
+add_executable (testEnvironment testEnvironment.cc)
+
+target_link_libraries (
+  testEnvironment
+  CORSIKAenvironment
+  CORSIKAthirdparty # for catch2
+  )
+
+add_test (NAME testGeometry COMMAND testGeometry)
+
diff --git a/Environment/Environment.h b/Environment/Environment.h
new file mode 100644
index 0000000000000000000000000000000000000000..5afe69c57000976ef529559696e6bcdca31da609
--- /dev/null
+++ b/Environment/Environment.h
@@ -0,0 +1,57 @@
+/**
+ * (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_Environment_h
+#define _include_Environment_h
+
+#include <corsika/environment/IMediumModel.h>
+#include <corsika/environment/VolumeTreeNode.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
+#include <corsika/setup/SetupEnvironment.h>
+
+namespace corsika::environment {
+  struct Universe : public corsika::geometry::Volume {
+    bool Contains(corsika::geometry::Point const&) const override { return true; }
+  };
+
+  // template <typename IEnvironmentModel>
+  class Environment {
+  public:
+    Environment()
+        : fUniverse(std::make_unique<VolumeTreeNode<IEnvironmentModel>>(
+              std::make_unique<Universe>())) {}
+
+    using IEnvironmentModel = corsika::setup::IEnvironmentModel;
+
+    auto& GetUniverse() { return fUniverse; }
+    auto const& GetUniverse() const { return fUniverse; }
+
+    auto const& GetCoordinateSystem() const {
+      return corsika::geometry::RootCoordinateSystem::GetInstance().GetRootCS();
+    }
+
+    // factory method for creation of VolumeTreeNodes
+    template <class VolumeType, typename... Args>
+    static auto CreateNode(Args&&... args) {
+      static_assert(std::is_base_of_v<corsika::geometry::Volume, VolumeType>,
+                    "unusable type provided, needs to be derived from "
+                    "\"corsika::geometry::Volume\"");
+
+      return std::make_unique<VolumeTreeNode<IEnvironmentModel>>(
+          std::make_unique<VolumeType>(std::forward<Args>(args)...));
+    }
+
+  private:
+    VolumeTreeNode<IEnvironmentModel>::VTNUPtr fUniverse;
+  };
+
+} // namespace corsika::environment
+
+#endif
diff --git a/Environment/HomogeneousMedium.h b/Environment/HomogeneousMedium.h
new file mode 100644
index 0000000000000000000000000000000000000000..a8fa5f7283cfe0e776383f24915b98748d54349f
--- /dev/null
+++ b/Environment/HomogeneousMedium.h
@@ -0,0 +1,59 @@
+/**
+ * (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
diff --git a/Environment/IMediumModel.h b/Environment/IMediumModel.h
new file mode 100644
index 0000000000000000000000000000000000000000..4a1eedd86562187c310d180a4516ee3425a05674
--- /dev/null
+++ b/Environment/IMediumModel.h
@@ -0,0 +1,34 @@
+#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
diff --git a/Environment/NuclearComposition.h b/Environment/NuclearComposition.h
new file mode 100644
index 0000000000000000000000000000000000000000..46a2bf5e859f10a7f5f13a4e4312d38cb1ae6086
--- /dev/null
+++ b/Environment/NuclearComposition.h
@@ -0,0 +1,36 @@
+#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
diff --git a/Environment/VolumeTreeNode.h b/Environment/VolumeTreeNode.h
new file mode 100644
index 0000000000000000000000000000000000000000..d92d77b398bd180d9d64dc5039f32e50ca9be487
--- /dev/null
+++ b/Environment/VolumeTreeNode.h
@@ -0,0 +1,119 @@
+/**
+ * (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
diff --git a/Environment/testEnvironment.cc b/Environment/testEnvironment.cc
new file mode 100644
index 0000000000000000000000000000000000000000..8a8e9c13b5a899ad983d0e852379e5621c9f1489
--- /dev/null
+++ b/Environment/testEnvironment.cc
@@ -0,0 +1,30 @@
+/**
+ * (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);
+}
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index e46960dc9c3d8e8c33254c2751a601ab72843a18..3fac21b57739e5a343baf9e826e7abd67b79984d 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -19,8 +19,6 @@
 
 #include <corsika/setup/SetupTrajectory.h>
 
-using namespace corsika::units::si;
-
 namespace corsika::cascade {
 
   template <typename Tracking, typename ProcessList, typename Stack>
@@ -34,12 +32,7 @@ namespace corsika::cascade {
     Cascade(Tracking& tr, ProcessList& pl, Stack& stack)
         : fTracking(tr)
         , fProcesseList(pl)
-        , fStack(stack) {
-      // static_assert(std::is_member_function_pointer<decltype(&ProcessList::DoDiscrete)>::value,
-      //"ProcessList has not function DoDiscrete.");
-      // static_assert(std::is_member_function_pointer<decltype(&ProcessList::DoContinuous)>::value,
-      //	    "ProcessList has not function DoContinuous.");
-    }
+        , fStack(stack) {}
 
     void Init() {
       fTracking.Init();
@@ -64,7 +57,8 @@ namespace corsika::cascade {
       fProcesseList.MinStepLength(particle, step);
 
       /// here the particle is actually moved along the trajectory to new position:
-      std::visit(corsika::setup::ParticleUpdate<Particle>{particle}, step);
+      // std::visit(corsika::setup::ParticleUpdate<Particle>{particle}, step);
+      particle.SetPosition(step.GetPosition(1));
 
       corsika::process::EProcessReturn status =
           fProcesseList.DoContinuous(particle, step, fStack);
diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc
index 6d0244555d7c9e66de6121e4f3d2a2e31bb4dda3..6bb0a10ddd39956df11e469739faf4dce3cfef7a 100644
--- a/Framework/Cascade/testCascade.cc
+++ b/Framework/Cascade/testCascade.cc
@@ -1,4 +1,3 @@
-
 /**
  * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
  *
@@ -9,6 +8,10 @@
  * the license.
  */
 
+#include <limits>
+
+#include <corsika/environment/Environment.h>
+
 #include <corsika/cascade/Cascade.h>
 
 #include <corsika/process/ProcessSequence.h>
@@ -18,6 +21,7 @@
 #include <corsika/stack/super_stupid/SuperStupidStack.h>
 
 #include <corsika/geometry/Point.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
 #include <corsika/geometry/Vector.h>
 
 #include <corsika/setup/SetupStack.h>
@@ -35,6 +39,7 @@ using namespace corsika::geometry;
 
 #include <iostream>
 using namespace std;
+using namespace corsika::units::si;
 
 static int fCount = 0;
 
@@ -74,8 +79,15 @@ private:
 };
 
 TEST_CASE("Cascade", "[Cascade]") {
-
-  tracking_line::TrackingLine<setup::Stack> tracking;
+  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;
diff --git a/Framework/Geometry/BaseTrajectory.h b/Framework/Geometry/BaseTrajectory.h
index 8289e0dca67fbdaa2865fd7c7799ccf51d619c35..feedace208c89eefc0ecd02cd066e0dc5eadf1a4 100644
--- a/Framework/Geometry/BaseTrajectory.h
+++ b/Framework/Geometry/BaseTrajectory.h
@@ -42,7 +42,10 @@ namespace corsika::geometry {
      * parameterized by \arg t1 and \arg t2. Requires \arg t2 > \arg t1.
      */
 
-    virtual LengthType GetDistance(corsika::units::si::TimeType 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(
diff --git a/Framework/Geometry/CMakeLists.txt b/Framework/Geometry/CMakeLists.txt
index 2b0c8d3af13a5ba00b5fbe6674febb0b21461f96..3e2a8e1d4298af618835e80a9f1f1751704bbe18 100644
--- a/Framework/Geometry/CMakeLists.txt
+++ b/Framework/Geometry/CMakeLists.txt
@@ -10,6 +10,7 @@ set (
   Point.h
   Line.h
   Sphere.h
+  Volume.h
   CoordinateSystem.h
   RootCoordinateSystem.h
   Helix.h
diff --git a/Framework/Geometry/CoordinateSystem.cc b/Framework/Geometry/CoordinateSystem.cc
index 9cf3a1a64125d4fe53236529bd821c6b2c5e57c5..cdc465c60842fb99d2f73539dd657d63131ef251 100644
--- a/Framework/Geometry/CoordinateSystem.cc
+++ b/Framework/Geometry/CoordinateSystem.cc
@@ -10,6 +10,7 @@
  */
 
 #include <corsika/geometry/CoordinateSystem.h>
+#include <stdexcept>
 
 using namespace corsika::geometry;
 
@@ -40,7 +41,7 @@ EigenTransform CoordinateSystem::GetTransformation(CoordinateSystem const& pFrom
     commonBase = a;
 
   } else {
-    throw std::string("no connection between coordinate systems found!");
+    throw std::runtime_error("no connection between coordinate systems found!");
   }
 
   EigenTransform t = EigenTransform::Identity();
diff --git a/Framework/Geometry/CoordinateSystem.h b/Framework/Geometry/CoordinateSystem.h
index 3a06dfbd8bc65b06a83a58184088d5bec42e8b03..6f0584f5e015aa521a0031919e0b0dda98d8d224 100644
--- a/Framework/Geometry/CoordinateSystem.h
+++ b/Framework/Geometry/CoordinateSystem.h
@@ -15,6 +15,7 @@
 #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;
@@ -60,7 +61,7 @@ namespace corsika::geometry {
 
     auto rotate(QuantityVector<phys::units::length_d> axis, double angle) const {
       if (axis.eVector.isZero()) {
-        throw std::string("null-vector given as axis parameter");
+        throw std::runtime_error("null-vector given as axis parameter");
       }
 
       EigenTransform const rotation{Eigen::AngleAxisd(angle, axis.eVector.normalized())};
@@ -71,7 +72,7 @@ namespace corsika::geometry {
     auto translateAndRotate(QuantityVector<phys::units::length_d> translation,
                             QuantityVector<phys::units::length_d> axis, double angle) {
       if (axis.eVector.isZero()) {
-        throw std::string("null-vector given as axis parameter");
+        throw std::runtime_error("null-vector given as axis parameter");
       }
 
       EigenTransform const transf{Eigen::AngleAxisd(angle, axis.eVector.normalized()) *
diff --git a/Framework/Geometry/Helix.h b/Framework/Geometry/Helix.h
index 8eb215820097fae23eef0e535631384d31db01c4..0fa07e64277f534a815e7e067c08e88e0901d783 100644
--- a/Framework/Geometry/Helix.h
+++ b/Framework/Geometry/Helix.h
@@ -58,10 +58,15 @@ namespace corsika::geometry {
 
     auto GetRadius() const { return radius; }
 
-    LengthType GetDistanceBetween(corsika::units::si::TimeType t1,
-                                  corsika::units::si::TimeType t2) const {
+    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
diff --git a/Framework/Geometry/Line.h b/Framework/Geometry/Line.h
index e093a0055fb82ceaeab2c656227b2c005f15958e..a79ba9ce85167ebdf5cef65622e060e844cc22b4 100644
--- a/Framework/Geometry/Line.h
+++ b/Framework/Geometry/Line.h
@@ -32,11 +32,18 @@ namespace corsika::geometry {
 
     Point GetPosition(corsika::units::si::TimeType t) const { return r0 + v0 * t; }
 
-    LengthType GetDistanceBetween(corsika::units::si::TimeType t1,
-                                  corsika::units::si::TimeType t2) const {
-      // assert(t2 >= t1);
+    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
diff --git a/Framework/Geometry/Point.h b/Framework/Geometry/Point.h
index c9c88d7bf599e3f41d33bd92dfe746706f5e5339..74c4b6f8da208064103c1a23018cc208b58b1659 100644
--- a/Framework/Geometry/Point.h
+++ b/Framework/Geometry/Point.h
@@ -26,27 +26,26 @@ namespace corsika::geometry {
    * A Point represents a point in position space. It is defined by its
    * coordinates with respect to some CoordinateSystem.
    */
-  class Point : public BaseVector<phys::units::length_d> {
+  class Point : public BaseVector<length_d> {
   public:
-    Point(CoordinateSystem const& pCS, QuantityVector<phys::units::length_d> pQVector)
-        : BaseVector<phys::units::length_d>(pCS, pQVector) {}
+    Point(CoordinateSystem const& pCS, QuantityVector<length_d> pQVector)
+        : BaseVector<length_d>(pCS, pQVector) {}
 
     Point(CoordinateSystem const& cs, LengthType x, LengthType y, LengthType z)
-        : BaseVector<phys::units::length_d>(cs, {x, y, 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<phys::units::length_d>::qVector; }
+    auto GetCoordinates() const { return BaseVector<length_d>::qVector; }
 
     /// this always returns a QuantityVector as triple
     auto GetCoordinates(CoordinateSystem const& pCS) const {
-      if (&pCS == BaseVector<phys::units::length_d>::cs) {
-        return BaseVector<phys::units::length_d>::qVector;
+      if (&pCS == BaseVector<length_d>::cs) {
+        return BaseVector<length_d>::qVector;
       } else {
-        return QuantityVector<phys::units::length_d>(
-            CoordinateSystem::GetTransformation(*BaseVector<phys::units::length_d>::cs,
-                                                pCS) *
-            BaseVector<phys::units::length_d>::qVector.eVector);
+        return QuantityVector<length_d>(
+            CoordinateSystem::GetTransformation(*BaseVector<length_d>::cs, pCS) *
+            BaseVector<length_d>::qVector.eVector);
       }
     }
 
@@ -55,22 +54,21 @@ namespace corsika::geometry {
      * coordinates interally
      */
     void rebase(CoordinateSystem const& pCS) {
-      BaseVector<phys::units::length_d>::qVector = GetCoordinates(pCS);
-      BaseVector<phys::units::length_d>::cs = &pCS;
+      BaseVector<length_d>::qVector = GetCoordinates(pCS);
+      BaseVector<length_d>::cs = &pCS;
     }
 
-    Point operator+(Vector<phys::units::length_d> const& pVec) const {
-      return Point(
-          *BaseVector<phys::units::length_d>::cs,
-          GetCoordinates() + pVec.GetComponents(*BaseVector<phys::units::length_d>::cs));
+    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<phys::units::length_d> operator-(Point const& pB) const {
-      auto& cs = *BaseVector<phys::units::length_d>::cs;
-      return Vector<phys::units::length_d>(cs, GetCoordinates() - pB.GetCoordinates(cs));
+    Vector<length_d> operator-(Point const& pB) const {
+      auto& cs = *BaseVector<length_d>::cs;
+      return Vector<length_d>(cs, GetCoordinates() - pB.GetCoordinates(cs));
     }
   };
 
diff --git a/Framework/Geometry/Sphere.h b/Framework/Geometry/Sphere.h
index 6ab7c8b9690a7d92e51a6cbb75fc3a73fc227f8c..3e5458b7e2328d062e61307f4c01589c8ee02f05 100644
--- a/Framework/Geometry/Sphere.h
+++ b/Framework/Geometry/Sphere.h
@@ -13,23 +13,27 @@
 #define _include_SPHERE_H_
 
 #include <corsika/geometry/Point.h>
+#include <corsika/geometry/Volume.h>
 #include <corsika/units/PhysicalUnits.h>
 
 namespace corsika::geometry {
 
-  class Sphere {
-    Point center;
-    LengthType const radius;
+  class Sphere : public Volume {
+    Point const fCenter;
+    LengthType const fRadius;
 
   public:
     Sphere(Point const& pCenter, LengthType const pRadius)
-        : center(pCenter)
-        , radius(pRadius) {}
+        : fCenter(pCenter)
+        , fRadius(pRadius) {}
 
-    //! returns true if the Point \a p is within the sphere
-    auto Contains(Point const& p) const {
-      return radius * radius > (center - p).squaredNorm();
+    //! returns true if the Point p is within the sphere
+    bool Contains(Point const& p) const override {
+      return fRadius * fRadius > (fCenter - p).squaredNorm();
     }
+
+    auto& GetCenter() const { return fCenter; }
+    auto GetRadius() const { return fRadius; }
   };
 
 } // namespace corsika::geometry
diff --git a/Framework/Geometry/Trajectory.h b/Framework/Geometry/Trajectory.h
index 686799756504138e349aacfa428b325166ef1e23..f3cabbac7886225846fd8cae1ca7079237812749 100644
--- a/Framework/Geometry/Trajectory.h
+++ b/Framework/Geometry/Trajectory.h
@@ -13,6 +13,7 @@
 #define _include_TRAJECTORY_H
 
 #include <corsika/units/PhysicalUnits.h>
+#include <corsika/geometry/Point.h>
 
 using corsika::units::si::LengthType;
 using corsika::units::si::TimeType;
@@ -25,7 +26,7 @@ namespace corsika::geometry {
     corsika::units::si::TimeType fTimeLength;
 
   public:
-    using T::GetDistanceBetween;
+    using T::ArcLength;
     using T::GetPosition;
 
     Trajectory(T const& theT, corsika::units::si::TimeType timeLength)
@@ -36,14 +37,14 @@ namespace corsika::geometry {
       return fTraj.GetPosition(t + fTStart);
       }*/
 
-    Point GetPosition(const double u) const { return T::GetPosition(fTimeLength * u); }
+    Point GetPosition(double u) const { return T::GetPosition(fTimeLength * u); }
 
     TimeType GetDuration() const { return fTimeLength; }
 
-    LengthType GetDistance(const corsika::units::si::TimeType t) const {
+    LengthType GetDistance(corsika::units::si::TimeType t) const {
       assert(t > fTimeLength);
       assert(t >= 0 * corsika::units::si::second);
-      return T::DistanceBetween(0, t);
+      return T::ArcLength(0, t);
     }
   };
 
diff --git a/Framework/Geometry/Vector.h b/Framework/Geometry/Vector.h
index cb4e4900544705c10a37ff16a0fd4ba00b0f9b29..94394b90dd859761d4f15d0a2389f173726ee34d 100644
--- a/Framework/Geometry/Vector.h
+++ b/Framework/Geometry/Vector.h
@@ -181,6 +181,17 @@ namespace corsika::geometry {
                                                              bareResult);
       }
     }
+
+    template <typename dim2>
+    auto dot(Vector<dim2> pV) const {
+      auto const c1 = GetComponents().eVector;
+      auto const c2 = pV.GetComponents(*BaseVector<dim>::cs).eVector;
+      auto const bareResult = c1.dot(c2);
+
+      using ProdQuantity = phys::units::detail::Product<dim, dim2, double, double>;
+
+      return ProdQuantity(phys::units::detail::magnitude_tag, bareResult);
+    }
   };
 
 } // namespace corsika::geometry
diff --git a/Framework/Geometry/Volume.h b/Framework/Geometry/Volume.h
new file mode 100644
index 0000000000000000000000000000000000000000..ac892e50b926758d76298675a3c0db80557ab9a4
--- /dev/null
+++ b/Framework/Geometry/Volume.h
@@ -0,0 +1,19 @@
+#ifndef _include_VOLUME_H_
+#define _include_VOLUME_H_
+
+#include <corsika/geometry/Point.h>
+
+namespace corsika::geometry {
+
+  class Volume {
+
+  public:
+    //! returns true if the Point p is within the volume
+    virtual bool Contains(Point const& p) const = 0;
+
+    virtual ~Volume() = default;
+  };
+
+} // namespace corsika::geometry
+
+#endif
diff --git a/Framework/Geometry/testGeometry.cc b/Framework/Geometry/testGeometry.cc
index 9e3be3310bbaf3ed6c9ae22bb6faa67701f20687..8130114aacb6c0bbac07b0f42f3d043d704b6c89 100644
--- a/Framework/Geometry/testGeometry.cc
+++ b/Framework/Geometry/testGeometry.cc
@@ -132,7 +132,15 @@ TEST_CASE("Sphere") {
   Point center(rootCS, {0_m, 3_m, 4_m});
   Sphere sphere(center, 5_m);
 
-  SECTION("isInside") {
+  SECTION("GetCenter") {
+    CHECK((sphere.GetCenter().GetCoordinates(rootCS) -
+           QuantityVector<length_d>(0_m, 3_m, 4_m))
+              .norm()
+              .magnitude() == Approx(0).margin(absMargin));
+    CHECK(sphere.GetRadius() / 5_m == Approx(1));
+  }
+
+  SECTION("Contains") {
     REQUIRE_FALSE(sphere.Contains(Point(rootCS, {100_m, 0_m, 0_m})));
     REQUIRE(sphere.Contains(Point(rootCS, {2_m, 3_m, 4_m})));
   }
@@ -152,11 +160,11 @@ TEST_CASE("Trajectories") {
             .norm()
             .magnitude() == Approx(0).margin(absMargin));
 
-    Trajectory<Line> base(line, 1_s);
-    CHECK(line.GetPosition(2_s).GetCoordinates() ==
-          base.GetPosition(2_s).GetCoordinates());
+    auto const t = 1_s;
+    Trajectory<Line> base(line, t);
+    CHECK(line.GetPosition(t).GetCoordinates() == base.GetPosition(1.).GetCoordinates());
 
-    CHECK(base.GetDistanceBetween(1_s, 2_s) / 1_m == Approx(1));
+    CHECK(base.ArcLength(1_s, 2_s) / 1_m == Approx(1));
   }
 
   SECTION("Helix") {
@@ -180,10 +188,10 @@ TEST_CASE("Trajectories") {
               .norm()
               .magnitude() == Approx(0).margin(absMargin));
 
-    Trajectory<Helix> const base(helix, 1_s);
-    CHECK(helix.GetPosition(1234_s).GetCoordinates() ==
-          base.GetPosition(1234_s).GetCoordinates());
+    auto const t = 1234_s;
+    Trajectory<Helix> const base(helix, t);
+    CHECK(helix.GetPosition(t).GetCoordinates() == base.GetPosition(1.).GetCoordinates());
 
-    CHECK(base.GetDistanceBetween(1_s, 2_s) / 1_m == Approx(5));
+    CHECK(base.ArcLength(1_s, 2_s) / 1_m == Approx(5));
   }
 }
diff --git a/Framework/ProcessSequence/BaseProcess.h b/Framework/ProcessSequence/BaseProcess.h
index 897497d83c82caee4dbc6c19f90996a5ee7f97fa..34dcb198e4d26f409837d253f3270cc0f488098b 100644
--- a/Framework/ProcessSequence/BaseProcess.h
+++ b/Framework/ProcessSequence/BaseProcess.h
@@ -29,8 +29,15 @@ namespace corsika::process {
   struct BaseProcess {
     derived& GetRef() { return static_cast<derived&>(*this); }
     const derived& GetRef() const { return static_cast<const derived&>(*this); }
+
+    template <typename Particle, typename Stack>
+    inline EProcessReturn DoDiscrete(Particle&, Stack&) const; // {}
+
+    template <typename Particle, typename Track, typename Stack>
+    inline EProcessReturn DoContinuous(Particle&, Track&, Stack&) const; // {}
   };
 
+  /*
   template <typename T>
   struct is_base {
     static const bool value = false;
@@ -39,6 +46,7 @@ namespace corsika::process {
   struct is_base<BaseProcess<T>> {
     static const bool value = true;
   };
+  */
 
 } // namespace corsika::process
 
diff --git a/Framework/ProcessSequence/ContinuousProcess.h b/Framework/ProcessSequence/ContinuousProcess.h
index c53365f1afcdfb64ac808580077dd0e284dcdbd2..a9bc21e5db1c6dcbb2ff8e13110ec957e6242e39 100644
--- a/Framework/ProcessSequence/ContinuousProcess.h
+++ b/Framework/ProcessSequence/ContinuousProcess.h
@@ -13,6 +13,7 @@
 #define _include_corsika_continuousprocess_h_
 
 #include <corsika/process/ProcessReturn.h> // for convenience
+//#include <corsika/setup/SetupTrajectory.h>
 
 namespace corsika::process {
 
@@ -29,6 +30,11 @@ namespace corsika::process {
   struct ContinuousProcess {
     derived& GetRef() { return static_cast<derived&>(*this); }
     const derived& GetRef() const { return static_cast<const derived&>(*this); }
+
+    // here starts the interface part
+    // -> enforce derived to implement DoContinuous...
+    template <typename D, typename T, typename S>
+    inline EProcessReturn DoContinuous(D&, T&, S&) const;
   };
 
 } // namespace corsika::process
diff --git a/Framework/ProcessSequence/DiscreteProcess.h b/Framework/ProcessSequence/DiscreteProcess.h
index 83064f5d8767c93619f6742d50c56aeae0863700..a540c92c783dcad931f0d2bf8278c057736bda79 100644
--- a/Framework/ProcessSequence/DiscreteProcess.h
+++ b/Framework/ProcessSequence/DiscreteProcess.h
@@ -13,7 +13,7 @@
 #define _include_corsika_discreteprocess_h_
 
 #include <corsika/process/ProcessReturn.h> // for convenience
-
+#include <corsika/setup/SetupTrajectory.h>
 #include <iostream> // debug
 
 namespace corsika::process {
@@ -30,25 +30,13 @@ namespace corsika::process {
   template <typename derived>
   struct DiscreteProcess {
 
-    // DiscreteProcess() {
-    // static_assert(mustProvide<derived>::mustProvide, "");
-    //}
-
     derived& GetRef() { return static_cast<derived&>(*this); }
     const derived& GetRef() const { return static_cast<const derived&>(*this); }
 
-    // here starts the interface part
+    /// here starts the interface-definition part
     // -> enforce derived to implement DoDiscrete...
     template <typename Particle, typename Stack>
     inline EProcessReturn DoDiscrete(Particle&, Stack&) const; // {}
-
-    //  private:
-    template <typename D, typename T, typename S>
-    inline EProcessReturn DoContinuous(D& d, T&, S&) const {
-      std::cout << "yeah" << std::endl;
-      return EProcessReturn::eOk;
-    } // find out how to make this FINAL
-    // void DoContinuous;
   };
 
 } // namespace corsika::process
diff --git a/Framework/ProcessSequence/ProcessSequence.h b/Framework/ProcessSequence/ProcessSequence.h
index a07e3ef84c91e6e4e54743c4a4e87c09d933bf2f..54ac90e7a18f7c40c1e5a53c82ebf854e5abe486 100644
--- a/Framework/ProcessSequence/ProcessSequence.h
+++ b/Framework/ProcessSequence/ProcessSequence.h
@@ -17,14 +17,12 @@
 #include <corsika/process/DiscreteProcess.h>
 #include <corsika/process/ProcessReturn.h>
 
-#include <corsika/setup/SetupTrajectory.h>
-
-#include <variant>
+//#include <corsika/setup/SetupTrajectory.h>
+// using corsika::setup::Trajectory;
+//#include <variant>
 
 //#include <type_traits> // still needed ?
 
-using corsika::setup::Trajectory;
-
 namespace corsika::process {
 
   /* namespace detail { */
@@ -149,31 +147,33 @@ namespace corsika::process {
     // example for a trait-based call:
     // void Hello() const  { detail::CallHello<T1,T2>::Call(A, B); }
 
-    template <typename Particle, typename Stack>
-    inline EProcessReturn DoContinuous(Particle& p, Trajectory& t, Stack& s) const {
+    template <typename Particle, typename Track, typename Stack>
+    inline EProcessReturn DoContinuous(Particle& p, Track& t, Stack& s) const {
       EProcessReturn ret = EProcessReturn::eOk;
       if constexpr (!std::is_base_of<DiscreteProcess<T1>, T1>::value) {
+        // A.DoContinuous(std::forward<Particle>(p), t, std::forward<Stack>(s));
         A.DoContinuous(p, t, s);
       }
       if constexpr (!std::is_base_of<DiscreteProcess<T2>, T2>::value) {
+        // B.DoContinuous(std::forward<Particle>(p), t, std::forward<Stack>(s));
         B.DoContinuous(p, t, s);
       }
       return ret;
     }
 
-    template <typename Particle>
-    inline void MinStepLength(Particle& p, Trajectory& step) const {
-      A.MinStepLength(p, step);
-      B.MinStepLength(p, step);
+    template <typename Particle, typename Track>
+    inline void MinStepLength(Particle& p, Track& track) const {
+      A.MinStepLength(p, track);
+      B.MinStepLength(p, track);
     }
 
     /*
-    template <typename Particle, typename Trajectory>
-    inline Trajectory Transport(Particle& p, double& length) const {
+    template <typename Particle, typename Track>
+    inline Track Transport(Particle& p, double& length) const {
       A.Transport(p, length); // todo: maybe check (?) if there is more than one Transport
                               // process implemented??
       return B.Transport(
-          p, length); // need to do this also to decide which Trajectory to return!!!!
+          p, length); // need to do this also to decide which Track to return!!!!
     }
     */
 
diff --git a/Framework/ProcessSequence/testProcessSequence.cc b/Framework/ProcessSequence/testProcessSequence.cc
index 9757298cdc2269a90fd875cdedfa68a16f19ab5a..7094cacb94e8da03bcfce1b9190160fb1b8c1898 100644
--- a/Framework/ProcessSequence/testProcessSequence.cc
+++ b/Framework/ProcessSequence/testProcessSequence.cc
@@ -19,76 +19,81 @@
 
 #include <corsika/process/ProcessSequence.h>
 
-#include <corsika/setup/SetupTrajectory.h> // TODO: maybe try to break this dependency later!
-using corsika::setup::Trajectory;
-#include <corsika/units/PhysicalUnits.h>
+using namespace corsika;
 using namespace corsika::units::si;
-
-using namespace std;
 using namespace corsika::process;
+using namespace std;
 
 static const int nData = 10;
 
+int globalCount = 0;
+
 class ContinuousProcess1 : public ContinuousProcess<ContinuousProcess1> {
+  int fV = 0;
+
 public:
-  ContinuousProcess1() {}
-  void Init() { cout << "ContinuousProcess1::Init" << endl; }
-  template <typename D, typename S>
-  inline EProcessReturn DoContinuous(D& d, Trajectory&, S&) const {
+  ContinuousProcess1(const int v)
+      : fV(v) {}
+  void Init() {
+    cout << "ContinuousProcess1::Init" << endl;
+    assert(globalCount == fV);
+    globalCount++;
+  }
+  template <typename D, typename T, typename S>
+  inline EProcessReturn DoContinuous(D& d, T&, S&) const {
     cout << "ContinuousProcess1::DoContinuous" << endl;
     for (int i = 0; i < nData; ++i) d.p[i] += 0.933;
     return EProcessReturn::eOk;
   }
-
-  template <typename Particle, typename Stack>
-  inline EProcessReturn DoDiscrete(Particle&, Stack&) const {
-    cout << "ContinuousProcess1::DoDiscrete" << endl;
-    return EProcessReturn::eOk;
-  }
 };
 
 class ContinuousProcess2 : public ContinuousProcess<ContinuousProcess2> {
+  int fV = 0;
+
 public:
-  ContinuousProcess2() {}
-  void Init() { cout << "ContinuousProcess2::Init" << endl; }
-  template <typename D, typename S>
-  inline EProcessReturn DoContinuous(D& d, Trajectory&, S&) const {
+  ContinuousProcess2(const int v)
+      : fV(v) {}
+  void Init() {
+    cout << "ContinuousProcess2::Init" << endl;
+    assert(globalCount == fV);
+    globalCount++;
+  }
+  template <typename D, typename T, typename S>
+  inline EProcessReturn DoContinuous(D& d, T&, S&) const {
     cout << "ContinuousProcess2::DoContinuous" << endl;
     for (int i = 0; i < nData; ++i) d.p[i] += 0.933;
     return EProcessReturn::eOk;
   }
-
-  template <typename Particle, typename Stack>
-  inline EProcessReturn DoDiscrete(Particle&, Stack&) const {
-    cout << "ContinuousProcess2::DoDiscrete" << endl;
-    return EProcessReturn::eOk;
-  }
 };
 
 class Process1 : public DiscreteProcess<Process1> {
 public:
-  Process1() {}
-  void Init() { cout << "Process1::Init" << endl; }
+  Process1(const int v)
+      : fV(v) {}
+  void Init() {
+    cout << "Process1::Init" << endl;
+    assert(globalCount == fV);
+    globalCount++;
+  }
   template <typename D, typename S>
-  inline EProcessReturn DoContinuous(D& d, Trajectory&, S&) const {
+  inline EProcessReturn DoDiscrete(D& d, S&) const {
     for (int i = 0; i < nData; ++i) d.p[i] += 1 + i;
     return EProcessReturn::eOk;
   }
-  template <typename Particle, typename Stack>
-  inline EProcessReturn DoDiscrete(Particle&, Stack&) const {
-    cout << "Process1::DoDiscrete" << endl;
-    return EProcessReturn::eOk;
-  }
+  // private:
+  int fV;
 };
 
 class Process2 : public DiscreteProcess<Process2> {
+  int fV = 0;
+
 public:
-  Process2() {}
-  void Init() { cout << "Process2::Init" << endl; }
-  template <typename D, typename S>
-  inline EProcessReturn DoContinuous(D& d, Trajectory&, S&) const {
-    for (int i = 0; i < nData; ++i) d.p[i] *= 0.7;
-    return EProcessReturn::eOk;
+  Process2(const int v)
+      : fV(v) {}
+  void Init() {
+    cout << "Process2::Init" << endl;
+    assert(globalCount == fV);
+    globalCount++;
   }
   template <typename Particle, typename Stack>
   inline EProcessReturn DoDiscrete(Particle&, Stack&) const {
@@ -98,13 +103,15 @@ public:
 };
 
 class Process3 : public DiscreteProcess<Process3> {
+  int fV = 0;
+
 public:
-  Process3() {}
-  void Init() { cout << "Process3::Init" << endl; }
-  template <typename D, typename S>
-  inline EProcessReturn DoContinuous(D& d, Trajectory&, S&) const {
-    for (int i = 0; i < nData; ++i) d.p[i] += 0.933;
-    return EProcessReturn::eOk;
+  Process3(const int v)
+      : fV(v) {}
+  void Init() {
+    cout << "Process3::Init" << endl;
+    assert(globalCount == fV);
+    globalCount++;
   }
   template <typename Particle, typename Stack>
   inline EProcessReturn DoDiscrete(Particle&, Stack&) const {
@@ -114,66 +121,84 @@ public:
 };
 
 class Process4 : public BaseProcess<Process4> {
+  int fV = 0;
+
 public:
-  Process4() {}
-  void Init() { cout << "Process4::Init" << endl; }
-  template <typename D, typename S>
-  inline EProcessReturn DoContinuous(D& d, Trajectory&, S&) const {
-    for (int i = 0; i < nData; ++i) d.p[i] /= 1.2;
+  Process4(const int v)
+      : fV(v) {}
+  void Init() {
+    cout << "Process4::Init" << endl;
+    assert(globalCount == fV);
+    globalCount++;
+  }
+  template <typename D, typename T, typename S>
+  inline EProcessReturn DoContinuous(D& d, T&, S&) const {
+    for (int i = 0; i < nData; ++i) { d.p[i] /= 1.2; }
     return EProcessReturn::eOk;
   }
   // inline double MinStepLength(D& d) {
-  // void DoDiscrete(Particle& p, Stack& s) const {
+  template <typename Particle, typename Stack>
+  EProcessReturn DoDiscrete(Particle&, Stack&) const {
+    return EProcessReturn::eOk;
+  }
 };
 
 struct DummyData {
   double p[nData] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
 };
 struct DummyStack {};
+struct DummyTrajectory {};
 
-TEST_CASE("Cascade", "[Cascade]") {
+TEST_CASE("Process Sequence", "[Process Sequence]") {
 
-  SECTION("sectionTwo") {
-
-    Process1 m1;
-    Process2 m2;
-    Process3 m3;
-    Process4 m4;
+  SECTION("Check init order") {
+    Process1 m1(0);
+    Process2 m2(1);
+    Process3 m3(2);
+    Process4 m4(3);
 
     const auto sequence = m1 + m2 + m3 + m4;
 
-    ContinuousProcess1 cp1;
-    ContinuousProcess2 cp2;
+    globalCount = 0;
+    sequence.Init();
+    // REQUIRE_NOTHROW( (sequence.Init()) );
+
+    // const auto sequence_wrong = m3 + m2 + m1 + m4;
+    // globalCount = 0;
+    // sequence_wrong.Init();
+    // REQUIRE_THROWS(sequence_wrong.Init());
+  }
+
+  SECTION("sectionTwo") {
+
+    ContinuousProcess1 cp1(0);
+    ContinuousProcess2 cp2(3);
+    Process2 m2(1);
+    Process3 m3(2);
 
     const auto sequence2 = cp1 + m2 + m3 + cp2;
 
     DummyData p;
     DummyStack s;
+    DummyTrajectory t;
 
-    cout << "-->init" << endl;
+    cout << "-->init sequence2" << endl;
+    globalCount = 0;
     sequence2.Init();
     cout << "-->docont" << endl;
 
-    // auto const root = corsika::geometry::CoordinateSystem::CreateRootCS();
-    // corsika::geometry::Point pos(root, {0_m, 0_m, 0_m});
-    // corsika::geometry::Vector<SpeedType::dimension_type> vec(root,
-    // {1_m/1_s,0_m/1_s,0_m/1_s}); corsika::geometry::Line traj(pos, vec);
-    Trajectory
-        t; //(corsika::geometry::Trajectory<corsika::geometry::Line>(traj, 0_s, 100_ns));
-
     sequence2.DoContinuous(p, t, s);
     cout << "-->dodisc" << endl;
     sequence2.DoDiscrete(p, s);
     cout << "-->done" << endl;
 
-    sequence.Init();
-
     const int nLoop = 5;
     cout << "Running loop with n=" << nLoop << endl;
-    for (int i = 0; i < nLoop; ++i) { sequence.DoContinuous(p, t, s); }
+    for (int i = 0; i < nLoop; ++i) {
+      sequence2.DoContinuous(p, t, s);
+      sequence2.DoDiscrete(p, s);
+    }
     for (int i = 0; i < nData; i++) { cout << "data[" << i << "]=" << p.p[i] << endl; }
     cout << "done" << endl;
   }
-
-  SECTION("sectionThree") {}
 }
diff --git a/Framework/Units/PhysicalUnits.h b/Framework/Units/PhysicalUnits.h
index 8975a6e9994cb316cf4c633c3ba6a0715c1d90c0..7083ce743811176e2c807f2d421dde2273fbb673 100644
--- a/Framework/Units/PhysicalUnits.h
+++ b/Framework/Units/PhysicalUnits.h
@@ -51,6 +51,8 @@ namespace corsika::units::si {
       phys::units::quantity<phys::units::electric_charge_d, double>;
   using EnergyType = phys::units::quantity<phys::units::energy_d, double>;
   using MassType = phys::units::quantity<phys::units::mass_d, double>;
+  using MassDensityType = phys::units::quantity<phys::units::mass_density_d, double>;
+  using GrammageType = phys::units::quantity<phys::units::dimensions<-2, 1, 0>, double>;
   using MomentumType = phys::units::quantity<momentum_d, double>;
   using CrossSectionType = phys::units::quantity<sigma_d, double>;
 
diff --git a/Processes/StackInspector/StackInspector.cc b/Processes/StackInspector/StackInspector.cc
index b0acefc4a86c72acba373feafc82b7434aa2b9a0..1986c83bb4471b1adee8b2b13653cb64c7ef5167 100644
--- a/Processes/StackInspector/StackInspector.cc
+++ b/Processes/StackInspector/StackInspector.cc
@@ -12,6 +12,7 @@
 #include <corsika/geometry/RootCoordinateSystem.h>
 #include <corsika/process/stack_inspector/StackInspector.h>
 #include <corsika/units/PhysicalUnits.h>
+#include <corsika/particles/ParticleProperties.h>
 
 #include <corsika/logging/Logger.h>
 
diff --git a/Processes/TrackingLine/CMakeLists.txt b/Processes/TrackingLine/CMakeLists.txt
index 2911c9b787787f43ba837e64fd67c89921c699b0..83ed893c150d6d041182bbe7ab9cd30046edb640 100644
--- a/Processes/TrackingLine/CMakeLists.txt
+++ b/Processes/TrackingLine/CMakeLists.txt
@@ -22,17 +22,16 @@ target_include_directories (
 
 install (FILES ${MODEL_HEADERS} DESTINATION include/${MODEL_NAMESPACE})
 
-
-
-# # --------------------
-# # code unit testing
-# add_executable (testStackInspector testStackInspector.cc)
-
-# target_link_libraries (
-#   testStackInspector
-#   CORSIKAunits
-#   CORSIKAthirdparty # for catch2
-#   )
-
-# add_test (NAME testStackInspector COMMAND testStackInspector)
-
+# #-- -- -- -- -- -- -- -- -- --
+# #code unit testing
+add_executable (testTrackingLine testTrackingLine.cc)
+
+target_link_libraries (
+   testTrackingLine
+   CORSIKAunits
+   CORSIKAenvironment
+   CORSIKAgeometry
+   CORSIKAthirdparty # for catch2
+)
+
+add_test (NAME testTrackingLine COMMAND testTrackingLine)
diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h
index e750342a1add1bff5e32bd270710e995d60ce4b0..b6b8864087b53fe577029c379996ca9d6ef29d35 100644
--- a/Processes/TrackingLine/TrackingLine.h
+++ b/Processes/TrackingLine/TrackingLine.h
@@ -1,14 +1,32 @@
-#ifndef _include_corsika_processes_TrackinLine_h_
-#define _include_corsika_processes_TrackinLine_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_processes_TrackingLine_h_
+#define _include_corsika_processes_TrackingLine_h_
 
 #include <corsika/geometry/Point.h>
+#include <corsika/geometry/Sphere.h>
 #include <corsika/geometry/Vector.h>
 
 #include <corsika/units/PhysicalUnits.h>
 
+#include <corsika/environment/Environment.h>
 #include <corsika/setup/SetupStack.h>
 #include <corsika/setup/SetupTrajectory.h>
 
+#include <algorithm>
+#include <iostream>
+#include <optional>
+#include <stdexcept>
+#include <utility>
+
 using namespace corsika;
 
 namespace corsika::process {
@@ -16,15 +34,96 @@ namespace corsika::process {
   namespace tracking_line {
 
     template <typename Stack>
-    class TrackingLine { // Newton-step, naja.. not yet
+    class TrackingLine { //
       typedef typename Stack::ParticleType Particle;
 
+      corsika::environment::Environment const& fEnvironment;
+
     public:
+      std::optional<std::pair<corsika::units::si::TimeType, corsika::units::si::TimeType>>
+      TimeOfIntersection(corsika::geometry::Line const& traj,
+                         geometry::Sphere const& sphere) {
+        using namespace corsika::units::si;
+        auto const& cs = fEnvironment.GetCoordinateSystem();
+        geometry::Point const origin(cs, 0_m, 0_m, 0_m);
+
+        auto const r0 = (traj.GetR0() - origin);
+        auto const v0 = traj.GetV0();
+        auto const c0 = (sphere.GetCenter() - origin);
+
+        auto const alpha = r0.dot(v0) - 2 * v0.dot(c0);
+        auto const beta = c0.squaredNorm() + r0.squaredNorm() + 2 * c0.dot(r0) -
+                          sphere.GetRadius() * sphere.GetRadius();
+
+        auto const discriminant = alpha * alpha - 4 * beta * v0.squaredNorm();
+
+        //~ std::cout << "discriminant: " << discriminant << std::endl;
+        //~ std::cout << "alpha: " << alpha << std::endl;
+        //~ std::cout << "beta: " << beta << std::endl;
+
+        if (discriminant.magnitude() > 0) {
+          (-alpha - sqrt(discriminant)) / (2 * v0.squaredNorm());
+          return std::make_pair((-alpha - sqrt(discriminant)) / (2 * v0.squaredNorm()),
+                                (-alpha + sqrt(discriminant)) / (2 * v0.squaredNorm()));
+        } else {
+          return {};
+        }
+      }
+
+      TrackingLine(corsika::environment::Environment const& pEnv)
+          : fEnvironment(pEnv) {}
       void Init() {}
-      setup::Trajectory GetTrack(Particle& p) {
-        geometry::Vector<SpeedType::dimension_type> v = p.GetDirection();
-        geometry::Line traj(p.GetPosition(), v);
-        return geometry::Trajectory<corsika::geometry::Line>(traj, 100_ns);
+
+      auto GetTrack(Particle& p) {
+        using namespace corsika::units::si;
+        geometry::Vector<SpeedType::dimension_type> const velocity =
+            p.GetMomentum() / p.GetEnergy() * corsika::units::si::constants::cSquared;
+
+        auto const currentPosition = p.GetPosition();
+        geometry::Line line(currentPosition, velocity);
+
+        auto const* currentVolumeNode =
+            fEnvironment.GetUniverse()->GetContainingNode(currentPosition);
+        auto const& children = currentVolumeNode->GetChildNodes();
+        auto const& excluded = currentVolumeNode->GetExcludedNodes();
+
+        std::vector<TimeType> intersectionTimes;
+
+        auto addIfIntersects = [&](auto& vtn) {
+          auto const& volume = vtn.GetVolume();
+          auto const& sphere = dynamic_cast<geometry::Sphere const&>(
+              volume); // for the moment we are a bit bold here and assume
+                       // everything is a sphere, crashes with exception if not
+
+          if (auto opt = TimeOfIntersection(line, sphere); opt.has_value()) {
+            auto const [t1, t2] = *opt;
+            if (t1.magnitude() >= 0)
+              intersectionTimes.push_back(t1);
+            else if (t2.magnitude() >= 0)
+              intersectionTimes.push_back(t2);
+          }
+        };
+
+        for (auto const& child : children) { addIfIntersects(*child); }
+
+        for (auto const* child : excluded) { addIfIntersects(*child); }
+
+        addIfIntersects(*currentVolumeNode);
+
+        auto const minIter =
+            std::min_element(intersectionTimes.cbegin(), intersectionTimes.cend());
+
+        TimeType min;
+
+        if (minIter == intersectionTimes.cend()) {
+          min = 1_s; // todo: do sth. more reasonable as soon as tracking is able
+                     // to handle the numerics properly
+          //~ throw std::runtime_error("no intersection with anything!");
+        } else {
+          min = *minIter;
+        }
+
+        return geometry::Trajectory<corsika::geometry::Line>(line, min);
       }
     };
 
diff --git a/Processes/TrackingLine/testTrackingLine.cc b/Processes/TrackingLine/testTrackingLine.cc
new file mode 100644
index 0000000000000000000000000000000000000000..ef44f2390a5071f645e383bff23373d8ec240993
--- /dev/null
+++ b/Processes/TrackingLine/testTrackingLine.cc
@@ -0,0 +1,98 @@
+/**
+ * (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/environment/Environment.h>
+
+#include <corsika/process/tracking_line/TrackingLine.h>
+
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/Vector.h>
+#include <corsika/geometry/Sphere.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;
+
+struct DummyParticle {
+    EnergyType fEnergy;
+    Vector<momentum_d> fMomentum;
+    Point fPosition;
+    
+    DummyParticle(EnergyType pEnergy, Vector<momentum_d> pMomentum, Point pPosition) : fEnergy(pEnergy), fMomentum(pMomentum), fPosition(pPosition) {}
+    
+    auto GetEnergy() const { return fEnergy; }
+    auto GetMomentum() const { return fMomentum; }
+    auto GetPosition() const { return fPosition; }
+};
+
+struct DummyStack {
+    using ParticleType = DummyParticle;
+};
+
+TEST_CASE("TrackingLine") {
+  corsika::environment::Environment env; // dummy environment
+  auto const& cs = env.GetCoordinateSystem();
+  
+  tracking_line::TrackingLine<DummyStack> tracking(env);
+
+  SECTION("intersection with sphere") {
+    Point const origin(cs, {0_m, 0_m, 0_m});
+    Point const center(cs, {0_m, 0_m, 10_m});
+    Sphere const sphere(center, 1_m);
+    Vector<corsika::units::si::SpeedType::dimension_type> v(cs, 0_m/second, 0_m/second, 1_m/second);
+    Line line(origin, v);
+    
+    geometry::Trajectory<Line> traj(line, 12345_s);
+    
+    auto const opt = tracking.TimeOfIntersection(traj, Sphere(Point(cs, {0_m, 0_m, 10_m}), 1_m));
+    REQUIRE(opt.has_value());
+    
+    auto [t1, t2] = opt.value();
+    REQUIRE(t1 / 9_s == Approx(1));
+    REQUIRE(t2 / 11_s == Approx(1));
+    
+    auto const optNoIntersection = tracking.TimeOfIntersection(traj, Sphere(Point(cs, {5_m, 0_m, 10_m}), 1_m));
+    REQUIRE_FALSE(optNoIntersection.has_value());
+  }
+  
+  SECTION("maximally possible propagation") {
+      auto& universe = *(env.GetUniverse());
+      
+      //~ std::cout << env.GetUniverse().get() << std::endl;
+      
+      DummyParticle p(1_J, Vector<momentum_d>(cs, 0*kilogram*meter/second, 0*kilogram*meter/second, 1*kilogram*meter/second), Point(cs, 0_m, 0_m,0_m));
+      
+      auto const radius = 20_m;
+      
+      auto theMedium = corsika::environment::Environment::CreateNode<Sphere>(
+      Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, radius);      
+      universe.AddChild(std::move(theMedium));
+      
+      Point const origin(cs, {0_m, 0_m, 0_m});
+      Vector<corsika::units::si::SpeedType::dimension_type> v(cs, 0_m/second, 0_m/second, 1_m/second);
+      Line line(origin, v);
+      
+      auto const traj = tracking.GetTrack(p);
+      
+      REQUIRE((traj.GetPosition(1.) - Point(cs, 0_m, 0_m, radius)).GetComponents(cs).norm().magnitude() == Approx(0).margin(1e-4));
+  }
+}
diff --git a/Setup/SetupEnvironment.h b/Setup/SetupEnvironment.h
index 1beb2748052770fc25b84991fdf3e651744c3dc8..e822b64ec0dd4b393fd7285360c96879ce2ccc0f 100644
--- a/Setup/SetupEnvironment.h
+++ b/Setup/SetupEnvironment.h
@@ -12,6 +12,10 @@
 #ifndef _include_corsika_setup_environment_h_
 #define _include_corsika_setup_environment_h_
 
-namespace corsika {}
+#include <corsika/environment/IMediumModel.h>
+
+namespace corsika::setup {
+    using IEnvironmentModel = corsika::environment::IMediumModel;
+}
 
 #endif
diff --git a/Setup/SetupTrajectory.h b/Setup/SetupTrajectory.h
index f1dfd497a6bec31eeb9e4d84967a8b80f380ebb0..acbbc373ac2c3b8454a0bcfb1fb42c84b33f49f5 100644
--- a/Setup/SetupTrajectory.h
+++ b/Setup/SetupTrajectory.h
@@ -26,9 +26,12 @@ namespace corsika::setup {
   using corsika::geometry::Line;
 
   /// definition of Trajectory base class, to be used in tracking and cascades
+  typedef corsika::geometry::Trajectory<Line> Trajectory;
+
+  /*
   typedef std::variant<std::monostate, corsika::geometry::Trajectory<Line>,
                        corsika::geometry::Trajectory<Helix>>
-      Trajectory;
+                       Trajectory;
 
   /// helper visitor to modify Particle by moving along Trajectory
   template <typename Particle>
@@ -58,7 +61,7 @@ namespace corsika::setup {
       return trajectory.GetDuration();
     }
   };
-
+  */
 } // namespace corsika::setup
 
 #endif
diff --git a/Stack/SuperStupidStack/SuperStupidStack.h b/Stack/SuperStupidStack/SuperStupidStack.h
index 064edffa2c4497743cfd4a5e5ef2caca6c86860a..d1c8bc2b38622698b359afe8805d22403d99a4ed 100644
--- a/Stack/SuperStupidStack/SuperStupidStack.h
+++ b/Stack/SuperStupidStack/SuperStupidStack.h
@@ -89,8 +89,11 @@ namespace corsika::stack {
       void Init() {}
 
       void Clear() {
-        fDataE.clear();
         fDataPID.clear();
+        fDataE.clear();
+        fMomentum.clear();
+        fPosition.clear();
+        fTime.clear();
       }
 
       int GetSize() const { return fDataPID.size(); }