CoREAS & ZHS advancements

corsika/modules/radio/CoREAS.hpp

@@ -15,13 +15,28 @@
 namespace corsika {
 
   /**
-   * 
+   * A concrete implementation of the Enpoints formalism.
    */
-  template <typename TPropagator = StraightPropagator>
-  class CoREAS final : public RadioProcess<CoREAS> {
+  template <typename TRadioDetector, typename TPropagator>
+  class CoREAS final : public RadioProcess<TRadioDetector, CoREAS<TRadioDetector, TPropagator>, TPropagator> {
 
+    using Base = RadioProcess<TRadioDetector, CoREAS<TRadioDetector, TPropagator>, TPropagator>;
+    using Base::detector_;
 
   public:
+    using ElectricFieldVector =
+    QuantityVector<ElectricFieldType::dimension_type>;
+    /**
+     * Construct a new CoREAS instance.
+     *
+     * This forwards the detector and other arguments to
+     * the RadioProcess parent.
+     *
+     */
+    template <typename... TArgs>
+    CoREAS(TRadioDetector& detector, TArgs&&... args)
+        : RadioProcess<TRadioDetector, CoREAS, TPropagator>(detector, args...){}
+
     /**
      * Simulate the radio emission from a particle across a track.
      *
@@ -32,20 +47,70 @@ namespace corsika {
      *
      */
     template <typename Particle, typename Track>
-    EProcessReturn simulate(Particle&, Track const&) const {
+    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()};
+
+      // beta is defined as velocity / speed of light
+      auto startBeta_ {track.getVelocity(0) / constants::c};
+      auto endBeta_ {track.getVelocity(1) / constants::c};
 
-        // loop over every antenna
-        for (auto& antenna : getAntennas()) {
+      // calculate gamma factor using beta (the proper way would be with energy over mass but not yet :( )
+      auto startGamma_ {1. / sqrt(1. - (startBeta_ * startBeta_))};
+      auto endGamma_ {1. / sqrt(1. - (endBeta_ * endBeta_))};
 
-            // use TPropagator to calculate Path from track to antenna
-            auto paths = ;/*  a collection of paths */
+      // we loop over each antenna in the collection
+      for (auto& antenna : detector_.getAntennas()) {
 
-            // do endpoint + ZHS formalism
+        // 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
+        // This is a SignalPathCollection
+        auto paths1{this->propagator_.propagate(startPoint, antenna.getLocation())};
 
-            // give the electric field and the direction to the antenna
-            antenna.receive(...); // << something
-        }
+        // 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.
+          // CoREAS calculation -> get ElectricFieldVector1
+          antenna.receive(/* startPoint_, receive vector, ElectricFieldVector1 */);
+
+        } // END: loop over paths
+
+        // get the Path from the track to the antenna
+        // This is a SignalPathCollection
+        auto paths2{this->propagator_.propagate(endPoint, antenna.getLocation())};
+
+        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.
+          //CoREAS calculation -> get ElectricFieldVector2
+          antenna.receive(/* endPoint_, receive vector, ElectricFieldVector2 */);
+
+        } // END: loop over paths
+
+      } // END: loop over antennas
     }
 
     /**
@@ -60,7 +125,16 @@ namespace corsika {
      */
     template <typename Particle, typename Track>
     LengthType MaxStepLength(Particle const& particle,
-                                        Track const& track) const;
+                             Track const& track) const {
+
+      // TODO : This is where we control the maximum step size
+      // of a particle track in order to maintain the accuracy
+      // of the particular formalism.
+      //
+      // This is part of the ZHS / CoReas formalisms and can
+      // be related from the magnetic field / acceleration, charge,
+      // etc. of the particle.
+    }
 
   }; // END: class RadioProcess
 

corsika/modules/radio/RadioProcess.hpp

@@ -7,13 +7,16 @@
  */
 #pragma once
 
-#include <corsika/framework/process/ContinuousProcess.hpp>
 #include <istream>
 #include <fstream>
 #include <iostream>
 #include <xtensor/xcsv.hpp>
 #include <xtensor/xtensor.hpp>
 #include <string>
+#include <corsika/framework/process/ContinuousProcess.hpp>
+#include <corsika/setup/SetupStack.hpp>
+#include <corsika/setup/SetupTrajectory.hpp>
+
 
 namespace corsika {
 
@@ -28,6 +31,14 @@ namespace corsika {
   class RadioProcess : public ContinuousProcess<
                            RadioProcess<TRadioDetector, TRadioImpl, TPropagator>> {
 
+    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_;
+
     /**
      * Get a reference to the underlying radio implementation.
      */
@@ -64,14 +75,38 @@ namespace corsika {
      */
     template <typename Particle, typename Track>
     ProcessReturn doContinuous(Particle& particle, Track const& track) const {
+      //we want the following particles:
+      // Code::Electron & Code::Positron & Code::Gamma
+
       // we wrap Simulate() in doContinuous as the plan is to add particle level
       // 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)
-      return this->implementation().simulate(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);
+    }
+
+    /**
        * TODO: This is placeholder so we can use text output while
        * we wait for the true output formatting to be ready.
        **/

corsika/modules/radio/ZHS.hpp

@@ -9,6 +9,8 @@
 
 #include <corsika/modules/radio/RadioProcess.hpp>
 #include <corsika/modules/radio/propagators/StraightPropagator.hpp>
+#include <corsika/framework/geometry/QuantityVector.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
 
 namespace corsika {
 
@@ -22,6 +24,8 @@ namespace corsika {
       using Base::detector_;
       
   public:
+    using ElectricFieldVector =
+    QuantityVector<ElectricFieldType::dimension_type>;
     /**
      * Construct a new ZHS instance.
      *
@@ -45,26 +49,27 @@ namespace corsika {
     template <typename Particle, typename Track>
     ProcessReturn simulate(Particle& particle, Track const& track) const {
 
-      auto global_time = particle.getTime(); // this is very shady at the moment...
+      auto globalTime_ = particle.getTime(); // this is very shady at the moment...
+
       //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 startTime_ = track.getDuration(0); // time at start point of track.
+//      auto endTime_ = track.getDuration(1); // time at end point of track.
 
       // we loop over each antenna in the collection
       for (auto& antenna : detector_.getAntennas()) {
 
         // auto start = /* TODO: get Point from Track */;
-        auto start = track.getStart(); // just another shady idea...
+        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(start, antenna.getLocation())};
+        auto paths{this->propagator_.propagate(startPoint, antenna.getLocation())};
 
         // 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_ + global_time;
+//          path.total_time_ + globalTime_;
 //          path.average_refractivity_;
 //          path.emit_;
 //          path.receive_;
@@ -73,6 +78,11 @@ namespace corsika {
           // combination along this path.
           // global time + time delay
           // and pass it to the antenna.
+
+//      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 */);
 
         } // END: loop over paths
@@ -94,7 +104,7 @@ namespace corsika {
     LengthType MaxStepLength(Particle const& particle,
                                         Track const& track) const {
 
-      // TODO : This is where te control the maximum step size
+      // TODO : This is where we control the maximum step size
       // of a particle track in order to maintain the accuracy
       // of the particular formalism.
       //