IAP GITLAB

Skip to content
Snippets Groups Projects
Commit bd5d5c8b authored by Alan Coleman's avatar Alan Coleman
Browse files

Merge branch 'expand_radio_test_coverage' into 'master' 

Expand radio test coverage, fix warnings

Closes #553

See merge request !481
parents 70014749 1c35790a
No related branches found
No related tags found
1 merge request!481Expand radio test coverage, fix warnings
Pipeline #9952 passed with warnings
Stage: quality
Stage: config
Stage: build_test_example
Stage: python
Stage: optional
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 864 additions and 477 deletions
corsika/detail/modules/radio/CoREAS.inl
+ 285
− 288
View file @ bd5d5c8b
......@@ -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;
}
......
corsika/detail/modules/radio/antennas/Antenna.inl
+ 1
− 1
View file @ bd5d5c8b
......@@ -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 {
......
corsika/detail/modules/radio/antennas/TimeDomainAntenna.inl
+ 34
− 19
View file @ bd5d5c8b
......@@ -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);
......
corsika/detail/modules/radio/propagators/SimplePropagator.inl
+ 7
− 5
View file @ bd5d5c8b
......@@ -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()
......
corsika/detail/modules/radio/propagators/StraightPropagator.inl
+ 24
− 18
View file @ bd5d5c8b
......@@ -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
corsika/modules/radio/CoREAS.hpp
+ 1
− 1
View file @ bd5d5c8b
......@@ -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.
......
corsika/modules/radio/ZHS.hpp
+ 1
− 1
View file @ bd5d5c8b
......@@ -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.
......
corsika/modules/radio/antennas/TimeDomainAntenna.hpp
+ 5
− 4
View file @ bd5d5c8b
......@@ -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.
......
tests/modules/testRadio.cpp
+ 506
− 140
View file @ bd5d5c8b
This diff is collapsed.
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

data privacy - impressum