added refractive index at source and fixed it in CoREAS - added TODOs that need to be reviewed

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

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_);

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) {}

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()
 

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();
 }

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