From fabf19c2603d09aac5e5440e343141499645accb Mon Sep 17 00:00:00 2001
From: Nikos <nckcoop@gmail.com>
Date: Wed, 31 Mar 2021 13:52:49 +0200
Subject: [PATCH] added refractive index at source and fixed it in CoREAS -
added TODOs that need to be reviewed
---
corsika/modules/radio/CoREAS.hpp | 38 ++++++++++---------
.../radio/antennas/TimeDomainAntenna.hpp | 4 +-
.../modules/radio/propagators/SignalPath.hpp | 8 ++--
.../radio/propagators/StraightPropagator.hpp | 3 +-
examples/radio_shower2.cpp | 1 +
tests/modules/testRadio.cpp | 8 ++--
6 files changed, 34 insertions(+), 28 deletions(-)
diff --git a/corsika/modules/radio/CoREAS.hpp b/corsika/modules/radio/CoREAS.hpp
index 1fc15e80c..07a4c0eef 100755
--- a/corsika/modules/radio/CoREAS.hpp
+++ b/corsika/modules/radio/CoREAS.hpp
@@ -66,6 +66,7 @@ namespace corsika {
auto startPoint_{track.getPosition(0)};
std::cout << "STARTPOINT : " << startPoint_ << std::endl;
auto endPoint_{track.getPosition(1)};
+// track.getVelocity(1); TODO: use this for velocity weight factors
std::cout << "ENDPOINT : " << endPoint_ << std::endl;
// beta is velocity / speed of light. Start & end should be the same!
@@ -85,27 +86,25 @@ 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_m)}; // 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()!!!!
// get the Path (path2) from the end "endpoint" to the antenna.
// This is a SignalPathCollection
auto paths2{this->propagator_.propagate(endPoint_, antenna.getLocation(), 1_m)};
// loop over both paths at once and directly compare 'start' and 'end' attributes
- for (size_t i = (paths1.size() == paths2.size()) ? 0 : paths1.size();
+ for (size_t i = (paths1.size() == paths2.size()) ? 0 : paths1.size(); // TODO: throw an exception if the sizes don't match
(i < paths1.size() && i < paths2.size()); i++) {
// First start with the 'start' point
// calculate preDoppler factor
- double preDoppler_{1. - paths1[i].average_refractive_index_ *
+ double preDoppler_{1. - paths1[i].refractive_index_source_ * // TODO: use the refractive index at source, not average!
beta_.dot(paths1[i].emit_)};
std::cout << "preDoppler: " << preDoppler_<< std::endl;
// calculate receive times for startpoint
- auto startPointReceiveTime_{paths1[i].total_time_ + startTime_};
- std::cout << "START RECEIVE TIME: " << startPointReceiveTime_ << std::endl;
+ auto startPointReceiveTime_{paths1[i].propagation_time_ + startTime_}; // TODO: total time -> propagation time
+ std::cout << "START RECEIVE TIME: " << startPointReceiveTime_ << std::endl; // TODO: time 0 is when the imaginary primary hits the ground
// get receive unit vector at 'start'
auto ReceiveVectorStart_ {paths1[i].receive_};
@@ -113,22 +112,21 @@ namespace corsika {
// calculate electric field vector for startpoint
ElectricFieldVector EV1_ =
- paths1[i]
- .receive_.cross(paths1[i].receive_.cross(beta_))
+ paths1[i].receive_.cross(paths1[i].receive_.cross(beta_))
.getComponents() /
(paths1[i].R_distance_ * preDoppler_) *
- (1 / (4 * M_PI * track.getDuration())) *
+ (1 / (4 * M_PI * track.getDuration())) * // TODO: divide by sample width not track.getDuration!
((1 / constants::epsilonZero) * (1 / constants::c)) * charge_;
// Now continue with the 'end' point
// calculate postDoppler factor
double postDoppler_{
- 1. - paths2[i].average_refractive_index_ *
+ 1. - paths2[i].refractive_index_source_ *
beta_.dot(paths2[i].emit_)}; // maybe this is path.receive_ (?)
std::cout << "postDoppler: " << postDoppler_<< std::endl;
// calculate receive times for endpoint
- auto endPointReceiveTime_{paths2[i].total_time_ + endTime_};
+ auto endPointReceiveTime_{paths2[i].propagation_time_ + endTime_};
std::cout << "END RECEIVE TIME: " << endPointReceiveTime_ << std::endl;
// get receive unit vector at 'end'
@@ -137,17 +135,18 @@ namespace corsika {
// calculate electric field vector for endpoint
ElectricFieldVector EV2_ =
- paths2[i]
- .receive_.cross(paths2[i].receive_.cross(beta_))
+ paths2[i].receive_.cross(paths2[i].receive_.cross(beta_))
.getComponents() /
(paths2[i].R_distance_ * postDoppler_) *
- ((-1) / (4 * M_PI * track.getDuration())) *
+ ((-1) / (4 * M_PI * track.getDuration())) * // TODO: divide by sample width not track.getDuration!
((1 / constants::epsilonZero) * (1 / constants::c)) * charge_;
//////////////////////////////////////////////////////////////////////////////
// start comparing stuff
if ((preDoppler_ < 1.e-9) || (postDoppler_ < 1.e-9)) {
+ std::cout << "Gets into if less than 1.e-9" << std::endl;
+
auto gridResolution_ {antenna.duration_};
auto deltaT_ { endPointReceiveTime_ - startPointReceiveTime_ };
@@ -218,6 +217,7 @@ namespace corsika {
} // End of if deltaT < gridresolution
} // End of if that checks small doppler factors
+ // TODO; fix this if with the one above, they should work together
// perform ZHS-like calculation close to Cherenkov angle
if (std::fabs(preDoppler_) <= approxThreshold_ || std::fabs(postDoppler_) <= approxThreshold_) {
@@ -229,7 +229,7 @@ namespace corsika {
auto midTime_{particle.getTime() - (track.getDuration() / 2)};
// get "mid" position of the track (that may not work properly)
- auto midPoint_{track.getPosition(0.5)};
+ auto midPoint_{track.getPosition(0.5)}; // TODO: get mid position geometrically
// get the Path (path3) from the middle "endpoint" to the antenna.
// This is a SignalPathCollection
@@ -238,7 +238,7 @@ namespace corsika {
// now loop over the paths for endpoint that we got above
for (auto const& path : paths3) {
- auto const midPointReceiveTime_{path.total_time_ + midTime_};
+ auto const midPointReceiveTime_{path.propagation_time_ + midTime_};
auto midDoppler_{1. - path.average_refractive_index_ * beta_.dot(path.emit_)};
// change the values of the receive unit vectors of start and end
@@ -346,8 +346,10 @@ namespace corsika {
std::cout << "RIGHT BEFORE RECEIVE INCIDENT :" << std::endl;
std::cout << "startTTimes_.at(index): " << startPointReceiveTime_ << std::endl;
- std::cout << "endTTimes_.at(index): " << endPointReceiveTime_ << std::endl;
+
antenna.receive(startPointReceiveTime_, ReceiveVectorStart_, EV1_);
+ std::cout << "SECOND RECEIVE ! ! !" << std::endl;
+ std::cout << "endTTimes_.at(index): " << endPointReceiveTime_ << std::endl;
antenna.receive(endPointReceiveTime_, ReceiveVectorEnd_, EV2_);
} // End of loop over both paths to get signal info
diff --git a/corsika/modules/radio/antennas/TimeDomainAntenna.hpp b/corsika/modules/radio/antennas/TimeDomainAntenna.hpp
index 9e5b539d0..a0c38711d 100644
--- a/corsika/modules/radio/antennas/TimeDomainAntenna.hpp
+++ b/corsika/modules/radio/antennas/TimeDomainAntenna.hpp
@@ -90,7 +90,7 @@ namespace corsika {
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.
@@ -119,7 +119,7 @@ namespace corsika {
xt::xtensor<double, 2> times_ (xt::zeros<double>({num_bins_, 1}));
for (int i = 0; i < num_bins_; i++) {
- times_.at(i,0) = static_cast<double>(start_time_ / 1_s + i * sample_rate_ * 1_s);
+ times_.at(i,0) = static_cast<double>(start_time_ / 1_s + i / (sample_rate_ * 1_s));
}
return std::make_pair(times_, waveformE_);
diff --git a/corsika/modules/radio/propagators/SignalPath.hpp b/corsika/modules/radio/propagators/SignalPath.hpp
index 741a90af4..f178ff627 100644
--- a/corsika/modules/radio/propagators/SignalPath.hpp
+++ b/corsika/modules/radio/propagators/SignalPath.hpp
@@ -23,8 +23,9 @@ namespace corsika {
using path = std::deque<Point>;
//TODO: discuss if we need average refractivity or average refractive index
- TimeType const total_time_; ///< The total propagation time.
+ TimeType const propagation_time_; ///< The total propagation time.
double const average_refractive_index_; ///< The average refractive index.
+ double const refractive_index_source_; ///< The refractive index at the source.
Vector<dimensionless_d> const emit_; ///< The (unit-length) emission vector.
Vector<dimensionless_d> const receive_; ///< The (unit-length) receive vector.
path const points_; ///< A collection of points that make up the geometrical path.
@@ -33,12 +34,13 @@ namespace corsika {
/**
* Create a new SignalPath instance.
*/
- SignalPath(TimeType const total_time, double const average_refractive_index,
+ SignalPath(TimeType const propagation_time, double const average_refractive_index, double const refractive_index_source,
Vector<dimensionless_d> const emit, Vector<dimensionless_d> const receive,
LengthType const R_distance, path const& points)
: Path(points)
- , total_time_(total_time)
+ , propagation_time_(propagation_time)
, average_refractive_index_(average_refractive_index)
+ , refractive_index_source_(refractive_index_source)
, emit_(emit)
, receive_(receive)
, R_distance_(R_distance) {}
diff --git a/corsika/modules/radio/propagators/StraightPropagator.hpp b/corsika/modules/radio/propagators/StraightPropagator.hpp
index 7bd009ca7..94df29496 100644
--- a/corsika/modules/radio/propagators/StraightPropagator.hpp
+++ b/corsika/modules/radio/propagators/StraightPropagator.hpp
@@ -98,6 +98,7 @@ namespace corsika {
auto const* node{universe->getContainingNode(destination)};
auto const refractive_index{node->getModelProperties().getRefractiveIndex(destination)};
// auto const refractive_index{1.000327};
+ auto const ri_source{refractive_index};
rindex.push_back(refractive_index);
points.push_back(destination);
@@ -129,7 +130,7 @@ namespace corsika {
// realize that emission and receive vector are 'direction' in this case.
//TODO: receive and emission vector should have opposite signs!
- return { SignalPath(time, averageRefractiveIndex_, direction , direction, distance_,points) };
+ return { SignalPath(time, averageRefractiveIndex_, ri_source, direction , direction, distance_,points) };
} // END: propagate()
diff --git a/examples/radio_shower2.cpp b/examples/radio_shower2.cpp
index 2135097ce..3f8d5353b 100644
--- a/examples/radio_shower2.cpp
+++ b/examples/radio_shower2.cpp
@@ -156,6 +156,7 @@ int main() {
Cascade EAS(env, tracking, sequence, stack);
EAS.run();
+ //TODO: this will run indefinetly due to no energy losses
// get radio output
coreas.writeOutput();
}
diff --git a/tests/modules/testRadio.cpp b/tests/modules/testRadio.cpp
index 84c6bf6cf..00c579b4c 100644
--- a/tests/modules/testRadio.cpp
+++ b/tests/modules/testRadio.cpp
@@ -163,8 +163,8 @@ TEST_CASE("Radio", "[processes]") {
// create times for the antenna
const TimeType t1{0_s}; // TODO: initialization of times to antennas! particle hits the observation level should be zero
- const TimeType t2{100_s};
- const InverseTimeType t3{1/1e+2_s};
+ const TimeType t2{10_s};
+ const InverseTimeType t3{1e+3_Hz};
const TimeType t4{11_s};
// check that I can create an antenna at (1, 2, 3)
@@ -198,7 +198,7 @@ TEST_CASE("Radio", "[processes]") {
auto const k{1_m * ((1_m) / ((1_s * 1_s) * 1_V))};
- auto const t = 10_s;
+ auto const t = 1_s;
LeapFrogTrajectory base(point4, v0, B0, k, t);
// create a new stack for each trial
@@ -263,7 +263,7 @@ TEST_CASE("Radio", "[processes]") {
coreas(detector, envCoREAS);
// check doContinuous and simulate methods
- coreas.doContinuous(particle1, base);
+ coreas.doContinuous(particle1, base, true);
// coreas1.simulate(particle1, base);
// check writeOutput method -> should produce 2 csv files for each antenna
--
GitLab