put formulas to CoREAS & improve ZHS & .writeOutput()

corsika/modules/radio/CoREAS.hpp

@@ -50,63 +50,64 @@ namespace corsika {
     ProcessReturn simulate(Particle& particle, Track const& track) const {
 
       //get global simulation time for that track. (This is my best guess for now)
-      auto startTime_ {particle.getTime()}; // time at start point of track.
-      auto endTime_ {particle.getTime() - track.getDuration()}; // time at end point of track.
-
-      // an alternative is the following which shouldn't work I think
-//      auto startTime_ {particle.getTime()};
-//      auto endTime_ {startTime_ + track.getDuration()};
+      auto startTime_ {particle.getTime()
+                      - track.getDuration()}; // time at start point of track.
+      auto endTime_ {particle.getTime()}; // time at end point of track.
 
       // beta is defined as velocity / speed of light
       auto startBeta_ {track.getVelocity(0) / constants::c};
       auto endBeta_ {track.getVelocity(1) / constants::c};
 
-      // calculate gamma factor using beta (the proper way would be with energy over mass but not yet :( )
+      // 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_))};
 
+      // get start and end position of the track
+      auto startPoint_ {track.getPosition(0)};
+      auto endPoint_ {track.getPosition(1)};
+
+      // get particle charge
+      auto const charge_ {get_charge(particle.getPID())};
+
       // we loop over each antenna in the collection
       for (auto& antenna : detector_.getAntennas()) {
 
-        // auto startPoint = /* TODO: get Point from Track */;
-        auto startPoint = track.getPosition(0); // this MIGHT work and also get it out of the for loop?
-        auto endPoint = track.getPosition(1); // I think we should get these guys out of the for loop
-        // get the Path from the track to the antenna
+        // get the Path (path1) from the start "endpoint" to the antenna.
         // This is a SignalPathCollection
-        auto paths1{this->propagator_.propagate(startPoint, antenna.getLocation())};
+        auto paths1{this->propagator_.propagate(startPoint_, antenna.getLocation())};
+        auto R1_ {(startPoint_ - antenna.getLocation()).getNorm()};
 
         // 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 : paths1) {
-//          auto startPoint_ {path.total_time_ + startTime_};
-//          path.average_refractivity_;
-//          path.emit_;
-//          path.receive_;
-
-          // combination along this path.
-          // global time + time delay
-          // and pass it to the antenna.
+          auto startPointReceiveTime_ {path.total_time_ + startTime_}; // might do it on the fly
+
           // CoREAS calculation -> get ElectricFieldVector1
-          antenna.receive(/* startPoint_, receive vector, ElectricFieldVector1 */);
+          ElectricFieldVector EV1_ {(charge_ / constants::c) *
+                                   path.receive_.cross(path.receive_.cross(startBeta_)) /
+                                        (R1_ * (1 - path.average_refractivity_ *
+                                                       startBeta_ * path.receive_))};
+
+          // pass it to the antenna
+          antenna.receive(startPointReceiveTime_, path.receive_, EV1_);
 
         } // END: loop over paths
 
-        // get the Path from the track to the antenna
+        // 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())};
+        auto R2_ {(endPoint_ - antenna.getLocation()).getNorm()};
 
         for (auto const& path : paths2) {
-//          auto endPoint_ {path.total_time + endTime_};
-//          path.average_refractivity_;
-//          path.emit_;
-//          path.receive_;
-
-          // combination along this path.
-          // global time + time delay
-          // and pass it to the antenna.
+          auto endPointReceiveTime_ {path.total_time + endTime_}; // might do it on the fly
+
           //CoREAS calculation -> get ElectricFieldVector2
-          antenna.receive(/* endPoint_, receive vector, ElectricFieldVector2 */);
+          ElectricFieldVector EV2_ {(charge_ / constants::c) *
+                                    path.receive_.cross(path.receive_.cross(endBeta_)) /
+                                    (R2_ * (1 - path.average_refractivity_ *
+                                                   endBeta_ * path.receive_))};
+
+          // pass it to the antenna
+          antenna.receive(endPointReceiveTime_, path.receive_, EV2_);
 
         } // END: loop over paths
 
@@ -134,6 +135,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;
     }
 
   }; // END: class RadioProcess

corsika/modules/radio/RadioProcess.hpp

@@ -31,13 +31,13 @@ namespace corsika {
   class RadioProcess : public ContinuousProcess<
                            RadioProcess<TRadioDetector, TRadioImpl, TPropagator>> {
 
-    using ParticleType = corsika::setup::Stack::particle_type;
-    using TrackType = corsika::LeapFrogTrajectory;
+//    using ParticleType = corsika::setup::Stack::particle_type;
+//    using TrackType = corsika::LeapFrogTrajectory;
 
     /**
      * A collection of filter objects for deciding on valid particles and tracks.
      */
-    std::vector<std::function<bool(ParticleType&, TrackType const&)>> filters_;
+    //std::vector<std::function<bool(ParticleType&, TrackType const&)>> filters_;
 
     /**
      * Get a reference to the underlying radio implementation.
@@ -82,29 +82,29 @@ namespace corsika {
       // filtering or thinning for calculation of the radio emission. This is
       // important for controlling the runtime of radio (by ignoring particles
       // that aren't going to contribute i.e. heavy hadrons)
-      if (valid(particle, track)) {
+      //if (valid(particle, track)) {
         return this->implementation().simulate(particle, track);
-      }
+      //}
     }
 
     /**
      * Decide whether this particle and track is valid for radio emission.
      */
-    template <typename Particle, typename Track>
-    auto valid(Particle& particle, Track const& track) const {
-
-      // loop over the filters in the our collection
-      for (auto& filter : filters_) {
-        // evaluate the filter. If the filter returns false,
-        // then this track is not valid for radio emission.
-        if (!filter(particle, track)) return false;
-      }
-    }
-
-    template <typename Particle, typename Track>
-    void addFilter(const std::function<bool(Particle&, Track const&)> filter) {
-      filters_.push_back(filter);
-    }
+//    template <typename Particle, typename Track>
+//    auto valid(Particle& particle, Track const& track) const {
+//
+//      // loop over the filters in the our collection
+//      for (auto& filter : filters_) {
+//        // evaluate the filter. If the filter returns false,
+//        // then this track is not valid for radio emission.
+//        if (!filter(particle, track)) return false;
+//      }
+//    }
+
+//    template <typename Particle, typename Track>
+//    void addFilter(const std::function<bool(Particle&, Track const&)> filter) {
+//      filters_.push_back(filter);
+//    }
 
     /**
        * TODO: This is placeholder so we can use text output while
@@ -114,10 +114,13 @@ namespace corsika {
 
         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, antenna.getWaveform());
+          xt::dump_csv(out_file, t, E);
           out_file.close();
           ++i;
+
         }
         // how this method should work:
         // 1. Loop over the antennas in the collection

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 <bits/stdc++.h>
 
 namespace corsika {
 
@@ -49,41 +50,46 @@ namespace corsika {
     template <typename Particle, typename Track>
     ProcessReturn simulate(Particle& particle, Track const& track) const {
 
-      auto globalTime_ = particle.getTime(); // this is very shady at the moment...
+      //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.
+      auto midTime_ {(startTime_ + endTime_) / 2};
 
-      //get global time for that track
-//      auto startTime_ = track.getDuration(0); // time at start point of track.
-//      auto endTime_ = track.getDuration(1); // time at end point of track.
+      auto startPoint_ = track.getPosition(0);
+
+      // track velocity
+      auto trackVelocity_ {(track.getVelocity(0) + track.getVelocity(1)) / 2};
+
+      // beta is defined as velocity / speed of light
+      auto beta_ { trackVelocity_ / constants::c};
+
+      // get particle charge
+      auto const charge_ {get_charge(particle.getPID())};
 
       // we loop over each antenna in the collection
       for (auto& antenna : detector_.getAntennas()) {
 
-        // auto start = /* TODO: get Point from Track */;
-        auto startPoint = track.getPosition(0); // this MIGHT work and also get it out of the for loop?
-
         // get the Path from the track to the antenna
         // This is a SignalPathCollection
-        auto paths{this->propagator_.propagate(startPoint, antenna.getLocation())};
+        auto paths{this->propagator_.propagate(startPoint_, antenna.getLocation())};
+        auto R_ {(startPoint_ - antenna.getLocation()).getNorm()};
 
         // 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) {
-//          path.total_time_ + globalTime_;
-//          path.average_refractivity_;
-//          path.emit_;
-//          path.receive_;
 
           // calculate the ZHS formalism for this particle-antenna
-          // combination along this path.
-          // global time + time delay
-          // and pass it to the 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_)));
 
-//      ElectricFieldVector EV_ = -(constants)*(V_perpendicular)*
-//      (delta(global time + time delay - (1 - average_refractivity*beta*costheta)t1)
-//       - delta(global time + time delay- (1 - average_refractivity*beta*costheta)t1))
-//      /(1 - average_refractivity*beta*costheta);
-          antenna.receive(/* global time + time delay, receive vector, ElectricFieldVector */);
+          // pass it to the antenna
+          antenna.receive(midTime_ + path.total_time_, path.receive_, EV_);
 
         } // END: loop over paths
 
@@ -111,6 +117,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 1000000000_m;
     }
 
   }; // END: class RadioProcess