corsika/media/UniformRefractiveIndex.hpp

@@ -19,7 +19,7 @@ namespace corsika {
    * for all evaluated locations.
    */
   template <typename T>
-  class UniformRefractiveIndex final : public T {
+  class UniformRefractiveIndex : public T {
 
     double n_; ///< The constant refractive index that we use.
 

corsika/modules/conex/CONEX_f.hpp

@@ -36,13 +36,18 @@ namespace conex {
 
   extern double double_rndm_interface();
 
-  extern "C" {}
-
   // the CONEX fortran interface
 
   extern "C" {
   extern struct { std::array<double, 16> dptl; } cxoptl_;
 
+  //! common block for atmosphere composition
+  extern struct {
+    std::array<double, 3> airz, aira, airw; //!< nuclear Z, A, composition fraction
+    double airavz, airava;                  //!< average Z, A
+    std::array<double, 3> airi;             //!< ionization potential, not used in cxroot
+  } cxair_;
+
   void cegs4_(int&, int&);
 
   void initconex_(int&, int*, int&, int&,

corsika/modules/radio/RadioProcess.hpp

@@ -7,11 +7,8 @@
  */
 #pragma once
 
-#include <istream>
-#include <fstream>
-#include <iostream>
-#include <string>
 #include <corsika/output/BaseOutput.hpp>
+#include <corsika/output/ParquetStreamer.hpp>
 #include <corsika/framework/process/ContinuousProcess.hpp>
 #include <corsika/framework/core/Step.hpp>
 #include <corsika/setup/SetupStack.hpp>
@@ -49,9 +46,11 @@ namespace corsika {
   protected:
     TAntennaCollection& antennas_; ///< The radio antennas we store into.
     TPropagator propagator_;       ///< The propagator implementation.
-    int event_{0};                 ///< The current event ID.
+    unsigned int showerId_{0};     ///< The current event ID.
+    ParquetStreamer output_;       //!< The parquet streamer for this process.
 
   public:
+    using axistype = std::vector<long double>;
     /**
      * Construct a new RadioProcess.
      */

corsika/modules/radio/antennas/Antenna.hpp

@@ -7,10 +7,10 @@
  */
 #pragma once
 
-#include <cnpy.hpp>
 #include <boost/filesystem.hpp>
 #include <corsika/framework/geometry/Point.hpp>
 #include <corsika/framework/core/PhysicalGeometry.hpp>
+#include <corsika/output/ParquetStreamer.hpp>
 
 namespace corsika {
 
@@ -28,7 +28,6 @@ namespace corsika {
     std::string const name_;                     ///< The name/identifier of this antenna.
     Point const location_;                       ///< The location of this antenna.
     CoordinateSystemPtr const coordinateSystem_; ///< The coordinate system of the antenna
-    std::string filename_ = ""; ///< The filename for the output file for this antenna.
 
   public:
     using axistype = std::vector<long double>;
@@ -102,19 +101,6 @@ namespace corsika {
      */
     std::vector<double> const& getWaveformZ() const;
 
-    /**
-     * Prepare for the start of the library.
-     */
-    void startOfLibrary(boost::filesystem::path const& directory,
-                        std::string const radioImplementation);
-
-    /**
-     * Flush the data from this shower to disk.
-     */
-    void endOfShower(int const event, std::string const& radioImplementation,
-                     double const sampleRate);
-
-  protected:
     /**
      * Get a reference to the underlying radio implementation.
      */

corsika/modules/radio/antennas/TimeDomainAntenna.hpp

@@ -33,10 +33,6 @@ namespace corsika {
 
   public:
     // import the methods from the antenna
-
-    // label this as a time-domain antenna.
-    static constexpr bool is_time_domain{true};
-
     using Antenna<TimeDomainAntenna>::getName;
     using Antenna<TimeDomainAntenna>::getLocation;
 
@@ -98,6 +94,14 @@ namespace corsika {
      */
     auto const& getWaveformZ() const;
 
+    /**
+     * Return a label that indicates that this is a time
+     * domain antenna
+     *
+     * This returns the string "Time".
+     */
+    std::string const getDomainLabel();
+
     /**
      * Creates time-units of each waveform.
      *
@@ -110,7 +114,7 @@ namespace corsika {
      *
      * This returns them in nanoseconds for ease of use.
      */
-    auto const& getAxis() const;
+    auto const getAxis() const;
 
     /**
      * Returns the sampling rate of the time domain antenna.

corsika/modules/writers/ParticleWriterParquet.hpp

@@ -71,7 +71,7 @@ namespace corsika {
     double countHadrons_ = 0; ///< count hadrons hitting plane
     double countMuons_ = 0;   ///< count muons hitting plane
     double countEM_ = 0;      ///< count EM particles hitting plane.
-    double countOthers_ = 0;  ///< count othe types of particles hitting plane
+    double countOthers_ = 0;  ///< count other types of particles hitting plane
 
     HEPEnergyType totalEnergy_; ///< energy absorbed in ground.
 

examples/clover_leaf.cpp

 /*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
  *
  * This software is distributed under the terms of the GNU General Public
  * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
@@ -37,8 +37,6 @@
 #include <corsika/modules/radio/detectors/AntennaCollection.hpp>
 #include <corsika/modules/radio/propagators/StraightPropagator.hpp>
 #include <corsika/modules/radio/propagators/SimplePropagator.hpp>
-#include <corsika/modules/radio/propagators/SignalPath.hpp>
-#include <corsika/modules/radio/propagators/RadioPropagator.hpp>
 #include <corsika/modules/TrackWriter.hpp>
 
 #include <corsika/modules/StackInspector.hpp>
@@ -47,13 +45,13 @@
 //#include <corsika/modules/TrackWriter.hpp>
 
 /*
-  NOTE, WARNING, ATTENTION
+ NOTE, WARNING, ATTENTION
 
-  The .../Random.hpppp implement the hooks of external modules to the C8 random
-  number generator. It has to occur excatly ONCE per linked
-  executable. If you include the header below multiple times and
-  link this together, it will fail.
- */
+ The .../Random.hpppp implement the hooks of external modules to the C8 random
+ number generator. It has to occur excatly ONCE per linked
+ executable. If you include the header below multiple times and
+ link this together, it will fail.
+*/
 #include <corsika/modules/sibyll/Random.hpp>
 #include <corsika/modules/urqmd/Random.hpp>
 
@@ -72,7 +70,7 @@ using namespace std;
 //
 int main() {
 
-  logging::set_level(logging::level::warn);
+  logging::set_level(logging::level::info);
   corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
 
   CORSIKA_LOG_INFO("Synchrotron radiation");
@@ -210,11 +208,10 @@ int main() {
   // assemble all processes into an ordered process list
   auto sequence = make_sequence(coreas, cut);
 
+  output.startOfLibrary();
   // define air shower object, run simulation
   Cascade EAS(env, tracking, sequence, output, stack);
-  output.startOfShower();
   EAS.run();
-  output.endOfShower();
 
   CORSIKA_LOG_INFO("|p| electron = {} and E electron = {}", plab.getNorm(), Elab);
   CORSIKA_LOG_INFO("period: {}", period);
@@ -225,4 +222,4 @@ int main() {
   CORSIKA_LOG_INFO("gamma: {}", gammaP);
 
   output.endOfLibrary();
-}
+}
\ No newline at end of file

examples/corsika.cpp

@@ -19,6 +19,7 @@
 #include <corsika/framework/geometry/PhysicalGeometry.hpp>
 #include <corsika/framework/geometry/Plane.hpp>
 #include <corsika/framework/geometry/Sphere.hpp>
+#include <corsika/framework/process/DynamicInteractionProcess.hpp>
 #include <corsika/framework/process/InteractionCounter.hpp>
 #include <corsika/framework/process/ProcessSequence.hpp>
 #include <corsika/framework/process/SwitchProcessSequence.hpp>
@@ -65,8 +66,10 @@
 #include <CLI/Config.hpp>
 #include <CLI/Formatter.hpp>
 
+#include <cstdlib>
 #include <iomanip>
 #include <limits>
+#include <string>
 #include <string_view>
 
 /*
@@ -165,6 +168,10 @@ int main(int argc, char** argv) {
       ->default_val("info")
       ->check(CLI::IsMember({"warn", "info", "debug", "trace"}))
       ->group("Misc.");
+  app.add_option("-M,--hadronModel", "High-energy hadronic interaction model")
+      ->default_val("SIBYLL-2.3d")
+      ->check(CLI::IsMember({"SIBYLL-2.3d", "QGSJet-II.04", "EPOS-LHC"}))
+      ->group("Misc.");
 
   // parse the command line options into the variables
   CLI11_PARSE(app, argc, argv);
@@ -287,8 +294,26 @@ int main(int argc, char** argv) {
   TrackWriter tracks;
   output.add("tracks", tracks);
 
-  corsika::sibyll::Interaction sibyll{env};
-  InteractionCounter sibyllCounted{sibyll};
+  DynamicInteractionProcess<setup::Stack<EnvType>> heModel;
+
+  // have SIBYLL always for PROPOSAL photo-hadronic interactions
+  auto sibyll = std::make_shared<corsika::sibyll::Interaction>(env);
+
+  if (auto const modelStr = app["--hadronModel"]->as<std::string_view>();
+      modelStr == "SIBYLL-2.3d") {
+    heModel = DynamicInteractionProcess<setup::Stack<EnvType>>{sibyll};
+  } else if (modelStr == "QGSJet-II.04") {
+    heModel = DynamicInteractionProcess<setup::Stack<EnvType>>{
+        std::make_shared<corsika::qgsjetII::Interaction>()};
+  } else if (modelStr == "EPOS-LHC") {
+    heModel = DynamicInteractionProcess<setup::Stack<EnvType>>{
+        std::make_shared<corsika::epos::Interaction>()};
+  } else {
+    CORSIKA_LOG_CRITICAL("invalid choice \"{}\"; also check argument parser", modelStr);
+    return EXIT_FAILURE;
+  }
+
+  InteractionCounter heCounted{heModel};
 
   corsika::pythia8::Decay decayPythia;
 
@@ -327,7 +352,7 @@ int main(int argc, char** argv) {
   HEPEnergyType heHadronModelThreshold = 63.1_GeV;
 
   corsika::proposal::Interaction emCascade(
-      env, sophia, sibyll.getHadronInteractionModel(), heHadronModelThreshold);
+      env, sophia, sibyll->getHadronInteractionModel(), heHadronModelThreshold);
 
   // use BetheBlochPDG for hadronic continuous losses, and proposal otherwise
   corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>> emContinuousProposal(
@@ -356,11 +381,9 @@ int main(int argc, char** argv) {
     bool operator()(const Particle& p) const { return (p.getKineticEnergy() < cutE_); }
   };
   auto hadronSequence =
-      make_select(EnergySwitch(heHadronModelThreshold), urqmdCounted, sibyllCounted);
+      make_select(EnergySwitch(heHadronModelThreshold), urqmdCounted, heCounted);
   auto decaySequence = make_sequence(decayPythia, decaySibyll);
 
-  TrackWriter trackWriter{tracks};
-
   // observation plane
   Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
   ObservationPlane<setup::Tracking, ParticleWriterParquet> observationLevel{
@@ -430,9 +453,7 @@ int main(int argc, char** argv) {
         Efinal / 1_GeV, dEdX.getEnergyLost() / 1_GeV,
         observationLevel.getEnergyGround() / 1_GeV, (Efinal / E0 - 1) * 100);
 
-    // auto const hists = heModelCounted.getHistogram() +
-    // urqmdCounted.getHistogram();
-    auto const hists = sibyllCounted.getHistogram() + urqmdCounted.getHistogram();
+    auto const hists = heCounted.getHistogram() + urqmdCounted.getHistogram();
 
     save_hist(hists.labHist(), labHist_file, true);
     save_hist(hists.CMSHist(), cMSHist_file, true);
@@ -440,4 +461,6 @@ int main(int argc, char** argv) {
 
   // and finalize the output on disk
   output.endOfLibrary();
+
+  return EXIT_SUCCESS;
 }

examples/em_shower.cpp

@@ -106,7 +106,7 @@ int main(int argc, char** argv) {
   // build a Linsley US Standard atmosphere into `env`
   create_5layer_atmosphere<EnvironmentInterface, MyExtraEnv>(
       env, AtmosphereId::LinsleyUSStd, center, Medium::AirDry1Atm,
-      MagneticFieldVector{rootCS, 20.4_uT, 0_T, 43.23_uT});
+      MagneticFieldVector{rootCS, 20.4_uT, 0_T, -43.23_uT});
 
   std::unordered_map<Code, HEPEnergyType> energy_resolution = {
       {Code::Electron, 2_MeV},

examples/radio_em_shower.cpp

@@ -51,8 +51,6 @@
 #include <corsika/modules/radio/detectors/AntennaCollection.hpp>
 #include <corsika/modules/radio/propagators/StraightPropagator.hpp>
 #include <corsika/modules/radio/propagators/SimplePropagator.hpp>
-#include <corsika/modules/radio/propagators/SignalPath.hpp>
-#include <corsika/modules/radio/propagators/RadioPropagator.hpp>
 
 #include <corsika/setup/SetupStack.hpp>
 #include <corsika/setup/SetupTrajectory.hpp>

examples/synchrotron_test_C8tracking.cpp

@@ -158,11 +158,10 @@ int main() {
   // assemble all processes into an ordered process list
   auto sequence = make_sequence(coreas, zhs, cut);
 
+  output.startOfLibrary();
   // define air shower object, run simulation
   Cascade EAS(env, tracking, sequence, output, stack);
-  output.startOfShower();
   EAS.run();
-  output.endOfShower();
 
   CORSIKA_LOG_INFO("|p| = {} and E = {}", plab.getNorm(), Elab);
   CORSIKA_LOG_INFO("period: {}", period);

examples/vertical_EAS.cpp

@@ -135,7 +135,7 @@ int main(int argc, char** argv) {
   // build a Linsley US Standard atmosphere into `env`
   create_5layer_atmosphere<EnvironmentInterface, MyExtraEnv>(
       env, AtmosphereId::LinsleyUSStd, center, Medium::AirDry1Atm,
-      MagneticFieldVector{rootCS, 20.4_uT, 0_T, 43.23_uT});
+      MagneticFieldVector{rootCS, 20.4_uT, 0_T, -43.23_uT});
 
   // Uncomment if you want to use PROPOSAL
   //  std::unordered_map<Code, HEPEnergyType> energy_resolution = {

conex @ c26ca0fa

-Subproject commit 344c875bed534590eb86aa5e12b0382a9868fc64
+Subproject commit c26ca0fada6ef182147d014bb7d413814f9cf202

tests/framework/CMakeLists.txt

@@ -4,6 +4,7 @@ set (test_framework_sources
   testFourVector.cpp
   testSaveBoostHistogram.cpp
   testClassTimer.cpp
+  testDynamicInteractionProcess.cpp;
   testLogging.cpp
   testParticles.cpp
   testStackInterface.cpp

tests/framework/testDynamicInteractionProcess.cpp

+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <corsika/framework/process/DynamicInteractionProcess.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/FourVector.hpp>
+#include <corsika/framework/geometry/CoordinateSystem.hpp>
+
+#include <catch2/catch.hpp>
+
+using namespace corsika;
+
+struct DummyStack {
+  using stack_view_type = int&;
+};
+
+struct DummyProcessA : public InteractionProcess<DummyProcessA> {
+  int id_;
+
+  DummyProcessA(int id)
+      : id_{id} {}
+
+  // mockup to "identify" the process via its ID
+  CrossSectionType getCrossSection(Code, Code, FourMomentum const&,
+                                   FourMomentum const&) const {
+    return id_ * 10_mb;
+  }
+
+  // mockup: set argument to our ID
+  template <typename TArgument>
+  void doInteraction(TArgument& argument, Code, Code, FourMomentum const&,
+                     FourMomentum const&) {
+    argument = id_ * 10;
+  }
+
+  // to make sure we can test proper handling of dangling references and
+  // stuff like that
+  ~DummyProcessA() { id_ = 0; }
+
+  // prevent copying
+  DummyProcessA(DummyProcessA const&) = delete;
+  DummyProcessA& operator=(DummyProcessA const&) = delete;
+
+  // allow moving
+  DummyProcessA(DummyProcessA&&) = default;
+  DummyProcessA& operator=(DummyProcessA&&) = default;
+};
+
+// same as above to have another class
+struct DummyProcessB : public InteractionProcess<DummyProcessB> {
+  int id_;
+
+  DummyProcessB(int id)
+      : id_{id} {}
+
+  CrossSectionType getCrossSection(Code, Code, FourMomentum const&,
+                                   FourMomentum const&) const {
+    return id_ * 1_mb;
+  }
+
+  template <typename TArgument>
+  void doInteraction(TArgument& argument, Code, Code, FourMomentum const&,
+                     FourMomentum const&) {
+    argument = id_;
+  }
+
+  // to make sure we can test proper handling of dangling references and
+  // stuff like that
+  ~DummyProcessB() { id_ = 0; }
+
+  // prevent copying
+  DummyProcessB(DummyProcessB const&) = delete;
+  DummyProcessB& operator=(DummyProcessB const&) = delete;
+
+  // allow moving
+  DummyProcessB(DummyProcessB&&) = default;
+  DummyProcessB& operator=(DummyProcessB&&) = default;
+};
+
+TEST_CASE("DynamicInteractionProcess", "[process]") {
+  auto const rootCS = get_root_CoordinateSystem();
+  FourMomentum const fourVec{1_GeV, {rootCS, 1_GeV, 1_GeV, 1_GeV}};
+
+  std::shared_ptr<DummyProcessA> a = std::make_shared<DummyProcessA>(1);
+  std::shared_ptr<DummyProcessB> b = std::make_shared<DummyProcessB>(2);
+
+  SECTION("getCrossSection") {
+    DynamicInteractionProcess<DummyStack> dynamic;
+
+    dynamic = a;
+    REQUIRE(dynamic.getCrossSection(Code::Electron, Code::Electron, fourVec, fourVec) /
+                10_mb ==
+            Approx(1).epsilon(1e-6));
+
+    dynamic = b;
+    REQUIRE(dynamic.getCrossSection(Code::Electron, Code::Electron, fourVec, fourVec) /
+                2_mb ==
+            Approx(1).epsilon(1e-6));
+
+    dynamic = DynamicInteractionProcess<DummyStack>{std::make_shared<DummyProcessA>(3)};
+    REQUIRE(dynamic.getCrossSection(Code::Electron, Code::Electron, fourVec, fourVec) /
+                30_mb ==
+            Approx(1).epsilon(1e-6));
+
+    dynamic = std::make_shared<DummyProcessB>(4);
+    REQUIRE(dynamic.getCrossSection(Code::Electron, Code::Electron, fourVec, fourVec) /
+                4_mb ==
+            Approx(1).epsilon(1e-6));
+
+    // test process going out of scope
+    { dynamic = std::make_shared<DummyProcessB>(5); }
+    REQUIRE(dynamic.getCrossSection(Code::Electron, Code::Electron, fourVec, fourVec) /
+                5_mb ==
+            Approx(1).epsilon(1e-6));
+  }
+
+  SECTION("doInteraction") {
+    int value{};
+
+    DynamicInteractionProcess<DummyStack> dynamic;
+
+    dynamic = a;
+    dynamic.doInteraction(value, Code::Electron, Code::Electron, fourVec, fourVec);
+    REQUIRE(value == 10);
+
+    dynamic = b;
+    dynamic.doInteraction(value, Code::Electron, Code::Electron, fourVec, fourVec);
+    REQUIRE(value == 2);
+
+    dynamic = DynamicInteractionProcess<DummyStack>{std::make_shared<DummyProcessA>(3)};
+    dynamic.doInteraction(value, Code::Electron, Code::Electron, fourVec, fourVec);
+    REQUIRE(value == 30);
+
+    dynamic = std::make_shared<DummyProcessB>(4);
+    dynamic.doInteraction(value, Code::Electron, Code::Electron, fourVec, fourVec);
+    REQUIRE(value == 4);
+  }
+}

tests/media/testRefractiveIndex.cpp

@@ -20,6 +20,7 @@
 #include <corsika/media/NuclearComposition.hpp>
 #include <corsika/media/UniformRefractiveIndex.hpp>
 #include <corsika/media/ExponentialRefractiveIndex.hpp>
+#include <corsika/media/GladstoneDaleRefractiveIndex.hpp>
 #include <corsika/media/CORSIKA7Atmospheres.hpp>
 
 #include <SetupTestTrajectory.hpp>
@@ -31,8 +32,11 @@ using namespace corsika;
 template <typename TInterface>
 using MyExtraEnv =
     ExponentialRefractiveIndex<MediumPropertyModel<UniformMagneticField<TInterface>>>;
+template <typename TInterface2>
+using MyExtraEnv2 =
+    GladstoneDaleRefractiveIndex<MediumPropertyModel<UniformMagneticField<TInterface2>>>;
 
-TEST_CASE("UniformRefractiveIndex w/ Homogeneous") {
+TEST_CASE("UniformRefractiveIndex w/ Homogeneous medium") {
 
   logging::set_level(logging::level::info);
 
@@ -40,7 +44,7 @@ TEST_CASE("UniformRefractiveIndex w/ Homogeneous") {
 
   Point const gOrigin(gCS, {0_m, 0_m, 0_m});
 
-  // setup our interface types
+  // set up our interface types
   using IModelInterface = IRefractiveIndexModel<IMediumModel>;
   using AtmModel = UniformRefractiveIndex<HomogeneousMedium<IModelInterface>>;
 
@@ -50,7 +54,7 @@ TEST_CASE("UniformRefractiveIndex w/ Homogeneous") {
   // the composition we use for the homogenous medium
   NuclearComposition const protonComposition({Code::Proton}, {1.});
 
-  // the refrative index that we use
+  // the refractive index that we use
   const double n{1.000327};
 
   // create the atmospheric model
@@ -213,3 +217,78 @@ TEST_CASE("ExponentialRefractiveIndex w/ 5-layered atmosphere") {
 
   CHECK(rIndex - n0 == Approx(0));
 }
+
+TEST_CASE("GladstoneDaleRefractiveIndex w/ Homogeneous medium") {
+
+  logging::set_level(logging::level::info);
+
+  // get a CS and a point
+  CoordinateSystemPtr const& gCS = get_root_CoordinateSystem();
+
+  Point const gOrigin(gCS, {0_m, 0_m, 0_m});
+
+  // setup interface types
+  using IModelInterface = IRefractiveIndexModel<IMediumModel>;
+  using AtmModel = GladstoneDaleRefractiveIndex<HomogeneousMedium<IModelInterface>>;
+
+  // the constant density
+  const auto density{19.2_g / cube(1_cm)};
+
+  // the composition we use for the homogenous medium
+  NuclearComposition const protonComposition({Code::Proton}, {1.});
+
+  // the refractive index at sea level
+  const double n0{1.000327};
+
+  // a point at the surface of the earth
+  Point const surface_{gCS, 0_m, 0_m, constants::EarthRadius::Mean};
+
+  // a random point in the atmosphere
+  Point const p1_{gCS, 1_km, 1_km, constants::EarthRadius::Mean + 10_km};
+
+  // create the atmospheric model and check refractive index
+  AtmModel medium(n0, surface_, density, protonComposition);
+
+  CHECK(n0 - medium.getRefractiveIndex(surface_) == Approx(0));
+  CHECK(n0 - medium.getRefractiveIndex(p1_) == Approx(0));
+}
+
+TEST_CASE("GladstoneDaleRefractiveIndex w/ 5-layered atmosphere") {
+
+  logging::set_level(logging::level::info);
+
+  // get a CS
+  CoordinateSystemPtr const& gCS = get_root_CoordinateSystem();
+
+  // the center of the earth
+  Point const center_{gCS, 0_m, 0_m, 0_m};
+
+  // a point at the surface of the earth
+  Point const surface_{gCS, 0_m, 0_m, constants::EarthRadius::Mean};
+
+  // the refractive index at sea level
+  const double n0{1.000327};
+
+  // a reference point to calculate the refractive index there
+  Point const ref_{gCS, 0_km, 0_km, constants::EarthRadius::Mean + 10_km};
+
+  // setup a 5-layered environment
+  using EnvironmentInterface =
+      IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
+  using EnvType = Environment<EnvironmentInterface>;
+  EnvType env;
+
+  create_5layer_atmosphere<EnvironmentInterface, MyExtraEnv2>(
+      env, AtmosphereId::LinsleyUSStd, center_, n0, surface_, Medium::AirDry1Atm,
+      MagneticFieldVector{gCS, 0_T, 50_uT, 0_T});
+
+  // get the universe for this environment
+  auto const* const universe{env.getUniverse().get()};
+  auto const* node{universe->getContainingNode(ref_)};
+  // get the refractive index
+  auto const rIndex1{node->getModelProperties().getRefractiveIndex(ref_)};
+  auto const rIndex2{node->getModelProperties().getRefractiveIndex(surface_)};
+
+  CHECK(rIndex1 - n0 == Approx(-0.0002591034));
+  CHECK(rIndex2 - n0 == Approx(0));
+}
\ No newline at end of file

tests/modules/testRadio.cpp

@@ -7,6 +7,7 @@
  */
 #include <catch2/catch.hpp>
 
+#include <corsika/modules/radio/RadioProcess.hpp>
 #include <corsika/modules/radio/ZHS.hpp>
 #include <corsika/modules/radio/CoREAS.hpp>
 #include <corsika/modules/radio/antennas/TimeDomainAntenna.hpp>
@@ -65,9 +66,10 @@ TEST_CASE("Radio", "[processes]") {
 
   SECTION("CoREAS process") {
 
-    // This serves as a compiler test for any changes in the CoREAS algorithm
-    // Environment
+    // This serves as a compiler test for any changes in the CoREAS algorithm and
+    // check the radio process output
 
+    // Environment
     using EnvironmentInterface =
         IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
 
@@ -178,6 +180,23 @@ TEST_CASE("Radio", "[processes]") {
       REQUIRE(total < (1e-12 * ant.getWaveformX().size()));
     }
 
+    // coreas output check
+    std::string const implemencoreas{"CoREAS"};
+
+    boost::filesystem::path const tempPathC{boost::filesystem::temp_directory_path() /
+                                            ("test_corsika_radio_" + implemencoreas)};
+    if (boost::filesystem::exists(tempPathC)) {
+      boost::filesystem::remove_all(tempPathC);
+    }
+
+    boost::filesystem::create_directory(tempPathC);
+    coreas.startOfLibrary(tempPathC);
+    auto const outputFileC = tempPathC / ("antennas.parquet");
+    CHECK(boost::filesystem::exists(outputFileC));
+    // run end of shower and make sure that something extra was added
+    auto const fileSizeC = boost::filesystem::file_size(outputFileC);
+    coreas.endOfShower(0);
+    CHECK(boost::filesystem::file_size(outputFileC) > fileSizeC);
     coreas.endOfLibrary();
 
   } // END: SECTION("CoREAS process")
@@ -374,6 +393,25 @@ TEST_CASE("Radio", "[processes]") {
     // check doContinuous and simulate methods
     zhs.doContinuous(step, true);
 
+    // zhs output check
+    std::string const implemenzhs{"ZHS"};
+
+    boost::filesystem::path const tempPathZ{boost::filesystem::temp_directory_path() /
+                                            ("test_corsika_radio_" + implemenzhs)};
+    if (boost::filesystem::exists(tempPathZ)) {
+      boost::filesystem::remove_all(tempPathZ);
+    }
+
+    boost::filesystem::create_directory(tempPathZ);
+    zhs.startOfLibrary(tempPathZ);
+    auto const outputFileZ = tempPathZ / ("antennas.parquet");
+    CHECK(boost::filesystem::exists(outputFileZ));
+    // run end of shower and make sure that something extra was added
+    auto const fileSizeZ = boost::filesystem::file_size(outputFileZ);
+    zhs.endOfShower(0);
+    CHECK(boost::filesystem::file_size(outputFileZ) > fileSizeZ);
+    zhs.endOfLibrary();
+
   } // END: SECTION("ZHS process")
 
   SECTION("Radio extreme cases") {
@@ -399,6 +437,7 @@ TEST_CASE("Radio", "[processes]") {
     Point const point_1(rootCSRadio, {1_m, 1_m, 0_m});
     Point const point_2(rootCSRadio, {2_km, 1_km, 0_m});
     Point const point_4(rootCSRadio, {0_m, 1_m, 0_m});
+    const auto point_b{Point(rootCSRadio, 30000_m, 0_m, 0_m)};
 
     // create times for the antenna
     const TimeType start{0_s};
@@ -407,17 +446,26 @@ TEST_CASE("Radio", "[processes]") {
     const TimeType duration_dummy{2_s};
     const InverseTimeType sample_dummy{1_Hz};
 
+    // create specific times for antenna to do timebin check
+    const TimeType start_b{0.994e-4_s};
+    const TimeType duration_b{1.07e-4_s - 0.994e-4_s};
+    const InverseTimeType sampleRate_b{5e+11_Hz};
+
     // check that I can create an antenna at (1, 2, 3)
     TimeDomainAntenna ant1("antenna_name", point1, rootCSRadio, start, duration, sample,
                            start);
     TimeDomainAntenna ant2("dummy", point1, rootCSRadio, start, duration_dummy,
                            sample_dummy, start);
+    TimeDomainAntenna ant_b("timebin", point_b, rootCSRadio, start_b, duration_b,
+                            sampleRate_b, start_b);
     // construct a radio detector instance to store our antennas
     AntennaCollection<TimeDomainAntenna> detector;
     AntennaCollection<TimeDomainAntenna> detector_dummy;
+    AntennaCollection<TimeDomainAntenna> detector_b;
     // add the antennas to the detector
     detector.addAntenna(ant1);
     detector_dummy.addAntenna(ant2);
+    detector_b.addAntenna(ant_b);
 
     // create a new stack for each trial
     setup::Stack<EnvType> stack;
@@ -428,7 +476,7 @@ TEST_CASE("Radio", "[processes]") {
     const auto pmass{get_mass(particle)};
     // construct an energy
     const HEPEnergyType E0{1_TeV};
-    // compute the necessary momentumn
+    // compute the necessary momentum
     const HEPMomentumType P0{sqrt(E0 * E0 - pmass * pmass)};
     // and create the momentum vector
     const auto plab{MomentumVector(rootCSRadio, {P0, 0_GeV, 0_GeV})};
@@ -455,7 +503,18 @@ TEST_CASE("Radio", "[processes]") {
     VelocityVector vh{(point_4 - point1) / th};
     Line lh{point1, vh};
     StraightTrajectory track_h{lh, th};
+    StraightTrajectory track_h_neg_time{lh, -th};
     Step step_h(particle_stack, track_h);
+    Step step_h_neg_time(particle_stack, track_h_neg_time);
+
+    // feed radio with an electron track that ends in a different antenna bin.
+    Point const point_start(rootCSRadio, {100_m, 0_m, 0_m});
+    Point const point_end(rootCSRadio, {100_m, 0.00628319_m, 0_m});
+    TimeType tb{(point_end - point_start).getNorm() / (0.999 * constants::c)};
+    VelocityVector vb{(point_end - point_start) / tb};
+    Line lb{point_start, vb};
+    StraightTrajectory track_b{lb, tb};
+    Step step_b(particle_stack, track_b);
 
     // create radio process instances
     RadioProcess<decltype(detector),
@@ -467,13 +526,12 @@ TEST_CASE("Radio", "[processes]") {
                  ZHS<decltype(detector), decltype(SimplePropagator(envRadio))>,
                  decltype(SimplePropagator(envRadio))>
         zhs(detector, envRadio);
-
     coreas.doContinuous(step_proton, true);
     zhs.doContinuous(step_proton, true);
     coreas.doContinuous(step_h, true);
     zhs.doContinuous(step_h, true);
+    zhs.doContinuous(step_h_neg_time, true);
 
-    // create radio processes with "dummy" antenna to trigger extreme time-binning
     RadioProcess<decltype(detector_dummy),
                  CoREAS<decltype(detector_dummy), decltype(SimplePropagator(envRadio))>,
                  decltype(SimplePropagator(envRadio))>
@@ -483,10 +541,24 @@ TEST_CASE("Radio", "[processes]") {
                  ZHS<decltype(detector_dummy), decltype(SimplePropagator(envRadio))>,
                  decltype(SimplePropagator(envRadio))>
         zhs_dummy(detector_dummy, envRadio);
-
+    coreas_dummy.doContinuous(step_proton, true);
+    zhs_dummy.doContinuous(step_proton, true);
     coreas_dummy.doContinuous(step_h, true);
     zhs_dummy.doContinuous(step_h, true);
 
+    // create radio process instances
+    RadioProcess<decltype(detector_b),
+                 CoREAS<decltype(detector_b), decltype(SimplePropagator(envRadio))>,
+                 decltype(SimplePropagator(envRadio))>
+        coreas_b(detector_b, envRadio);
+
+    RadioProcess<decltype(detector_b),
+                 ZHS<decltype(detector_b), decltype(SimplePropagator(envRadio))>,
+                 decltype(SimplePropagator(envRadio))>
+        zhs_b(detector_b, envRadio);
+    coreas_b.doContinuous(step_b, true);
+    zhs_b.doContinuous(step_b, true);
+
   } // END: SECTION("Radio extreme cases")
 
   SECTION("Process Library") {
@@ -574,7 +646,7 @@ TEST_CASE("Radio", "[processes]") {
     CHECK(config["antennas"]["antenna_name2"]["location"][0].as<double>() == 4);
     CHECK(config["antennas"]["antenna_name2"]["location"][1].as<double>() == 80);
     CHECK(config["antennas"]["antenna_name2"]["location"][2].as<double>() == 6);
-  }
+  } // END: SECTION("Process Library")
 
 } // END: TEST_CASE("Radio", "[processes]")
 
@@ -827,7 +899,7 @@ TEST_CASE("Antennas") {
 
   } // END: SECTION("TimeDomainAntenna AntennaCollection")
 
-  SECTION("TimeDomainAntenna Library") {
+  SECTION("TimeDomainAntenna Config File") {
     // Runs checks that the file readers are working properly
     Environment<IRefractiveIndexModel<IMediumModel>> env;
     const auto rootCS = env.getCoordinateSystem();
@@ -837,36 +909,24 @@ TEST_CASE("Antennas") {
     InverseTimeType const sampleRate(1_GHz);
     TimeType const groundHitTime(1e3_ns);
 
-    TimeDomainAntenna antenna("test_antenna", antPos, rootCS, tStart, duration,
-                              sampleRate, groundHitTime);
-
-    // Run the start of lib and end of shower functions on the antenna
-    std::vector<std::string> implementations{"CoREAS", "ZHS"};
-    for (auto const implemen : implementations) {
-      // For each implementation type, save a file
-      boost::filesystem::path const tempPath{boost::filesystem::temp_directory_path() /
-                                             ("test_corsika_radio_" + implemen)};
-      if (boost::filesystem::exists(tempPath)) {
-        boost::filesystem::remove_all(tempPath);
-      }
-      boost::filesystem::create_directory(tempPath);
-      antenna.startOfLibrary(tempPath, implemen);
-      auto const outputFile = tempPath / (antenna.getName() + ".npz");
-      CHECK(boost::filesystem::exists(outputFile));
-
-      // run end of shower and make sure that something extra was added
-      auto const fileSize = boost::filesystem::file_size(outputFile);
-      antenna.endOfShower(0, implemen, sampleRate / 1_Hz);
-      CHECK(boost::filesystem::file_size(outputFile) > fileSize);
-    }
+    TimeDomainAntenna antennaC("test_antennaCoREAS", antPos, rootCS, tStart, duration,
+                               sampleRate, groundHitTime);
+    TimeDomainAntenna antennaZ("test_antennaZHS", antPos, rootCS, tStart, duration,
+                               sampleRate, groundHitTime);
 
     // Check the YAML file output
-    auto const config = antenna.getConfig();
-    CHECK(config["type"].as<std::string>() == "TimeDomainAntenna");
-    CHECK(config["start_time"].as<double>() == tStart / 1_ns);
-    CHECK(config["duration"].as<double>() == duration / 1_ns);
-    CHECK(config["sample_rate"].as<double>() == sampleRate / 1_GHz);
-  } // END: SECTION("TimeDomainAntenna Library")
+    auto const configC = antennaC.getConfig();
+    CHECK(configC["type"].as<std::string>() == "TimeDomainAntenna");
+    CHECK(configC["start_time"].as<double>() == tStart / 1_ns);
+    CHECK(configC["duration"].as<double>() == duration / 1_ns);
+    CHECK(configC["sample_rate"].as<double>() == sampleRate / 1_GHz);
+
+    auto const configZ = antennaZ.getConfig();
+    CHECK(configZ["type"].as<std::string>() == "TimeDomainAntenna");
+    CHECK(configZ["start_time"].as<double>() == tStart / 1_ns);
+    CHECK(configZ["duration"].as<double>() == duration / 1_ns);
+    CHECK(configZ["sample_rate"].as<double>() == sampleRate / 1_GHz);
+  } // END: SECTION("TimeDomainAntenna Config File")
 
 } // END: TEST_CASE("Antennas")
 
@@ -1156,7 +1216,7 @@ TEST_CASE("Propagators") {
             Approx(0).margin(absMargin));
       CHECK(path.average_refractive_index_ == Approx(0.210275935));
       CHECK(path.refractive_index_source_ == Approx(2));
-      // CHECK(path.refractive_index_destination_ == Approx(0.0000000041));
+      CHECK(path.refractive_index_destination_ == Approx(4.12231e-09));
       CHECK(path.emit_.getComponents() == vvv1.getComponents());
       CHECK(path.receive_.getComponents() == vvv2.getComponents());
       CHECK(path.R_distance_ == 10_m);