diff --git a/corsika/detail/modules/radio/CoREAS.inl b/corsika/detail/modules/radio/CoREAS.inl
index 9f25ae97265f97cddf955f19ad2eee0c199c01b8..f79a0f8d77f0d23a3e78a4165bbd3ef2485f06d3 100644
--- a/corsika/detail/modules/radio/CoREAS.inl
+++ b/corsika/detail/modules/radio/CoREAS.inl
@@ -23,11 +23,9 @@ namespace corsika {
// should use something similar to fCoreHitTime (?)
auto const endTime_{step.getTimePost()}; // time at end point of track.
- // LCOV_EXCL_START
if (startTime_ == endTime_) {
CORSIKA_LOG_ERROR("Time at the start and end of the track coincides! - radio");
return ProcessReturn::Ok;
- // LCOV_EXCL_STOP
} else {
// get start and end position of the track
@@ -66,281 +64,191 @@ namespace corsika {
// get the SignalPathCollection (path2) from the end "endpoint" to the antenna.
auto paths2{this->propagator_.propagate(endPoint_, antenna.getLocation(), 1_m)};
- // throw an exception if path sizes don't match
- try {
- // loop over both paths at once and directly compare 'start' and 'end'
- // attributes
- for (size_t i = (paths1.size() == paths2.size()) ? 0 : throw(i = paths1.size());
- (i < paths1.size() && i < paths2.size()); i++) {
-
- // calculate preDoppler factor
- double preDoppler_{1.0 - paths1[i].refractive_index_source_ *
- beta_.dot(paths1[i].emit_)};
-
- // check if preDoppler has become zero in case of refractive index of unity
- // because of numerical limitations here you might need std::fabs(preDoppler)
- // in the if statement - same with post & mid
- // LCOV_EXCL_START
- if (preDoppler_ == 0) {
- CORSIKA_LOG_WARN("preDoppler factor numerically zero in COREAS");
- // redo calculation with higher precision
- auto const& beta_components_{beta_.getComponents(cs_)};
- auto const& emit_components_{paths1[i].emit_.getComponents(cs_)};
- long double const indexL_{paths1[i].refractive_index_source_};
- long double const betaX_{static_cast<double>(beta_components_.getX())};
- long double const betaY_{static_cast<double>(beta_components_.getY())};
- long double const betaZ_{static_cast<double>(beta_components_.getZ())};
- long double const startX_{static_cast<double>(emit_components_.getX())};
- long double const startY_{static_cast<double>(emit_components_.getY())};
- long double const startZ_{static_cast<double>(emit_components_.getZ())};
- long double const doppler =
- 1.0l -
- indexL_ * (betaX_ * startX_ + betaY_ * startY_ + betaZ_ * startZ_);
- preDoppler_ = doppler;
- }
- // LCOV_EXCL_STOP
-
- // calculate postDoppler factor
- double postDoppler_{1.0 - paths2[i].refractive_index_source_ *
- beta_.dot(paths2[i].emit_)};
-
- // check if postDoppler has become zero in case of refractive index of unity
- // because of numerical limitations
- // LCOV_EXCL_START
- if (postDoppler_ == 0) {
- CORSIKA_LOG_WARN("postDoppler factor numerically zero in CoREAS");
- // redo calculation with higher precision
- auto const& beta_components_{beta_.getComponents(cs_)};
- auto const& emit_components_{paths2[i].emit_.getComponents(cs_)};
- long double const indexL_{paths2[i].refractive_index_source_};
- long double const betaX_{static_cast<double>(beta_components_.getX())};
- long double const betaY_{static_cast<double>(beta_components_.getY())};
- long double const betaZ_{static_cast<double>(beta_components_.getZ())};
- long double const endX_{static_cast<double>(emit_components_.getX())};
- long double const endY_{static_cast<double>(emit_components_.getY())};
- long double const endZ_{static_cast<double>(emit_components_.getZ())};
- long double const doppler =
- 1.0l - indexL_ * (betaX_ * endX_ + betaY_ * endY_ + betaZ_ * endZ_);
- postDoppler_ = doppler;
- }
- // LCOV_EXCL_STOP
-
- // calculate receive time for startpoint (aka time delay)
- auto startPointReceiveTime_{
- startTime_ +
- paths1[i].propagation_time_}; // TODO: time 0 is when the imaginary
- // primary hits the ground
-
- // calculate receive time for endpoint
- auto endPointReceiveTime_{endTime_ + paths2[i].propagation_time_};
-
- // get unit vector for startpoint at antenna location
- auto ReceiveVectorStart_{paths1[i].receive_};
-
- // get unit vector for endpoint at antenna location
- auto ReceiveVectorEnd_{paths2[i].receive_};
-
- // perform ZHS-like calculation close to Cherenkov angle and for refractive
- // index at antenna location greater than 1
- if ((paths1[i].refractive_index_destination_ > 1) &&
- ((std::fabs(preDoppler_) < approxThreshold_) ||
- (std::fabs(postDoppler_) < approxThreshold_))) {
- CORSIKA_LOG_DEBUG("Used ZHS-like approximation in CoREAS - radio");
-
- // clear the existing paths for this particle and track, since we don't need
- // them anymore
- paths1.clear();
- paths2.clear();
-
- auto const halfVector_{(startPoint_ - endPoint_) * 0.5};
- auto const midPoint_{endPoint_ + halfVector_};
-
- // get global simulation time for the middle point of that track.
- TimeType const midTime_{(startTime_ + endTime_) * 0.5};
-
- // get the SignalPathCollection (path3) from the middle "endpoint" to the
- // antenna.
- auto paths3{
- this->propagator_.propagate(midPoint_, antenna.getLocation(), 1_m)};
-
- // now loop over the paths for endpoint that we got above
- for (auto const& path : paths3) {
-
- auto const midPointReceiveTime_{midTime_ + path.propagation_time_};
- double midDoppler_{1.0 -
- path.refractive_index_source_ * beta_.dot(path.emit_)};
-
- // check if midDoppler has become zero because of numerical limitations
- // LCOV_EXCL_START
- if (midDoppler_ == 0) {
- CORSIKA_LOG_WARN("midDoppler factor numerically zero in COREAS");
- // redo calculation with higher precision
- auto const& beta_components_{beta_.getComponents(cs_)};
- auto const& emit_components_{path.emit_.getComponents(cs_)};
- long double const indexL_{path.refractive_index_source_};
- long double const betaX_{static_cast<double>(beta_components_.getX())};
- long double const betaY_{static_cast<double>(beta_components_.getY())};
- long double const betaZ_{static_cast<double>(beta_components_.getZ())};
- long double const midX_{static_cast<double>(emit_components_.getX())};
- long double const midY_{static_cast<double>(emit_components_.getY())};
- long double const midZ_{static_cast<double>(emit_components_.getZ())};
- long double const doppler =
- 1.0l - indexL_ * (betaX_ * midX_ + betaY_ * midY_ + betaZ_ * midZ_);
- midDoppler_ = doppler;
- }
- // LCOV_EXCL_STOP
-
- // change the values of the receive unit vectors of start and end
- ReceiveVectorStart_ = path.receive_;
- ReceiveVectorEnd_ = path.receive_;
-
- // CoREAS calculation -> get ElectricFieldVector for "midPoint"
- ElectricFieldVector EVmid_ = (path.emit_.cross(path.emit_.cross(beta_))) /
- midDoppler_ / path.R_distance_ * constants_ *
- antenna.getSampleRate();
-
- ElectricFieldVector EV1_{EVmid_};
- ElectricFieldVector EV2_{EVmid_ * (-1.0)};
-
- TimeType deltaT_{tracklength_ / (constants::c * corrBetaValue) *
- std::fabs(midDoppler_)}; // TODO: Caution with this!
-
- if (startPointReceiveTime_ <
- endPointReceiveTime_) // EVstart_ arrives earlier
- {
- startPointReceiveTime_ = midPointReceiveTime_ - 0.5 * deltaT_;
- endPointReceiveTime_ = midPointReceiveTime_ + 0.5 * deltaT_;
- } else // EVend_ arrives earlier
- {
- startPointReceiveTime_ = midPointReceiveTime_ + 0.5 * deltaT_;
- endPointReceiveTime_ = midPointReceiveTime_ - 0.5 * deltaT_;
- }
-
- TimeType const gridResolution_{1 / antenna.getSampleRate()};
- deltaT_ = endPointReceiveTime_ - startPointReceiveTime_;
-
- // redistribute contributions over time scale defined by the observation
- // time resolution
- if (abs(deltaT_) < (gridResolution_)) {
-
- EV1_ *= std::fabs((deltaT_ / gridResolution_));
- EV2_ *= std::fabs((deltaT_ / gridResolution_));
-
- // ToDO: be careful with times in C8!!! where is the zero (time). Is it
- // close-by?
- long const startBin = static_cast<long>(
- std::floor(startPointReceiveTime_ / gridResolution_ + 0.5l));
- long const endBin = static_cast<long>(
- std::floor(endPointReceiveTime_ / gridResolution_ + 0.5l));
- double const startBinFraction =
- (startPointReceiveTime_ / gridResolution_) -
- std::floor(startPointReceiveTime_ / gridResolution_);
- double const endBinFraction =
- (endPointReceiveTime_ / gridResolution_) -
- std::floor(endPointReceiveTime_ / gridResolution_);
-
- // only do timing modification if contributions would land in same bin
- if (startBin == endBin) {
-
- // if startE arrives before endE
- if ((deltaT_) >= 0_s) {
- if ((startBinFraction >= 0.5) &&
- (endBinFraction >= 0.5)) // both points left of bin center
- {
- startPointReceiveTime_ -=
- gridResolution_; // shift EV1_ to previous gridpoint
- } else if ((startBinFraction < 0.5) &&
- (endBinFraction <
- 0.5)) // both points right of bin center
- {
- endPointReceiveTime_ +=
- gridResolution_; // shift EV2_ to next gridpoint
- } else // points on both sides of bin center
- {
- double const leftDist = 1.0 - startBinFraction;
- double const rightDist = endBinFraction;
- // check if asymmetry to right or left
- if (rightDist >= leftDist) {
- endPointReceiveTime_ +=
- gridResolution_; // shift EV2_ to next gridpoint
- } else {
- startPointReceiveTime_ -=
- gridResolution_; // shift EV1_ to previous gridpoint
- }
- }
- } else // if endE arrives before startE
- {
- if ((startBinFraction >= 0.5) &&
- (endBinFraction >= 0.5)) // both points left of bin center
- {
- endPointReceiveTime_ -=
- gridResolution_; // shift EV2_ to previous gridpoint
- } else if ((startBinFraction < 0.5) &&
- (endBinFraction <
- 0.5)) // both points right of bin center
- {
- startPointReceiveTime_ +=
- gridResolution_; // shift EV1_ to next gridpoint
- } else // points on both sides of bin center
- {
- double const leftDist = 1.0 - endBinFraction;
- double const rightDist = startBinFraction;
- // check if asymmetry to right or left
- if (rightDist >= leftDist) {
- startPointReceiveTime_ +=
- gridResolution_; // shift EV1_ to next gridpoint
- } else {
- endPointReceiveTime_ -=
- gridResolution_; // shift EV2_ to previous gridpoint
- }
- }
- } // End of else statement
- } // End of if for startbin == endbin
- } // End of if deltaT < gridresolution
-
- // TODO: Be very careful with this. Maybe the EVs should be fed after the
- // for loop of paths3
- antenna.receive(startPointReceiveTime_, ReceiveVectorStart_, EV1_);
- antenna.receive(endPointReceiveTime_, ReceiveVectorEnd_, EV2_);
- } // End of looping over paths3
-
- } // end of ZHS-like approximation
- else {
-
- // calculate electric field vector for startpoint
- ElectricFieldVector EV1_ =
- (paths1[i].emit_.cross(paths1[i].emit_.cross(beta_))) / preDoppler_ /
- paths1[i].R_distance_ * constants_ * antenna.getSampleRate();
-
- // calculate electric field vector for endpoint
- ElectricFieldVector EV2_ =
- (paths2[i].emit_.cross(paths2[i].emit_.cross(beta_))) / postDoppler_ /
- paths2[i].R_distance_ * constants_ * (-1.0) * antenna.getSampleRate();
-
- if ((preDoppler_ < 1.e-9) || (postDoppler_ < 1.e-9)) {
-
- TimeType const gridResolution_{1 / antenna.getSampleRate()};
- TimeType deltaT_{endPointReceiveTime_ - startPointReceiveTime_};
-
- if (abs(deltaT_) < (gridResolution_)) {
-
- EV1_ *= std::fabs(deltaT_ / gridResolution_); // Todo: rename EV1 and 2
- EV2_ *= std::fabs(deltaT_ / gridResolution_);
-
- long const startBin = static_cast<long>(
- std::floor(startPointReceiveTime_ / gridResolution_ + 0.5l));
- long const endBin = static_cast<long>(
- std::floor(endPointReceiveTime_ / gridResolution_ + 0.5l));
- double const startBinFraction =
- (startPointReceiveTime_ / gridResolution_) -
- std::floor(startPointReceiveTime_ / gridResolution_);
- double const endBinFraction =
- (endPointReceiveTime_ / gridResolution_) -
- std::floor(endPointReceiveTime_ / gridResolution_);
-
- // only do timing modification if contributions would land in same bin
- if (startBin == endBin) {
+ // LCOV_EXCL_START
+ // This should never happen unless someone implements a bad propagator
+ // good to catch this, hard to test without writing a custom and broken propagator
+ if (paths1.size() != paths2.size()) {
+ CORSIKA_LOG_CRITICAL(
+ "Signal Paths do not have the same size in CoREAS. Path starts: {} and "
+ "path ends: {}",
+ paths1.size(), paths2.size());
+ }
+ // LCOV_EXCL_STOP
+ // loop over both paths at once and directly compare 'start' and 'end'
+ // attributes
+ for (size_t i = 0; (i < paths1.size() && i < paths2.size()); i++) {
+
+ // calculate preDoppler factor
+ double preDoppler_{1.0 - paths1[i].refractive_index_source_ *
+ beta_.dot(paths1[i].emit_)};
+
+ // check if preDoppler has become zero in case of refractive index of unity
+ // because of numerical limitations here you might need std::fabs(preDoppler)
+ // in the if statement - same with post & mid
+ // LCOV_EXCL_START
+ if (preDoppler_ == 0) {
+ CORSIKA_LOG_WARN("preDoppler factor numerically zero in COREAS");
+ // redo calculation with higher precision
+ auto const& beta_components_{beta_.getComponents(cs_)};
+ auto const& emit_components_{paths1[i].emit_.getComponents(cs_)};
+ long double const indexL_{paths1[i].refractive_index_source_};
+ long double const betaX_{static_cast<double>(beta_components_.getX())};
+ long double const betaY_{static_cast<double>(beta_components_.getY())};
+ long double const betaZ_{static_cast<double>(beta_components_.getZ())};
+ long double const startX_{static_cast<double>(emit_components_.getX())};
+ long double const startY_{static_cast<double>(emit_components_.getY())};
+ long double const startZ_{static_cast<double>(emit_components_.getZ())};
+ long double const doppler =
+ 1.0l - indexL_ * (betaX_ * startX_ + betaY_ * startY_ + betaZ_ * startZ_);
+ preDoppler_ = doppler;
+ }
+ // LCOV_EXCL_STOP
+
+ // calculate postDoppler factor
+ double postDoppler_{1.0 - paths2[i].refractive_index_source_ *
+ beta_.dot(paths2[i].emit_)};
+
+ // check if postDoppler has become zero in case of refractive index of unity
+ // because of numerical limitations
+ // LCOV_EXCL_START
+ if (postDoppler_ == 0) {
+ CORSIKA_LOG_WARN("postDoppler factor numerically zero in CoREAS");
+ // redo calculation with higher precision
+ auto const& beta_components_{beta_.getComponents(cs_)};
+ auto const& emit_components_{paths2[i].emit_.getComponents(cs_)};
+ long double const indexL_{paths2[i].refractive_index_source_};
+ long double const betaX_{static_cast<double>(beta_components_.getX())};
+ long double const betaY_{static_cast<double>(beta_components_.getY())};
+ long double const betaZ_{static_cast<double>(beta_components_.getZ())};
+ long double const endX_{static_cast<double>(emit_components_.getX())};
+ long double const endY_{static_cast<double>(emit_components_.getY())};
+ long double const endZ_{static_cast<double>(emit_components_.getZ())};
+ long double const doppler =
+ 1.0l - indexL_ * (betaX_ * endX_ + betaY_ * endY_ + betaZ_ * endZ_);
+ postDoppler_ = doppler;
+ }
+ // LCOV_EXCL_STOP
+
+ // calculate receive time for startpoint (aka time delay)
+ auto startPointReceiveTime_{
+ startTime_ +
+ paths1[i].propagation_time_}; // TODO: time 0 is when the imaginary
+ // primary hits the ground
+
+ // calculate receive time for endpoint
+ auto endPointReceiveTime_{endTime_ + paths2[i].propagation_time_};
+
+ // get unit vector for startpoint at antenna location
+ auto ReceiveVectorStart_{paths1[i].receive_};
+
+ // get unit vector for endpoint at antenna location
+ auto ReceiveVectorEnd_{paths2[i].receive_};
+
+ // perform ZHS-like calculation close to Cherenkov angle and for refractive
+ // index at antenna location greater than 1
+ if ((paths1[i].refractive_index_destination_ > 1) &&
+ ((std::fabs(preDoppler_) < approxThreshold_) ||
+ (std::fabs(postDoppler_) < approxThreshold_))) {
+ CORSIKA_LOG_DEBUG("Used ZHS-like approximation in CoREAS - radio");
+
+ // clear the existing paths for this particle and track, since we don't need
+ // them anymore
+ paths1.clear();
+ paths2.clear();
+
+ auto const halfVector_{(startPoint_ - endPoint_) * 0.5};
+ auto const midPoint_{endPoint_ + halfVector_};
+
+ // get global simulation time for the middle point of that track.
+ TimeType const midTime_{(startTime_ + endTime_) * 0.5};
+
+ // get the SignalPathCollection (path3) from the middle "endpoint" to the
+ // antenna.
+ auto paths3{
+ this->propagator_.propagate(midPoint_, antenna.getLocation(), 1_m)};
+
+ // now loop over the paths for endpoint that we got above
+ for (auto const& path : paths3) {
+
+ auto const midPointReceiveTime_{midTime_ + path.propagation_time_};
+ double midDoppler_{1.0 -
+ path.refractive_index_source_ * beta_.dot(path.emit_)};
+
+ // check if midDoppler has become zero because of numerical limitations
+ // LCOV_EXCL_START
+ if (midDoppler_ == 0) {
+ CORSIKA_LOG_WARN("midDoppler factor numerically zero in COREAS");
+ // redo calculation with higher precision
+ auto const& beta_components_{beta_.getComponents(cs_)};
+ auto const& emit_components_{path.emit_.getComponents(cs_)};
+ long double const indexL_{path.refractive_index_source_};
+ long double const betaX_{static_cast<double>(beta_components_.getX())};
+ long double const betaY_{static_cast<double>(beta_components_.getY())};
+ long double const betaZ_{static_cast<double>(beta_components_.getZ())};
+ long double const midX_{static_cast<double>(emit_components_.getX())};
+ long double const midY_{static_cast<double>(emit_components_.getY())};
+ long double const midZ_{static_cast<double>(emit_components_.getZ())};
+ long double const doppler =
+ 1.0l - indexL_ * (betaX_ * midX_ + betaY_ * midY_ + betaZ_ * midZ_);
+ midDoppler_ = doppler;
+ }
+ // LCOV_EXCL_STOP
+
+ // change the values of the receive unit vectors of start and end
+ ReceiveVectorStart_ = path.receive_;
+ ReceiveVectorEnd_ = path.receive_;
+
+ // CoREAS calculation -> get ElectricFieldVector for "midPoint"
+ ElectricFieldVector EVmid_ = (path.emit_.cross(path.emit_.cross(beta_))) /
+ midDoppler_ / path.R_distance_ * constants_ *
+ antenna.getSampleRate();
+
+ ElectricFieldVector EV1_{EVmid_};
+ ElectricFieldVector EV2_{EVmid_ * (-1.0)};
+
+ TimeType deltaT_{tracklength_ / (constants::c * corrBetaValue) *
+ std::fabs(midDoppler_)}; // TODO: Caution with this!
+
+ if (startPointReceiveTime_ <
+ endPointReceiveTime_) // EVstart_ arrives earlier
+ {
+ startPointReceiveTime_ = midPointReceiveTime_ - 0.5 * deltaT_;
+ endPointReceiveTime_ = midPointReceiveTime_ + 0.5 * deltaT_;
+ } else // EVend_ arrives earlier
+ {
+ startPointReceiveTime_ = midPointReceiveTime_ + 0.5 * deltaT_;
+ endPointReceiveTime_ = midPointReceiveTime_ - 0.5 * deltaT_;
+ }
+
+ TimeType const gridResolution_{1 / antenna.getSampleRate()};
+ deltaT_ = endPointReceiveTime_ - startPointReceiveTime_;
+
+ // redistribute contributions over time scale defined by the observation
+ // time resolution
+ if (abs(deltaT_) < (gridResolution_)) {
+
+ EV1_ *= std::fabs((deltaT_ / gridResolution_));
+ EV2_ *= std::fabs((deltaT_ / gridResolution_));
+
+ // ToDO: be careful with times in C8!!! where is the zero (time). Is it
+ // close-by?
+ long const startBin = static_cast<long>(
+ std::floor(startPointReceiveTime_ / gridResolution_ + 0.5l));
+ long const endBin = static_cast<long>(
+ std::floor(endPointReceiveTime_ / gridResolution_ + 0.5l));
+ double const startBinFraction =
+ (startPointReceiveTime_ / gridResolution_) -
+ std::floor(startPointReceiveTime_ / gridResolution_);
+ double const endBinFraction =
+ (endPointReceiveTime_ / gridResolution_) -
+ std::floor(endPointReceiveTime_ / gridResolution_);
+
+ // only do timing modification if contributions would land in same bin
+ if (startBin == endBin) {
+
+ // if startE arrives before endE
+ if ((deltaT_) >= 0_s) {
if ((startBinFraction >= 0.5) &&
(endBinFraction >= 0.5)) // both points left of bin center
{
@@ -364,22 +272,111 @@ namespace corsika {
gridResolution_; // shift EV1_ to previous gridpoint
}
}
+ } else // if endE arrives before startE
+ {
+ if ((startBinFraction >= 0.5) &&
+ (endBinFraction >= 0.5)) // both points left of bin center
+ {
+ endPointReceiveTime_ -=
+ gridResolution_; // shift EV2_ to previous gridpoint
+ } else if ((startBinFraction < 0.5) &&
+ (endBinFraction < 0.5)) // both points right of bin center
+ {
+ startPointReceiveTime_ +=
+ gridResolution_; // shift EV1_ to next gridpoint
+ } else // points on both sides of bin center
+ {
+ double const leftDist = 1.0 - endBinFraction;
+ double const rightDist = startBinFraction;
+ // check if asymmetry to right or left
+ if (rightDist >= leftDist) {
+ startPointReceiveTime_ +=
+ gridResolution_; // shift EV1_ to next gridpoint
+ } else {
+ endPointReceiveTime_ -=
+ gridResolution_; // shift EV2_ to previous gridpoint
+ }
+ }
+ } // End of else statement
+ } // End of if for startbin == endbin
+ } // End of if deltaT < gridresolution
- } // End of if for startbin == endbin
- } // End of if deltaT < gridresolution
- } // End of if that checks small doppler factors
+ // TODO: Be very careful with this. Maybe the EVs should be fed after the
+ // for loop of paths3
antenna.receive(startPointReceiveTime_, ReceiveVectorStart_, EV1_);
antenna.receive(endPointReceiveTime_, ReceiveVectorEnd_, EV2_);
- } // End of else that does not perform ZHS-like approximation
+ } // End of looping over paths3
+
+ } // end of ZHS-like approximation
+ else {
+
+ // calculate electric field vector for startpoint
+ ElectricFieldVector EV1_ =
+ (paths1[i].emit_.cross(paths1[i].emit_.cross(beta_))) / preDoppler_ /
+ paths1[i].R_distance_ * constants_ * antenna.getSampleRate();
+
+ // calculate electric field vector for endpoint
+ ElectricFieldVector EV2_ =
+ (paths2[i].emit_.cross(paths2[i].emit_.cross(beta_))) / postDoppler_ /
+ paths2[i].R_distance_ * constants_ * (-1.0) * antenna.getSampleRate();
+
+ if ((preDoppler_ < 1.e-9) || (postDoppler_ < 1.e-9)) {
+
+ TimeType const gridResolution_{1 / antenna.getSampleRate()};
+ TimeType deltaT_{endPointReceiveTime_ - startPointReceiveTime_};
+
+ if (abs(deltaT_) < (gridResolution_)) {
+
+ EV1_ *= std::fabs(deltaT_ / gridResolution_); // Todo: rename EV1 and 2
+ EV2_ *= std::fabs(deltaT_ / gridResolution_);
+
+ long const startBin = static_cast<long>(
+ std::floor(startPointReceiveTime_ / gridResolution_ + 0.5l));
+ long const endBin = static_cast<long>(
+ std::floor(endPointReceiveTime_ / gridResolution_ + 0.5l));
+ double const startBinFraction =
+ (startPointReceiveTime_ / gridResolution_) -
+ std::floor(startPointReceiveTime_ / gridResolution_);
+ double const endBinFraction =
+ (endPointReceiveTime_ / gridResolution_) -
+ std::floor(endPointReceiveTime_ / gridResolution_);
+
+ // only do timing modification if contributions would land in same bin
+ if (startBin == endBin) {
+
+ if ((startBinFraction >= 0.5) &&
+ (endBinFraction >= 0.5)) // both points left of bin center
+ {
+ startPointReceiveTime_ -=
+ gridResolution_; // shift EV1_ to previous gridpoint
+ } else if ((startBinFraction < 0.5) &&
+ (endBinFraction < 0.5)) // both points right of bin center
+ {
+ endPointReceiveTime_ +=
+ gridResolution_; // shift EV2_ to next gridpoint
+ } else // points on both sides of bin center
+ {
+ double const leftDist = 1.0 - startBinFraction;
+ double const rightDist = endBinFraction;
+ // check if asymmetry to right or left
+ if (rightDist >= leftDist) {
+ endPointReceiveTime_ +=
+ gridResolution_; // shift EV2_ to next gridpoint
+ } else {
+ startPointReceiveTime_ -=
+ gridResolution_; // shift EV1_ to previous gridpoint
+ }
+ }
- } // End of loop over both paths to get signal info
- } // End of try block
- // LCOV_EXCL_START
- catch (size_t i) {
- CORSIKA_LOG_ERROR("Signal Paths do not have the same size in CoREAS");
- }
- // LCOV_EXCL_STOP
- } // End of looping over antennas
+ } // End of if for startbin == endbin
+ } // End of if deltaT < gridresolution
+ } // End of if that checks small doppler factors
+ antenna.receive(startPointReceiveTime_, ReceiveVectorStart_, EV1_);
+ antenna.receive(endPointReceiveTime_, ReceiveVectorEnd_, EV2_);
+ } // End of else that does not perform ZHS-like approximation
+
+ } // End of loop over both paths to get signal info
+ } // End of looping over antennas
return ProcessReturn::Ok;
}
diff --git a/corsika/detail/modules/radio/antennas/Antenna.inl b/corsika/detail/modules/radio/antennas/Antenna.inl
index 0b97560aaefc129779fdf1e39dfdfe8f68e9e61b..89d511257c5d9b9069df9036a243956e17d9e06c 100644
--- a/corsika/detail/modules/radio/antennas/Antenna.inl
+++ b/corsika/detail/modules/radio/antennas/Antenna.inl
@@ -16,7 +16,7 @@ namespace corsika {
CoordinateSystemPtr const& coordinateSystem)
: name_(name)
, location_(location)
- , coordinateSystem_(coordinateSystem){};
+ , coordinateSystem_(coordinateSystem) {}
template <typename TAntennaImpl>
inline Point const& Antenna<TAntennaImpl>::getLocation() const {
diff --git a/corsika/detail/modules/radio/antennas/TimeDomainAntenna.inl b/corsika/detail/modules/radio/antennas/TimeDomainAntenna.inl
index 920b303b9d63281e708ad59b278b8089a4a98cf8..c32f51657d745f2cd7f4ee0f5ace92979f98ca88 100644
--- a/corsika/detail/modules/radio/antennas/TimeDomainAntenna.inl
+++ b/corsika/detail/modules/radio/antennas/TimeDomainAntenna.inl
@@ -20,56 +20,71 @@ namespace corsika {
TimeType const ground_hit_time)
: Antenna(name, location, coordinateSystem)
, start_time_(start_time)
- , duration_(duration)
- , sample_rate_(sample_rate)
+ , duration_(std::abs(duration / 1_s) * 1_s)
+ , sample_rate_(std::abs(sample_rate / 1_Hz) * 1_Hz)
, ground_hit_time_(ground_hit_time)
- , num_bins_(static_cast<std::size_t>(duration * sample_rate + 1.5l))
+ , num_bins_(static_cast<std::size_t>(duration_ * sample_rate_ + 1.5l))
, waveformEX_(num_bins_, 0)
, waveformEY_(num_bins_, 0)
, waveformEZ_(num_bins_, 0)
- , time_axis_(createTimeAxis()){};
+ , time_axis_(createTimeAxis()) {
+ if (0_s == duration_) {
+ CORSIKA_LOG_WARN(
+ "Antenna: \"{}\" has a duration of zero. Nothing will be injected into it",
+ name);
+ } else if (duration_ != duration) {
+ CORSIKA_LOG_WARN(
+ "Antenna: \"{}\" was given a negative duration. Set to absolute value.", name);
+ }
+
+ if (sample_rate_ != sample_rate) {
+ CORSIKA_LOG_WARN(
+ "Antenna: \"{}\" was given a negative sampling rate. Set to absolute value.",
+ name);
+ }
+ }
- inline void TimeDomainAntenna::receive(const TimeType time,
- const Vector<dimensionless_d>& receive_vector,
- const ElectricFieldVector& efield) {
+ inline void TimeDomainAntenna::receive(
+ const TimeType time, [[maybe_unused]] const Vector<dimensionless_d>& receive_vector,
+ const ElectricFieldVector& efield) {
if (time < start_time_ || time > (start_time_ + duration_)) {
return;
} else {
// figure out the correct timebin to store the E-field value.
// NOTE: static cast is implicitly flooring
- auto timebin_{static_cast<std::size_t>(
+ auto timebin{static_cast<std::size_t>(
std::floor((time - start_time_) * sample_rate_ + 0.5l))};
// store the x,y,z electric field components.
auto const& Electric_field_components{efield.getComponents(coordinateSystem_)};
- waveformEX_.at(timebin_) += (Electric_field_components.getX() * (1_m / 1_V));
- waveformEY_.at(timebin_) += (Electric_field_components.getY() * (1_m / 1_V));
- waveformEZ_.at(timebin_) += (Electric_field_components.getZ() * (1_m / 1_V));
+ waveformEX_.at(timebin) += (Electric_field_components.getX() * (1_m / 1_V));
+ waveformEY_.at(timebin) += (Electric_field_components.getY() * (1_m / 1_V));
+ waveformEZ_.at(timebin) += (Electric_field_components.getZ() * (1_m / 1_V));
// TODO: Check how they are stored in memory, row-wise or column-wise? Probably use
// a 3D object
}
}
- inline void TimeDomainAntenna::receive(const TimeType time,
- const Vector<dimensionless_d>& receive_vector,
- const VectorPotential& vectorP) {
+ inline void TimeDomainAntenna::receive(
+ const TimeType time, [[maybe_unused]] const Vector<dimensionless_d>& receive_vector,
+ const VectorPotential& vectorP) {
if (time < start_time_ || time > (start_time_ + duration_)) {
return;
} else {
// figure out the correct timebin to store the E-field value.
// NOTE: static cast is implicitly flooring
- auto timebin_{static_cast<std::size_t>(
+ auto timebin{static_cast<std::size_t>(
std::floor((time - start_time_) * sample_rate_ + 0.5l))};
// store the x,y,z electric field components.
auto const& Vector_potential_components{vectorP.getComponents(coordinateSystem_)};
- waveformEX_.at(timebin_) +=
+ waveformEX_.at(timebin) +=
(Vector_potential_components.getX() * (1_m / (1_V * 1_s)));
- waveformEY_.at(timebin_) +=
+ waveformEY_.at(timebin) +=
(Vector_potential_components.getY() * (1_m / (1_V * 1_s)));
- waveformEZ_.at(timebin_) +=
+ waveformEZ_.at(timebin) +=
(Vector_potential_components.getZ() * (1_m / (1_V * 1_s)));
// TODO: Check how they are stored in memory, row-wise or column-wise? Probably use
// a 3D object
@@ -91,7 +106,7 @@ namespace corsika {
auto sample_period{1 / sample_rate_};
// fill in every time-value
- for (std::size_t i = 0; i < num_bins_; i++) {
+ for (uint64_t i = 0; i < num_bins_; i++) {
// create the current time in nanoseconds
times.at(i) = static_cast<long double>(
((start_time_ - ground_hit_time_) + i * sample_period) / 1_ns);
diff --git a/corsika/detail/modules/radio/propagators/SimplePropagator.inl b/corsika/detail/modules/radio/propagators/SimplePropagator.inl
index 4101878983cb7f3060e0743392a08caca9689508..36cf165e155cf85012ff837f2f5354e8775da4c7 100644
--- a/corsika/detail/modules/radio/propagators/SimplePropagator.inl
+++ b/corsika/detail/modules/radio/propagators/SimplePropagator.inl
@@ -13,12 +13,13 @@ namespace corsika {
template <typename TEnvironment>
inline SimplePropagator<TEnvironment>::SimplePropagator(TEnvironment const& env)
- : RadioPropagator<SimplePropagator, TEnvironment>(env){};
+ : RadioPropagator<SimplePropagator, TEnvironment>(env) {}
template <typename TEnvironment>
inline typename SimplePropagator<TEnvironment>::SignalPathCollection
- SimplePropagator<TEnvironment>::propagate(Point const& source, Point const& destination,
- LengthType const stepsize) const {
+ SimplePropagator<TEnvironment>::propagate(
+ Point const& source, Point const& destination,
+ [[maybe_unused]] LengthType const stepsize) const {
/**
* This is the simplest case of straight propagator
@@ -62,8 +63,9 @@ namespace corsika {
// compute the total time delay.
TimeType const time = averageRefractiveIndex_ * (distance_ / constants::c);
- return {SignalPath(time, averageRefractiveIndex_, ri_source, ri_destination, emit_,
- receive_, distance_, points)};
+ return std::vector<SignalPath>(
+ 1, SignalPath(time, averageRefractiveIndex_, ri_source, ri_destination, emit_,
+ receive_, distance_, points));
} // END: propagate()
diff --git a/corsika/detail/modules/radio/propagators/StraightPropagator.inl b/corsika/detail/modules/radio/propagators/StraightPropagator.inl
index 9962fc133003d855ab682aa15b12c275e225c035..59b5c2c8e6b6cd85ca1aad199a18ed0ac0fd1505 100644
--- a/corsika/detail/modules/radio/propagators/StraightPropagator.inl
+++ b/corsika/detail/modules/radio/propagators/StraightPropagator.inl
@@ -14,7 +14,7 @@ namespace corsika {
template <typename TEnvironment>
// TODO: maybe the constructor doesn't take any arguments for the environment (?)
inline StraightPropagator<TEnvironment>::StraightPropagator(TEnvironment const& env)
- : RadioPropagator<StraightPropagator, TEnvironment>(env){};
+ : RadioPropagator<StraightPropagator, TEnvironment>(env) {}
template <typename TEnvironment>
inline typename StraightPropagator<TEnvironment>::SignalPathCollection
@@ -30,17 +30,17 @@ namespace corsika {
*/
// these are used for the direction of emission and reception of signal at the antenna
- auto const emit_{(destination - source).normalized()};
- auto const receive_{-emit_};
+ auto const emit{(destination - source).normalized()};
+ auto const receive{-emit};
// the distance from the point of emission to an observer
- auto const distance_{(destination - source).getNorm()};
+ auto const distance{(destination - source).getNorm()};
try {
- if (stepsize <= 0.5 * distance_) {
+ if (stepsize <= 0.5 * distance) {
// "step" is the direction vector with length `stepsize`
- auto const step{emit_ * stepsize};
+ auto const step{emit * stepsize};
// calculate the number of points (roughly) for the numerical integration
auto const n_points{(destination - source).getNorm() / stepsize};
@@ -93,7 +93,7 @@ namespace corsika {
auto N = rindex.size();
std::size_t index = 0;
double sum = rindex.at(index);
- auto refra_ = rindex.at(index);
+ auto refra = rindex.at(index);
TimeType time{0_s};
if ((N - 1) % 2 == 0) {
@@ -102,15 +102,15 @@ namespace corsika {
for (std::size_t index = 1; index < (N - 1); index += 2) {
sum += 4 * rindex.at(index);
- refra_ += rindex.at(index);
+ refra += rindex.at(index);
}
for (std::size_t index = 2; index < (N - 1); index += 2) {
sum += 2 * rindex.at(index);
- refra_ += rindex.at(index);
+ refra += rindex.at(index);
}
index = N - 1;
sum = sum + rindex.at(index);
- refra_ += rindex.at(index);
+ refra += rindex.at(index);
// compute the total time delay.
time = sum * (h / (3 * constants::c));
@@ -133,15 +133,15 @@ namespace corsika {
auto const h = ((extrapoint2_ - source).getNorm()) / (N - 1);
for (std::size_t index = 1; index < (N - 1); index += 2) {
sum += 4 * rindex.at(index);
- refra_ += rindex.at(index);
+ refra += rindex.at(index);
}
for (std::size_t index = 2; index < (N - 1); index += 2) {
sum += 2 * rindex.at(index);
- refra_ += rindex.at(index);
+ refra += rindex.at(index);
}
index = N - 1;
sum = sum + rindex.at(index);
- refra_ += rindex.at(index);
+ refra += rindex.at(index);
// compute the total time delay including the correction
time =
@@ -150,13 +150,14 @@ namespace corsika {
}
// uncomment the following if you want to skip the integration for fast tests
- // TimeType time = ri_destination * (distance_ / constants::c);
+ // TimeType time = ri_destination * (distance / constants::c);
// compute the average refractive index.
- auto const averageRefractiveIndex_ = refra_ / N;
+ auto const averageRefractiveIndex = refra / N;
- return {SignalPath(time, averageRefractiveIndex_, ri_source, ri_destination,
- emit_, receive_, distance_, points)};
+ return std::vector<SignalPath>(
+ 1, SignalPath(time, averageRefractiveIndex, ri_source, ri_destination, emit,
+ receive, distance, points));
} else {
throw stepsize;
}
@@ -164,6 +165,11 @@ namespace corsika {
CORSIKA_LOG_ERROR("Please choose a smaller stepsize for the numerical integration");
}
- } // END: propagate()
+ std::deque<Point> const defaultPoints;
+ return std::vector<SignalPath>(
+ 1, SignalPath(0_s, 0, 0, 0, emit, receive, distance,
+ defaultPoints)); // Dummy return that is never called (combat
+ // compile warnings)
+ } // END: propagate()
} // namespace corsika
diff --git a/corsika/modules/radio/CoREAS.hpp b/corsika/modules/radio/CoREAS.hpp
index 0756a78af2b932093312fd3416780929c319379b..058126aa7fb3070d2dac8e3c6de26989ad211e75 100644
--- a/corsika/modules/radio/CoREAS.hpp
+++ b/corsika/modules/radio/CoREAS.hpp
@@ -34,7 +34,7 @@ namespace corsika {
*/
template <typename... TArgs>
CoREAS(TRadioDetector& detector, TArgs&&... args)
- : RadioProcess<TRadioDetector, CoREAS, TPropagator>(detector, args...){};
+ : RadioProcess<TRadioDetector, CoREAS, TPropagator>(detector, args...) {}
/**
* Simulate the radio emission from a particle across a track.
diff --git a/corsika/modules/radio/ZHS.hpp b/corsika/modules/radio/ZHS.hpp
index 5ec38895b8a8c75a9cb6974b996253a430de1f3e..49153ba25c5ac2707e5c58adfc98153928846c26 100644
--- a/corsika/modules/radio/ZHS.hpp
+++ b/corsika/modules/radio/ZHS.hpp
@@ -36,7 +36,7 @@ namespace corsika {
*/
template <typename... TArgs>
ZHS(TRadioDetector& detector, TArgs&&... args)
- : RadioProcess<TRadioDetector, ZHS, TPropagator>(detector, args...){};
+ : RadioProcess<TRadioDetector, ZHS, TPropagator>(detector, args...) {}
/**
* Simulate the radio emission from a particle across a track.
diff --git a/corsika/modules/radio/antennas/TimeDomainAntenna.hpp b/corsika/modules/radio/antennas/TimeDomainAntenna.hpp
index 801c529b87a035a021ab1c3aeb52cd86adf9e14d..85fd9dd100743f674fef937f0aff8c69570b5d9c 100644
--- a/corsika/modules/radio/antennas/TimeDomainAntenna.hpp
+++ b/corsika/modules/radio/antennas/TimeDomainAntenna.hpp
@@ -8,6 +8,7 @@
#pragma once
#include <corsika/modules/radio/antennas/Antenna.hpp>
+#include <yaml-cpp/yaml.h>
#include <vector>
namespace corsika {
@@ -22,11 +23,11 @@ namespace corsika {
TimeType const start_time_; ///< The start time of this waveform.
TimeType const duration_; ///< The duration of this waveform.
InverseTimeType const sample_rate_; ///< The sampling rate of this antenna.
- int const num_bins_; ///< The number of bins used.
- std::vector<double> waveformEX_; ///< EX polarization.
- std::vector<double> waveformEY_; ///< EY polarization.
- std::vector<double> waveformEZ_; ///< EZ polarization.
TimeType const ground_hit_time_; ///< The time the primary particle hits the ground.
+ uint64_t const num_bins_; ///< The number of bins used.
+ std::vector<double> waveformEX_; ///< EX polarization.
+ std::vector<double> waveformEY_; ///< EY polarization.
+ std::vector<double> waveformEZ_; ///< EZ polarization.
std::vector<long double> const
time_axis_; ///< The time axis corresponding to the electric field.
diff --git a/tests/modules/testRadio.cpp b/tests/modules/testRadio.cpp
index 2260107939952a45b0f8d7bdb122f67fc1de082d..e31434fd4be10d403b5bab43ac0d7373d9e71739 100644
--- a/tests/modules/testRadio.cpp
+++ b/tests/modules/testRadio.cpp
@@ -16,8 +16,10 @@
#include <corsika/modules/radio/propagators/SignalPath.hpp>
#include <corsika/modules/radio/propagators/RadioPropagator.hpp>
+#include <boost/filesystem.hpp>
#include <vector>
#include <istream>
+#include <filesystem>
#include <fstream>
#include <iostream>
@@ -65,60 +67,49 @@ TEST_CASE("Radio", "[processes]") {
// This serves as a compiler test for any changes in the CoREAS algorithm
// Environment
+
using EnvironmentInterface =
IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
- // using EnvType = setup::Environment;
using EnvType = Environment<EnvironmentInterface>;
+ // using EnvType = setup::Environment;
EnvType envCoREAS;
- CoordinateSystemPtr const& rootCSCoREAS = envCoREAS.getCoordinateSystem();
- Point const center{rootCSCoREAS, 0_m, 0_m, 0_m};
+ CoordinateSystemPtr const& rootCS = envCoREAS.getCoordinateSystem();
+ Point const center{rootCS, 0_m, 0_m, 0_m};
- // 1.000327,
create_5layer_atmosphere<EnvironmentInterface, MyExtraEnv>(
envCoREAS, AtmosphereId::LinsleyUSStd, center, 1.000327, Medium::AirDry1Atm,
- MagneticFieldVector{rootCSCoREAS, 0_T, 50_uT, 0_T});
+ MagneticFieldVector{rootCS, 0_T, 50_uT, 0_T});
- // now create antennas and detectors
- // the antennas location
- const auto point1{Point(envCoREAS.getCoordinateSystem(), 100_m, 2_m, 3_m)};
- const auto point2{Point(envCoREAS.getCoordinateSystem(), 4_m, 80_m, 6_m)};
- const auto point3{Point(envCoREAS.getCoordinateSystem(), 7_m, 8_m, 9_m)};
- const auto point4{Point(envCoREAS.getCoordinateSystem(), 5_m, 5_m, 10_m)};
-
- // create times for the antenna
+ // create the detector
+ const auto ant1Loc{Point(rootCS, 100_m, 2_m, 3_m)};
+ const auto ant2Loc{Point(rootCS, 4_m, 80_m, 6_m)};
const TimeType t1{0_s};
const TimeType t2{10_s};
const InverseTimeType t3{1e+3_Hz};
- const TimeType t4{11_s};
-
- // check that I can create an antenna at (1, 2, 3)
- TimeDomainAntenna ant1("antenna_name", point1, rootCSCoREAS, t1, t2, t3, t1);
- TimeDomainAntenna ant2("antenna_name2", point2, rootCSCoREAS, t1, t2, t3, t1);
-
- // construct a radio detector instance to store our antennas
+ TimeDomainAntenna ant1("antenna_name", ant1Loc, rootCS, t1, t2, t3, t1);
+ TimeDomainAntenna ant2("antenna_name2", ant2Loc, rootCS, t1, t2, t3, t1);
AntennaCollection<TimeDomainAntenna> detector;
-
- // add the antennas to the detector
detector.addAntenna(ant1);
detector.addAntenna(ant2);
- // create a particle
- const Code particle{Code::Electron};
- // const Code particle{Code::Proton};
- const auto pmass{get_mass(particle)};
+ const auto trackStart{Point(rootCS, 7_m, 8_m, 9_m)};
+ const auto trackEnd{Point(rootCS, 5_m, 5_m, 10_m)};
+
+ // create an electron
+ const Code electron{Code::Electron};
+ const auto pmass{get_mass(electron)};
- VelocityVector v0(rootCSCoREAS, {5e+2_m / second, 5e+2_m / second, 5e+2_m / second});
+ VelocityVector v0(rootCS, {5e+2_m / second, 5e+2_m / second, 5e+2_m / second});
- Vector B0(rootCSCoREAS, 5_T, 5_T, 5_T);
+ Vector B0(rootCS, 5_T, 5_T, 5_T);
- Line const line(point3, v0);
+ Line const line(trackStart, v0);
auto const k{1_m * ((1_m) / ((1_s * 1_s) * 1_V))};
auto const t = 1e-12_s;
- LeapFrogTrajectory base(point4, v0, B0, k, t);
- // std::cout << "Leap Frog Trajectory is: " << base << std::endl;
+ LeapFrogTrajectory base(trackEnd, v0, B0, k, t);
// create a new stack for each trial
setup::Stack<EnvType> stack;
@@ -128,20 +119,16 @@ TEST_CASE("Radio", "[processes]") {
// compute the necessary momentumn
const HEPMomentumType P0{sqrt(E0 * E0 - pmass * pmass)};
-
- // and create the momentum vector
- const auto plab{MomentumVector(rootCSCoREAS, {0_GeV, 0_GeV, P0})};
+ const auto plab{MomentumVector(rootCS, {0_GeV, 0_GeV, P0})};
// and create the location of the particle in this coordinate system
- const Point pos(rootCSCoREAS, 50_m, 10_m, 80_m);
+ const Point pos(rootCS, 50_m, 10_m, 80_m);
// add the particle to the stack
auto const particle1{stack.addParticle(std::make_tuple(
- particle, calculate_kinetic_energy(plab.getNorm(), get_mass(particle)),
+ electron, calculate_kinetic_energy(plab.getNorm(), get_mass(electron)),
plab.normalized(), pos, 0_ns))};
- auto const charge_{get_charge(particle1.getPID())};
-
// create a radio process instance using CoREAS
RadioProcess<decltype(detector),
CoREAS<decltype(detector), decltype(StraightPropagator(envCoREAS))>,
@@ -150,9 +137,149 @@ TEST_CASE("Radio", "[processes]") {
Step step(particle1, base);
// check doContinuous and simulate methods
- coreas.doContinuous(step, true);
+ auto const result = coreas.doContinuous(step, true);
+ REQUIRE(ProcessReturn::Ok == result);
+
+ for (auto const& ant : detector.getAntennas()) {
+ // make sure something was put into the antenna
+ auto totalX = ant.getWaveformX()[0];
+ auto totalY = ant.getWaveformY()[0];
+ auto totalZ = ant.getWaveformZ()[0];
+ for (size_t i = 0; i < ant.getWaveformX().size(); i++) {
+ totalX += ant.getWaveformX()[i];
+ totalY += ant.getWaveformY()[i];
+ totalZ += ant.getWaveformZ()[i];
+ }
+ REQUIRE((totalX + totalY + totalZ) > (totalX * 0));
+ }
+
+ //////////////////////////////////////
+ // reset everything for a new particle
+ //////////////////////////////////////
+ ant1.reset();
+ ant2.reset();
+ stack.purge();
+
+ // add the particle to the stack that is VERY late
+ auto const particle2{stack.addParticle(std::make_tuple(
+ electron, calculate_kinetic_energy(plab.getNorm(), get_mass(electron)),
+ plab.normalized(), pos, t1 + t2 * 100000))};
+ Step step2(particle2, base);
+ auto const result2 = coreas.doContinuous(step2, true);
+ REQUIRE(ProcessReturn::Ok == result2);
+ for (auto const& ant : detector.getAntennas()) {
+ // make sure something was put into the antenna
+ auto total = ant.getWaveformX()[0];
+ for (size_t i = 0; i < ant.getWaveformX().size(); i++) {
+ total += ant.getWaveformX()[i] * ant.getWaveformX()[i];
+ total += ant.getWaveformY()[i] * ant.getWaveformY()[i];
+ total += ant.getWaveformZ()[i] * ant.getWaveformZ()[i];
+ }
+ REQUIRE(total < (1e-12 * ant.getWaveformX().size()));
+ }
+
+ coreas.endOfLibrary();
+
} // END: SECTION("CoREAS process")
+ SECTION("CoREAS Edge Cases") {
+ using IModelInterface =
+ IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
+ using AtmModel = UniformRefractiveIndex<
+ MediumPropertyModel<UniformMagneticField<HomogeneousMedium<IModelInterface>>>>;
+ using EnvType = Environment<AtmModel>;
+ EnvType envCoREAS;
+ CoordinateSystemPtr const& rootCS = envCoREAS.getCoordinateSystem();
+ Point const center{rootCS, 0_m, 0_m, 0_m};
+
+ Vector B1(rootCS, 0_T, 0_T, 1_T);
+ NuclearComposition const protonComposition({Code::Proton}, {1.});
+ const double refractiveIndex{1.000327};
+ const auto density{1_g / cube(1_cm)};
+ auto Medium = EnvType::createNode<Sphere>(
+ center, 10_km * std::numeric_limits<double>::infinity());
+ auto const props = Medium->setModelProperties<AtmModel>(
+ refractiveIndex, Medium::AirDry1Atm, B1, density, protonComposition);
+ envCoREAS.getUniverse()->addChild(std::move(Medium));
+
+ // create the detector
+ const auto ant1Loc{Point(rootCS, 100_m, 2_m, 3_m)};
+ const auto ant2Loc{Point(rootCS, 4_m, 80_m, 6_m)};
+ const TimeType t1{0_s};
+ const TimeType t2{10_s};
+ const InverseTimeType t3{1e+3_Hz};
+ TimeDomainAntenna ant1("antenna_name", ant1Loc, rootCS, t1, t2, t3, t1);
+ TimeDomainAntenna ant2("antenna_name2", ant2Loc, rootCS, t1, t2, t3, t1);
+ AntennaCollection<TimeDomainAntenna> detector;
+ detector.addAntenna(ant1);
+ detector.addAntenna(ant2);
+
+ const auto trackStart{Point(rootCS, 7_m, 8_m, 9_m)};
+ const auto trackEnd{Point(rootCS, 5_m, 5_m, 10_m)};
+
+ // create an electron
+ const Code electron{Code::Electron};
+ const auto pmass{get_mass(electron)};
+
+ VelocityVector v0(rootCS, {1_m / second, 0_m / second, 0_m / second});
+
+ Vector B0(rootCS, 5_T, 5_T, 5_T);
+
+ Line const line(trackStart, v0);
+
+ // create a new stack for each trial
+ setup::Stack<EnvType> stack;
+
+ // construct an energy
+ const HEPEnergyType E0{1_TeV};
+
+ // compute the necessary momentumn
+ const HEPMomentumType P0{sqrt(E0 * E0 - pmass * pmass)};
+ const auto plab{MomentumVector(rootCS, {0_GeV, 0_GeV, P0})};
+
+ // and create the location of the particle in this coordinate system
+ const Point pos(rootCS, 50_m, 10_m, 80_m);
+
+ // add the particle to the stack
+ auto const particle1{stack.addParticle(std::make_tuple(
+ electron, calculate_kinetic_energy(plab.getNorm(), get_mass(electron)),
+ plab.normalized(), pos, 0_ns))};
+
+ // create a radio process instance using CoREAS
+ RadioProcess<decltype(detector),
+ CoREAS<decltype(detector), decltype(StraightPropagator(envCoREAS))>,
+ decltype(StraightPropagator(envCoREAS))>
+ coreas(detector, envCoREAS);
+
+ auto result = coreas.doContinuous(
+ Step(particle1, StraightTrajectory(line, 0_ns, 0_ns, v0, v0)), true);
+ REQUIRE(ProcessReturn::Ok == result);
+ result = coreas.doContinuous(
+ Step(particle1, StraightTrajectory(line, 0_ns, 1_ns, v0, v0)), true);
+ REQUIRE(ProcessReturn::Ok == result);
+ result = coreas.doContinuous(
+ Step(particle1, StraightTrajectory(line, 0_ns, -1_ns, v0, v0)), true);
+ REQUIRE(ProcessReturn::Ok == result);
+ result = coreas.doContinuous(
+ Step(particle1, StraightTrajectory(line, 1_ns, -1_ns, v0, v0)), true);
+ REQUIRE(ProcessReturn::Ok == result);
+ result = coreas.doContinuous(
+ Step(particle1, StraightTrajectory(line, -1_ns, 1_ns, v0, v0)), true);
+ REQUIRE(ProcessReturn::Ok == result);
+ result = coreas.doContinuous(
+ Step(particle1, StraightTrajectory(line, -1_ns, 1_ns, v0, -v0)), true);
+ REQUIRE(ProcessReturn::Ok == result);
+
+ // Use ZHS-like loop
+ auto const vParallel =
+ VelocityVector(rootCS, {0_m / second, 1_m / second, 0_m / second});
+ result = coreas.doContinuous(
+ Step(particle1, StraightTrajectory(Line(trackStart, vParallel), 0_ns, 1_ns,
+ vParallel, vParallel)),
+ true);
+ REQUIRE(ProcessReturn::Ok == result);
+ }
+
SECTION("ZHS process") {
// This section serves as a compiler test for any changes in the ZHS algorithm
@@ -163,9 +290,9 @@ TEST_CASE("Radio", "[processes]") {
MediumPropertyModel<UniformMagneticField<HomogeneousMedium<IModelInterface>>>>;
using EnvType = Environment<AtmModel>;
EnvType envZHS;
- CoordinateSystemPtr const& rootCSZHS = envZHS.getCoordinateSystem();
+ CoordinateSystemPtr const& rootCS = envZHS.getCoordinateSystem();
// get the center point
- Point const center{rootCSZHS, 0_m, 0_m, 0_m};
+ Point const center{rootCS, 0_m, 0_m, 0_m};
// a refractive index
const double ri_{1.000327};
@@ -176,7 +303,7 @@ TEST_CASE("Radio", "[processes]") {
NuclearComposition const protonComposition({Code::Proton}, {1.});
// create magnetic field vector
- Vector B1(rootCSZHS, 0_T, 0_T, 1_T);
+ Vector B1(rootCS, 0_T, 0_T, 1_T);
auto Medium = EnvType::createNode<Sphere>(
center, 1_km * std::numeric_limits<double>::infinity());
@@ -186,25 +313,18 @@ TEST_CASE("Radio", "[processes]") {
envZHS.getUniverse()->addChild(std::move(Medium));
// the antennas location
- const auto point1{Point(envZHS.getCoordinateSystem(), 100_m, 2_m, 3_m)};
- const auto point2{Point(envZHS.getCoordinateSystem(), 4_m, 80_m, 6_m)};
- const auto point3{Point(envZHS.getCoordinateSystem(), 7_m, 8_m, 9_m)};
- const auto point4{Point(envZHS.getCoordinateSystem(), 5_m, 5_m, 10_m)};
+ const auto trackStart{Point(envZHS.getCoordinateSystem(), 7_m, 8_m, 9_m)};
+ const auto trackEnd{Point(envZHS.getCoordinateSystem(), 5_m, 5_m, 10_m)};
- // create times for the antenna
+ // create the detector
+ const auto ant1Loc{Point(rootCS, 100_m, 2_m, 3_m)};
+ const auto ant2Loc{Point(rootCS, 4_m, 80_m, 6_m)};
const TimeType t1{0_s};
const TimeType t2{10_s};
const InverseTimeType t3{1e+3_Hz};
- const TimeType t4{11_s};
-
- // check that I can create an antenna at (1, 2, 3)
- TimeDomainAntenna ant1("antenna_zhs", point1, rootCSZHS, t1, t2, t3, t1);
- TimeDomainAntenna ant2("antenna_zhs2", point2, rootCSZHS, t1, t2, t3, t1);
-
- // construct a radio detector instance to store our antennas
+ TimeDomainAntenna ant1("antenna_name", ant1Loc, rootCS, t1, t2, t3, t1);
+ TimeDomainAntenna ant2("antenna_name2", ant2Loc, rootCS, t1, t2, t3, t1);
AntennaCollection<TimeDomainAntenna> detector;
-
- // add the antennas to the detector
detector.addAntenna(ant1);
detector.addAntenna(ant2);
@@ -212,17 +332,16 @@ TEST_CASE("Radio", "[processes]") {
auto const particle{Code::Electron};
const auto pmass{get_mass(particle)};
- VelocityVector v0(rootCSZHS, {5e+2_m / second, 5e+2_m / second, 5e+2_m / second});
+ VelocityVector v0(rootCS, {5e+2_m / second, 5e+2_m / second, 5e+2_m / second});
- Vector B0(rootCSZHS, 5_T, 5_T, 5_T);
+ Vector B0(rootCS, 5_T, 5_T, 5_T);
- Line const line(point3, v0);
+ Line const line(trackStart, v0);
auto const k{1_m * ((1_m) / ((1_s * 1_s) * 1_V))};
auto const t = 1e-12_s;
- LeapFrogTrajectory base(point4, v0, B0, k, t);
- // std::cout << "Leap Frog Trajectory is: " << base << std::endl;
+ LeapFrogTrajectory base(trackEnd, v0, B0, k, t);
// create a new stack for each trial
setup::Stack<EnvType> stack;
@@ -232,12 +351,10 @@ TEST_CASE("Radio", "[processes]") {
// compute the necessary momentum
const HEPMomentumType P0{sqrt(E0 * E0 - pmass * pmass)};
-
- // and create the momentum vector
- const auto plab{MomentumVector(rootCSZHS, {0_GeV, 0_GeV, P0})};
+ const auto plab{MomentumVector(rootCS, {0_GeV, 0_GeV, P0})};
// and create the location of the particle in this coordinate system
- const Point pos(rootCSZHS, 50_m, 10_m, 80_m);
+ const Point pos(rootCS, 50_m, 10_m, 80_m);
// add the particle to the stack
auto const particle1{stack.addParticle(std::make_tuple(
@@ -309,7 +426,6 @@ TEST_CASE("Radio", "[processes]") {
const Code particle2{Code::Proton};
const auto pmass{get_mass(particle)};
- const auto pmass2{get_mass(particle2)};
// construct an energy
const HEPEnergyType E0{1_TeV};
// compute the necessary momentumn
@@ -373,11 +489,271 @@ TEST_CASE("Radio", "[processes]") {
} // END: SECTION("Radio extreme cases")
+ SECTION("Process Library") {
+ using IModelInterface =
+ IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
+ using AtmModel = UniformRefractiveIndex<
+ MediumPropertyModel<UniformMagneticField<HomogeneousMedium<IModelInterface>>>>;
+ using EnvType = Environment<AtmModel>;
+ EnvType envCoREAS;
+ CoordinateSystemPtr const& rootCS = envCoREAS.getCoordinateSystem();
+ Point const center{rootCS, 0_m, 0_m, 0_m};
+
+ Vector B1(rootCS, 0_T, 0_T, 1_T);
+ NuclearComposition const protonComposition({Code::Proton}, {1.});
+ const double refractiveIndex{1.000327};
+ const auto density{1_g / cube(1_cm)};
+ auto Medium = EnvType::createNode<Sphere>(
+ center, 10_km * std::numeric_limits<double>::infinity());
+ auto const props = Medium->setModelProperties<AtmModel>(
+ refractiveIndex, Medium::AirDry1Atm, B1, density, protonComposition);
+ envCoREAS.getUniverse()->addChild(std::move(Medium));
+
+ // create the detector
+ const auto ant1Loc{Point(rootCS, 100_m, 2_m, 3_m)};
+ const auto ant2Loc{Point(rootCS, 4_m, 80_m, 6_m)};
+ const TimeType t1{0_s};
+ const TimeType t2{10_s};
+ const InverseTimeType t3{1e+3_Hz};
+ TimeDomainAntenna ant1("antenna_name", ant1Loc, rootCS, t1, t2, t3, t1);
+ TimeDomainAntenna ant2("antenna_name2", ant2Loc, rootCS, t1, t2, t3, t1);
+ AntennaCollection<TimeDomainAntenna> detector;
+ detector.addAntenna(ant1);
+ detector.addAntenna(ant2);
+
+ const auto trackStart{Point(rootCS, 7_m, 8_m, 9_m)};
+ const auto trackEnd{Point(rootCS, 5_m, 5_m, 10_m)};
+
+ // create an electron
+ const Code electron{Code::Electron};
+ const auto pmass{get_mass(electron)};
+
+ VelocityVector v0(rootCS, {1_m / second, 0_m / second, 0_m / second});
+
+ Vector B0(rootCS, 5_T, 5_T, 5_T);
+
+ Line const line(trackStart, v0);
+
+ // create a new stack for each trial
+ setup::Stack<EnvType> stack;
+
+ // construct an energy
+ const HEPEnergyType E0{1_TeV};
+
+ // compute the necessary momentumn
+ const HEPMomentumType P0{sqrt(E0 * E0 - pmass * pmass)};
+ const auto plab{MomentumVector(rootCS, {0_GeV, 0_GeV, P0})};
+
+ // and create the location of the particle in this coordinate system
+ const Point pos(rootCS, 50_m, 10_m, 80_m);
+
+ // add the particle to the stack
+ auto const particle1{stack.addParticle(std::make_tuple(
+ electron, calculate_kinetic_energy(plab.getNorm(), get_mass(electron)),
+ plab.normalized(), pos, 0_ns))};
+
+ // create a radio process instance using CoREAS
+ RadioProcess<decltype(detector),
+ CoREAS<decltype(detector), decltype(StraightPropagator(envCoREAS))>,
+ decltype(StraightPropagator(envCoREAS))>
+ coreas(detector, envCoREAS);
+
+ const auto config = coreas.getConfig();
+
+ CHECK(config["type"].as<std::string>() == "RadioProcess");
+ CHECK(config["algorithm"].as<std::string>() == "CoREAS");
+ CHECK(config["units"]["time"].as<std::string>() == "ns");
+ CHECK(config["units"]["frequency"].as<std::string>() == "GHz");
+ CHECK(config["units"]["electric field"].as<std::string>() == "V/m");
+ CHECK(config["units"]["distance"].as<std::string>() == "m");
+
+ CHECK(config["antennas"]["antenna_name"]["location"][0].as<double>() == 100);
+ CHECK(config["antennas"]["antenna_name"]["location"][1].as<double>() == 2);
+ CHECK(config["antennas"]["antenna_name"]["location"][2].as<double>() == 3);
+
+ 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: TEST_CASE("Radio", "[processes]")
TEST_CASE("Antennas") {
- SECTION("TimeDomainAntenna") {
+ SECTION("TimeDomainAntenna Constructor") {
+ Environment<IRefractiveIndexModel<IMediumModel>> env;
+ const auto rootCS = env.getCoordinateSystem();
+ auto const antPos = Point(rootCS, {0_m, 0_m, 0_m});
+ TimeType const tStart(0_s);
+ TimeType const duration(10_ns);
+ InverseTimeType const sampleRate(1_GHz);
+ TimeType const groundHitTime(1e3_ns);
+
+ TimeDomainAntenna const antenna("antenna", antPos, rootCS, tStart, duration,
+ sampleRate, groundHitTime);
+
+ // All waveforms are of equal non-zero size
+ CHECK(antenna.getWaveformX().size() == antenna.getWaveformY().size());
+ CHECK(antenna.getWaveformX().size() == antenna.getWaveformZ().size());
+ CHECK(antenna.getWaveformX().size() > 0);
+
+ // All waveform values are initialized to zero
+ for (auto const& val : antenna.getWaveformX()) { CHECK(val * 0 == val); }
+ for (auto const& val : antenna.getWaveformY()) { CHECK(val * 0 == val); }
+ for (auto const& val : antenna.getWaveformZ()) { CHECK(val * 0 == val); }
+
+ // check that variables were set properly
+ CHECK("antenna" == antenna.getName());
+ CHECK(sampleRate == antenna.getSampleRate());
+ CHECK(tStart == antenna.getStartTime());
+
+ // and check that the antenna is at the right location
+ CHECK((antenna.getLocation() - antPos).getNorm() < 1e-12 * 1_m);
+ } // END: SECTION("TimeDomainAntenna Constructor")
+
+ SECTION("TimeDomainAntenna Bad Constructor") {
+ Environment<IRefractiveIndexModel<IMediumModel>> env;
+ const auto rootCS = env.getCoordinateSystem();
+ auto const antPos = Point(rootCS, {0_m, 0_m, 0_m});
+ TimeType const tStart(0_s);
+ TimeType const duration(1e3_ns);
+ InverseTimeType const sampleRate(1_GHz);
+ TimeType const groundHitTime(10_ns);
+
+ // Giving zero or negative values for sampling rate and duration
+ TimeDomainAntenna const antenna_bad1("bad_antenna", antPos, rootCS, tStart, -13_ns,
+ sampleRate, groundHitTime);
+ TimeDomainAntenna const antenna_bad2("bad_antenna", antPos, rootCS, tStart, 0_ns,
+ sampleRate, groundHitTime);
+ TimeDomainAntenna const antenna_bad3("bad_antenna", antPos, rootCS, tStart, duration,
+ -1_GHz, groundHitTime);
+ } // END: SECTION("TimeDomainAntenna Bad Constructor")
+
+ SECTION("TimeDomainAntenna Receive Efield") {
+ // Checks that the basic functionality of the receive function is working properly
+
+ using EnvType = Environment<IRefractiveIndexModel<IMediumModel>>;
+ EnvType env;
+
+ const auto rootCS = env.getCoordinateSystem();
+
+ auto const point1 = Point(rootCS, {1_m, 2_m, 3_m});
+ auto const point2 = Point(rootCS, {4_m, 5_m, 6_m});
+
+ // create times for the antenna
+ const TimeType t1{10_s};
+ const TimeType t2{10_s};
+ const InverseTimeType t3{1 / 1_s};
+ const TimeType t4{11_s};
+
+ // make the two antennas with different start times
+ TimeDomainAntenna ant1("antenna_name", point1, rootCS, t1, t2, t3, t1);
+ TimeDomainAntenna ant2("antenna_name", point2, rootCS, t4, t2, t3, t4);
+
+ Vector<dimensionless_d> receiveVec1(rootCS, {0, 0, 1});
+ Vector<dimensionless_d> receiveVec2(rootCS, {0, 1, 0});
+
+ Vector<ElectricFieldType::dimension_type> eField1(
+ rootCS, {10_V / 1_m, 10_V / 1_m, 10_V / 1_m});
+ Vector<ElectricFieldType::dimension_type> eField2(
+ rootCS, {20_V / 1_m, 20_V / 1_m, 20_V / 1_m});
+
+ // inject efield into ant1
+ ant1.receive(15_s, receiveVec1, eField1);
+ REQUIRE(ant1.getWaveformX()[5] - 10 == 0);
+ REQUIRE(ant1.getWaveformX()[5] == ant1.getWaveformY()[5]);
+ REQUIRE(ant1.getWaveformX()[5] == ant1.getWaveformZ()[5]);
+
+ // inject efield but with different receive vector into ant2
+ ant2.receive(16_s, receiveVec2, eField1);
+ REQUIRE(ant1.getWaveformX()[5] ==
+ ant2.getWaveformX()[5]); // Currently receive vector does nothing
+ ant2.reset();
+ REQUIRE(ant2.getWaveformX()[5] == 0); // reset was successful
+
+ // inject the other eField into ant2
+ ant2.receive(16_s, receiveVec2, eField2);
+ REQUIRE(ant2.getWaveformX()[5] - 20 == 0);
+ REQUIRE(ant2.getWaveformX()[5] == ant2.getWaveformY()[5]);
+ REQUIRE(ant2.getWaveformX()[5] == ant2.getWaveformZ()[5]);
+
+ // make sure the next one is empty before filling it
+ REQUIRE(ant2.getWaveformX()[6] == 0);
+ ant2.receive(17_s, receiveVec2, eField2);
+ REQUIRE(ant2.getWaveformX()[6] - 20 == 0);
+
+ // reset ant1 and then put values in out of range
+ ant1.reset();
+ ant1.receive(-1000_s, receiveVec1, eField1);
+ for (auto const& val : ant1.getWaveformX()) { CHECK(val * 0 == val); }
+ ant1.reset();
+ ant1.receive(t1 + t2 + 1_s, receiveVec1, eField1);
+ for (auto const& val : ant1.getWaveformX()) { CHECK(val * 0 == val); }
+ } // END: SECTION("TimeDomainAntenna Receive EField")
+
+ SECTION("TimeDomainAntenna Receive Vector Potential") {
+ // Checks that the basic functionality of the receive function is working properly
+
+ using EnvType = Environment<IRefractiveIndexModel<IMediumModel>>;
+ EnvType env;
+
+ const auto rootCS = env.getCoordinateSystem();
+
+ auto const point1 = Point(rootCS, {1_m, 2_m, 3_m});
+ auto const point2 = Point(rootCS, {4_m, 5_m, 6_m});
+
+ // create times for the antenna
+ const TimeType t1{10_s};
+ const TimeType t2{10_s};
+ const InverseTimeType t3{1 / 1_s};
+ const TimeType t4{11_s};
+
+ // make the two antennas with different start times
+ TimeDomainAntenna ant1("antenna_name", point1, rootCS, t1, t2, t3, t1);
+ TimeDomainAntenna ant2("antenna_name", point2, rootCS, t4, t2, t3, t4);
+
+ Vector<dimensionless_d> receiveVec1(rootCS, {0, 0, 1});
+ Vector<dimensionless_d> receiveVec2(rootCS, {0, 1, 0});
+
+ Vector<VectorPotentialType::dimension_type> vectorPotential1(
+ rootCS, {10_V * 1_s / 1_m, 10_V * 1_s / 1_m, 10_V * 1_s / 1_m});
+ Vector<VectorPotentialType::dimension_type> vectorPotential2(
+ rootCS, {20_V * 1_s / 1_m, 20_V * 1_s / 1_m, 20_V * 1_s / 1_m});
+
+ // inject efield into ant1
+ ant1.receive(15_s, receiveVec1, vectorPotential1);
+ REQUIRE(ant1.getWaveformX()[5] - 10 == 0);
+ REQUIRE(ant1.getWaveformX()[5] == ant1.getWaveformY()[5]);
+ REQUIRE(ant1.getWaveformX()[5] == ant1.getWaveformZ()[5]);
+
+ // inject efield but with different receive vector into ant2
+ ant2.receive(16_s, receiveVec2, vectorPotential1);
+ REQUIRE(ant1.getWaveformX()[5] ==
+ ant2.getWaveformX()[5]); // Currently receive vector does nothing
+ ant2.reset();
+ REQUIRE(ant2.getWaveformX()[5] == 0); // reset was successful
+
+ // inject the other eField into ant2
+ ant2.receive(16_s, receiveVec2, vectorPotential2);
+ REQUIRE(ant2.getWaveformX()[5] - 20 == 0);
+ REQUIRE(ant2.getWaveformX()[5] == ant2.getWaveformY()[5]);
+ REQUIRE(ant2.getWaveformX()[5] == ant2.getWaveformZ()[5]);
+
+ // make sure the next one is empty before filling it
+ REQUIRE(ant2.getWaveformX()[6] == 0);
+ ant2.receive(17_s, receiveVec2, vectorPotential2);
+ REQUIRE(ant2.getWaveformX()[6] - 20 == 0);
+
+ // reset ant1 and then put values in out of range
+ ant1.reset();
+ ant1.receive(-1000_s, receiveVec1, vectorPotential1);
+ for (auto const& val : ant1.getWaveformX()) { CHECK(val * 0 == val); }
+ ant1.reset();
+ ant1.receive(t1 + t2 + 1_s, receiveVec1, vectorPotential1);
+ for (auto const& val : ant1.getWaveformX()) { CHECK(val * 0 == val); }
+ } // END: SECTION("TimeDomainAntenna Receive Vector Potential")
+
+ SECTION("TimeDomainAntenna AntennaCollection") {
// create an environment so we can get a coordinate system
using EnvType = Environment<IRefractiveIndexModel<IMediumModel>>;
@@ -391,10 +767,10 @@ TEST_CASE("Antennas") {
const auto point2{Point(env6.getCoordinateSystem(), 4_m, 5_m, 6_m)};
// get a coordinate system
- const CoordinateSystemPtr rootCS6 = env6.getCoordinateSystem();
+ const CoordinateSystemPtr rootCS = env6.getCoordinateSystem();
auto Medium6 = EnvType::createNode<Sphere>(
- Point{rootCS6, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
+ Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
auto const props6 = Medium6->setModelProperties<UniRIndex>(
1, 1_kg / (1_m * 1_m * 1_m), NuclearComposition({Code::Nitrogen}, {1.}));
@@ -407,100 +783,90 @@ TEST_CASE("Antennas") {
const InverseTimeType t3{1 / 1_s};
const TimeType t4{11_s};
- // check that I can create an antenna at (1, 2, 3)
- TimeDomainAntenna ant1("antenna_name", point1, rootCS6, t1, t2, t3, t1);
- TimeDomainAntenna ant2("antenna_name2", point2, rootCS6, t4, t2, t3, t4);
-
- // assert that the antenna name is correct
- REQUIRE(ant1.getName() == "antenna_name");
- REQUIRE(ant2.getName() == "antenna_name2");
-
- // and check that the antenna is at the right location
- REQUIRE((ant1.getLocation() - point1).getNorm() < 1e-12 * 1_m);
- REQUIRE((ant2.getLocation() - point2).getNorm() < 1e-12 * 1_m);
-
// construct a radio detector instance to store our antennas
AntennaCollection<TimeDomainAntenna> detector;
- // add this antenna to the process
- detector.addAntenna(ant1);
- detector.addAntenna(ant2);
- CHECK(detector.size() == 2);
-
- // get a unit vector
- Vector<dimensionless_d> v1(rootCS6, {0, 0, 1});
- Vector<ElectricFieldType::dimension_type> v11(rootCS6,
- {10_V / 1_m, 10_V / 1_m, 10_V / 1_m});
-
- Vector<dimensionless_d> v2(rootCS6, {0, 1, 0});
- Vector<ElectricFieldType::dimension_type> v22(rootCS6,
- {20_V / 1_m, 20_V / 1_m, 20_V / 1_m});
-
- // use receive methods
- ant1.receive(15_s, v1, v11);
- ant2.receive(16_s, v2, v22);
-
- // use getDataX,Y,Z() and getAxis() methods
- auto Ex = ant1.getWaveformX();
- CHECK(Ex[5] - 10 == 0);
- auto tx = ant1.getAxis();
- CHECK(tx[5] - 5 * 1_s / 1_ns == Approx(0.0));
- auto Ey = ant1.getWaveformY();
- CHECK(Ey[5] - 10 == 0);
- auto Ez = ant1.getWaveformZ();
- CHECK(Ez[5] - 10 == 0);
- auto ty = ant1.getAxis();
- auto tz = ant1.getAxis();
- CHECK(tx[5] - ty[5] == 0);
- CHECK(ty[5] - tz[5] == 0);
- auto Ex2 = ant2.getWaveformX();
- CHECK(Ex2[5] - 20 == 0);
- auto Ey2 = ant2.getWaveformY();
- CHECK(Ey2[5] - 20 == 0);
- auto Ez2 = ant2.getWaveformZ();
- CHECK(Ez2[5] - 20 == 0);
-
// the following creates a star-shaped pattern of antennas in the ground
- AntennaCollection<TimeDomainAntenna> detector__;
const auto point11{Point(env6.getCoordinateSystem(), 1000_m, 20_m, 30_m)};
const TimeType t2222{1e-6_s};
const InverseTimeType t3333{1e+9_Hz};
std::vector<std::string> antenna_names;
std::vector<Point> antenna_locations;
- for (auto radius_ = 100_m; radius_ <= 200_m; radius_ += 100_m) {
- for (auto phi_ = 0; phi_ <= 315; phi_ += 45) {
- auto phiRad_ = phi_ / 180. * M_PI;
- auto const point_{Point(env6.getCoordinateSystem(), radius_ * cos(phiRad_),
- radius_ * sin(phiRad_), 0_m)};
- antenna_locations.push_back(point_);
- auto time__{(point11 - point_).getNorm() / constants::c};
- const int rr_ = static_cast<int>(radius_ / 1_m);
- std::string var_ = "antenna_R=" + std::to_string(rr_) +
- "_m-Phi=" + std::to_string(phi_) + "degrees";
- antenna_names.push_back(var_);
- TimeDomainAntenna ant111(var_, point_, rootCS6, time__, t2222, t3333, time__);
- detector__.addAntenna(ant111);
+ for (auto radius = 100_m; radius <= 200_m; radius += 100_m) {
+ for (auto phi = 0; phi <= 315; phi += 45) {
+ auto phiRad = phi / 180. * M_PI;
+ auto const point{Point(env6.getCoordinateSystem(), radius * cos(phiRad),
+ radius * sin(phiRad), 0_m)};
+ antenna_locations.push_back(point);
+ auto time__{(point11 - point).getNorm() / constants::c};
+ const int rr_ = static_cast<int>(radius / 1_m);
+ std::string name = "antenna_R=" + std::to_string(rr_) +
+ "_m-Phi=" + std::to_string(phi) + "degrees";
+ antenna_names.push_back(name);
+ TimeDomainAntenna ant(name, point, rootCS, time__, t2222, t3333, time__);
+ detector.addAntenna(ant);
}
}
- CHECK(detector__.size() == 16);
- CHECK(detector__.getAntennas().size() == 16);
+ CHECK(detector.size() == 16);
+ CHECK(detector.getAntennas().size() == 16);
int i = 0;
// this prints out the antenna names and locations
- for (auto const antenna : detector__.getAntennas()) {
+ for (auto const& antenna : detector.getAntennas()) {
CHECK(antenna.getName() == antenna_names[i]);
CHECK(distance(antenna.getLocation(), antenna_locations[i]) / 1_m == 0);
i++;
}
// Check the .at() method for radio detectors
- for (size_t i = 0; i <= (detector__.size() - 1); i++) {
- CHECK(detector__.at(i).getName() == antenna_names[i]);
- CHECK(distance(detector__.at(i).getLocation(), antenna_locations[i]) / 1_m == 0);
+ for (int i = 0; i <= (detector.size() - 1); i++) {
+ CHECK(detector.at(i).getName() == antenna_names[i]);
+ CHECK(distance(detector.at(i).getLocation(), antenna_locations[i]) / 1_m == 0);
+ }
+
+ } // END: SECTION("TimeDomainAntenna AntennaCollection")
+
+ SECTION("TimeDomainAntenna Library") {
+ // Runs checks that the file readers are working properly
+ Environment<IRefractiveIndexModel<IMediumModel>> env;
+ const auto rootCS = env.getCoordinateSystem();
+ auto const antPos = Point(rootCS, {0_m, 0_m, 0_m});
+ TimeType const tStart(0_s);
+ TimeType const duration(10_ns);
+ 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);
}
- } // END: SECTION("TimeDomainAntenna")
+ // 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")
} // END: TEST_CASE("Antennas")