From 35158b337d246ed61ee050f5c3172e928c36683f Mon Sep 17 00:00:00 2001
From: Nikos Karastathis <n.karastathis@kit.edu>
Date: Fri, 26 Feb 2021 15:05:41 +0100
Subject: [PATCH] Antennas are tested and work. TODO: fix the foor loop for
antennas.
---
corsika/modules/radio/CoREAS.hpp | 414 +++++++++++++-----
corsika/modules/radio/RadioProcess.hpp | 5 +-
corsika/modules/radio/ZHS.hpp | 25 +-
corsika/modules/radio/antennas/Antenna.hpp | 2 +-
.../radio/antennas/TimeDomainAntenna.hpp | 40 +-
.../modules/radio/detectors/RadioDetector.hpp | 2 +-
tests/modules/testRadio.cpp | 281 +++++++++++-
7 files changed, 591 insertions(+), 178 deletions(-)
diff --git a/corsika/modules/radio/CoREAS.hpp b/corsika/modules/radio/CoREAS.hpp
index a9f8f9c05..e0a4ea1f7 100755
--- a/corsika/modules/radio/CoREAS.hpp
+++ b/corsika/modules/radio/CoREAS.hpp
@@ -11,6 +11,7 @@
#include <corsika/modules/radio/propagators/StraightPropagator.hpp>
#include <corsika/framework/geometry/QuantityVector.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/modules/radio/propagators/SignalPath.hpp>
#include <cmath>
namespace corsika {
@@ -50,152 +51,325 @@ namespace corsika {
template <typename Particle, typename Track>
ProcessReturn simulate(Particle& particle, Track const& track) const {
- // set threshold for application of ZHS-like approximation
- const double approxThreshold_ {1.0e-3};
-
- //get global simulation time for that track. (This is my best guess for now)
- auto startTime_ {particle.getTime()
- - track.getDuration()}; // time at start point of track.
- auto endTime_ {particle.getTime()}; // time at end point of track.
+ // get the global simulation time for that track. (best guess for now)
+ auto startTime_ {particle.getTime() - track.getDuration()}; // time at the start point of the track hopefully.
+ auto endTime_ {particle.getTime()};
- // beta is defined as velocity / speed of light
- auto startBeta_ {track.getVelocity(0) / constants::c};
- auto endBeta_ {track.getVelocity(1) / constants::c};
- //TODO: check if they are different!!! They shouldn't!!!
-// auto startBeta_ {(track.getPosition(0) / (particle.getTime()
-// - track.getDuration()) ) / constants::c};
-// auto endBeta_ {(track.getPosition(1) / particle.getTime() ) / constants::c};
-
- // calculate gamma factor using beta (the proper way would be with energy over mass)
- auto startGamma_ {1. / sqrt(1. - (startBeta_ * startBeta_))};
- auto endGamma_ {1. / sqrt(1. - (endBeta_ * endBeta_))};
+ // gamma factor is calculated using beta
+// auto startGamma_ {1. / sqrt(1. - (startBeta_ * startBeta_))};
+// auto endGamma_ {1. / sqrt(1. - (endBeta_ * endBeta_))};
// get start and end position of the track
auto startPoint_ {track.getPosition(0)};
auto endPoint_ {track.getPosition(1)};
+ // beta is velocity / speed of light. Start & end should be the same!
+ auto beta_ {((endPoint_ - startPoint_) / (endTime_ - startTime_)).normalized()};
+
// get particle charge
auto const charge_ {get_charge(particle.getPID())};
- // we loop over each antenna in the collection
- for (auto& antenna : detector_.getAntennas()) {
+ // set threshold for application of ZHS-like approximation.
+ const double approxThreshold_ {1.0e-3};
- ElectricFieldVector EV1_ {0_V / 0_m};
- ElectricFieldVector EV2_ {0_V / 0_m};
- ElectricFieldVector EV3_ {0_V / 0_m};
+ // we loop over each antenna in the collection.
+ for (auto& antenna : detector_.getAntennas()) {
- auto startPointReceiveTime_ {0_ns};
- auto endPointReceiveTime_ {0_ns};
+ std::vector<ElectricFieldVector> EVstart_;
+ std::vector<ElectricFieldVector> EVend_;
+ std::vector<TimeType> startTTimes_;
+ std::vector<TimeType> endTTimes_;
+ std::vector<double> preDoppler;
+ std::vector<double> postDoppler;
+ std::vector<QuantityVector<dimensionless_d>> ReceiveVectorsStart_;
+ std::vector<QuantityVector<dimensionless_d>> ReceiveVectorsEnd_;
// get the Path (path1) from the start "endpoint" to the antenna.
- // This is a SignalPathCollection
- auto paths1{this->propagator_.propagate(startPoint_, antenna.getLocation())};
+ // This is a Signal Path Collection
+ auto paths1 {this->propagator_.propagate(startPoint_, antenna.getLocation(), 1_nm)}; // TODO: Need to add the stepsize to .propagate()!!!!
- // set postDoppler to approximate threshold in advance and check after loop
- // if we need to perform ZHS-like approximation
- auto preDoppler_ {approxThreshold_};
// now loop over the paths for startpoint that we got above
for (auto const& path : paths1) {
- preDoppler_ = 1. - path.average_refractive_index_ *
- startBeta_ * path.emit_;
+ // calculate preDoppler factor
+ double preDoppler_{1. - path.average_refractive_index_ *
+ beta_ * path.emit_};
+
+ // store it to the preDoppler std::vector for later comparisons
+ preDoppler.push_back(preDoppler_);
+
+ // calculate receive times for startpoint
+ auto startPointReceiveTime_ {path.total_time_ + startTime_};
+
+ // store it to startTTimes_ std::vector for later use
+ startTTimes_.push_back(startPointReceiveTime_);
- if (preDoppler_ > approxThreshold_) {
- startPointReceiveTime_ = path.total_time_ +
- startTime_ ; // might do it on the fly
+ // store the receive unit vector
+ ReceiveVectorsStart_.push_back(path.receive_);
- // CoREAS calculation -> get ElectricFieldVector1
- EV1_= (charge_ / constants::c) *
- path.receive_.cross(path.receive_.cross(startBeta_)) /
- (path.R_distance_ * preDoppler_);
+ // calculate electric field vector for startpoint
+ auto EV1_= (charge_ / constants::c) *
+ path.receive_.cross(path.receive_.cross(beta_)) /
+ (path.R_distance_ * preDoppler_);
- // pass it to the antenna
- antenna.receive(startPointReceiveTime_, path.receive_, EV1_);
- } else {
- continue;
- } // END preDoppler check
+ // store it to EVstart_ std::vector for later use
+ EVstart_.push_back(EV1_);
- } // END: loop over paths for startpoint
+ } // End of looping over paths1
// get the Path (path2) from the end "endpoint" to the antenna.
// This is a SignalPathCollection
- auto paths2{this->propagator_.propagate(endPoint_, antenna.getLocation())};
+ auto paths2 {this->propagator_.propagate(endPoint_, antenna.getLocation())};
- // set postDoppler to approximate threshold in advance and check after loop
- // if we need to perform ZHS-like approximation
- auto postDoppler_ {approxThreshold_};
// now loop over the paths for endpoint that we got above
for (auto const& path : paths2) {
- postDoppler_ = 1. - path.average_refractive_index_ *
- endBeta_ * path.emit_;
-
- if (preDoppler_ > approxThreshold_) {
- auto endPointReceiveTime_ = path.total_time + endTime_ ; // might do it on the fly
-
- //CoREAS calculation -> get ElectricFieldVector2
- EV2_ = (charge_ / constants::c) *
- path.receive_.cross(path.receive_.cross(endBeta_)) /
- (path.R_distance_ * postDoppler_);
-
- // pass it to the antenna
- antenna.receive(endPointReceiveTime_, path.receive_, EV2_);
- } else {
- continue;
- } // END postDoppler check
-
- } // END: loop over paths for endpoint
-
- // perform ZHS-like calculation close to Cherenkov angle
- if (fabs(preDoppler_) <= approxThreshold_ || fabs(postDoppler_) <= approxThreshold_) {
- // get global simulation time for the middle point of that track. (This is my best guess for now)
- auto midTime_{particle.getTime() - (track.getDuration() / 2)};
-
- // beta is defined as velocity / speed of light
- auto midBeta_{track.getVelocity(0.5) / constants::c};
- // auto midBeta_ {(track.getPosition(0.5) / (particle.getTime()
- // - (track.getDuration() / 2) ) / constants::c};
-
- // calculate gamma factor using beta (the proper way would be with energy over mass)
- // auto midGamma_ {1. / sqrt(1. - (midBeta_ * midBeta_))};
-
- // get start and end position of the track
- auto midPoint_{track.getPosition(0.5)};
-
- // get the Path (path3) from the middle "endpoint" to the antenna.
- // This is a SignalPathCollection
- auto paths3{this->propagator_.propagate(midPoint_, antenna.getLocation())};
-
- // now loop over the paths for endpoint that we got above
- for (auto const& path : paths3) {
- midDoppler_ = 1. - path.average_refractive_index_ * midBeta_ * path.emit_;
-
- auto const midPointReceiveTime_{path.total_time_ +
- midTime_}; // might do it on the fly
-
- // CoREAS calculation -> get ElectricFieldVector3 for "midPoint"
- EV3_ = (charge_ / constants::c) *
- path.receive_.cross(path.receive_.cross(midBeta_)) /
- (path.R_distance_ * midDoppler_);
-
- // pass it to the antenna but first check if "start" or "end" are double counted!!!
- if (fabs(preDoppler_) > approxThreshold_ &&
- fabs(postDoppler_) <= approxThreshold_) {
-
- antenna.receive(startPointReceiveTime_, path.receive_, -EV1_);
- antenna.receive(midPointReceiveTime_, path.receive_, EV3_);
- } else if (fabs(preDoppler_) <= approxThreshold_ &&
- fabs(postDoppler_) > approxThreshold_) {
- antenna.receive(endPointReceiveTime_, path.receive_, -EV2_);
- antenna.receive(midPointReceiveTime_, path.receive_, EV3_);
-
- } else {
- antenna.receive(midPointReceiveTime_, path.receive_, EV3_);
- } // end deleting double-counted values
- }
- } // end of ZHS-like approximation
-
- } // END: loop over antennas
- }
+
+ double postDoppler_{1. - path.average_refractive_index_ *
+ beta_ * path.emit_}; // maybe this is path.receive_ (?)
+
+ // store it to the postDoppler std::vector for later comparisons
+ postDoppler.push_back(postDoppler_);
+
+ // calculate receive times for endpoint
+ auto endPointReceiveTime_ {path.total_time_ + endTime_};
+
+ // store it to endTTimes_ std::vector for later use
+ endTTimes_.push_back(endPointReceiveTime_);
+
+ // store the receive unit vector
+ ReceiveVectorsEnd_.push_back(path.receive_);
+
+ // calculate electric field vector for endpoint
+ auto EV2_= (charge_ / constants::c) *
+ path.receive_.cross(path.receive_.cross(beta_)) /
+ (path.R_distance_ * postDoppler_);
+
+ // store it to EVstart_ std::vector for later use
+ EVend_.push_back(EV2_);
+
+ } // End of looping over paths2
+
+ // start doing comparisons for preDoppler and postDoppler
+ // first check that start and end paths have the same number of paths
+ if (EVstart_.size() == EVend_.size()) {
+
+ // use this to access different elements of std::vectors
+ std::size_t index = 0;
+ for (auto& preDoppler__ : preDoppler) {
+
+ // redistribute contributions over time scale defined by the observation time resolution
+ // this is to make sure that "start" and "end" won't end up in the same bin (xtensor)!!
+ if ((preDoppler__ < 1.e-9) || (postDoppler.at(index) < 1.e-9)) {
+
+ auto gridResolution_ {antenna.duration_};
+ auto deltaT_ { endTTimes_.at(index) - startTTimes_.at(index) };
+
+ if (fabs(deltaT_) < gridResolution_) {
+
+ EVstart_.at(index) = EVstart_.at(index) * fabs(deltaT_ / gridResolution_);
+ EVend_.at(index) = EVend_.at(index) * fabs(deltaT_ / gridResolution_);
+
+ const long startBin = static_cast<long>(floor(startTTimes_.at(index)/gridResolution_+0.5l));
+ const long endBin = static_cast<long>(floor(endTTimes_.at(index)/gridResolution_+0.5l));
+ const double startBinFraction = (startTTimes_.at(index)/gridResolution_)-floor(startTTimes_.at(index)/gridResolution_);
+ const double endBinFraction = (endTTimes_.at(index)/gridResolution_)-floor(endTTimes_.at(index)/gridResolution_);
+
+ // only do timing modification if contributions would land in same bin
+ if (startBin == endBin) {
+
+ // if startE arrives before endE
+ if (deltaT_ >= 0) {
+ if ((startBinFraction >= 0.5) && (endBinFraction >= 0.5)) // both points left of bin center
+ {
+ startTTimes_.at(index) = startTTimes_.at(index) - gridResolution_; // shift EV1_ to previous gridpoint
+ }
+ else if ((startBinFraction < 0.5) && (endBinFraction < 0.5)) // both points right of bin center
+ {
+ endTTimes_.at(index) = endTTimes_.at(index) + gridResolution_; // shift EV2_ to next gridpoint
+ }
+ else // points on both sides of bin center
+ {
+ const double leftDist = 1.0-startBinFraction;
+ const double rightDist = endBinFraction;
+ // check if asymmetry to right or left
+ if (rightDist >= leftDist)
+ {
+ endTTimes_.at(index) = endTTimes_.at(index) + gridResolution_; // shift EV2_ to next gridpoint
+ }
+ else
+ {
+ startTTimes_.at(index) = startTTimes_.at(index) - 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
+ {
+ endTTimes_.at(index) = endTTimes_.at(index) - gridResolution_; // shift EV2_ to previous gridpoint
+ }
+ else if ((startBinFraction < 0.5) && (endBinFraction < 0.5)) // both points right of bin center
+ {
+ startTTimes_.at(index) = startTTimes_.at(index) + gridResolution_; // shift EV1_ to next gridpoint
+ }
+ else // points on both sides of bin center
+ {
+ const double leftDist = 1.0-endBinFraction;
+ const double rightDist = startBinFraction;
+ // check if asymmetry to right or left
+ if (rightDist >= leftDist)
+ {
+ startTTimes_.at(index) = startTTimes_.at(index) + gridResolution_; // shift EV1_ to next gridpoint
+ }
+ else
+ {
+ endTTimes_.at(index) = endTTimes_.at(index) - gridResolution_; // shift EV2_ to previous gridpoint
+ }
+ }
+ } // End of else statement
+ } // End of if for startbin == endbin
+ } // End of if deltaT < gridresolution
+ } // End of checking for very small doppler factors
+
+ // perform ZHS-like calculation close to Cherenkov angle
+ if (fabs(preDoppler__) <= approxThreshold_ || fabs(postDoppler.at(index)) <= approxThreshold_) {
+
+ // get global simulation time for the middle point of that track. (This is my best guess for now)
+ auto midTime_{particle.getTime() - (track.getDuration() / 2)};
+
+ // get "mid" position of the track (that may not work properly)
+ auto midPoint_{track.getPosition(0.5)};
+
+ // get the Path (path3) from the middle "endpoint" to the antenna.
+ // This is a SignalPathCollection
+ auto paths3{this->propagator_.propagate(midPoint_, antenna.getLocation())};
+
+// std::size_t j_index {0}; // this will be useful for multiple paths (aka curved propagators)
+ // now loop over the paths for endpoint that we got above
+ for (auto const& path : paths3) {
+
+// EVstart_.erase(EVstart_.begin() + index + j_index); // this should work for curved + curved propagators
+// EVend_.erase(EVend_.begin() + index + j_index); // for now just use one index and not j_index since at the moment you are working with StraightPropagator
+
+ auto const midPointReceiveTime_{path.total_time_ + midTime_};
+ auto midDoppler_{1. - path.average_refractive_index_ * beta_ * path.emit_};
+
+ // change the values of the receive unit vectors of start and end
+ ReceiveVectorsStart_.at(index) = path.receive_;
+ ReceiveVectorsEnd_.at(index) = path.receive_;
+
+ // CoREAS calculation -> get ElectricFieldVector3 for "midPoint"
+ ElectricFieldVector EVmid_ = (charge_ / constants::c) *
+ path.receive_.cross(path.receive_.cross(beta_)) /
+ (path.R_distance_ * midDoppler_);
+
+// EVstart_.insert(EVstart_.begin() + index + j_index, EVmid_); // this should work for curved + curved propagators
+// EVend_.insert(EVend_.begin() + index + j_index, - EVmid_); // for now just use one index and not j_index since at the moment you are working with StraightPropagator
+ EVstart_.at(index) = EVmid_;
+ EVend_.at(index) = - EVmid_;
+
+ auto deltaT_{midPoint_.getNorm() / (constants::c * beta_ * fabs(midDoppler_))};
+
+ if (startTTimes_.at(index) < endTTimes_.at(index)) // EVstart_ arrives earlier
+ {
+ startTTimes_.at(index) = midPointReceiveTime_ - 0.5 * deltaT_;
+ endTTimes_.at(index) = midPointReceiveTime_ + 0.5 * deltaT_;
+ }
+ else // EVend_ arrives earlier
+ {
+ startTTimes_.at(index) = midPointReceiveTime_ + 0.5 * deltaT_;
+ endTTimes_.at(index) = midPointReceiveTime_ - 0.5 * deltaT_;
+ }
+
+ const long double gridResolution_{antenna.duration_};
+ deltaT_ = endTTimes_.at(index) - startTTimes_.at(index);
+
+ // redistribute contributions over time scale defined by the observation time resolution
+ if (fabs(deltaT_) < gridResolution_) {
+
+ EVstart_.at(index) = EVstart_.at(index) * fabs(deltaT_ / gridResolution_);
+ EVend_.at(index) = EVend_.at(index) * fabs(deltaT_ / gridResolution_);
+
+ const long startBin = static_cast<long>(floor(startTTimes_.at(index)/gridResolution_+0.5l));
+ const long endBin = static_cast<long>(floor(endTTimes_.at(index)/gridResolution_+0.5l));
+ const double startBinFraction = (startTTimes_.at(index)/gridResolution_)-floor(startTTimes_.at(index)/gridResolution_);
+ const double endBinFraction = (endTTimes_.at(index)/gridResolution_)-floor(endTTimes_.at(index)/gridResolution_);
+
+ // only do timing modification if contributions would land in same bin
+ if (startBin == endBin) {
+
+ // if startE arrives before endE
+ if (deltaT_ >= 0) {
+ if ((startBinFraction >= 0.5) && (endBinFraction >= 0.5)) // both points left of bin center
+ {
+ startTTimes_.at(index) = startTTimes_.at(index) - gridResolution_; // shift EV1_ to previous gridpoint
+ }
+ else if ((startBinFraction < 0.5) && (endBinFraction < 0.5)) // both points right of bin center
+ {
+ endTTimes_.at(index) = endTTimes_.at(index) + gridResolution_; // shift EV2_ to next gridpoint
+ }
+ else // points on both sides of bin center
+ {
+ const double leftDist = 1.0-startBinFraction;
+ const double rightDist = endBinFraction;
+ // check if asymmetry to right or left
+ if (rightDist >= leftDist)
+ {
+ endTTimes_.at(index) = endTTimes_.at(index) + gridResolution_; // shift EV2_ to next gridpoint
+ }
+ else
+ {
+ startTTimes_.at(index) = startTTimes_.at(index) - 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
+ {
+ endTTimes_.at(index) = endTTimes_.at(index) - gridResolution_; // shift EV2_ to previous gridpoint
+ }
+ else if ((startBinFraction < 0.5) && (endBinFraction < 0.5)) // both points right of bin center
+ {
+ startTTimes_.at(index) = startTTimes_.at(index) + gridResolution_; // shift EV1_ to next gridpoint
+ }
+ else // points on both sides of bin center
+ {
+ const double leftDist = 1.0-endBinFraction;
+ const double rightDist = startBinFraction;
+ // check if asymmetry to right or left
+ if (rightDist >= leftDist)
+ {
+ startTTimes_.at(index) = startTTimes_.at(index) + gridResolution_; // shift EV1_ to next gridpoint
+ }
+ else
+ {
+ endTTimes_.at(index) = endTTimes_.at(index) - gridResolution_; // shift EV2_ to previous gridpoint
+ }
+ }
+ } // End of else statement
+ } // End of if for startbin == endbin
+ } // End of if deltaT < gridresolution
+
+ } // End of looping over paths3
+
+ } // end of ZHS-like approximation
+
+ // Feed start and end to the antenna
+ antenna.receive(startTTimes_.at(index), ReceiveVectorsStart_.at(index), EVstart_.at(index));
+ antenna.receive(endTTimes_.at(index), ReceiveVectorsEnd_.at(index), EVend_.at(index));
+
+ // update index
+ index = index + 1;
+
+ } // End of for loop for preDoppler factor (this includes checking for postDoppler factors)
+
+ } // End of checking of vector sizes
+
+ } // End of looping over the antennas.
+
+ } // End of simulate method.
+
/**
* Return the maximum step length for this particle and track.
@@ -218,7 +392,7 @@ namespace corsika {
// This is part of the ZHS / CoReas formalisms and can
// be related from the magnetic field / acceleration, charge,
// etc. of the particle.
- return 1000000000000;
+ return 1000000000000_m;
}
}; // END: class RadioProcess
diff --git a/corsika/modules/radio/RadioProcess.hpp b/corsika/modules/radio/RadioProcess.hpp
index 93cfb6911..4fecfc1b6 100644
--- a/corsika/modules/radio/RadioProcess.hpp
+++ b/corsika/modules/radio/RadioProcess.hpp
@@ -111,13 +111,14 @@ namespace corsika {
* we wait for the true output formatting to be ready.
**/
bool writeOutput() const {
-
+ // this for loop still has some issues
int i = 1;
for (auto& antenna : detector_.getAntennas()) {
auto [t,E] = antenna.getWaveform();
std::ofstream out_file ("antenna" + to_string(i) + "_output.csv");
- xt::dump_csv(out_file, t, E);
+ xt::dump_csv(out_file, t);
+ xt::dump_csv(out_file, E);
out_file.close();
++i;
diff --git a/corsika/modules/radio/ZHS.hpp b/corsika/modules/radio/ZHS.hpp
index 0aa1e0617..f9f25b25b 100644
--- a/corsika/modules/radio/ZHS.hpp
+++ b/corsika/modules/radio/ZHS.hpp
@@ -11,6 +11,7 @@
#include <corsika/modules/radio/propagators/StraightPropagator.hpp>
#include <corsika/framework/geometry/QuantityVector.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/modules/radio/propagators/SignalPath.hpp>
#include <bits/stdc++.h>
namespace corsika {
@@ -72,24 +73,26 @@ namespace corsika {
// get the Path from the track to the antenna
// This is a SignalPathCollection
- auto paths{this->propagator_.propagate(startPoint_, antenna.getLocation())};
- auto R_ {(startPoint_ - antenna.getLocation()).getNorm()};
+ auto paths{this->propagator_.propagate(startPoint_, antenna.getLocation(), 1_nm)};
// now loop over the paths that we got above
// Note: for the StraightPropagator, there will only be a single
// path but other propagators may return more than one.
for (auto const& path : paths) {
+ auto t1 = path.total_time_;
+
+ QuantityVector<ElectricFieldType::dimension_type> v11{10_V / 1_m, 10_V / 1_m, 10_V / 1_m};
// calculate the ZHS formalism for this particle-antenna
- ElectricFieldVector EV_ = ((- constants) * trackVelocity_.dot(path.emit)) *
- ((midTime_ + path.total_time_ - (1 - path.average_refractivity_ * beta_ *
- acos(track.getDirection(0).dot(path.emit_))) * startTime_)
- - (midTime_ + path.total_time_ - (1 - path.average_refractivity_ * beta_ *
- acos(track.getDirection(0).dot(path.emit_))) * endTime_))
- / (1 - path.average_refractivity_ * beta_ * acos(track.getDirection(0).dot(path.emit_)));
-
- // pass it to the antenna
- antenna.receive(midTime_ + path.total_time_, path.receive_, EV_);
+// ElectricFieldVector EV_ = ((- constants) * trackVelocity_.dot(path.emit)) *
+// ((midTime_ + path.total_time_ - (1 - path.average_refractivity_ * beta_ *
+// acos(track.getDirection(0).dot(path.emit_))) * startTime_)
+// - (midTime_ + path.total_time_ - (1 - path.average_refractivity_ * beta_ *
+// acos(track.getDirection(0).dot(path.emit_))) * endTime_))
+// / (1 - path.average_refractivity_ * beta_ * acos(track.getDirection(0).dot(path.emit_)));
+//
+// // pass it to the antenna
+// antenna.receive(midTime_ + path.total_time_, path.receive_, EV_);
} // END: loop over paths
diff --git a/corsika/modules/radio/antennas/Antenna.hpp b/corsika/modules/radio/antennas/Antenna.hpp
index 9aaabc365..41ebf33b0 100644
--- a/corsika/modules/radio/antennas/Antenna.hpp
+++ b/corsika/modules/radio/antennas/Antenna.hpp
@@ -34,7 +34,7 @@ namespace corsika {
// QuantityVector<MagneticFieldType::dimension_type>;
// a dimensionless vector used for the incident direction
- using Vector = QuantityVector<dimensionless_d>;
+// using Vector = QuantityVector<dimensionless_d>;
/**
* \brief Construct a base antenna instance.
diff --git a/corsika/modules/radio/antennas/TimeDomainAntenna.hpp b/corsika/modules/radio/antennas/TimeDomainAntenna.hpp
index 1462b23e1..4d095f87d 100644
--- a/corsika/modules/radio/antennas/TimeDomainAntenna.hpp
+++ b/corsika/modules/radio/antennas/TimeDomainAntenna.hpp
@@ -18,26 +18,25 @@ namespace corsika {
/**
* An implementation of a time-domain antenna that has a customized
- * start time, sampling period, and waveform duration.
+ * start time, sampling rate, and waveform duration.
*
*/
- class TimeDomainAntenna final : public Antenna<TimeDomainAntenna> {
+ class TimeDomainAntenna : public Antenna<TimeDomainAntenna> {
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 num_bins_; ///< The number of bins used.
- xt::xtensor<double,3> waveformE_; ///< The waveform stored by this antenna.
- std::pair<xt::xtensor<double, 1>,
- xt::xtensor<double,3>> waveform_; ///< useful for .getWaveform()
+ xt::xtensor<double,2> waveformE_; ///< The waveform stored by this antenna.
+ std::pair<xt::xtensor<double, 2>,
+ xt::xtensor<double,2>> waveform_; ///< useful for .getWaveform()
protected:
// expose the CRTP interfaces constructor
public:
// import the methods from the antenna
- using Antenna<TimeDomainAntenna>::Antenna;
- using Antenna<TimeDomainAntenna>::receive;
+
using Antenna<TimeDomainAntenna>::getName;
using Antenna<TimeDomainAntenna>::getLocation;
@@ -62,7 +61,7 @@ namespace corsika {
, duration_(duration)
, sample_rate_(sample_rate)
, num_bins_ (static_cast<int>(duration * sample_rate))
- , waveformE_ (xt::zeros<double>({3, num_bins_}))
+ , waveformE_ (xt::zeros<double>({num_bins_, 3}))
{};
/**
@@ -76,20 +75,20 @@ namespace corsika {
* @param field The incident electric field vector.
*
*/
- void receive(TimeType const time, Vector const& receive_vector,
+ void receive(TimeType const time, Vector<dimensionless_d> const& receive_vector,
ElectricFieldVector const& efield) {
- if (time < start_time_ || time > start_time_ + duration_) {
+ if (time < start_time_ || time >= start_time_ + duration_) {
return;
} else {
// figure out the correct timebin to store the E-field value.
auto timebin_ {static_cast<std::size_t>((time - start_time_) * sample_rate_)};
// store the x,y,z electric field components.
- waveformE_.at(0, timebin_) += (efield.getX().magnitude());
- waveformE_.at(1, timebin_) += (efield.getY().magnitude());
- waveformE_.at(2, timebin_) += (efield.getZ().magnitude());
-
+ waveformE_.at(timebin_, 0) += efield.getX().magnitude();
+ waveformE_.at(timebin_, 1) += efield.getY().magnitude();
+ waveformE_.at(timebin_, 2) += efield.getZ().magnitude();
+ //TODO: Check how they are stored in memory, row-wise or column-wise?
}
}
@@ -100,17 +99,17 @@ namespace corsika {
*
* @returns A pair of the sample times, and the field
*/
- std::pair<xt::xtensor<double, 1>, xt::xtensor<double,3>> getWaveform() const {
+ std::pair<xt::xtensor<double, 2>, xt::xtensor<double,2>> getWaveform() const {
+ // TODO: divide by Δt for CoREAS ONLY!
// create a 1-D xtensor to store time values so we can print them later.
- xt::xtensor<double, 1> times_ {xt::zeros<double>({num_bins_, 1})};
+ xt::xtensor<double, 2> times_ (xt::zeros<double>({num_bins_, 1}));
- for (auto i = 0; i <= num_bins_; i++) {
- times_.at(i) = static_cast<double>(start_time_ / 1_ns) +
- static_cast<double>(i * sample_rate_ * 1_ns);
+ for (int i = 0; i < num_bins_; i++) {
+ times_.at(i,0) = static_cast<double>(start_time_ / 1_s + i * sample_rate_ * 1_s);
}
- return waveform_;
+ return std::make_pair(times_, waveformE_);
};
/**
@@ -118,6 +117,7 @@ namespace corsika {
*/
void reset() {
waveformE_ = xt::zeros_like(waveformE_);
+// times_ = xt::zeros_like(times_);
};
}; // END: class Antenna final
diff --git a/corsika/modules/radio/detectors/RadioDetector.hpp b/corsika/modules/radio/detectors/RadioDetector.hpp
index 1d27a8ddd..804b8ccb2 100644
--- a/corsika/modules/radio/detectors/RadioDetector.hpp
+++ b/corsika/modules/radio/detectors/RadioDetector.hpp
@@ -37,7 +37,7 @@ namespace corsika {
*
* @returns An iterable mutable reference to the antennas.
*/
- std::vector<TAntennaImpl> const& getAntennas() { return antennas_; }
+ std::vector<TAntennaImpl> const& getAntennas() { return antennas_; } // maybe the const& here is an issue
/**
* Reset all the antenna waveforms.
diff --git a/tests/modules/testRadio.cpp b/tests/modules/testRadio.cpp
index 50f033e3f..afba63cba 100644
--- a/tests/modules/testRadio.cpp
+++ b/tests/modules/testRadio.cpp
@@ -8,12 +8,22 @@
#include <catch2/catch.hpp>
#include <corsika/modules/radio/ZHS.hpp>
+//#include <corsika/modules/radio/CoREAS.hpp>
#include <corsika/modules/radio/antennas/TimeDomainAntenna.hpp>
-#include <corsika/modules/radio/detectors/TimeDomainDetector.hpp>
+#include <corsika/modules/radio/detectors/RadioDetector.hpp>
#include <corsika/modules/radio/propagators/StraightPropagator.hpp>
#include <corsika/modules/radio/propagators/SignalPath.hpp>
#include <corsika/modules/radio/propagators/RadioPropagator.hpp>
+#include <vector>
+#include <xtensor/xtensor.hpp>
+#include <xtensor/xbuilder.hpp>
+#include <xtensor/xio.hpp>
+#include <xtensor/xcsv.hpp>
+#include <istream>
+#include <fstream>
+#include <iostream>
+
#include <corsika/media/Environment.hpp>
#include <corsika/media/HomogeneousMedium.hpp>
@@ -39,45 +49,158 @@ double constexpr absMargin = 1.0e-7;
TEST_CASE("Radio", "[processes]") {
-// logging::set_level(logging::level::info);
-// corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
+ logging::set_level(logging::level::info);
+ corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
// check that I can create and reset a TimeDomain process
- SECTION("TimeDomainDetector") {
-
- // construct a time domain detector
- TimeDomainDetector<TimeDomainAntenna> detector;
-
- // // and construct some time domain radio process
- // TimeDomain<ZHS<>, TimeDomainDetector<TimeDomainAntenna>> process(detector);
- // TimeDomain<ZHS<CPU>, TimeDomainDetector<TimeDomainAntenna>> process2(detector);
- // TimeDomain<ZHS<CPU>, decltype(detector)> process3(detector);
- }
+// SECTION("TimeDomainDetector") {
+//
+// // construct a time domain detector
+// AntennaCollection<TimeDomainAntenna> detector;
+//
+// // // and construct some time domain radio process
+// // TimeDomain<ZHS<>, TimeDomainDetector<TimeDomainAntenna>> process(detector);
+// // TimeDomain<ZHS<CPU>, TimeDomainDetector<TimeDomainAntenna>> process2(detector);
+// // TimeDomain<ZHS<CPU>, decltype(detector)> process3(detector);
+// }
// check that I can create time domain antennas
- SECTION("TimeDomainDetector") {
+ SECTION("TimeDomainAntenna") {
// create an environment so we can get a coordinate system
Environment<IMediumModel> env;
// the antenna location
- const auto point{Point(env.getCoordinateSystem(), 1_m, 2_m, 3_m)};
+ const auto point1{Point(env.getCoordinateSystem(), 1_m, 2_m, 3_m)};
+ const auto point2{Point(env.getCoordinateSystem(), 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};
// check that I can create an antenna at (1, 2, 3)
- const auto ant{TimeDomainAntenna("antenna_name", point)};
+ TimeDomainAntenna ant1("antenna_name", point1, t1, t2, t3);
+ TimeDomainAntenna ant2("antenna_name2", point2, 11_s, t2, t3);
// assert that the antenna name is correct
- REQUIRE(ant.getName() == "antenna_name");
+ REQUIRE(ant1.getName() == "antenna_name");
// and check that the antenna is at the right location
- REQUIRE((ant.getLocation() - point).getNorm() < 1e-12 * 1_m);
+ REQUIRE((ant1.getLocation() - point1).getNorm() < 1e-12 * 1_m);
// construct a radio detector instance to store our antennas
- TimeDomainDetector<TimeDomainAntenna> detector;
+ AntennaCollection<TimeDomainAntenna> detector;
// add this antenna to the process
- detector.addAntenna(ant);
+ detector.addAntenna(ant1);
+ detector.addAntenna(ant2);
+
+ // create an environment with uniform refractive index of 1
+ using UniRIndex =
+ UniformRefractiveIndex<HomogeneousMedium<IRefractiveIndexModel<IMediumModel>>>;
+
+ using EnvType = Environment<IRefractiveIndexModel<IMediumModel>>;
+ EnvType env6;
+
+ // get a coordinate system
+ const CoordinateSystemPtr rootCS6 = env6.getCoordinateSystem();
+
+ auto Medium6 = EnvType::createNode<Sphere>(
+ Point{rootCS6, 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(
+ std::vector<Code>{Code::Nitrogen},
+ std::vector<float>{1.f}));
+
+ env6.getUniverse()->addChild(std::move(Medium6));
+
+
+ // get a unit vector
+ Vector<dimensionless_d> v1(rootCS6, {0, 0, 1});
+ QuantityVector<ElectricFieldType::dimension_type> v11{10_V / 1_m, 10_V / 1_m, 10_V / 1_m};
+
+ Vector<dimensionless_d> v2(rootCS6, {0, 1, 0});
+ QuantityVector<ElectricFieldType::dimension_type> v22{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 getWaveform() method
+ auto [t11, E1] = ant1.getWaveform();
+ CHECK(E1(5,0) - 10 == 0);
+
+ auto [t222, E2] = ant2.getWaveform();
+ CHECK(E2(5,0) -20 == 0);
+
+
+ // the rest is just random ideas
+ //////////////////////////////////////////////////////////////////////////////////////
+
+// for (auto& kostas : detector.getAntennas()) {
+// std::cout << kostas.getName() << std::endl;
+// kostas.receive(15_s, v1, v11);
+// std::cout << kostas.waveformE_;
+// auto [t, E] = kostas.getWaveform();
+// std::cout << E(5,0) << std::endl;
+// kostas.receive(16_s, v2, v22);
+// std::cout << E(5,0) << std::endl;
+// }
+
+
+
+ // auto t33{static_cast<int>(t1 / t2)};
+// std::cout << t33 << std::endl;
+// xt::xtensor<double,2> waveformE_ = xt::zeros<double>({2, 3});
+// xt::xtensor<double,2> res = xt::view(waveformE_);
+// std::cout << ant1.waveformE_ << std::endl;
+// std::cout << " " << std::endl;
+// std::cout << ant2.waveformE_ << std::endl;
+// std::cout << " " << std::endl;
+// std::cout << ant1.times_ << std::endl;
+// std::cout << " " << std::endl;
+// std::cout << ant2.times_ << std::endl;
+// auto [t, E] = ant2.getWaveform();
+//
+
+//
+//
+// std::ofstream out_file("antenna_output.csv");
+// xt::dump_csv(out_file, t);
+// xt::dump_csv(out_file, E);
+// out_file.close();
+
+
+// for (auto& antenna : detector.getAntennas()) {
+// auto [t, E] = antenna.getWaveform();
+// }
+
+// int i = 1;
+// auto x = detector.getAntennas();
+// for (auto& antenna : x) {
+//
+// auto [t, E] = antenna.getWaveform();
+// std::ofstream out_file("antenna" + to_string(i) + "_output.csv");
+// std::cout << E(5,0) << std::endl;
+// xt::dump_csv(out_file, t);
+// xt::dump_csv(out_file, E);
+// out_file.close();
+// ++i;
+// }
+
+
+// detector.writeOutput();
+// ant1.reset();
+// ant2.reset();
+//
+// std::cout << ant1.waveformE_ << std::endl;
+// std::cout << ant2.waveformE_ << std::endl;
+
+ //////////////////////////////////////////////////////////////////////////////////////
}
// check that I can create working Straight Propagators in different environments
@@ -133,9 +256,10 @@ TEST_CASE("Radio", "[processes]") {
for (auto const& path : paths_) {
CHECK((path.total_time_ / 1_s) - ((34_m / (3 * constants::c)) / 1_s) ==
Approx(0).margin(absMargin));
- CHECK(path.average_refractivity_ == Approx(1));
+ CHECK(path.average_refractive_index_ == Approx(1));
CHECK(path.emit_.getComponents() == v1.getComponents());
CHECK(path.receive_.getComponents() == v1.getComponents());
+ CHECK(path.R_distance_ == 10_m);
// CHECK(std::equal(P1.begin(), P1.end(), path.points_.begin(),[]
// (Point a, Point b) { return (a - b).norm() / 1_m < 1e-5;}));
//TODO:THINK ABOUT THE POINTS IN THE SIGNALPATH.H
@@ -198,9 +322,10 @@ TEST_CASE("Radio", "[processes]") {
for (auto const& path :paths1_) {
CHECK( (path.total_time_ / 1_s) - ((34_m / (3 * constants::c)) / 1_s)
== Approx(0).margin(absMargin) );
- CHECK( path.average_refractivity_ == Approx(1) );
+ CHECK( path.average_refractive_index_ == Approx(1) );
CHECK( path.emit_.getComponents() == vv1.getComponents() );
CHECK( path.receive_.getComponents() == vv1.getComponents() );
+ CHECK( path.R_distance_ == 10_m );
}
CHECK( paths1_.size() == 1 );
@@ -249,9 +374,10 @@ TEST_CASE("Radio", "[processes]") {
for (auto const& path :paths2_) {
CHECK( (path.total_time_ / 1_s) - ((3.177511688_m / (3 * constants::c)) / 1_s)
== Approx(0).margin(absMargin) );
- CHECK( path.average_refractivity_ == Approx(0.210275935) );
+ CHECK( path.average_refractive_index_ == Approx(0.210275935) );
CHECK( path.emit_.getComponents() == vvv1.getComponents() );
CHECK( path.receive_.getComponents() == vvv1.getComponents() );
+ CHECK( path.R_distance_ == 10_m );
}
CHECK( paths2_.size() == 1 );
@@ -262,6 +388,115 @@ TEST_CASE("Radio", "[processes]") {
//
// // TODO: construct the environment for the propagator
+ //////////////////////////////////////////////////////////////////////////////////////
+// // useful information
+// // create an environment so we can get a coordinate system
+// environment::Environment<environment::IMediumModel> env;
+//
+// // the antenna location
+// const auto point{geometry::Point(env.GetCoordinateSystem(), 1_m, 2_m, 3_m)};
+//
+// // check that I can create an antenna at (1, 2, 3)
+// const auto ant{TimeDomainAntenna("antenna_name", point)};
+//
+// // assert that the antenna name is correct
+// REQUIRE(ant.GetName() == "antenna_name");
+//
+// // and check that the antenna is at the right location
+// REQUIRE((ant.GetLocation() - point).norm() < 1e-12 * 1_m);
+//
+// // construct a radio detector instance to store our antennas
+// TimeDomainDetector<TimeDomainAntenna> detector;
+//
+// // add this antenna to the process
+// detector.AddAntenna(ant);
+//
+// // create a TimeDomain process
+// auto zhs{TimeDomain<ZHS<CPU>, decltype(detector)>(detector)};
+//
+// // get the antenna back from the process and check the properties
+// REQUIRE(detector.GetAntennas().cbegin()->GetName() == "antenna_name");
+//
+// // and check the location is the same
+// REQUIRE((detector.GetAntennas().cbegin()->GetLocation() - point).norm() <
+// 1e-12 * 1_m);
+//
+// // and ANOTHER antenna
+// const auto ant2{TimeDomainAntenna("antenna_name", point)};
+// detector.AddAntenna(ant2);
+// // and check that the number of antennas visible to the process has increased
+// REQUIRE(zhs.GetDetector().GetAntennas().size() == 2);
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+// using VelocityVec = Vector<corsika::units::si::SpeedType::dimension_type>;
+//
+// // setup environment, geometry
+// environment::Environment<environment::IMediumModel> env;
+//
+// // and get the coordinate systm
+// geometry::CoordinateSystem const& cs{env.GetCoordinateSystem()};
+//
+// // a test particle
+// const auto pcode{particles::Code::Electron};
+// const auto mass{particles::Electron::GetMass()};
+//
+// // create a new stack for each trial
+// setup::Stack stack;
+//
+// // construct an energy
+// const HEPEnergyType E0{1_TeV};
+//
+// // compute the necessary momentumn
+// const HEPMomentumType P0{sqrt(E0 * E0 - mass * mass)};
+//
+// // and create the momentum vector
+// const auto plab{corsika::stack::MomentumVector(cs, {0_GeV, 0_GeV, P0})};
+//
+// // and create the location of the particle in this coordinate system
+// const geometry::Point pos(cs, 5_m, 1_m, 8_m);
+//
+// // add the particle to the stack
+// auto particle{stack.AddParticle(std::make_tuple(pcode, E0, plab, pos, 0_ns))};
+//
+// // get a view into the stack
+// corsika::stack::SecondaryView stackview(particle);
+//
+// // create a trajectory
+// const auto& track{
+// Trajectory(Line(pos, VelocityVec(cs, 0_m / 1_s, 0_m / 1_s, 0.1_m / 1_ns)), 1_ns)};
+//
+// // and create a view vector
+// const auto& view{Vector(cs, 0_m, 0_m, 1_m)};
+//
+// // construct a radio detector instance to store our antennas
+// TimeDomainDetector<TimeDomainAntenna> detector;
+//
+// // the antenna location
+// const auto point{geometry::Point(cs, 1_m, 2_m, 3_m)};
+//
+// // check that I can create an antenna at (1, 2, 3)
+// const auto ant{TimeDomainAntenna("antenna_name", point)};
+//
+// // add this antenna to the process
+// detector.AddAntenna(ant);
+//
+// // create a TimeDomain process
+// auto zhs{TimeDomain<ZHS<CPU>, decltype(detector)>(detector)};
+//
+// // get the projectile
+// auto projectile{stackview.GetProjectile()};
+//
+// // call ZHS over the track
+// zhs.DoContinuous(projectile, track);
+// zhs.Simulate(projectile, track);
+//
+// // try and call the ZHS implementation directly for a given view direction
+// ZHS<CPU>::Emit(projectile, track, view);
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+
// create an environment with uniform refractive index of 1
// using UniRIndex =
// UniformRefractiveIndex<HomogeneousMedium<IRefractiveIndexModel<IMediumModel>>>;
--
GitLab