diff --git a/corsika/modules/radio/CoREAS.hpp b/corsika/modules/radio/CoREAS.hpp
index acf22484b58e6eef01a784914e2b4ad5f3429a01..ad15692f63b7168b21c6b5dd860fbddaaa0301d7 100755
--- a/corsika/modules/radio/CoREAS.hpp
+++ b/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
diff --git a/corsika/modules/radio/RadioProcess.hpp b/corsika/modules/radio/RadioProcess.hpp
index 96f7043332081edb8aad006ce138004726f0226b..93cfb69110121599f40e7bcb5fea69071560bba1 100644
--- a/corsika/modules/radio/RadioProcess.hpp
+++ b/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
diff --git a/corsika/modules/radio/ZHS.hpp b/corsika/modules/radio/ZHS.hpp
index eb5437191d7913c9a160fa6dba2f4c50addee291..0aa1e0617fe39381ce8e75fa512ec31924e11570 100644
--- a/corsika/modules/radio/ZHS.hpp
+++ b/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