TabulatedFlatAtmospherePropagator now handles all cases of particles coming from the shower.

corsika/detail/modules/radio/propagators/TabulatedFlatAtmospherePropagator.inl

 /*
- * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ * (c) Copyright 2023 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
@@ -21,9 +21,10 @@ namespace corsika {
       , upperLimit_(upperLimit)
       , lowerLimit_(lowerLimit)
       , step_(step)
-      , inverseStep_(1 / step) {
-    auto const maxHeight_ = upperLimit_.getCoordinates().getZ();
-    auto const minHeight_ = lowerLimit_.getCoordinates().getZ();
+      , inverseStep_(1 / step)
+      , maxHeight_(upperLimit_.getCoordinates().getZ())
+      , minHeight_(lowerLimit_.getCoordinates().getZ() - 1_km)
+  {
     auto const minX_ = lowerLimit_.getCoordinates().getX();
     auto const minY_ = lowerLimit_.getCoordinates().getY();
     std::size_t const nBins_ = (maxHeight_ - minHeight_) * inverseStep_ + 1;
@@ -49,23 +50,30 @@ namespace corsika {
     auto const intRerfZero_ = refractivityTable_.at(0) * stepOverMeter_;
     integratedRefractivityTable_.push_back(intRerfZero_);
     for (std::size_t i = 1; i < nBins_; i++) {
-      auto const intRefrI_ =
-          integratedRefractivityTable_[i - 1] + refractivityTable_[i] * stepOverMeter_;
+      auto const intRefrI_ = integratedRefractivityTable_[i-1] + refractivityTable_[i] * stepOverMeter_;
       integratedRefractivityTable_.push_back(intRefrI_);
     }
+
+    // this is used to interpolate refractivity and integrated refractivity for particles below the lower limit
+    slopeRefrLower_ = (refractivityTable_.at(10) - refractivityTable_.front()) / 10.;
+    slopeIntRefrLower_ = (integratedRefractivityTable_.at(10) - integratedRefractivityTable_.front()) / 10.;
+
+    // this is used to interpolate refractivity and integrated refractivity for particles above the upper limit
+    lastElement_ = integratedRefractivityTable_.size() - 1;
+    slopeRefrUpper_ = (refractivityTable_.at(lastElement_) - refractivityTable_.at(lastElement_ - 10)) / 10.;
+    slopeIntRefrUpper_ = (integratedRefractivityTable_.at(lastElement_) - integratedRefractivityTable_.at(lastElement_ - 10)) / 10.;
   };
 
   template <typename TEnvironment>
   template <typename Particle>
   inline typename TabulatedFlatAtmospherePropagator<TEnvironment>::SignalPathCollection
-  TabulatedFlatAtmospherePropagator<TEnvironment>::propagate(Particle const& particle, Point const& source,
-                                                             Point const& destination) {
+  TabulatedFlatAtmospherePropagator<TEnvironment>::propagate( [[maybe_unused]] Particle const& particle, Point const& source, Point const& destination) {
 
     /**
-     * This is a simple case of straight propagator where
-     * tabulated values of refractive index are called assuming
-     * a flat atmosphere.
-     *
+    * This is a simple case of straight propagator where
+    * tabulated values of refractive index are called assuming
+    * a flat atmosphere.
+    *
      */
 
     // these are used for the direction of emission and reception of signal at the antenna
@@ -78,33 +86,89 @@ namespace corsika {
     // clear the refractive index vector and points deque for this signal propagation.
     rindex.clear();
     points.clear();
+    auto const sourceZ_{source.getCoordinates().getZ()};
+    auto const checkMaxZ_{(maxHeight_ - sourceZ_) / 1_m};
+    auto ri_source{1.};
+    auto intRerfr_source{1.};
+    auto ri_destination{1.};
+    auto intRerfr_destination{1.};
+    auto height_{1.};
+
+    double const sourceHeight_{(sourceZ_ - heightTable_.front()) * inverseStep_};
+    double const destinationHeight_{(destination.getCoordinates().getZ() - heightTable_.front()) * inverseStep_};
+
+    if ((sourceHeight_ + 0.5) >= lastElement_) { // source particle is above maximum height
+      ri_source = refractivityTable_.back() + slopeRefrUpper_ * std::abs((sourceZ_ - heightTable_.back()) * inverseStep_) + 1;
+      intRerfr_source = integratedRefractivityTable_.back() + slopeIntRefrUpper_ * std::abs(checkMaxZ_);
+      rindex.push_back(ri_source);
+      points.push_back(source);
+
+      // add the refractive index of last point 'destination' and store it.
+      std::size_t const indexDestination_{static_cast<std::size_t>(destinationHeight_ + 0.5)};
+      ri_destination = refractivityTable_.at(indexDestination_) + 1;
+      intRerfr_destination = integratedRefractivityTable_.at(indexDestination_);
+      rindex.push_back(ri_destination);
+      points.push_back(destination);
+
+      height_ = sourceHeight_ - destinationHeight_;
+
+    } else if ((sourceHeight_ + 0.5 < lastElement_) && sourceHeight_ > 0) { // source particle in the table
+      // get and store the refractive index of the first point 'source'.
+      std::size_t const indexSource_{static_cast<std::size_t>(sourceHeight_ + 0.5)};
+      ri_source = refractivityTable_.at(indexSource_) + 1;
+      intRerfr_source = integratedRefractivityTable_.at(indexSource_);
+      rindex.push_back(ri_source);
+      points.push_back(source);
+
+      // add the refractive index of last point 'destination' and store it.
+      std::size_t const indexDestination_{static_cast<std::size_t>(destinationHeight_ + 0.5)};
+      ri_destination = refractivityTable_.at(indexDestination_) + 1;
+      intRerfr_destination = integratedRefractivityTable_.at(indexDestination_);
+      rindex.push_back(ri_destination);
+      points.push_back(destination);
+
+      height_ = (heightTable_.at(indexSource_) - heightTable_.at(indexDestination_)) / 1_m;
+      if (height_ == 0) {
+        height_ = 1.;
+      }
+    } else if (sourceHeight_ == 0) { // source particle is exactly at the lower edge of the table
+      // get and store the refractive index of the first point 'source'.
+      std::size_t const indexSource_{static_cast<std::size_t>(sourceHeight_ + 0.5)};
+      auto const ri_source{refractivityTable_.at(indexSource_) + 1};
+      intRerfr_source = integratedRefractivityTable_.at(indexSource_);
+      rindex.push_back(ri_source);
+      points.push_back(source);
+
+      // add the refractive index of last point 'destination' and store it.
+      std::size_t const indexDestination_{static_cast<std::size_t>(destinationHeight_ + 0.5)};
+      auto const ri_destination{refractivityTable_.at(indexDestination_) + 1};
+      intRerfr_destination = integratedRefractivityTable_.at(indexDestination_);
+      rindex.push_back(ri_destination);
+      points.push_back(destination);
+
+      height_ = destinationHeight_ - sourceHeight_;
+    } else if (sourceHeight_ < 0) { // source particle is in the ground.
+      // get and store the refractive index of the first point 'source'.
+      ri_source = refractivityTable_.front() + slopeRefrLower_ * std::abs(sourceHeight_) + 1;
+      intRerfr_source = integratedRefractivityTable_.front() + slopeIntRefrLower_ * std::abs(sourceHeight_);
+      rindex.push_back(ri_source);
+      points.push_back(source);
+
+      // add the refractive index of last point 'destination' and store it.
+      std::size_t const indexDestination_{static_cast<std::size_t>(destinationHeight_ + 0.5)};
+      auto const ri_destination{refractivityTable_.at(indexDestination_) + 1};
+      intRerfr_destination = integratedRefractivityTable_.at(indexDestination_);
+      rindex.push_back(ri_destination);
+      points.push_back(destination);
+
+      height_ = destinationHeight_ - std::abs(sourceHeight_);
+    }
 
-    // get and store the refractive index of the first point 'source'.
-    std::size_t const indexSource_{static_cast<std::size_t>(
-        (source.getCoordinates().getZ() - heightTable_.front()) * inverseStep_ +
-        0.5)}; // ToDo: this does no interpolation for particles in ground, it just stops
-               // on the surface.
-    auto const ri_source{refractivityTable_.at(indexSource_) + 1};
-    rindex.push_back(ri_source);
-    points.push_back(source);
-
-    // add the refractive index of last point 'destination' and store it.
-    std::size_t const indexDestination_{static_cast<std::size_t>(
-        (destination.getCoordinates().getZ() - heightTable_.front()) * inverseStep_ +
-        0.5)};
-    auto const ri_destination{refractivityTable_.at(indexDestination_) + 1};
-    rindex.push_back(ri_destination);
-    points.push_back(destination);
-
-    auto const height_ =
-        (heightTable_.at(indexSource_) - heightTable_.at(indexDestination_)) / 1_m;
     // compute the average refractive index.
     auto const averageRefractiveIndex_ = (ri_source + ri_destination) * 0.5;
 
     // compute the total time delay.
-    TimeType const time = (1 + (integratedRefractivityTable_.at(indexSource_) -
-                                integratedRefractivityTable_.at(indexDestination_)) /
-                                   height_) *
+    TimeType const time = (1 + (intRerfr_source - intRerfr_destination) / height_) *
                           (distance_ / constants::c);
 
     return {SignalPath(time, averageRefractiveIndex_, ri_source, ri_destination, emit_,
@@ -113,3 +177,4 @@ namespace corsika {
   } // END: propagate()
 
 } // namespace corsika
+

corsika/modules/radio/propagators/TabulatedFlatAtmospherePropagator.hpp

 /*
- * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ * (c) Copyright 2023 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
@@ -20,7 +20,8 @@ namespace corsika {
    * This class implements a tabulated propagator that approximates
    * the Earth's atmosphere as flat. Signal propagation is rectilinear
    * and this is intended to be used for vertical showers
-   * (<60 degrees zenith angle) for fast simulations.
+   * (<60 degrees zenith angle) for fast simulations. The table needs
+   * to have at least 10 elements.
    *
    */
   template <typename TEnvironment>
@@ -50,24 +51,29 @@ namespace corsika {
     SignalPathCollection propagate(Particle const& particle, Point const& source, Point const& destination);
 
   private:
-    Point const upperLimit_;
-    Point const lowerLimit_;
-    LengthType const step_;
-    InverseLengthType const inverseStep_;
-    std::vector<double> refractivityTable_;
-    std::vector<double> integratedRefractivityTable_;
-    std::vector<LengthType> heightTable_;
-    std::deque<Point> points;
-    std::vector<double> rindex;
+    Point const upperLimit_; ///< the upper point of the table.
+    Point const lowerLimit_; ///< the lowest point of the table (ideally earth's surface).
+    LengthType const minHeight_; ///< z coordinate of lower limit (minimum height) - 1_km for safety reasons.
+    LengthType const maxHeight_; ///< z coordinate of upper limit (maximum height).
+    LengthType const step_; ///< the tabulation step.
+    InverseLengthType const inverseStep_; ///< inverse of the step used to speed up calculations.
+    std::vector<double> refractivityTable_; ///< the table that stores refractivity.
+    std::vector<double> integratedRefractivityTable_; ///< the table that stores integrated refractivity.
+    std::vector<LengthType> heightTable_; ///< the table that stores the height using the step above.
+    std::deque<Point> points; ///< the points that the signal has propagated through.
+    std::vector<double> rindex; ///< the refractive index values along the signal propagation path.
+    double slopeRefrLower_; ///< used to interpolate refractivity for particles below the lowest point.
+    double slopeIntRefrLower_; ///< used to interpolate integrated refractivity for particles below the lowest point.
+    double slopeRefrUpper_; ///< used to interpolate refractivity for particles above the highest point.
+    double slopeIntRefrUpper_; ///< used to interpolate integrated refractivity for particles above the highest point.
+    double lastElement_ ; ///< last index of tables.
 
   }; // End: FlatEarthPropagator
 
   template <typename TEnvironment>
   TabulatedFlatAtmospherePropagator<TEnvironment>
-  make_tabulated_flat_atmosphere_radio_propagator(TEnvironment const& env,
-                                                  Point const& upperLimit,
-                                                  Point const& lowerLimit,
-                                                  LengthType const step) {
+  make_tabulated_flat_atmosphere_radio_propagator(TEnvironment const& env, Point const& upperLimit, Point const& lowerLimit,
+                                                  LengthType const step){
     return TabulatedFlatAtmospherePropagator<TEnvironment>(env, upperLimit, lowerLimit,
                                                            step);
   }