added ZHS-like approximation to CoREAS + general updates

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 <cmath>
 
 namespace corsika {
 
@@ -49,6 +50,9 @@ 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.
@@ -57,6 +61,10 @@ namespace corsika {
       // 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_))};
@@ -72,44 +80,119 @@ namespace corsika {
       // we loop over each antenna in the collection
       for (auto& antenna : detector_.getAntennas()) {
 
+        ElectricFieldVector EV1_ {0_V / 0_m};
+        ElectricFieldVector EV2_ {0_V / 0_m};
+        ElectricFieldVector EV3_ {0_V / 0_m};
+
+        auto startPointReceiveTime_ {0_ns};
+        auto endPointReceiveTime_ {0_ns};
+
         // get the Path (path1) from the start "endpoint" to the antenna.
         // This is a SignalPathCollection
         auto paths1{this->propagator_.propagate(startPoint_, antenna.getLocation())};
-        auto R1_ {(startPoint_ - antenna.getLocation()).getNorm()};
 
-        // now loop over the paths that we got above
+        // 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) {
-          auto startPointReceiveTime_ {path.total_time_ + startTime_}; // might do it on the fly
 
-          // CoREAS calculation -> get ElectricFieldVector1
-          ElectricFieldVector EV1_ {(charge_ / constants::c) *
-                                   path.receive_.cross(path.receive_.cross(startBeta_)) /
-                                        (R1_ * (1 - path.average_refractivity_ *
-                                                       startBeta_ * path.receive_))};
+          preDoppler_ = 1. - path.average_refractive_index_ *
+                               startBeta_ * path.emit_;
 
-          // pass it to the antenna
-          antenna.receive(startPointReceiveTime_, path.receive_, EV1_);
+          if (preDoppler_ > approxThreshold_) {
+            startPointReceiveTime_ = path.total_time_ +
+                                              startTime_ ; // might do it on the fly
 
-        } // END: loop over paths
+            // CoREAS calculation -> get ElectricFieldVector1
+            EV1_= (charge_ / constants::c) *
+                path.receive_.cross(path.receive_.cross(startBeta_)) /
+                (path.R_distance_ * preDoppler_);
+
+            // pass it to the antenna
+            antenna.receive(startPointReceiveTime_, path.receive_, EV1_);
+          } else {
+            continue;
+          } // END preDoppler check
+
+        } // END: loop over paths for startpoint
 
         // get the Path (path2) from the end "endpoint" to the antenna.
         // This is a SignalPathCollection
         auto paths2{this->propagator_.propagate(endPoint_, antenna.getLocation())};
-        auto R2_ {(endPoint_ - antenna.getLocation()).getNorm()};
 
+        // 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) {
-          auto endPointReceiveTime_ {path.total_time + endTime_}; // might do it on the fly
+          postDoppler_ = 1. - path.average_refractive_index_ *
+                                endBeta_ * path.emit_;
 
-          //CoREAS calculation -> get ElectricFieldVector2
-          ElectricFieldVector EV2_ {(charge_ / constants::c) *
-                                    path.receive_.cross(path.receive_.cross(endBeta_)) /
-                                    (R2_ * (1 - path.average_refractivity_ *
-                                                   endBeta_ * path.receive_))};
+          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_);
-
-        } // END: loop over paths
+          } 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
     }