diff --git a/corsika/detail/modules/radio/propagators/TabulatedFlatAtmospherePropagator.inl b/corsika/detail/modules/radio/propagators/TabulatedFlatAtmospherePropagator.inl
index f80550dea62444bb3be572d0f3bec975a9e02ea5..0c213daf91e19d36fa201253ed7a83534b91a3a5 100644
--- a/corsika/detail/modules/radio/propagators/TabulatedFlatAtmospherePropagator.inl
+++ b/corsika/detail/modules/radio/propagators/TabulatedFlatAtmospherePropagator.inl
@@ -1,5 +1,5 @@
/*
- * (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
+
diff --git a/corsika/modules/radio/propagators/TabulatedFlatAtmospherePropagator.hpp b/corsika/modules/radio/propagators/TabulatedFlatAtmospherePropagator.hpp
index 4675daee3c0f43a1fe12581170ec2263eb446092..a9d78cc6566eadc0c51d74909985fd4dc2b48c0e 100644
--- a/corsika/modules/radio/propagators/TabulatedFlatAtmospherePropagator.hpp
+++ b/corsika/modules/radio/propagators/TabulatedFlatAtmospherePropagator.hpp
@@ -1,5 +1,5 @@
/*
- * (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);
}