diff --git a/corsika/modules/radio/CoREAS.hpp b/corsika/modules/radio/CoREAS.hpp
index 7730d07db31b04f775b5333f979901d1532078b4..627ec20631adf9be7bcd88c1bc46ed99ec98ef1f 100755
--- a/corsika/modules/radio/CoREAS.hpp
+++ b/corsika/modules/radio/CoREAS.hpp
@@ -52,8 +52,11 @@ namespace corsika {
ProcessReturn simulate(Particle& particle, Track const& 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()};
+ auto startTime_ {particle.getTime()}; // time at the start point of the track hopefully.
+ std::cout << "startTime_: " << startTime_ << std::endl;
+ auto endTime_ {particle.getTime() + track.getDuration()};
+ std::cout << "endTime_: " << endTime_ << std::endl;
+ std::cout << "track.getDuration(): " << track.getDuration() << std::endl;
// gamma factor is calculated using beta
// auto startGamma_ {1. / sqrt(1. - (startBeta_ * startBeta_))};
@@ -61,12 +64,15 @@ namespace corsika {
// get start and end position of the track
auto startPoint_ {track.getPosition(0)};
+ std::cout << "EDO EINAI I ARXI: " << startPoint_ << std::endl;
auto endPoint_ {track.getPosition(1)};
+ std::cout << "EDO EINAI TO TELOS: " << endPoint_ << std::endl;
// beta is velocity / speed of light. Start & end should be the same!
// auto beta_ {static_cast<double>((endPoint_.distance_to(startPoint_) / 1_m ) / (endTime_/ 1_s - startTime_/ 1_s) )};
// auto beta_ {((endPoint_.distance_to(startPoint_)) / (endTime_ - startTime_)).normalized()};
auto beta_ {((endPoint_ - startPoint_) / (constants::c * (endTime_ - startTime_))).normalized()};
+ std::cout << "BETA_: " << beta_ << std::endl;
// get particle charge
auto const charge_ {get_charge(particle.getPID())};
@@ -89,19 +95,27 @@ namespace corsika {
// get the Path (path1) from the start "endpoint" to the antenna.
// This is a Signal Path Collection
- auto paths1 {this->propagator_.propagate(startPoint_, antenna.getLocation(), 1_nm)}; // TODO: Need to add the stepsize to .propagate()!!!!
+ auto paths1 {this->propagator_.propagate(startPoint_, antenna.getLocation(), 1_m)}; // TODO: Need to add the stepsize to .propagate()!!!!
// now loop over the paths for startpoint that we got above
for (auto const& path : paths1) {
// calculate preDoppler factor
- auto preDoppler_{1. - path.average_refractive_index_ * beta_.dot(path.emit_)};
+ double preDoppler_{1. - path.average_refractive_index_ * beta_.dot(path.emit_)};
+
+ std::cout << "Path Emit PRE: " << beta_.dot(path.emit_) << std::endl;
+
+ std::cout << "preDoppler: " << preDoppler_<< std::endl;
+ long double preDoppler_1{1. - path.average_refractive_index_ * beta_.dot(path.emit_)};
+ std::cout << "preDoppler1: " << preDoppler_1<< std::endl;
+
// 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_};
+ std::cout << "START RECEIVE TIME: " << startPointReceiveTime_ << std::endl;
// store it to startTTimes_ std::vector for later use
startTTimes_.push_back(startPointReceiveTime_);
@@ -112,7 +126,7 @@ namespace corsika {
// calculate electric field vector for startpoint
ElectricFieldVector EV1_=
path.receive_.cross(path.receive_.cross(beta_)).getComponents() /
- (path.R_distance_ * preDoppler_) * (1 / (4 * M_PI * 1_s)) * ((1 / constants::epsilonZero) * (1 / constants::c)) * charge_;
+ (path.R_distance_ * preDoppler_) * (1 / (4 * M_PI * track.getDuration())) * ((1 / constants::epsilonZero) * (1 / constants::c)) * charge_;
// store it to EVstart_ std::vector for later use
EVstart_.push_back(EV1_);
@@ -121,7 +135,7 @@ namespace corsika {
// get the Path (path2) from the end "endpoint" to the antenna.
// This is a SignalPathCollection
- auto paths2 {this->propagator_.propagate(endPoint_, antenna.getLocation(), 1_nm)};
+ auto paths2 {this->propagator_.propagate(endPoint_, antenna.getLocation(), 1_m)};
// now loop over the paths for endpoint that we got above
for (auto const& path : paths2) {
@@ -129,11 +143,20 @@ namespace corsika {
double postDoppler_{1. - path.average_refractive_index_ *
beta_.dot(path.emit_)}; // maybe this is path.receive_ (?)
+ std::cout << "Path Emit POST: " << beta_.dot(path.emit_) << std::endl;
+
+ std::cout << "postDoppler: " << postDoppler_<< std::endl;
+ long double postDoppler_1{1. - path.average_refractive_index_ *
+ beta_.dot(path.emit_)}; // maybe this is path.receive_ (?)
+ std::cout << "postDoppler1: " << postDoppler_1 << std::endl;
+
// 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_};
+ std::cout << "END RECEIVE TIME: " << endPointReceiveTime_ << std::endl;
+
// store it to endTTimes_ std::vector for later use
endTTimes_.push_back(endPointReceiveTime_);
@@ -144,7 +167,7 @@ namespace corsika {
// calculate electric field vector for endpoint
ElectricFieldVector EV2_=
path.receive_.cross(path.receive_.cross(beta_)).getComponents() /
- (path.R_distance_ * postDoppler_) * (1 / (4 * M_PI * 1_s)) * ((1 / constants::epsilonZero) * (1 / constants::c)) * charge_;
+ (path.R_distance_ * postDoppler_) * ((-1) / (4 * M_PI * track.getDuration())) * ((1 / constants::epsilonZero) * (1 / constants::c)) * charge_;
// store it to EVstart_ std::vector for later use
EVend_.push_back(EV2_);
@@ -158,7 +181,10 @@ namespace corsika {
// use this to access different elements of std::vectors
std::size_t index = 0;
for (auto& preDoppler__ : preDoppler) {
+ std::cout << "preDoppler__: " << preDoppler__ << std::endl;
+ std::cout << "postDoppler.at(index): " << postDoppler.at(index) << std::endl;
+ // TODO: We need to check this. Something funny is happening here.
// 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)) {
@@ -244,7 +270,7 @@ namespace corsika {
// get the Path (path3) from the middle "endpoint" to the antenna.
// This is a SignalPathCollection
- auto paths3{this->propagator_.propagate(midPoint_, antenna.getLocation(), 1_nm)};
+ auto paths3{this->propagator_.propagate(midPoint_, antenna.getLocation(), 1_m)};
// 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
@@ -263,11 +289,7 @@ namespace corsika {
// CoREAS calculation -> get ElectricFieldVector3 for "midPoint"
ElectricFieldVector EVmid_ =
path.receive_.cross(path.receive_.cross(beta_)).getComponents() /
- (path.R_distance_ * midDoppler_) * (1 / (4 * M_PI * 1_s)) * ((1 / constants::epsilonZero) * (1 / constants::c)) * charge_;
-
-// ElectricFieldVector EVmid2_ = (- charge_ / constants::c) *
-// path.receive_.cross(path.receive_.cross(beta_)) /
-// (path.R_distance_ * midDoppler_);
+ (path.R_distance_ * midDoppler_) * (1 / (4 * M_PI * track.getDuration())) * ((1 / constants::epsilonZero) * (1 / constants::c)) * charge_;
// 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
@@ -362,6 +384,9 @@ namespace corsika {
} // end of ZHS-like approximation
// Feed start and end to the antenna
+ std::cout << "LIGO PRIN TO RECEIVE :" << std::endl;
+ std::cout << "startTTimes_.at(index): " << startTTimes_.at(index) << std::endl;
+ std::cout << "endTTimes_.at(index): " << endTTimes_.at(index) << std::endl;
antenna.receive(startTTimes_.at(index), ReceiveVectorsStart_.at(index), EVstart_.at(index));
antenna.receive(endTTimes_.at(index), ReceiveVectorsEnd_.at(index), EVend_.at(index));
@@ -370,6 +395,7 @@ namespace corsika {
} // End of for loop for preDoppler factor (this includes checking for postDoppler factors)
+ index = 0;
} // End of checking of vector sizes
} // End of looping over the antennas.
diff --git a/corsika/modules/radio/antennas/TimeDomainAntenna.hpp b/corsika/modules/radio/antennas/TimeDomainAntenna.hpp
index 873a3903b7db8bd0998e7c00468c0da84fe2aa0d..23647895710ad3f4a7691993fd6c34d291fe4cd5 100644
--- a/corsika/modules/radio/antennas/TimeDomainAntenna.hpp
+++ b/corsika/modules/radio/antennas/TimeDomainAntenna.hpp
@@ -94,6 +94,7 @@ namespace corsika {
} else {
// figure out the correct timebin to store the E-field value.
auto timebin_ {static_cast<std::size_t>((time - start_time_) * sample_rate_)};
+ std::cout << "TIMEBIN IS: " << timebin_ << std::endl;
// store the x,y,z electric field components.
waveformE_.at(timebin_, 0) += efield.getX().magnitude();
diff --git a/corsika/modules/radio/detectors/RadioDetector.hpp b/corsika/modules/radio/detectors/RadioDetector.hpp
index 3a8f9e6c74d44be93e860fefcd032c5999f1e4e4..38173fcd6ca308c9dd463a7fe45fcf3df246067a 100644
--- a/corsika/modules/radio/detectors/RadioDetector.hpp
+++ b/corsika/modules/radio/detectors/RadioDetector.hpp
@@ -30,7 +30,7 @@ namespace corsika {
*
* @param antenna The antenna to add
*/
- void addAntenna(TAntennaImpl const antenna) { antennas_.push_back(antenna); }
+ void addAntenna(TAntennaImpl const& antenna) { antennas_.push_back(antenna); }
/**
* Get the specific antenna at that place in the collection
diff --git a/examples/radio_shower.cpp b/examples/radio_shower.cpp
index e9612b4e36c2c25b712f8c120ebdccecbb8a9c64..f38aa4f3ce68116fa8eb0ebd9c72d8ab30f0625e 100644
--- a/examples/radio_shower.cpp
+++ b/examples/radio_shower.cpp
@@ -51,6 +51,7 @@
#include <corsika/modules/radio/RadioProcess.hpp>
#include <corsika/modules/radio/CoREAS.hpp>
+#include <corsika/modules/radio/antennas/Antenna.hpp>
#include <corsika/modules/radio/antennas/TimeDomainAntenna.hpp>
#include <corsika/modules/radio/detectors/RadioDetector.hpp>
#include <corsika/modules/radio/propagators/StraightPropagator.hpp>
@@ -248,7 +249,8 @@ int main(int argc, char** argv) {
corsika::proposal::ContinuousProcess emContinuous(env);
InteractionCounter emCascadeCounted(emCascade);
// put radio here
- CoREAS<decltype(detector), decltype(StraightPropagator(env))> coreas(detector, env);
+ RadioProcess<decltype(detector), CoREAS<decltype(detector), decltype(StraightPropagator(envCoREAS))>, decltype(StraightPropagator(envCoREAS))>
+ coreas(detector, envCoREAS);
OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
TrackWriter trackWriter("tracks.dat");
diff --git a/tests/modules/testRadio.cpp b/tests/modules/testRadio.cpp
index 6064495d4ce11c5f1a716425c24f0bae03691e38..19bfa885caf9b5aa802b705b18782f05c3805e9d 100644
--- a/tests/modules/testRadio.cpp
+++ b/tests/modules/testRadio.cpp
@@ -89,34 +89,37 @@ TEST_CASE("Radio", "[processes]") {
// now create antennas and detectors
// the antennas location
- const auto point1{Point(envCoREAS.getCoordinateSystem(), 1_m, 2_m, 3_m)};
- const auto point2{Point(envCoREAS.getCoordinateSystem(), 4_m, 5_m, 6_m)};
+ 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
- const TimeType t1{10_s};
- const TimeType t2{10_s};
+ const TimeType t1{0_s};
+ const TimeType t2{100_s};
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, t1, t2, t3);
- TimeDomainAntenna ant2("antenna_name2", point2, t4, t2, t3);
+ TimeDomainAntenna ant2("antenna_name2", point2, t1, t2, t3);
// construct a radio detector instance to store our antennas
AntennaCollection<TimeDomainAntenna> detector;
-// std::vector<TimeDomainAntenna> detector;
+ // add the antennas to the detector
detector.addAntenna(ant1);
-// detector.push_back(ant1);
+ detector.addAntenna(ant2);
+
// create a particle
auto const particle{Code::Electron};
+// auto const particle{Code::Gamma};
const auto pmass{get_mass(particle)};
- VelocityVector v0(rootCSCoREAS, {5_m / second, 5_m / second, 5_m / second});
+ VelocityVector v0(rootCSCoREAS, {5e+2_m / second, 5e+2_m / second, 5e+2_m / second});
Vector B0(rootCSCoREAS, 5_T, 5_T, 5_T);
@@ -124,8 +127,8 @@ TEST_CASE("Radio", "[processes]") {
auto const k{1_m * ((1_m) / ((1_s * 1_s) * 1_V))};
- auto const t = 10_s;
- LeapFrogTrajectory base(point1, v0, B0, k, t);
+ auto const t = 10000_ms;
+ LeapFrogTrajectory base(point4, v0, B0, k, t);
// create a new stack for each trial
setup::Stack stack;
@@ -146,351 +149,362 @@ TEST_CASE("Radio", "[processes]") {
auto const particle1{stack.addParticle(std::make_tuple(particle, E0, plab, pos, 0_ns))};
auto const charge_ {get_charge(particle1.getPID())};
- std::cout << "charge: " << charge_ << std::endl;
- std::cout << "1 / c: " << 1. / constants::c << std::endl;
+// std::cout << "charge: " << charge_ << std::endl;
+// std::cout << "1 / c: " << 1. / constants::c << std::endl;
// set up a track object
// setup::Tracking tracking;
- auto startPoint_ {base.getPosition(0)};
- auto midPoint_ {base.getPosition(0.5)};
- auto endPoint_ {base.getPosition(1)};
- std::cout << "startPoint_: " << startPoint_ << std::endl;
- std::cout << "midPoint_: " << midPoint_ << std::endl;
- std::cout << "endPoint_: " << endPoint_ << std::endl;
-
- auto velo_ {base.getVelocity(0)};
- std::cout << "velocity: " << velo_ << std::endl;
+// auto startPoint_ {base.getPosition(0)};
+// auto midPoint_ {base.getPosition(0.5)};
+// auto endPoint_ {base.getPosition(1)};
+// std::cout << "startPoint_: " << startPoint_ << std::endl;
+// std::cout << "midPoint_: " << midPoint_ << std::endl;
+// std::cout << "endPoint_: " << endPoint_ << std::endl;
- auto startTime_ {particle1.getTime() - base.getDuration()}; // time at the start point of the track hopefully.
- auto endTime_ {particle1.getTime()};
- std::cout << "startTime_: " << startTime_ << std::endl;
- std::cout << "endTime_: " << endTime_ << std::endl;
+// auto velo_ {base.getVelocity(0)};
+// std::cout << "velocity: " << velo_ << std::endl;
- auto beta_ {((endPoint_ - startPoint_) / (constants::c * (endTime_ - startTime_))).normalized()};
- std::cout << "beta_: " << beta_ << std::endl;
+// auto startTime_ {particle1.getTime() - base.getDuration()}; // time at the start point of the track hopefully.
+// auto endTime_ {particle1.getTime()};
+// std::cout << "startTime_: " << startTime_ << std::endl;
+// std::cout << "endTime_: " << endTime_ << std::endl;
+//
+// auto beta_ {((endPoint_ - startPoint_) / (constants::c * (endTime_ - startTime_))).normalized()};
+// std::cout << "beta_: " << beta_ << std::endl;
- Vector<dimensionless_d> v1(rootCSCoREAS, {0, 0, 1});
- std::cout << "v1: " << v1.getComponents() << std::endl;
+// Vector<dimensionless_d> v1(rootCSCoREAS, {0, 0, 1});
+// std::cout << "v1: " << v1.getComponents() << std::endl;
+//
+// std::cout << "beta_.dot(v1): " << beta_.dot(v1) << std::endl;
+//
+// std::cout << "Pi: " << 1/M_PI << std::endl;
+//
+// std::cout << "speed of light: " << constants::c << std::endl;
+//
+// std::cout << "vacuum permitivity: " << constants::epsilonZero << std::endl;
- std::cout << "beta_.dot(v1): " << beta_.dot(v1) << std::endl;
+ // Create a CoREAS instance
+// CoREAS<decltype(detector), decltype(StraightPropagator(envCoREAS))> coreas1(detector, envCoREAS);
- // Create a ZHS instance
- CoREAS<decltype(detector), decltype(StraightPropagator(envCoREAS))> coreas(detector, envCoREAS);
+ // create a radio process instance using CoREAS
+ RadioProcess<decltype(detector), CoREAS<decltype(detector), decltype(StraightPropagator(envCoREAS))>, decltype(StraightPropagator(envCoREAS))>
+ coreas(detector, envCoREAS);
- // call CoREAS over the track
-// coreas.doContinuous(particle1, base);
-// coreas.simulate(particle1, base);
+ // check doContinuous and simulate methods
+ coreas.doContinuous(particle1, base);
+// coreas1.simulate(particle1, base);
-// coreas.writeOutput();
+ // check writeOutput method -> should produce 2 csv files for each antenna
+ coreas.writeOutput();
}
- SECTION("TimeDomainAntenna") {
-
- // create an environment so we can get a coordinate system
- using EnvType = Environment<IRefractiveIndexModel<IMediumModel>>;
- EnvType env6;
-
- using UniRIndex =
- UniformRefractiveIndex<HomogeneousMedium<IRefractiveIndexModel<IMediumModel>>>;
-
-
- // the antenna location
- const auto point1{Point(env6.getCoordinateSystem(), 1_m, 2_m, 3_m)};
- const auto point2{Point(env6.getCoordinateSystem(), 4_m, 5_m, 6_m)};
-
-
- // 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));
-
- // 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};
-
- // check that I can create an antenna at (1, 2, 3)
- TimeDomainAntenna ant1("antenna_name", point1, t1, t2, t3);
- TimeDomainAntenna ant2("antenna_name2", point2, t4, t2, t3);
-
- // 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});
- 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 [t111, E1] = ant1.getWaveform();
- CHECK(E1(5,0) - 10 == 0);
- auto [t222, E2] = ant2.getWaveform();
- CHECK(E2(5,0) -20 == 0);
-
- // use the receive method in a for loop. It works now!
- for (auto& xx : detector.getAntennas()) {
- xx.receive(15_s, v1, v11);
- }
-
- // t & E are correct! uncomment to see them
+// SECTION("TimeDomainAntenna") {
+//
+// // create an environment so we can get a coordinate system
+// using EnvType = Environment<IRefractiveIndexModel<IMediumModel>>;
+// EnvType env6;
+//
+// using UniRIndex =
+// UniformRefractiveIndex<HomogeneousMedium<IRefractiveIndexModel<IMediumModel>>>;
+//
+//
+// // the antenna location
+// const auto point1{Point(env6.getCoordinateSystem(), 1_m, 2_m, 3_m)};
+// const auto point2{Point(env6.getCoordinateSystem(), 4_m, 5_m, 6_m)};
+//
+//
+// // 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));
+//
+// // 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};
+//
+// // check that I can create an antenna at (1, 2, 3)
+// TimeDomainAntenna ant1("antenna_name", point1, t1, t2, t3);
+// TimeDomainAntenna ant2("antenna_name2", point2, t4, t2, t3);
+//
+// // 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});
+// 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 [t111, E1] = ant1.getWaveform();
+// CHECK(E1(5,0) - 10 == 0);
+// auto [t222, E2] = ant2.getWaveform();
+// CHECK(E2(5,0) -20 == 0);
+//
+// // use the receive method in a for loop. It works now!
// for (auto& xx : detector.getAntennas()) {
-// auto [t,E] = xx.getWaveform();
-// std::cout << t << std::endl;
-// std::cout << " ... "<< std::endl;
-// std::cout << E << std::endl;
-// std::cout << " ... "<< std::endl;
+// xx.receive(15_s, v1, v11);
// }
-
- // check output files. It works, uncomment to see.
-// int i = 1;
-// for (auto& antenna : detector.getAntennas()) {
//
-// auto [t,E] = antenna.getWaveform();
-// auto c0 = xt::hstack(xt::xtuple(t,E));
-// std::ofstream out_file ("antenna" + to_string(i) + "_output.csv");
-// xt::dump_csv(out_file, c0);
-// out_file.close();
-// ++i;
+// // t & E are correct! uncomment to see them
+//// for (auto& xx : detector.getAntennas()) {
+//// auto [t,E] = xx.getWaveform();
+//// std::cout << t << std::endl;
+//// std::cout << " ... "<< std::endl;
+//// std::cout << E << std::endl;
+//// std::cout << " ... "<< std::endl;
+//// }
+//
+// // check output files. It works, uncomment to see.
+//// int i = 1;
+//// for (auto& antenna : detector.getAntennas()) {
+////
+//// auto [t,E] = antenna.getWaveform();
+//// auto c0 = xt::hstack(xt::xtuple(t,E));
+//// std::ofstream out_file ("antenna" + to_string(i) + "_output.csv");
+//// xt::dump_csv(out_file, c0);
+//// out_file.close();
+//// ++i;
+////
+//// }
+//
+// // check reset method for antennas. Uncomment to see they are zero
+//// ant1.reset();
+//// ant2.reset();
+////
+//// std::cout << ant1.waveformE_ << std::endl;
+//// std::cout << ant2.waveformE_ << std::endl;
+////
+//// std::cout << " ... "<< std::endl;
+//// std::cout << " ... "<< std::endl;
+//
+// // check reset method for antenna collection. Uncomment to see they are zero
+//// detector.reset();
+//// for (auto& xx : detector.getAntennas()) {
+//// std::cout << xx.waveformE_ << std::endl;
+//// std::cout << " ... "<< std::endl;
+//// }
+//
+//
+// }
+//
+// // check that I can create working Straight Propagators in different environments
+// SECTION("Straight Propagator w/ Uniform Refractive Index") {
+//
+// // create an environment with uniform refractive index of 1
+// using UniRIndex =
+// UniformRefractiveIndex<HomogeneousMedium<IRefractiveIndexModel<IMediumModel>>>;
+//
+// using EnvType = Environment<IRefractiveIndexModel<IMediumModel>>;
+// EnvType env;
+//
+// // get a coordinate system
+// const CoordinateSystemPtr rootCS = env.getCoordinateSystem();
+//
+// auto Medium = EnvType::createNode<Sphere>(
+// Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
+//
+// auto const props = Medium->setModelProperties<UniRIndex>(
+// 1, 1_kg / (1_m * 1_m * 1_m),
+// NuclearComposition(
+// std::vector<Code>{Code::Nitrogen},
+// std::vector<float>{1.f}));
+//
+// env.getUniverse()->addChild(std::move(Medium));
+//
+// // get some points
+// Point p0(rootCS, {0_m, 0_m, 0_m});
+// // Point p1(rootCS, {0_m, 0_m, 1_m});
+// // Point p2(rootCS, {0_m, 0_m, 2_m});
+// // Point p3(rootCS, {0_m, 0_m, 3_m});
+// // Point p4(rootCS, {0_m, 0_m, 4_m});
+// // Point p5(rootCS, {0_m, 0_m, 5_m});
+// // Point p6(rootCS, {0_m, 0_m, 6_m});
+// // Point p7(rootCS, {0_m, 0_m, 7_m});
+// // Point p8(rootCS, {0_m, 0_m, 8_m});
+// // Point p9(rootCS, {0_m, 0_m, 9_m});
+// Point p10(rootCS, {0_m, 0_m, 10_m});
+//
+// // get a unit vector
+// Vector<dimensionless_d> v1(rootCS, {0, 0, 1});
+//
+// // // get a geometrical path of points
+// // Path P1({p0,p1,p2,p3,p4,p5,p6,p7,p8,p9,p10});
+//
+// // construct a Straight Propagator given the uniform refractive index environment
+// StraightPropagator SP(env);
+//
+// // store the outcome of the Propagate method to paths_
+// auto const paths_ = SP.propagate(p0, p10, 1_m);
+//
+// // perform checks to paths_ components
+// 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_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
+//
+//// std::cout << "path.total_time_: " << path.total_time_ << std::endl;
+//// std::cout << "path.average_refractive_index_: " << path.average_refractive_index_ << std::endl;
+//// std::cout << "path.emit_: " << path.emit_.getComponents() << std::endl;
+//// std::cout << "path.R_distance_: " << path.R_distance_ << std::endl;
+//
//
// }
-
- // check reset method for antennas. Uncomment to see they are zero
-// ant1.reset();
-// ant2.reset();
//
-// std::cout << ant1.waveformE_ << std::endl;
-// std::cout << ant2.waveformE_ << std::endl;
+// CHECK(paths_.size() == 1);
+// }
+//
+// SECTION("Straight Propagator w/ Exponential Refractive Index") {
+//
+//// logging::set_level(logging::level::info);
+//// corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
+//
+// // create an environment with exponential refractive index (n_0 = 1 & lambda = 0)
+// using ExpoRIndex = ExponentialRefractiveIndex<HomogeneousMedium
+// <IRefractiveIndexModel<IMediumModel>>>;
+//
+// using EnvType = Environment<IRefractiveIndexModel<IMediumModel>>;
+// EnvType env1;
+//
+// //get another coordinate system
+// const CoordinateSystemPtr rootCS1 = env1.getCoordinateSystem();
+//
+// auto Medium1 = EnvType::createNode<Sphere>(
+// Point{rootCS1, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
+//
+// auto const props1 =
+// Medium1
+// ->setModelProperties<ExpoRIndex>( 1, 0 / 1_m,
+// 1_kg / (1_m * 1_m * 1_m),
+// NuclearComposition(
+// std::vector<Code>{Code::Nitrogen},
+// std::vector<float>{1.f}));
+//
+// env1.getUniverse()->addChild(std::move(Medium1));
+//
+// // get some points
+// Point pp0(rootCS1, {0_m, 0_m, 0_m});
+//// Point pp1(rootCS1, {0_m, 0_m, 1_m});
+//// Point pp2(rootCS1, {0_m, 0_m, 2_m});
+//// Point pp3(rootCS1, {0_m, 0_m, 3_m});
+//// Point pp4(rootCS1, {0_m, 0_m, 4_m});
+//// Point pp5(rootCS1, {0_m, 0_m, 5_m});
+//// Point pp6(rootCS1, {0_m, 0_m, 6_m});
+//// Point pp7(rootCS1, {0_m, 0_m, 7_m});
+//// Point pp8(rootCS1, {0_m, 0_m, 8_m});
+//// Point pp9(rootCS1, {0_m, 0_m, 9_m});
+// Point pp10(rootCS1, {0_m, 0_m, 10_m});
+//
+// // get a unit vector
+// Vector<dimensionless_d> vv1(rootCS1, {0, 0, 1});
+//
+// // construct a Straight Propagator given the exponential refractive index environment
+// StraightPropagator SP1(env1);
+//
+// // store the outcome of Propagate method to paths1_
+// auto const paths1_ = SP1.propagate(pp0, pp10, 1_m);
+//
+// // perform checks to paths1_ components (this is just a sketch for now)
+// 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_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 );
+//
+// /*
+// * A second environment with another exponential refractive index
+// */
+//
+// // create an environment with exponential refractive index (n_0 = 2 & lambda = 2)
+// using ExpoRIndex = ExponentialRefractiveIndex<HomogeneousMedium
+// <IRefractiveIndexModel<IMediumModel>>>;
+//
+// using EnvType = Environment<IRefractiveIndexModel<IMediumModel>>;
+// EnvType env2;
+//
+// //get another coordinate system
+// const CoordinateSystemPtr rootCS2 = env2.getCoordinateSystem();
+//
+// auto Medium2 = EnvType::createNode<Sphere>(
+// Point{rootCS2, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
+//
+// auto const props2 =
+// Medium2
+// ->setModelProperties<ExpoRIndex>( 2, 2 / 1_m,
+// 1_kg / (1_m * 1_m * 1_m),
+// NuclearComposition(
+// std::vector<Code>{Code::Nitrogen},
+// std::vector<float>{1.f}));
+//
+// env2.getUniverse()->addChild(std::move(Medium2));
+//
+// // get some points
+// Point ppp0(rootCS2, {0_m, 0_m, 0_m});
+// Point ppp10(rootCS2, {0_m, 0_m, 10_m});
+//
+// // get a unit vector
+// Vector<dimensionless_d> vvv1(rootCS2, {0, 0, 1});
+//
+// // construct a Straight Propagator given the exponential refractive index environment
+// StraightPropagator SP2(env2);
+//
+// // store the outcome of Propagate method to paths1_
+// auto const paths2_ = SP2.propagate(ppp0, ppp10, 1_m);
+//
+// // perform checks to paths1_ components (this is just a sketch for now)
+// 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_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 );
//
-// std::cout << " ... "<< std::endl;
-// std::cout << " ... "<< std::endl;
-
- // check reset method for antenna collection. Uncomment to see they are zero
-// detector.reset();
-// for (auto& xx : detector.getAntennas()) {
-// std::cout << xx.waveformE_ << std::endl;
-// std::cout << " ... "<< std::endl;
// }
-
-
- }
-
- // check that I can create working Straight Propagators in different environments
- SECTION("Straight Propagator w/ Uniform Refractive Index") {
-
- // create an environment with uniform refractive index of 1
- using UniRIndex =
- UniformRefractiveIndex<HomogeneousMedium<IRefractiveIndexModel<IMediumModel>>>;
-
- using EnvType = Environment<IRefractiveIndexModel<IMediumModel>>;
- EnvType env;
-
- // get a coordinate system
- const CoordinateSystemPtr rootCS = env.getCoordinateSystem();
-
- auto Medium = EnvType::createNode<Sphere>(
- Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
-
- auto const props = Medium->setModelProperties<UniRIndex>(
- 1, 1_kg / (1_m * 1_m * 1_m),
- NuclearComposition(
- std::vector<Code>{Code::Nitrogen},
- std::vector<float>{1.f}));
-
- env.getUniverse()->addChild(std::move(Medium));
-
- // get some points
- Point p0(rootCS, {0_m, 0_m, 0_m});
- // Point p1(rootCS, {0_m, 0_m, 1_m});
- // Point p2(rootCS, {0_m, 0_m, 2_m});
- // Point p3(rootCS, {0_m, 0_m, 3_m});
- // Point p4(rootCS, {0_m, 0_m, 4_m});
- // Point p5(rootCS, {0_m, 0_m, 5_m});
- // Point p6(rootCS, {0_m, 0_m, 6_m});
- // Point p7(rootCS, {0_m, 0_m, 7_m});
- // Point p8(rootCS, {0_m, 0_m, 8_m});
- // Point p9(rootCS, {0_m, 0_m, 9_m});
- Point p10(rootCS, {0_m, 0_m, 10_m});
-
- // get a unit vector
- Vector<dimensionless_d> v1(rootCS, {0, 0, 1});
-
- // // get a geometrical path of points
- // Path P1({p0,p1,p2,p3,p4,p5,p6,p7,p8,p9,p10});
-
- // construct a Straight Propagator given the uniform refractive index environment
- StraightPropagator SP(env);
-
- // store the outcome of the Propagate method to paths_
- auto const paths_ = SP.propagate(p0, p10, 1_m);
-
- // perform checks to paths_ components
- 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_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
-
- std::cout << "path.total_time_: " << path.total_time_ << std::endl;
- std::cout << "path.average_refractive_index_: " << path.average_refractive_index_ << std::endl;
- std::cout << "path.emit_: " << path.emit_.getComponents() << std::endl;
- std::cout << "path.R_distance_: " << path.R_distance_ << std::endl;
-
-
- }
-
- CHECK(paths_.size() == 1);
- }
-
- SECTION("Straight Propagator w/ Exponential Refractive Index") {
-
-// logging::set_level(logging::level::info);
-// corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
-
- // create an environment with exponential refractive index (n_0 = 1 & lambda = 0)
- using ExpoRIndex = ExponentialRefractiveIndex<HomogeneousMedium
- <IRefractiveIndexModel<IMediumModel>>>;
-
- using EnvType = Environment<IRefractiveIndexModel<IMediumModel>>;
- EnvType env1;
-
- //get another coordinate system
- const CoordinateSystemPtr rootCS1 = env1.getCoordinateSystem();
-
- auto Medium1 = EnvType::createNode<Sphere>(
- Point{rootCS1, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
-
- auto const props1 =
- Medium1
- ->setModelProperties<ExpoRIndex>( 1, 0 / 1_m,
- 1_kg / (1_m * 1_m * 1_m),
- NuclearComposition(
- std::vector<Code>{Code::Nitrogen},
- std::vector<float>{1.f}));
-
- env1.getUniverse()->addChild(std::move(Medium1));
-
- // get some points
- Point pp0(rootCS1, {0_m, 0_m, 0_m});
-// Point pp1(rootCS1, {0_m, 0_m, 1_m});
-// Point pp2(rootCS1, {0_m, 0_m, 2_m});
-// Point pp3(rootCS1, {0_m, 0_m, 3_m});
-// Point pp4(rootCS1, {0_m, 0_m, 4_m});
-// Point pp5(rootCS1, {0_m, 0_m, 5_m});
-// Point pp6(rootCS1, {0_m, 0_m, 6_m});
-// Point pp7(rootCS1, {0_m, 0_m, 7_m});
-// Point pp8(rootCS1, {0_m, 0_m, 8_m});
-// Point pp9(rootCS1, {0_m, 0_m, 9_m});
- Point pp10(rootCS1, {0_m, 0_m, 10_m});
-
- // get a unit vector
- Vector<dimensionless_d> vv1(rootCS1, {0, 0, 1});
-
- // construct a Straight Propagator given the exponential refractive index environment
- StraightPropagator SP1(env1);
-
- // store the outcome of Propagate method to paths1_
- auto const paths1_ = SP1.propagate(pp0, pp10, 1_m);
-
- // perform checks to paths1_ components (this is just a sketch for now)
- 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_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 );
-
- /*
- * A second environment with another exponential refractive index
- */
-
- // create an environment with exponential refractive index (n_0 = 2 & lambda = 2)
- using ExpoRIndex = ExponentialRefractiveIndex<HomogeneousMedium
- <IRefractiveIndexModel<IMediumModel>>>;
-
- using EnvType = Environment<IRefractiveIndexModel<IMediumModel>>;
- EnvType env2;
-
- //get another coordinate system
- const CoordinateSystemPtr rootCS2 = env2.getCoordinateSystem();
-
- auto Medium2 = EnvType::createNode<Sphere>(
- Point{rootCS2, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
-
- auto const props2 =
- Medium2
- ->setModelProperties<ExpoRIndex>( 2, 2 / 1_m,
- 1_kg / (1_m * 1_m * 1_m),
- NuclearComposition(
- std::vector<Code>{Code::Nitrogen},
- std::vector<float>{1.f}));
-
- env2.getUniverse()->addChild(std::move(Medium2));
-
- // get some points
- Point ppp0(rootCS2, {0_m, 0_m, 0_m});
- Point ppp10(rootCS2, {0_m, 0_m, 10_m});
-
- // get a unit vector
- Vector<dimensionless_d> vvv1(rootCS2, {0, 0, 1});
-
- // construct a Straight Propagator given the exponential refractive index environment
- StraightPropagator SP2(env2);
-
- // store the outcome of Propagate method to paths1_
- auto const paths2_ = SP2.propagate(ppp0, ppp10, 1_m);
-
- // perform checks to paths1_ components (this is just a sketch for now)
- 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_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 );
-
- }
} // END: TEST_CASE("Radio", "[processes]")