From f945e4ea55171ce97ea01a45ac3fef54abbcd943 Mon Sep 17 00:00:00 2001
From: Remy Prechelt <prechelt@hawaii.edu>
Date: Mon, 10 May 2021 16:07:36 -1000
Subject: [PATCH] Ported radio branch to new output architecture.
---
corsika/modules/radio/CoREAS.hpp | 11 +-
corsika/modules/radio/RadioProcess.hpp | 190 +++++++++-----
corsika/modules/radio/ZHS.hpp | 242 +++++++++---------
corsika/modules/radio/antennas/Antenna.hpp | 87 +++++--
.../radio/antennas/TimeDomainAntenna.hpp | 240 ++++++++---------
tests/modules/testRadio.cpp | 12 +-
6 files changed, 465 insertions(+), 317 deletions(-)
diff --git a/corsika/modules/radio/CoREAS.hpp b/corsika/modules/radio/CoREAS.hpp
index 67af0a70b..a7933810e 100755
--- a/corsika/modules/radio/CoREAS.hpp
+++ b/corsika/modules/radio/CoREAS.hpp
@@ -23,7 +23,7 @@ namespace corsika {
class CoREAS final : public RadioProcess<TRadioDetector, CoREAS<TRadioDetector, TPropagator>, TPropagator> {
using Base = RadioProcess<TRadioDetector, CoREAS<TRadioDetector, TPropagator>, TPropagator>;
- using Base::detector_;
+ using Base::antennas_;
public:
using ElectricFieldVector =
@@ -36,8 +36,9 @@ namespace corsika {
*
*/
template <typename... TArgs>
- CoREAS(TRadioDetector& detector, TArgs&&... args)
- : RadioProcess<TRadioDetector, CoREAS, TPropagator>(detector, args...) {}
+ CoREAS(std::string const& name,
+ TRadioDetector& detector, TArgs&&... args)
+ : RadioProcess<TRadioDetector, CoREAS, TPropagator>(name, detector, args...) {}
/**
* Simulate the radio emission from a particle across a track.
@@ -77,7 +78,7 @@ namespace corsika {
const double approxThreshold_{1.0e-3};
// loop over each antenna in the antenna collection (detector)
- for (auto& antenna : detector_.getAntennas()) {
+ for (auto& antenna : antennas_.getAntennas()) {
// check with which antenna we work in this loop
std::cout << "ANTENNA: " << antenna.getName() << std::endl;
@@ -382,4 +383,4 @@ namespace corsika {
}; // END: class CoREAS
-} // namespace corsika
\ No newline at end of file
+} // namespace corsika
diff --git a/corsika/modules/radio/RadioProcess.hpp b/corsika/modules/radio/RadioProcess.hpp
index ecc0943e2..527837657 100644
--- a/corsika/modules/radio/RadioProcess.hpp
+++ b/corsika/modules/radio/RadioProcess.hpp
@@ -10,34 +10,34 @@
#include <istream>
#include <fstream>
#include <iostream>
-#include <xtensor/xcsv.hpp>
#include <xtensor/xtensor.hpp>
#include <string>
+#include <corsika/output/BaseOutput.hpp>
#include <corsika/framework/process/ContinuousProcess.hpp>
#include <corsika/setup/SetupStack.hpp>
#include <corsika/setup/SetupTrajectory.hpp>
-
namespace corsika {
/**
* The base interface for radio emission processes.
*
* TRadioImpl is the concrete implementation of the radio algorithm.
- * TRadioDetector is the detector instance that stores antennas
+ * TAntennaCollection is the detector instance that stores antennas
* and is responsible for managing the output writing.
*/
- template <typename TRadioDetector, typename TRadioImpl, typename TPropagator>
+ template <typename TAntennaCollection, typename TRadioImpl, typename TPropagator>
class RadioProcess : public ContinuousProcess<
- RadioProcess<TRadioDetector, TRadioImpl, TPropagator>> {
+ RadioProcess<TAntennaCollection, TRadioImpl, TPropagator>>,
+ public BaseOutput {
-// 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.
@@ -52,16 +52,19 @@ namespace corsika {
}
protected:
- TRadioDetector& detector_; ///< The radio detector we store into.
- TPropagator propagator_; ///< The propagator implementation.
+ std::string const name_; ///< The name of this radio process.
+ TAntennaCollection& antennas_; ///< The radio antennas we store into.
+ TPropagator propagator_; ///< The propagator implementation.
+ int event_{0}; ///< The current event ID.
public:
/**
* Construct a new RadioProcess.
*/
template <typename... TArgs>
- RadioProcess(TRadioDetector& detector, TArgs&&... args)
- : detector_(detector)
+ RadioProcess(std::string const& name, TAntennaCollection& antennas, TArgs&&... args)
+ : name_(name)
+ , antennas_(antennas)
, propagator_(args...) {}
/**
@@ -75,14 +78,14 @@ namespace corsika {
*/
template <typename Particle, typename Track>
ProcessReturn doContinuous(Particle& particle, Track const& track, bool const) {
- //we want the following particles:
+ // 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)
- //if (valid(particle, track)) {
+ // if (valid(particle, track)) {
if (particle.getPID() == Code::Electron || particle.getPID() == Code::Positron) {
return this->implementation().simulate(particle, track);
}
@@ -92,62 +95,131 @@ namespace corsika {
/**
* 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.
- **/
- bool writeOutput() const {
- // this for loop still has some issues
- int i = 1;
- for (auto& antenna : detector_.getAntennas()) {
-
- auto [t,E] = antenna.getWaveform();
- auto c = xt::hstack(xt::xtuple(t,E));
- std::ofstream out_file ("antenna" + to_string(i) + "_output.csv");
- xt::dump_csv(out_file, c);
- out_file.close();
- ++i;
-
- }
+ // 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.
+ // **/
+ // bool writeOutput() const {
+ // // this for loop still has some issues
+ // int i = 1;
+ // for (auto& antenna : antennas_.getAntennas()) {
+
+ // auto [t, E] = antenna.getWaveform();
+ // auto c = xt::hstack(xt::xtuple(t, E));
+ // std::ofstream out_file("antenna" + to_string(i) + "_output.csv");
+ // xt::dump_csv(out_file, c);
+ // out_file.close();
+ // ++i;
+ // }
// how this method should work:
// 1. Loop over the antennas in the collection
// 2. Get their waveforms
// 3. Create a text file for each antenna
// 4. and write out two columns, time and field.
-
- }
+ // }
/**
- * Return the maximum step length for this particle and track.
- *
- * This must be provided by the TRadioImpl.
- *
- * @param particle The current particle.
- * @param track The current track.
- *
- * @returns The maximum length of this track.
- */
+ * Return the maximum step length for this particle and track.
+ *
+ * This must be provided by the TRadioImpl.
+ *
+ * @param particle The current particle.
+ * @param track The current track.
+ *
+ * @returns The maximum length of this track.
+ */
LengthType getMaxStepLength(setup::Stack::particle_type const& vParticle,
setup::Trajectory const& vTrack) const {
return meter * std::numeric_limits<double>::infinity();
}
+ /**
+ * Called at the start of each library.
+ */
+ void startOfLibrary(boost::filesystem::path const& directory) final override {
+
+ // loop over every antenna and set the initial path
+ // this also writes the time-bins to disk.
+ for (auto& antenna : antennas_.getAntennas()) {
+ antenna.startOfLibrary(directory);
+ }
+ }
+
+ /**
+ * Called at the end of each shower.
+ */
+ void endOfShower() final override {
+
+ // loop over every antenna and instruct them to
+ // flush data to disk, and then reset the antenna
+ // before the next event
+ for (auto& antenna : antennas_.getAntennas()) {
+ antenna.endOfShower(event_);
+ antenna.reset();
+ }
+
+ // increment our event counter
+ event_++;
+
+ }
+
+ /**
+ * Called at the end of each library.
+ *
+ * This must also increment the run number since we override
+ * the default behaviour of BaseOutput.
+ */
+ void endOfLibrary() final override;
+
+ /**
+ * Get the configuration of this output.
+ */
+ YAML::Node getConfig() const final {
+
+ // top-level YAML node
+ YAML::Node config;
+
+ // fill in some basics
+ config["type"] = "RadioProcess";
+ config["units"]["time"] = "ns";
+ config["units"]["frequency"] = "GHz";
+ config["units"]["electric field"] = "V/m";
+ config["units"]["distance"] = "m";
+
+ for (auto& antenna : antennas_.getAntennas()) {
+ // get the name/location of this antenna
+ auto name = antenna.getName();
+ auto location = antenna.getLocation().getCoordinates();
+
+ // get the antennas config
+ config["antennas"][name] = antenna.getConfig();
+
+ // write the location of this antenna
+ config["antennas"][name]["location"].push_back(location.getX() / 1_m);
+ config["antennas"][name]["location"].push_back(location.getY() / 1_m);
+ config["antennas"][name]["location"].push_back(location.getZ() / 1_m);
+
+ }
+
+ return config;
+ }
+
}; // END: class RadioProcess
-} // namespace corsika
\ No newline at end of file
+} // namespace corsika
diff --git a/corsika/modules/radio/ZHS.hpp b/corsika/modules/radio/ZHS.hpp
index c4a7ab282..f5bbb0049 100644
--- a/corsika/modules/radio/ZHS.hpp
+++ b/corsika/modules/radio/ZHS.hpp
@@ -20,13 +20,15 @@ namespace corsika {
* A concrete implementation of the ZHS algorithm.
*/
template <typename TRadioDetector, typename TPropagator>
- class ZHS final : public RadioProcess<TRadioDetector, ZHS<TRadioDetector, TPropagator>, TPropagator> {
+ class ZHS final : public RadioProcess<TRadioDetector, ZHS<TRadioDetector, TPropagator>,
+ TPropagator> {
- using Base = RadioProcess<TRadioDetector, ZHS<TRadioDetector, TPropagator>, TPropagator>;
- using Base::detector_;
+ using Base =
+ RadioProcess<TRadioDetector, ZHS<TRadioDetector, TPropagator>, TPropagator>;
+ using Base::antennas_;
public:
-// using PotentialVector = QuantityVector<PotentialVectorType::dimension_type>;
+ // using PotentialVector = QuantityVector<PotentialVectorType::dimension_type>;
using ElectricFieldVector = QuantityVector<ElectricFieldType::dimension_type>;
/**
* Construct a new ZHS instance.
@@ -36,8 +38,8 @@ namespace corsika {
*
*/
template <typename... TArgs>
- ZHS(TRadioDetector& detector, TArgs&&... args)
- : RadioProcess<TRadioDetector, ZHS, TPropagator>(detector, args...){}
+ ZHS(std::string const& name, TRadioDetector& detector, TArgs&&... args)
+ : RadioProcess<TRadioDetector, ZHS, TPropagator>(name, detector, args...) {}
/**
* Simulate the radio emission from a particle across a track.
@@ -51,129 +53,141 @@ namespace corsika {
template <typename Particle, typename Track>
ProcessReturn simulate(Particle& particle, Track const& track) const {
- // TODO: think if we reuse these variables for the case of not being in the Fraunhoffer approx.
- auto const startTime_ {particle.getTime()
- - track.getDuration()}; // time at start point of track.
- auto const endTime_ {particle.getTime()}; // time at end point of track.
+ // TODO: think if we reuse these variables for the case of not being in the
+ // Fraunhoffer approx.
+ auto const startTime_{particle.getTime() -
+ track.getDuration()}; // time at start point of track.
+ auto const endTime_{particle.getTime()}; // time at end point of track.
auto const startPoint_{track.getPosition(0)};
auto const endPoint_{track.getPosition(1)};
- LengthType trackLength_ {(startPoint_ - endPoint_).getNorm()};
+ LengthType trackLength_{(startPoint_ - endPoint_).getNorm()};
// track velocity
- auto const trackVelocity_ {(track.getVelocity(0) + track.getVelocity(1)) / 2};
+ auto const trackVelocity_{(track.getVelocity(0) + track.getVelocity(1)) / 2};
// beta is defined as velocity / speed of light
- auto const beta_ { trackVelocity_ / constants::c};
+ auto const beta_{trackVelocity_ / constants::c};
// get particle charge
- auto const charge_ {get_charge(particle.getPID())};
+ auto const charge_{get_charge(particle.getPID())};
// get "mid" position of the track geometrically
- auto const midVector_ { (startPoint_ - endPoint_) / 2 };
- auto const midPoint_ {Point(midVector_.getCoordinateSystem(), midVector_.getComponents().getX(),
- midVector_.getComponents().getY(), midVector_.getComponents().getZ())};
- //changed if deltaT1_ > deltaT2_, so not const
- auto constants {charge_/ (4 * M_PI)/ (constants::epsilonZero) / constants::c};
+ auto const midVector_{(startPoint_ - endPoint_) / 2};
+ auto const midPoint_{
+ Point(midVector_.getCoordinateSystem(), midVector_.getComponents().getX(),
+ midVector_.getComponents().getY(), midVector_.getComponents().getZ())};
+ // changed if deltaT1_ > deltaT2_, so not const
+ auto constants{charge_ / (4 * M_PI) / (constants::epsilonZero) / constants::c};
// we loop over each antenna in the collection
- for (auto& antenna : detector_.getAntennas()) {
+ for (auto& antenna : antennas_.getAntennas()) {
// Probably will be reused in the case of not being in Fraunhoffer.
auto midPaths{this->propagator_.propagate(midPoint_, antenna.getLocation(), 1_m)};
- //Loop over midPaths to first check Fraunhoffer limit
- for (size_t i{0}; i < midPaths.size(); i++){
- //Maybe I can recalculate fraunhLimit without Midpaths to avoid calculating this third
- //path which is something really slow, and only use path1 and path2.
- double const betaTimesK {beta_.dot(midPaths[i].emit_)/beta_.getNorm()};
- double const slitSize_ {1. - betaTimesK * betaTimesK * trackLength_ * trackLength_ /1_m /1_m};
- //Parameter that determines the limit for the Fraunhoffer limit (probably related to antenna sampling rate
- double const lambda {0.1};
- double const fraunhLimit {slitSize_/midPaths[i].R_distance_/lambda * 1_m};
- if (fraunhLimit > 1.) // Checks if we are in fraunhoffer domain (maybe it should be less?)
- {
- ///code for dividing track and calculating field.
- std::cout << "This code hasnt been implemented!" << std::endl;
+ // Loop over midPaths to first check Fraunhoffer limit
+ for (size_t i{0}; i < midPaths.size(); i++) {
+ // Maybe I can recalculate fraunhLimit without Midpaths to avoid calculating
+ // this third path which is something really slow, and only use path1 and path2.
+ double const betaTimesK{beta_.dot(midPaths[i].emit_) / beta_.getNorm()};
+ double const slitSize_{1. - betaTimesK * betaTimesK * trackLength_ *
+ trackLength_ / 1_m / 1_m};
+ // Parameter that determines the limit for the Fraunhoffer limit (probably
+ // related to antenna sampling rate
+ double const lambda{0.1};
+ double const fraunhLimit{slitSize_ / midPaths[i].R_distance_ / lambda * 1_m};
+ if (fraunhLimit >
+ 1.) // Checks if we are in fraunhoffer domain (maybe it should be less?)
+ {
+ /// code for dividing track and calculating field.
+ std::cout << "This code hasnt been implemented!" << std::endl;
+ } else // Calculate vector potential of whole track
+ {
+ // paths from start of track to antenna
+ auto paths1{
+ this->propagator_.propagate(startPoint_, antenna.getLocation(), 1_m)};
+
+ // paths from end of track to antenna
+ auto paths2{
+ this->propagator_.propagate(endPoint_, antenna.getLocation(), 1_m)};
+ // First implementation, need to check if there are the same number of paths,
+ // but once its working.
+
+ // modified later if detectionTime1_ > detectionTime2_
+ TimeType detectionTime1_{startTime_ + paths1[i].propagation_time_};
+ TimeType detectionTime2_{endTime_ + paths2[i].propagation_time_};
+
+ // make detectionTime1_ be the smallest time => changes step function order so
+ // constants is changed to account for it
+ if (detectionTime1_ > detectionTime2_) {
+ detectionTime1_ = endTime_ + paths2[i].propagation_time_;
+ detectionTime2_ = startTime_ + paths1[i].propagation_time_;
+ constants = -constants;
}
- else // Calculate vector potential of whole track
- {
- // paths from start of track to antenna
- auto paths1{this->propagator_.propagate(startPoint_, antenna.getLocation(), 1_m)};
-
- // paths from end of track to antenna
- auto paths2{this->propagator_.propagate(endPoint_, antenna.getLocation(), 1_m)};
- // First implementation, need to check if there are the same number of paths, but once its working.
-
- //modified later if detectionTime1_ > detectionTime2_
- TimeType detectionTime1_ {startTime_ + paths1[i].propagation_time_};
- TimeType detectionTime2_ {endTime_ + paths2[i].propagation_time_};
-
- //make detectionTime1_ be the smallest time => changes step function order so
- // constants is changed to account for it
- if (detectionTime1_ > detectionTime2_)
- {
- detectionTime1_ = endTime_ + paths2[i].propagation_time_;
- detectionTime2_ = startTime_ + paths1[i].propagation_time_;
- constants = - constants;
- }
-
- double const startBin {std::floor(detectionTime1_ * antenna.sample_rate_)};
- double const endBin {std::floor(detectionTime2_ * antenna.sample_rate_)};
-
- auto const betaPerp_ {midPaths[i].emit_.cross(beta_.cross(midPaths[i].emit_))};
- double const denominator {1 - midPaths[i].refractive_index_source_ *
- beta_.dot(midPaths[i].emit_)};
-
-
- //IMPORTANT!!!!! Using ElectricFieldVector instead of PotentialVector until antenna is adapted
- if (startBin == endBin) {
- //track contained in bin
- //if not in Cerenkov angle then
- if (std::fabs(denominator)>1.e-15L) {
- double const f {std::fabs((detectionTime2_ * antenna.sample_rate_ - detectionTime1_ * antenna.sample_rate_) )};
- //should be PotentialVector const Vp_ = betaPerp_.getComponents()/denominator/
- // midPaths[i].R_distance_ * constants * f;
- // but to make it compile until antenna is adapted it stays like that to test it.
- ElectricFieldVector const Vp_ = betaPerp_.getComponents()/denominator/
- midPaths[i].R_distance_ * constants * f / 1_s;
- antenna.receive(detectionTime2_, betaPerp_, Vp_);
- }
- else {//If emission in Cerenkov angle => approximation
- double const f {(detectionTime2_ - detectionTime1_) * antenna.sample_rate_};
- ElectricFieldVector const Vp_ = betaPerp_.getComponents()/
- midPaths[i].R_distance_ * constants * f / 1_s;
- antenna.receive(detectionTime2_, betaPerp_, Vp_);
- }//end if Cerenkov angle approx
- } else {
- /*Track is contained in more than one bin*/
- int const numberOfBins {static_cast<int>(endBin - startBin)};
- //TODO: should we check for Cerenkov angle?
- //first contribution
- double f{std::fabs(startBin + 1. - detectionTime1_ * antenna.sample_rate_)};
- ElectricFieldVector Vp_ = betaPerp_.getComponents() * f *
- constants / denominator /
- midPaths[i].R_distance_ / 1_s;
- antenna.receive(detectionTime1_, betaPerp_, Vp_);
- //intermidiate contributions
- for (int it{1}; it < numberOfBins; ++it) {
- Vp_ = betaPerp_.getComponents() * constants / denominator /
- midPaths[i].R_distance_ / 1_s;
- antenna.receive(detectionTime1_ + static_cast<double>(it) / antenna.sample_rate_,
- betaPerp_, Vp_);
- }// end loop over bins in which potential vector is not zero
- //final contribution
- f = std::fabs(detectionTime2_ * antenna.sample_rate_ - endBin);
- Vp_ = betaPerp_.getComponents() * f * constants / denominator /
- midPaths[i].R_distance_ / 1_s;
- antenna.receive(detectionTime2_, betaPerp_, Vp_);
- }//end if statement for track in multiple bins
-
- }//finish if statement of track in fraunhoffer or not
-
- }//end loop over mid paths
-
- } // END: loop over antennas
- return ProcessReturn::Ok;
- } //end simulate
+
+ double const startBin{std::floor(detectionTime1_ * antenna.sample_rate_)};
+ double const endBin{std::floor(detectionTime2_ * antenna.sample_rate_)};
+
+ auto const betaPerp_{midPaths[i].emit_.cross(beta_.cross(midPaths[i].emit_))};
+ double const denominator{1 - midPaths[i].refractive_index_source_ *
+ beta_.dot(midPaths[i].emit_)};
+
+ // IMPORTANT!!!!! Using ElectricFieldVector instead of PotentialVector until
+ // antenna is adapted
+ if (startBin == endBin) {
+ // track contained in bin
+ // if not in Cerenkov angle then
+ if (std::fabs(denominator) > 1.e-15L) {
+ double const f{std::fabs((detectionTime2_ * antenna.sample_rate_ -
+ detectionTime1_ * antenna.sample_rate_))};
+ // should be PotentialVector const Vp_ =
+ // betaPerp_.getComponents()/denominator/
+ // midPaths[i].R_distance_ * constants
+ // * f;
+ // but to make it compile until antenna is adapted it stays like that to
+ // test it.
+ ElectricFieldVector const Vp_ = betaPerp_.getComponents() / denominator /
+ midPaths[i].R_distance_ * constants * f /
+ 1_s;
+ antenna.receive(detectionTime2_, betaPerp_, Vp_);
+ } else { // If emission in Cerenkov angle => approximation
+ double const f{(detectionTime2_ - detectionTime1_) *
+ antenna.sample_rate_};
+ ElectricFieldVector const Vp_ = betaPerp_.getComponents() /
+ midPaths[i].R_distance_ * constants * f /
+ 1_s;
+ antenna.receive(detectionTime2_, betaPerp_, Vp_);
+ } // end if Cerenkov angle approx
+ } else {
+ /*Track is contained in more than one bin*/
+ int const numberOfBins{static_cast<int>(endBin - startBin)};
+ // TODO: should we check for Cerenkov angle?
+ // first contribution
+ double f{std::fabs(startBin + 1. - detectionTime1_ * antenna.sample_rate_)};
+ ElectricFieldVector Vp_ = betaPerp_.getComponents() * f * constants /
+ denominator / midPaths[i].R_distance_ / 1_s;
+ antenna.receive(detectionTime1_, betaPerp_, Vp_);
+ // intermidiate contributions
+ for (int it{1}; it < numberOfBins; ++it) {
+ Vp_ = betaPerp_.getComponents() * constants / denominator /
+ midPaths[i].R_distance_ / 1_s;
+ antenna.receive(
+ detectionTime1_ + static_cast<double>(it) / antenna.sample_rate_,
+ betaPerp_, Vp_);
+ } // end loop over bins in which potential vector is not zero
+ // final contribution
+ f = std::fabs(detectionTime2_ * antenna.sample_rate_ - endBin);
+ Vp_ = betaPerp_.getComponents() * f * constants / denominator /
+ midPaths[i].R_distance_ / 1_s;
+ antenna.receive(detectionTime2_, betaPerp_, Vp_);
+ } // end if statement for track in multiple bins
+
+ } // finish if statement of track in fraunhoffer or not
+
+ } // end loop over mid paths
+
+ } // END: loop over antennas
+ return ProcessReturn::Ok;
+ } // end simulate
}; // END: class ZHS
-} // namespace corsika
\ No newline at end of file
+} // namespace corsika
diff --git a/corsika/modules/radio/antennas/Antenna.hpp b/corsika/modules/radio/antennas/Antenna.hpp
index 8fe910cc8..2acc44b94 100644
--- a/corsika/modules/radio/antennas/Antenna.hpp
+++ b/corsika/modules/radio/antennas/Antenna.hpp
@@ -9,6 +9,9 @@
*/
#pragma once
+#include <cnpy.hpp>
+#include <xtensor/xtensor.hpp>
+#include <boost/filesystem.hpp>
#include <corsika/framework/geometry/Point.hpp>
namespace corsika {
@@ -20,23 +23,22 @@ namespace corsika {
* type Antenna<T> where T is a concrete antenna implementation.
*
*/
- template <typename AntennaImpl>
+ template <typename TAntennaImpl>
class Antenna {
-
+ protected:
+ /**
+ * Get a reference to the underlying radio implementation.
+ */
+ TAntennaImpl& implementation() { return static_cast<TAntennaImpl&>(*this); }
public:
- std::string const name_; ///< The name/identifier of this antenna.
- Point const location_; ///< The location of this antenna.
+ std::string const name_; ///< The name/identifier of this antenna.
+ Point const location_; ///< The location of this antenna.
+ std::string filename_ = ""; ///< The filename for the output file for this antenna.
// this stores the polarization vector of an electric field
- using ElectricFieldVector =
- QuantityVector<ElectricFieldType::dimension_type>;
-// using MagneticFieldVector =
-// QuantityVector<MagneticFieldType::dimension_type>;
-
- // a dimensionless vector used for the incident direction
-// using Vector = QuantityVector<dimensionless_d>;
+ using ElectricFieldVector = QuantityVector<ElectricFieldType::dimension_type>;
/**
* \brief Construct a base antenna instance.
@@ -49,11 +51,6 @@ namespace corsika {
: name_(name)
, location_(location){};
- // copy constructor
- Antenna(const Antenna& Ant)
- : name_(Ant.name_)
- , location_(Ant.location_){};
-
/**
* Receive a signal at this antenna.
*
@@ -81,6 +78,64 @@ namespace corsika {
*/
void reset();
+ /**
+ * Return a reference to the x-axis labels (i.e. time or frequency).
+ *
+ * This should be an xtensor-convertible type with
+ * a ->data() method that converts to a raw pointer.
+ */
+ xt::xtensor<double, 1> getAxis() const;
+
+ /**
+ * Return a reference to the underlying data.
+ *
+ * This is used when writing the antenna information to disk
+ * and will be converted to a 32-bit float before writing.
+ */
+ xt::xtensor<double, 2>& getData() const;
+
+ /**
+ * Prepare for the start of the library.
+ */
+ void startOfLibrary(boost::filesystem::path const& directory) {
+
+ // calculate and save our filename
+ filename_ = (directory / this->getName()).string() + ".npz";
+
+ // get the axis labels for this antenna and write the first row.
+ xt::xtensor<float, 1> axis = xt::cast<float>(this->implementation().getAxis());
+
+ // check for the axis name
+ std::string label = "Unknown";
+ if constexpr (TAntennaImpl::is_time_domain) {
+ label = "Time";
+ } else if constexpr (TAntennaImpl::is_freq_domain) {
+ label = "Frequency";
+ }
+
+ // explicitly convert the arrays to the needed type for cnpy
+ float const* raw_data = axis.data();
+ std::vector<size_t> N = {axis.size()}; // cnpy needs a vector here
+
+ // write the labels to the first row of the NumPy file
+ cnpy::npz_save(filename_, label, raw_data, N, "w");
+ }
+
+ /**
+ * Flush the data from this shower to disk.
+ */
+ void endOfShower(int const event) {
+
+ // get the copy of the waveform data for this event
+ auto data{this->implementation().getData()};
+
+ // cnpy needs a vector for the shape
+ std::vector<size_t> shape = {data.shape()[0], data.shape()[1]};
+
+ // and write this event to the .npz archive
+ cnpy::npz_save(filename_, std::to_string(event), data.data(), shape, "a");
+ }
+
}; // END: class Antenna final
} // namespace corsika
diff --git a/corsika/modules/radio/antennas/TimeDomainAntenna.hpp b/corsika/modules/radio/antennas/TimeDomainAntenna.hpp
index e3cf1a745..879432f12 100644
--- a/corsika/modules/radio/antennas/TimeDomainAntenna.hpp
+++ b/corsika/modules/radio/antennas/TimeDomainAntenna.hpp
@@ -16,125 +16,131 @@
namespace corsika {
+ /**
+ * An implementation of a time-domain antenna that has a customized
+ * start time, sampling rate, and waveform duration.
+ *
+ */
+ class TimeDomainAntenna : public Antenna<TimeDomainAntenna> {
+
+ public:
+ // import the methods from the antenna
+
+ // label this as a time-domain antenna.
+ static constexpr bool is_time_domain{true};
+
+ TimeType const start_time_; ///< The start time of this waveform.
+ TimeType const duration_; ///< The duration of this waveform.
+ InverseTimeType const sample_rate_; ///< The sampling rate of this antenna.
+ int num_bins_; ///< The number of bins used.
+ xt::xtensor<double, 2> waveformE_; ///< The waveform stored by this antenna.
+
+ using Antenna<TimeDomainAntenna>::getName;
+ using Antenna<TimeDomainAntenna>::getLocation;
+
+
+ /**
+ * Construct a new TimeDomainAntenna.
+ *
+ * @param name The name of this antenna.
+ * @param location The location of this antenna.
+ * @param start_time The starting time of this waveform.
+ * @param duration The duration of this waveform.
+ * @param sample_rate The sample rate of this waveform.
+ * @param num_bins_ The number of timebins to store E-field.
+ * @param waveformE_ The xtensor initialized to zero for E-field.
+ *
+ */
+ TimeDomainAntenna(std::string const& name, Point const& location,
+ TimeType const& start_time, TimeType const& duration,
+ InverseTimeType const& sample_rate)
+ : Antenna(name, location)
+ , start_time_(start_time)
+ , duration_(duration)
+ , sample_rate_(sample_rate)
+ , num_bins_(static_cast<int>(duration * sample_rate))
+ , waveformE_(xt::zeros<double>({num_bins_, 3})){};
+
/**
- * An implementation of a time-domain antenna that has a customized
- * start time, sampling rate, and waveform duration.
+ * Receive an electric field at this antenna.
+ *
+ * This assumes that the antenna will receive
+ * an *instantaneous* electric field modeled as a delta function (or timebin).
+ *
+ * @param time The (global) time at which this signal is received.
+ * @param receive_vector The incident unit vector. (not used at the moment)
+ * @param field The incident electric field vector.
*
*/
- class TimeDomainAntenna : public Antenna<TimeDomainAntenna> {
-
-
-// protected:
-// // expose the CRTP interfaces constructor
-
- public:
- // import the methods from the antenna
-
- TimeType const start_time_; ///< The start time of this waveform.
- TimeType const duration_; ///< The duration of this waveform.
- InverseTimeType const sample_rate_; ///< The sampling rate of this antenna.
- int num_bins_; ///< The number of bins used.
- xt::xtensor<double,2> waveformE_; ///< The waveform stored by this antenna.
- std::pair<xt::xtensor<double, 2>,
- xt::xtensor<double,2>> waveform_; ///< useful for .getWaveform()
-
- using Antenna<TimeDomainAntenna>::getName;
- using Antenna<TimeDomainAntenna>::getLocation;
-
- /**
- * Construct a new TimeDomainAntenna.
- *
- * @param name The name of this antenna.
- * @param location The location of this antenna.
- * @param start_time The starting time of this waveform.
- * @param duration The duration of this waveform.
- * @param sample_rate The sample rate of this waveform.
- * @param num_bins_ The number of timebins to store E-field.
- * @param waveformE_ The xtensor initialized to zero for E-field.
- *
- */
- TimeDomainAntenna(std::string const& name, Point const& location,
- TimeType const& start_time,
- TimeType const& duration,
- InverseTimeType const& sample_rate)
- : Antenna(name, location)
- , start_time_(start_time)
- , duration_(duration)
- , sample_rate_(sample_rate)
- , num_bins_ (static_cast<int>(duration * sample_rate))
- , waveformE_ (xt::zeros<double>({num_bins_, 3}))
- {};
-
- // copy constructor
- TimeDomainAntenna(const TimeDomainAntenna& Tant)
- : Antenna(Tant.name_, Tant.location_)
- , start_time_(Tant.start_time_)
- , duration_(Tant.duration_)
- , sample_rate_(Tant.sample_rate_)
- , num_bins_(Tant.num_bins_)
- , waveformE_(Tant.waveformE_)
- {};
-
- /**
- * Receive an electric field at this antenna.
- *
- * This assumes that the antenna will receive
- * an *instantaneous* electric field modeled as a delta function (or timebin).
- *
- * @param time The (global) time at which this signal is received.
- * @param receive_vector The incident unit vector. (not used at the moment)
- * @param field The incident electric field vector.
- *
- */
- // TODO: rethink this method a bit. If the endpoint is at the end of the antenna resolution then you get the startpoint signal but you lose the endpoint signal!
- void receive(TimeType const time, Vector<dimensionless_d> const& receive_vector,
- ElectricFieldVector const& efield) {
-
- if (time < start_time_ || time > start_time_ + duration_) {
- return;
- } else {
- // figure out the correct timebin to store the E-field value.
- // NOTE: static cast is implicitly flooring
- auto timebin_ {static_cast<std::size_t>((time - start_time_) * sample_rate_)};
- std::cout << "TIMEBIN IS: " << timebin_ << std::endl;
-
- // store the x,y,z electric field components.
- waveformE_.at(timebin_, 0) += (efield.getX() / (1_V/1_m));
- waveformE_.at(timebin_, 1) += (efield.getY() / (1_V/1_m));
- waveformE_.at(timebin_, 2) += (efield.getZ() / (1_V/1_m));
- //TODO: Check how they are stored in memory, row-wise or column-wise?
- }
- }
-
- /**
- * Get the current waveform for this antenna.
- *
- * NOTE: Currently returns ns and V/m but this should be UNITful
- *
- * @returns A pair of the sample times, and the field
- */
- std::pair<xt::xtensor<double, 2>, xt::xtensor<double,2>> getWaveform() const {
- // TODO: divide by Δt for CoREAS ONLY!
-
- // create a 1-D xtensor to store time values so we can print them later.
- xt::xtensor<double, 2> times_ (xt::zeros<double>({num_bins_, 1}));
-
- for (int i = 0; i < num_bins_; i++) {
- // copy here waveformE_ (maybe that solves the segmentation error)
- times_.at(i,0) = static_cast<double>(start_time_ / 1_s + i / (sample_rate_ * 1_s));
- }
-
- return std::make_pair(times_, waveformE_);
- };
-
- /**
- * Reset the antenna before starting a new simulation.
- */
- void reset() {
- waveformE_ = xt::zeros_like(waveformE_);
-// times_ = xt::zeros_like(times_);
- };
-
- }; // END: class Antenna final
+ // TODO: rethink this method a bit. If the endpoint is at the end of the antenna
+ // resolution then you get the startpoint signal but you lose the endpoint signal!
+ void receive(TimeType const time, Vector<dimensionless_d> const& receive_vector,
+ ElectricFieldVector const& efield) {
+
+ if (time < start_time_ || time > start_time_ + duration_) {
+ return;
+ } else {
+ // figure out the correct timebin to store the E-field value.
+ // NOTE: static cast is implicitly flooring
+ auto timebin_{static_cast<std::size_t>((time - start_time_) * sample_rate_)};
+ std::cout << "TIMEBIN IS: " << timebin_ << std::endl;
+
+ // store the x,y,z electric field components.
+ waveformE_.at(timebin_, 0) += (efield.getX() / (1_V / 1_m));
+ waveformE_.at(timebin_, 1) += (efield.getY() / (1_V / 1_m));
+ waveformE_.at(timebin_, 2) += (efield.getZ() / (1_V / 1_m));
+ // TODO: Check how they are stored in memory, row-wise or column-wise?
+ }
+ }
+
+ /**
+ * Return the time-units of each waveform.
+ *
+ * This returns them in nanoseconds for ease of use.
+ */
+ auto& getData() const { return waveformE_; }
+
+ /**
+ * Return the time-units of each waveform.
+ *
+ * This returns them in nanoseconds for ease of use.
+ */
+ auto getAxis() const {
+
+ // create a 1-D xtensor to store time values so we can print them later.
+ xt::xtensor<double, 1> times(xt::zeros<double>({num_bins_}));
+
+ for (int i = 0; i < num_bins_; i++) {
+ // create the current time in nanoseconds
+ times.at(i) = static_cast<double>(start_time_ / 1_ns + i / (sample_rate_ * 1_ns));
+ }
+
+ return times;
+ }
+
+ auto getWaveform() const { return std::make_pair(getAxis(), waveformE_); }
+
+ /**
+ * Reset the antenna before starting a new simulation.
+ */
+ void reset() { waveformE_ = xt::zeros_like(waveformE_); };
+
+ /**
+ * Return a YAML configuration for this antenna.
+ */
+ YAML::Node getConfig() const {
+
+ // top-level config
+ YAML::Node config;
+
+ config["type"] = "TimeDomainAntenna";
+ config["start_time"] = start_time_ / 1_ns;
+ config["duration"] = duration_ / 1_ns;
+ config["sample_rate"] = sample_rate_ / 1_GHz;
+
+ return config;
+ }
+
+ }; // END: class Antenna final
} // namespace corsika
diff --git a/tests/modules/testRadio.cpp b/tests/modules/testRadio.cpp
index a533b3357..c05f2a6b3 100644
--- a/tests/modules/testRadio.cpp
+++ b/tests/modules/testRadio.cpp
@@ -258,7 +258,7 @@ TEST_CASE("Radio", "[processes]") {
// create a radio process instance using CoREAS
RadioProcess<decltype(detector), CoREAS<decltype(detector), decltype(StraightPropagator(envCoREAS))>, decltype(StraightPropagator(envCoREAS))>
- coreas(detector, envCoREAS);
+ coreas("CoREAS", detector, envCoREAS);
// check doContinuous and simulate methods
coreas.doContinuous(particle1, base, true);
@@ -363,7 +363,7 @@ TEST_CASE("Radio", "[processes]") {
// create a radio process instance using CoREAS
RadioProcess<decltype(detector), ZHS<decltype(detector), decltype(StraightPropagator(envZHS))>, decltype(StraightPropagator(envZHS))>
- zhs(detector, envZHS);
+ zhs("ZHS", detector, envZHS);
// check doContinuous and simulate methods
zhs.doContinuous(particle1, base, true);
@@ -895,7 +895,7 @@ TEST_CASE("Radio", "[processes]") {
// create a radio process instance using CoREAS
RadioProcess<decltype(detector), CoREAS<decltype(detector), decltype(StraightPropagator(env))>, decltype(StraightPropagator(env))>
- coreas(detector, env);
+ coreas("CoREAS", detector, env);
TimeType timeCounter {0._s};
@@ -925,7 +925,7 @@ TEST_CASE("Radio", "[processes]") {
coreas.doContinuous(particle1,track,true);
// get the output
- coreas.writeOutput();
+ // coreas.writeOutput();
}
@@ -1006,7 +1006,7 @@ const HEPEnergyType E0{11.4_MeV};
// create a radio process instance using CoREAS
RadioProcess<decltype(detector), CoREAS<decltype(detector), decltype(StraightPropagator(env))>, decltype(StraightPropagator(env))>
-coreas(detector, env);
+coreas("CoREAS", detector, env);
// loop over all the tracks except the last one
int const n_points {100000};
@@ -1108,7 +1108,7 @@ const HEPEnergyType E0{11.4_MeV};
// create a radio process instance using CoREAS or ZHS
RadioProcess<decltype(detector), CoREAS<decltype(detector), decltype(StraightPropagator(env))>, decltype(StraightPropagator(env))>
-coreas(detector, env);
+coreas("CoREAS", detector, env);
// loop over all the tracks except the last one
int const n_points {60000};
--
GitLab