diff --git a/corsika/detail/modules/TimeCut.inl b/corsika/detail/modules/TimeCut.inl
new file mode 100644
index 0000000000000000000000000000000000000000..3dbdecb7f28dd0c7804308f2841b5e7ae78ebed1
--- /dev/null
+++ b/corsika/detail/modules/TimeCut.inl
@@ -0,0 +1,28 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <corsika/modules/TimeCut.hpp>
+
+namespace corsika {
+
+ inline TimeCut::TimeCut(const TimeType time)
+ : time_(time) {}
+
+ template <typename Particle>
+ inline ProcessReturn TimeCut::doContinuous(Step<Particle> const& step, bool const) {
+ CORSIKA_LOG_TRACE("TimeCut::doContinuous");
+ if (step.getTimePost() >= time_) {
+ CORSIKA_LOG_TRACE("stopping continuous process");
+ return ProcessReturn::ParticleAbsorbed;
+ }
+ return ProcessReturn::Ok;
+ }
+
+} // namespace corsika
diff --git a/corsika/detail/modules/proposal/ContinuousProcess.inl b/corsika/detail/modules/proposal/ContinuousProcess.inl
index 7d0cb7eb56876bcff01b675a9b0430835fa8b82b..9ffb5d8fde6798f9adfcb2d21bc7a21edb2ea840 100644
--- a/corsika/detail/modules/proposal/ContinuousProcess.inl
+++ b/corsika/detail/modules/proposal/ContinuousProcess.inl
@@ -35,7 +35,7 @@ namespace corsika::proposal {
// Use higland multiple scattering and deactivate stochastic deflection by
// passing an empty vector
- static constexpr auto ms_type = PROPOSAL::MultipleScatteringType::Highland;
+ static constexpr auto ms_type = PROPOSAL::MultipleScatteringType::Moliere;
auto s_type = std::vector<PROPOSAL::InteractionType>();
// Build displacement integral and scattering object and interpolate them too and
diff --git a/corsika/detail/modules/radio/CoREAS.inl b/corsika/detail/modules/radio/CoREAS.inl
new file mode 100644
index 0000000000000000000000000000000000000000..9f25ae97265f97cddf955f19ad2eee0c199c01b8
--- /dev/null
+++ b/corsika/detail/modules/radio/CoREAS.inl
@@ -0,0 +1,388 @@
+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <corsika/modules/radio/CoREAS.hpp>
+
+namespace corsika {
+
+ template <typename TRadioDetector, typename TPropagator>
+ template <typename Particle>
+ inline ProcessReturn CoREAS<TRadioDetector, TPropagator>::simulate(
+ Step<Particle> const& step) {
+
+ // get the global simulation time for that track.
+ auto const startTime_{
+ step.getTimePre()}; // time at the start point of the track. I
+ // should use something similar to fCoreHitTime (?)
+ auto const endTime_{step.getTimePost()}; // time at end point of track.
+
+ // LCOV_EXCL_START
+ if (startTime_ == endTime_) {
+ CORSIKA_LOG_ERROR("Time at the start and end of the track coincides! - radio");
+ return ProcessReturn::Ok;
+ // LCOV_EXCL_STOP
+ } else {
+
+ // get start and end position of the track
+ Point const startPoint_{step.getPositionPre()};
+ Point const endPoint_{step.getPositionPost()};
+ // get the coordinate system of the startpoint and hence the track
+ auto const cs_{startPoint_.getCoordinateSystem()};
+
+ auto const currDirection{(endPoint_ - startPoint_).normalized()};
+ // calculate the track length
+ auto const tracklength_{(endPoint_ - startPoint_).getNorm()};
+
+ // beta is velocity / speed of light. Start & end should be the same in endpoints!
+ auto const corrBetaValue{(endPoint_ - startPoint_).getNorm() /
+ (constants::c * (endTime_ - startTime_))};
+ auto const beta_{currDirection * corrBetaValue};
+
+ // get particle charge
+ auto const charge_{get_charge(step.getParticlePre().getPID())};
+
+ // constants for electric field vector calculation
+ auto const constants_{charge_ / (4 * M_PI) / (constants::epsilonZero) /
+ constants::c};
+
+ // set threshold for application of ZHS-like approximation.
+ const double approxThreshold_{1.0e-3};
+
+ // loop over each antenna in the antenna collection (detector)
+ for (auto& antenna : antennas_.getAntennas()) {
+
+ // get the SignalPathCollection (path1) from the start "endpoint" to the antenna.
+ auto paths1{this->propagator_.propagate(
+ startPoint_, antenna.getLocation(),
+ 1_m)}; // TODO: Add the stepsize to .propagate() at some point
+
+ // get the SignalPathCollection (path2) from the end "endpoint" to the antenna.
+ auto paths2{this->propagator_.propagate(endPoint_, antenna.getLocation(), 1_m)};
+
+ // throw an exception if path sizes don't match
+ try {
+ // loop over both paths at once and directly compare 'start' and 'end'
+ // attributes
+ for (size_t i = (paths1.size() == paths2.size()) ? 0 : throw(i = paths1.size());
+ (i < paths1.size() && i < paths2.size()); i++) {
+
+ // calculate preDoppler factor
+ double preDoppler_{1.0 - paths1[i].refractive_index_source_ *
+ beta_.dot(paths1[i].emit_)};
+
+ // check if preDoppler has become zero in case of refractive index of unity
+ // because of numerical limitations here you might need std::fabs(preDoppler)
+ // in the if statement - same with post & mid
+ // LCOV_EXCL_START
+ if (preDoppler_ == 0) {
+ CORSIKA_LOG_WARN("preDoppler factor numerically zero in COREAS");
+ // redo calculation with higher precision
+ auto const& beta_components_{beta_.getComponents(cs_)};
+ auto const& emit_components_{paths1[i].emit_.getComponents(cs_)};
+ long double const indexL_{paths1[i].refractive_index_source_};
+ long double const betaX_{static_cast<double>(beta_components_.getX())};
+ long double const betaY_{static_cast<double>(beta_components_.getY())};
+ long double const betaZ_{static_cast<double>(beta_components_.getZ())};
+ long double const startX_{static_cast<double>(emit_components_.getX())};
+ long double const startY_{static_cast<double>(emit_components_.getY())};
+ long double const startZ_{static_cast<double>(emit_components_.getZ())};
+ long double const doppler =
+ 1.0l -
+ indexL_ * (betaX_ * startX_ + betaY_ * startY_ + betaZ_ * startZ_);
+ preDoppler_ = doppler;
+ }
+ // LCOV_EXCL_STOP
+
+ // calculate postDoppler factor
+ double postDoppler_{1.0 - paths2[i].refractive_index_source_ *
+ beta_.dot(paths2[i].emit_)};
+
+ // check if postDoppler has become zero in case of refractive index of unity
+ // because of numerical limitations
+ // LCOV_EXCL_START
+ if (postDoppler_ == 0) {
+ CORSIKA_LOG_WARN("postDoppler factor numerically zero in CoREAS");
+ // redo calculation with higher precision
+ auto const& beta_components_{beta_.getComponents(cs_)};
+ auto const& emit_components_{paths2[i].emit_.getComponents(cs_)};
+ long double const indexL_{paths2[i].refractive_index_source_};
+ long double const betaX_{static_cast<double>(beta_components_.getX())};
+ long double const betaY_{static_cast<double>(beta_components_.getY())};
+ long double const betaZ_{static_cast<double>(beta_components_.getZ())};
+ long double const endX_{static_cast<double>(emit_components_.getX())};
+ long double const endY_{static_cast<double>(emit_components_.getY())};
+ long double const endZ_{static_cast<double>(emit_components_.getZ())};
+ long double const doppler =
+ 1.0l - indexL_ * (betaX_ * endX_ + betaY_ * endY_ + betaZ_ * endZ_);
+ postDoppler_ = doppler;
+ }
+ // LCOV_EXCL_STOP
+
+ // calculate receive time for startpoint (aka time delay)
+ auto startPointReceiveTime_{
+ startTime_ +
+ paths1[i].propagation_time_}; // TODO: time 0 is when the imaginary
+ // primary hits the ground
+
+ // calculate receive time for endpoint
+ auto endPointReceiveTime_{endTime_ + paths2[i].propagation_time_};
+
+ // get unit vector for startpoint at antenna location
+ auto ReceiveVectorStart_{paths1[i].receive_};
+
+ // get unit vector for endpoint at antenna location
+ auto ReceiveVectorEnd_{paths2[i].receive_};
+
+ // perform ZHS-like calculation close to Cherenkov angle and for refractive
+ // index at antenna location greater than 1
+ if ((paths1[i].refractive_index_destination_ > 1) &&
+ ((std::fabs(preDoppler_) < approxThreshold_) ||
+ (std::fabs(postDoppler_) < approxThreshold_))) {
+ CORSIKA_LOG_DEBUG("Used ZHS-like approximation in CoREAS - radio");
+
+ // clear the existing paths for this particle and track, since we don't need
+ // them anymore
+ paths1.clear();
+ paths2.clear();
+
+ auto const halfVector_{(startPoint_ - endPoint_) * 0.5};
+ auto const midPoint_{endPoint_ + halfVector_};
+
+ // get global simulation time for the middle point of that track.
+ TimeType const midTime_{(startTime_ + endTime_) * 0.5};
+
+ // get the SignalPathCollection (path3) from the middle "endpoint" to the
+ // antenna.
+ auto paths3{
+ this->propagator_.propagate(midPoint_, antenna.getLocation(), 1_m)};
+
+ // now loop over the paths for endpoint that we got above
+ for (auto const& path : paths3) {
+
+ auto const midPointReceiveTime_{midTime_ + path.propagation_time_};
+ double midDoppler_{1.0 -
+ path.refractive_index_source_ * beta_.dot(path.emit_)};
+
+ // check if midDoppler has become zero because of numerical limitations
+ // LCOV_EXCL_START
+ if (midDoppler_ == 0) {
+ CORSIKA_LOG_WARN("midDoppler factor numerically zero in COREAS");
+ // redo calculation with higher precision
+ auto const& beta_components_{beta_.getComponents(cs_)};
+ auto const& emit_components_{path.emit_.getComponents(cs_)};
+ long double const indexL_{path.refractive_index_source_};
+ long double const betaX_{static_cast<double>(beta_components_.getX())};
+ long double const betaY_{static_cast<double>(beta_components_.getY())};
+ long double const betaZ_{static_cast<double>(beta_components_.getZ())};
+ long double const midX_{static_cast<double>(emit_components_.getX())};
+ long double const midY_{static_cast<double>(emit_components_.getY())};
+ long double const midZ_{static_cast<double>(emit_components_.getZ())};
+ long double const doppler =
+ 1.0l - indexL_ * (betaX_ * midX_ + betaY_ * midY_ + betaZ_ * midZ_);
+ midDoppler_ = doppler;
+ }
+ // LCOV_EXCL_STOP
+
+ // change the values of the receive unit vectors of start and end
+ ReceiveVectorStart_ = path.receive_;
+ ReceiveVectorEnd_ = path.receive_;
+
+ // CoREAS calculation -> get ElectricFieldVector for "midPoint"
+ ElectricFieldVector EVmid_ = (path.emit_.cross(path.emit_.cross(beta_))) /
+ midDoppler_ / path.R_distance_ * constants_ *
+ antenna.getSampleRate();
+
+ ElectricFieldVector EV1_{EVmid_};
+ ElectricFieldVector EV2_{EVmid_ * (-1.0)};
+
+ TimeType deltaT_{tracklength_ / (constants::c * corrBetaValue) *
+ std::fabs(midDoppler_)}; // TODO: Caution with this!
+
+ if (startPointReceiveTime_ <
+ endPointReceiveTime_) // EVstart_ arrives earlier
+ {
+ startPointReceiveTime_ = midPointReceiveTime_ - 0.5 * deltaT_;
+ endPointReceiveTime_ = midPointReceiveTime_ + 0.5 * deltaT_;
+ } else // EVend_ arrives earlier
+ {
+ startPointReceiveTime_ = midPointReceiveTime_ + 0.5 * deltaT_;
+ endPointReceiveTime_ = midPointReceiveTime_ - 0.5 * deltaT_;
+ }
+
+ TimeType const gridResolution_{1 / antenna.getSampleRate()};
+ deltaT_ = endPointReceiveTime_ - startPointReceiveTime_;
+
+ // redistribute contributions over time scale defined by the observation
+ // time resolution
+ if (abs(deltaT_) < (gridResolution_)) {
+
+ EV1_ *= std::fabs((deltaT_ / gridResolution_));
+ EV2_ *= std::fabs((deltaT_ / gridResolution_));
+
+ // ToDO: be careful with times in C8!!! where is the zero (time). Is it
+ // close-by?
+ long const startBin = static_cast<long>(
+ std::floor(startPointReceiveTime_ / gridResolution_ + 0.5l));
+ long const endBin = static_cast<long>(
+ std::floor(endPointReceiveTime_ / gridResolution_ + 0.5l));
+ double const startBinFraction =
+ (startPointReceiveTime_ / gridResolution_) -
+ std::floor(startPointReceiveTime_ / gridResolution_);
+ double const endBinFraction =
+ (endPointReceiveTime_ / gridResolution_) -
+ std::floor(endPointReceiveTime_ / gridResolution_);
+
+ // only do timing modification if contributions would land in same bin
+ if (startBin == endBin) {
+
+ // if startE arrives before endE
+ if ((deltaT_) >= 0_s) {
+ if ((startBinFraction >= 0.5) &&
+ (endBinFraction >= 0.5)) // both points left of bin center
+ {
+ startPointReceiveTime_ -=
+ gridResolution_; // shift EV1_ to previous gridpoint
+ } else if ((startBinFraction < 0.5) &&
+ (endBinFraction <
+ 0.5)) // both points right of bin center
+ {
+ endPointReceiveTime_ +=
+ gridResolution_; // shift EV2_ to next gridpoint
+ } else // points on both sides of bin center
+ {
+ double const leftDist = 1.0 - startBinFraction;
+ double const rightDist = endBinFraction;
+ // check if asymmetry to right or left
+ if (rightDist >= leftDist) {
+ endPointReceiveTime_ +=
+ gridResolution_; // shift EV2_ to next gridpoint
+ } else {
+ startPointReceiveTime_ -=
+ gridResolution_; // shift EV1_ to previous gridpoint
+ }
+ }
+ } else // if endE arrives before startE
+ {
+ if ((startBinFraction >= 0.5) &&
+ (endBinFraction >= 0.5)) // both points left of bin center
+ {
+ endPointReceiveTime_ -=
+ gridResolution_; // shift EV2_ to previous gridpoint
+ } else if ((startBinFraction < 0.5) &&
+ (endBinFraction <
+ 0.5)) // both points right of bin center
+ {
+ startPointReceiveTime_ +=
+ gridResolution_; // shift EV1_ to next gridpoint
+ } else // points on both sides of bin center
+ {
+ double const leftDist = 1.0 - endBinFraction;
+ double const rightDist = startBinFraction;
+ // check if asymmetry to right or left
+ if (rightDist >= leftDist) {
+ startPointReceiveTime_ +=
+ gridResolution_; // shift EV1_ to next gridpoint
+ } else {
+ endPointReceiveTime_ -=
+ gridResolution_; // shift EV2_ to previous gridpoint
+ }
+ }
+ } // End of else statement
+ } // End of if for startbin == endbin
+ } // End of if deltaT < gridresolution
+
+ // TODO: Be very careful with this. Maybe the EVs should be fed after the
+ // for loop of paths3
+ antenna.receive(startPointReceiveTime_, ReceiveVectorStart_, EV1_);
+ antenna.receive(endPointReceiveTime_, ReceiveVectorEnd_, EV2_);
+ } // End of looping over paths3
+
+ } // end of ZHS-like approximation
+ else {
+
+ // calculate electric field vector for startpoint
+ ElectricFieldVector EV1_ =
+ (paths1[i].emit_.cross(paths1[i].emit_.cross(beta_))) / preDoppler_ /
+ paths1[i].R_distance_ * constants_ * antenna.getSampleRate();
+
+ // calculate electric field vector for endpoint
+ ElectricFieldVector EV2_ =
+ (paths2[i].emit_.cross(paths2[i].emit_.cross(beta_))) / postDoppler_ /
+ paths2[i].R_distance_ * constants_ * (-1.0) * antenna.getSampleRate();
+
+ if ((preDoppler_ < 1.e-9) || (postDoppler_ < 1.e-9)) {
+
+ TimeType const gridResolution_{1 / antenna.getSampleRate()};
+ TimeType deltaT_{endPointReceiveTime_ - startPointReceiveTime_};
+
+ if (abs(deltaT_) < (gridResolution_)) {
+
+ EV1_ *= std::fabs(deltaT_ / gridResolution_); // Todo: rename EV1 and 2
+ EV2_ *= std::fabs(deltaT_ / gridResolution_);
+
+ long const startBin = static_cast<long>(
+ std::floor(startPointReceiveTime_ / gridResolution_ + 0.5l));
+ long const endBin = static_cast<long>(
+ std::floor(endPointReceiveTime_ / gridResolution_ + 0.5l));
+ double const startBinFraction =
+ (startPointReceiveTime_ / gridResolution_) -
+ std::floor(startPointReceiveTime_ / gridResolution_);
+ double const endBinFraction =
+ (endPointReceiveTime_ / gridResolution_) -
+ std::floor(endPointReceiveTime_ / gridResolution_);
+
+ // only do timing modification if contributions would land in same bin
+ if (startBin == endBin) {
+
+ if ((startBinFraction >= 0.5) &&
+ (endBinFraction >= 0.5)) // both points left of bin center
+ {
+ startPointReceiveTime_ -=
+ gridResolution_; // shift EV1_ to previous gridpoint
+ } else if ((startBinFraction < 0.5) &&
+ (endBinFraction < 0.5)) // both points right of bin center
+ {
+ endPointReceiveTime_ +=
+ gridResolution_; // shift EV2_ to next gridpoint
+ } else // points on both sides of bin center
+ {
+ double const leftDist = 1.0 - startBinFraction;
+ double const rightDist = endBinFraction;
+ // check if asymmetry to right or left
+ if (rightDist >= leftDist) {
+ endPointReceiveTime_ +=
+ gridResolution_; // shift EV2_ to next gridpoint
+ } else {
+ startPointReceiveTime_ -=
+ gridResolution_; // shift EV1_ to previous gridpoint
+ }
+ }
+
+ } // End of if for startbin == endbin
+ } // End of if deltaT < gridresolution
+ } // End of if that checks small doppler factors
+ antenna.receive(startPointReceiveTime_, ReceiveVectorStart_, EV1_);
+ antenna.receive(endPointReceiveTime_, ReceiveVectorEnd_, EV2_);
+ } // End of else that does not perform ZHS-like approximation
+
+ } // End of loop over both paths to get signal info
+ } // End of try block
+ // LCOV_EXCL_START
+ catch (size_t i) {
+ CORSIKA_LOG_ERROR("Signal Paths do not have the same size in CoREAS");
+ }
+ // LCOV_EXCL_STOP
+ } // End of looping over antennas
+
+ return ProcessReturn::Ok;
+ }
+ } // End of simulate method
+
+} // namespace corsika
\ No newline at end of file
diff --git a/corsika/detail/modules/radio/RadioProcess.inl b/corsika/detail/modules/radio/RadioProcess.inl
new file mode 100644
index 0000000000000000000000000000000000000000..0548b71a4226e554ddb20b59ed4fa13b0b69ff2a
--- /dev/null
+++ b/corsika/detail/modules/radio/RadioProcess.inl
@@ -0,0 +1,128 @@
+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+#pragma once
+
+#include <corsika/modules/radio/RadioProcess.hpp>
+
+namespace corsika {
+
+ template <typename TAntennaCollection, typename TRadioImpl, typename TPropagator>
+ inline TRadioImpl&
+ RadioProcess<TAntennaCollection, TRadioImpl, TPropagator>::implementation() {
+ return static_cast<TRadioImpl&>(*this);
+ }
+
+ template <typename TAntennaCollection, typename TRadioImpl, typename TPropagator>
+ inline TRadioImpl const&
+ RadioProcess<TAntennaCollection, TRadioImpl, TPropagator>::implementation() const {
+ return static_cast<TRadioImpl const&>(*this);
+ }
+
+ template <typename TAntennaCollection, typename TRadioImpl, typename TPropagator>
+ template <typename... TArgs>
+ inline RadioProcess<TAntennaCollection, TRadioImpl, TPropagator>::RadioProcess(
+ TAntennaCollection& antennas, TArgs&&... args)
+ : antennas_(antennas)
+ , propagator_(args...) {}
+
+ template <typename TAntennaCollection, typename TRadioImpl, typename TPropagator>
+ template <typename Particle>
+ inline ProcessReturn RadioProcess<TAntennaCollection, TRadioImpl,
+ TPropagator>::doContinuous(const Step<Particle>& step,
+ const bool) {
+ // 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)) {
+ auto const particleID_{step.getParticlePre().getPID()};
+ if ((particleID_ == Code::Electron) || (particleID_ == Code::Positron)) {
+ CORSIKA_LOG_DEBUG("Particle for radio calculation: {} ", particleID_);
+ return this->implementation().simulate(step);
+ } else {
+ CORSIKA_LOG_DEBUG("Particle {} is irrelevant for radio", particleID_);
+ return ProcessReturn::Ok;
+ }
+ //}
+ }
+
+ template <typename TAntennaCollection, typename TRadioImpl, typename TPropagator>
+ template <typename Particle, typename Track>
+ inline LengthType
+ RadioProcess<TAntennaCollection, TRadioImpl, TPropagator>::getMaxStepLength(
+ const Particle& vParticle, const Track& vTrack) const {
+ return meter * std::numeric_limits<double>::infinity();
+ }
+
+ // this should all be moved at a separate radio output function
+ // LCOV_EXCL_START
+ template <typename TAntennaCollection, typename TRadioImpl, typename TPropagator>
+ inline void RadioProcess<TAntennaCollection, TRadioImpl, TPropagator>::startOfLibrary(
+ const boost::filesystem::path& directory) {
+
+ // 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, this->implementation().algorithm);
+ }
+ }
+
+ template <typename TAntennaCollection, typename TRadioImpl, typename TPropagator>
+ inline void RadioProcess<TAntennaCollection, TRadioImpl, TPropagator>::endOfShower(
+ const unsigned int) {
+
+ // 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_, this->implementation().algorithm,
+ antenna.getSampleRate() * 1_s);
+ antenna.reset();
+ }
+
+ // increment our event counter
+ event_++;
+ }
+
+ template <typename TAntennaCollection, typename TRadioImpl, typename TPropagator>
+ inline YAML::Node RadioProcess<TAntennaCollection, TRadioImpl, TPropagator>::getConfig()
+ const {
+
+ // top-level YAML node
+ YAML::Node config;
+
+ // fill in some basics
+ config["type"] = "RadioProcess";
+ config["algorithm"] = this->implementation().algorithm;
+ 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;
+ }
+ // LCOV_EXCL_STOP
+
+} // namespace corsika
\ No newline at end of file
diff --git a/corsika/detail/modules/radio/ZHS.inl b/corsika/detail/modules/radio/ZHS.inl
new file mode 100644
index 0000000000000000000000000000000000000000..5c9e4aacc824b06499b65b988883d9bc1f3cbb00
--- /dev/null
+++ b/corsika/detail/modules/radio/ZHS.inl
@@ -0,0 +1,242 @@
+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <corsika/modules/radio/ZHS.hpp>
+
+namespace corsika {
+
+ template <typename TRadioDetector, typename TPropagator>
+ template <typename Particle>
+ inline ProcessReturn ZHS<TRadioDetector, TPropagator>::simulate(
+ Step<Particle> const& step) const {
+ auto const startTime{step.getTimePre()};
+ auto const endTime{step.getTimePost()};
+
+ // LCOV_EXCL_START
+ if (startTime == endTime) {
+ CORSIKA_LOG_ERROR("Time at the start and end of the track coincides! - radio");
+ return ProcessReturn::Ok;
+ // LCOV_EXCL_STOP
+ } else {
+
+ auto const startPoint{step.getPositionPre()};
+ auto const endPoint{step.getPositionPost()};
+ LengthType const trackLength{(startPoint - endPoint).getNorm()};
+
+ auto const betaModule{(endPoint - startPoint).getNorm() /
+ (constants::c * (endTime - startTime))};
+ auto const beta{(endPoint - startPoint).normalized() * betaModule};
+
+ auto const charge{get_charge(step.getParticlePre().getPID())};
+
+ // // get "mid" position of the track geometrically
+
+ auto const halfVector{(startPoint - endPoint) / 2};
+ auto const midPoint{endPoint + halfVector};
+
+ auto const constants{charge / (4 * M_PI) / (constants::epsilonZero) / constants::c};
+
+ // we loop over each antenna in the collection
+ for (auto& antenna : antennas_.getAntennas()) {
+ auto midPaths{this->propagator_.propagate(midPoint, antenna.getLocation(), 1_m)};
+ // Loop over midPaths, first check Fraunhoffer limit
+ for (size_t i{0}; i < midPaths.size(); i++) {
+ double const uTimesK{beta.dot(midPaths[i].emit_) / betaModule};
+ double const sinTheta2{1. - uTimesK * uTimesK};
+ LengthType const lambda{constants::c / antenna.getSampleRate()};
+ double const fraunhLimit{sinTheta2 * trackLength * trackLength /
+ midPaths[i].R_distance_ / lambda * 2 * M_PI};
+ // Checks if we are in fraunhoffer domain
+ if (fraunhLimit > 1.0) {
+ /// code for dividing track and calculating field.
+ double const nSubTracks{sqrt(fraunhLimit) + 1};
+ auto const step_{(endPoint - startPoint) / nSubTracks};
+ TimeType const timeStep{(endTime - startTime) / nSubTracks};
+ // energy should be divided up when it is possible to get the energy at end of
+ // track!!!!
+ auto point1{startPoint};
+ TimeType time1{startTime};
+ for (int j{0}; j < nSubTracks; j++) {
+ auto const point2{point1 + step_};
+ TimeType const time2{time1 + timeStep};
+ auto const newHalfVector{(point1 - point2) / 2.};
+ auto const newMidPoint{point2 + newHalfVector};
+ auto const newMidPaths{
+ this->propagator_.propagate(newMidPoint, antenna.getLocation(), 1_m)};
+ // A function for calculating the field should be made since it is repeated
+ // later
+ for (size_t k{0}; k < newMidPaths.size(); k++) {
+
+ double const n_source{newMidPaths[k].refractive_index_source_};
+
+ double const betaTimesK{beta.dot(newMidPaths[k].emit_)};
+ TimeType const midTime{(time1 + time2) / 2.};
+ TimeType detectionTime1{time1 + newMidPaths[k].propagation_time_ -
+ n_source * betaTimesK * (time1 - midTime)};
+ TimeType detectionTime2{time2 + newMidPaths[k].propagation_time_ -
+ n_source * betaTimesK * (time2 - midTime)};
+
+ // make detectionTime1_ be the smallest time =>
+ // changes step function order so sign is changed to account for it
+ double sign = 1.;
+ if (detectionTime1 > detectionTime2) {
+ detectionTime1 = time2 + newMidPaths[k].propagation_time_ -
+ n_source * betaTimesK * (time2 - midTime);
+ detectionTime2 = time1 + newMidPaths[k].propagation_time_ -
+ n_source * betaTimesK * (time1 - midTime);
+ sign = -1.;
+ } // end if statement for time structure
+
+ double const startBin{std::floor(
+ (detectionTime1 - antenna.getStartTime()) * antenna.getSampleRate() +
+ 0.5)};
+ double const endBin{std::floor((detectionTime2 - antenna.getStartTime()) *
+ antenna.getSampleRate() +
+ 0.5)};
+
+ auto const betaPerp{
+ newMidPaths[k].emit_.cross(beta.cross(newMidPaths[k].emit_))};
+ double const denominator{1. - n_source * betaTimesK};
+
+ if (startBin == endBin) {
+ // track contained in bin
+ // if not in Cerenkov angle then
+ if (std::fabs(denominator) > 1.e-15) {
+ double const f{std::fabs((detectionTime2 * antenna.getSampleRate() -
+ detectionTime1 * antenna.getSampleRate()))};
+ VectorPotential const Vp = betaPerp * sign * constants * f /
+ denominator / newMidPaths[k].R_distance_;
+ antenna.receive(detectionTime2, betaPerp, Vp);
+ } else { // If emission in Cerenkov angle => approximation
+ double const f{time2 * antenna.getSampleRate() -
+ time1 * antenna.getSampleRate()};
+ VectorPotential const Vp =
+ betaPerp * sign * constants * f / newMidPaths[k].R_distance_;
+ 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)};
+ // first contribution/ plus 1 bin minus 0.5 from new antenna ruonding
+ double f{std::fabs(startBin + 0.5 -
+ (detectionTime1 - antenna.getStartTime()) *
+ antenna.getSampleRate())};
+ VectorPotential Vp = betaPerp * sign * constants * f / denominator /
+ newMidPaths[k].R_distance_;
+ antenna.receive(detectionTime1, betaPerp, Vp);
+ // intermidiate contributions
+ for (int it{1}; it < numberOfBins; ++it) {
+ Vp = betaPerp * constants / denominator / newMidPaths[k].R_distance_;
+ antenna.receive(detectionTime1 +
+ static_cast<double>(it) / antenna.getSampleRate(),
+ betaPerp, Vp);
+ } // end loop over bins in which potential vector is not zero
+ // final contribution// f +0.5 from new antenna rounding
+ f = std::fabs((detectionTime2 - antenna.getStartTime()) *
+ antenna.getSampleRate() +
+ 0.5 - endBin);
+ Vp = betaPerp * sign * constants * f / denominator /
+ newMidPaths[k].R_distance_;
+ antenna.receive(detectionTime2, betaPerp, Vp);
+ } // end if statement for track in multiple bins
+
+ } // end of loop over newMidPaths
+ // update points for next sub track
+ point1 = point1 + step_;
+ time1 = time1 + timeStep;
+ }
+
+ } else // Calculate vector potential of whole track
+ {
+ double const n_source{midPaths[i].refractive_index_source_};
+
+ double const betaTimesK{beta.dot(midPaths[i].emit_)};
+ TimeType const midTime{(startTime + endTime) / 2};
+ TimeType detectionTime1{startTime + midPaths[i].propagation_time_ -
+ n_source * betaTimesK * (startTime - midTime)};
+ TimeType detectionTime2{endTime + midPaths[i].propagation_time_ -
+ n_source * betaTimesK * (endTime - midTime)};
+
+ // make detectionTime1_ be the smallest time =>
+ // changes step function order so sign is changed to account for it
+ double sign = 1.;
+ if (detectionTime1 > detectionTime2) {
+ detectionTime1 = endTime + midPaths[i].propagation_time_ -
+ n_source * betaTimesK * (endTime - midTime);
+ detectionTime2 = startTime + midPaths[i].propagation_time_ -
+ n_source * betaTimesK * (startTime - midTime);
+ sign = -1.;
+ } // end if statement for time structure
+
+ double const startBin{std::floor((detectionTime1 - antenna.getStartTime()) *
+ antenna.getSampleRate() +
+ 0.5)};
+ double const endBin{std::floor((detectionTime2 - antenna.getStartTime()) *
+ antenna.getSampleRate() +
+ 0.5)};
+
+ auto const betaPerp{midPaths[i].emit_.cross(beta.cross(midPaths[i].emit_))};
+ double const denominator{1. -
+ midPaths[i].refractive_index_source_ * betaTimesK};
+
+ if (startBin == endBin) {
+ // track contained in bin
+ // if not in Cerenkov angle then
+ if (std::fabs(denominator) > 1.e-15) {
+ double const f{std::fabs((detectionTime2 * antenna.getSampleRate() -
+ detectionTime1 * antenna.getSampleRate()))};
+
+ VectorPotential const Vp = betaPerp * sign * constants * f / denominator /
+ midPaths[i].R_distance_;
+ antenna.receive(detectionTime2, betaPerp, Vp);
+ } else { // If emission in Cerenkov angle => approximation
+ double const f{endTime * antenna.getSampleRate() -
+ startTime * antenna.getSampleRate()};
+ VectorPotential const Vp =
+ betaPerp * sign * constants * f / midPaths[i].R_distance_;
+ 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 + 0.5 -
+ (detectionTime1 - antenna.getStartTime()) *
+ antenna.getSampleRate())};
+ VectorPotential Vp =
+ betaPerp * sign * constants * f / denominator / midPaths[i].R_distance_;
+ antenna.receive(detectionTime1, betaPerp, Vp);
+ // intermidiate contributions
+ for (int it{1}; it < numberOfBins; ++it) {
+ Vp = betaPerp * sign * constants / denominator / midPaths[i].R_distance_;
+ antenna.receive(
+ detectionTime1 + static_cast<double>(it) / antenna.getSampleRate(),
+ betaPerp, Vp);
+ } // end loop over bins in which potential vector is not zero
+ // final contribution
+ f = std::fabs((detectionTime2 - antenna.getStartTime()) *
+ antenna.getSampleRate() +
+ 0.5 - endBin);
+ Vp =
+ betaPerp * sign * constants * f / denominator / midPaths[i].R_distance_;
+ 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
+
+} // namespace corsika
\ No newline at end of file
diff --git a/corsika/detail/modules/radio/antennas/Antenna.inl b/corsika/detail/modules/radio/antennas/Antenna.inl
new file mode 100644
index 0000000000000000000000000000000000000000..0b97560aaefc129779fdf1e39dfdfe8f68e9e61b
--- /dev/null
+++ b/corsika/detail/modules/radio/antennas/Antenna.inl
@@ -0,0 +1,115 @@
+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+#pragma once
+
+#include <corsika/modules/radio/antennas/Antenna.hpp>
+
+namespace corsika {
+
+ template <typename TAntennaImpl>
+ inline Antenna<TAntennaImpl>::Antenna(std::string const& name, Point const& location,
+ CoordinateSystemPtr const& coordinateSystem)
+ : name_(name)
+ , location_(location)
+ , coordinateSystem_(coordinateSystem){};
+
+ template <typename TAntennaImpl>
+ inline Point const& Antenna<TAntennaImpl>::getLocation() const {
+ return location_;
+ }
+
+ template <typename TAntennaImpl>
+ inline std::string const& Antenna<TAntennaImpl>::getName() const {
+ return name_;
+ }
+
+ template <typename TAntennaImpl>
+ inline void Antenna<TAntennaImpl>::startOfLibrary(
+ const boost::filesystem::path& directory, const std::string radioImplementation) {
+
+ // calculate and save our filename
+ filename_ = (directory / this->getName()).string() + ".npz";
+
+ // get the axis labels for this antenna and write the first row.
+ axistype axis = 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";
+ }
+
+ if (radioImplementation == "ZHS" && TAntennaImpl::is_time_domain) {
+ for (size_t i = 0; i < axis.size() - 1; i++) {
+ axis.at(i) = (axis.at(i + 1) + axis.at(i)) / 2.;
+ }
+ // explicitly convert the arrays to the needed type for cnpy
+ axis.pop_back();
+ long double const* raw_data = axis.data();
+ std::vector<size_t> N = {
+ axis.size()}; // cnpy needs a vector here -- this should be axis.size() - 1 --
+ // write the labels to the first row of the NumPy file
+ cnpy::npz_save(filename_, label, raw_data, N, "w");
+ } else {
+ // explicitly convert the arrays to the needed type for cnpy
+ long double 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");
+ }
+ }
+
+ template <typename TAntennaImpl>
+ inline void Antenna<TAntennaImpl>::endOfShower(const int event,
+ std::string const& radioImplementation,
+ const double sampleRate) {
+
+ // get the copy of the waveform data for this event
+ // we transpose it so that we can match dimensions with the
+ // time array that is already in the output file
+ std::vector<double> const& dataX = this->implementation().getWaveformX();
+ std::vector<double> const& dataY = this->implementation().getWaveformY();
+ std::vector<double> const& dataZ = this->implementation().getWaveformZ();
+
+ if (radioImplementation == "ZHS") {
+ std::vector<double> electricFieldX(dataX.size() - 1,
+ 0); // num_bins_, std::vector<double>(3, 0)
+ std::vector<double> electricFieldY(dataY.size() - 1, 0);
+ std::vector<double> electricFieldZ(dataZ.size() - 1, 0);
+ for (size_t i = 0; i < electricFieldX.size(); i++) {
+ electricFieldX.at(i) = -(dataX.at(i + 1) - dataX.at(i)) * sampleRate;
+ electricFieldY.at(i) = -(dataY.at(i + 1) - dataY.at(i)) * sampleRate;
+ electricFieldZ.at(i) = -(dataZ.at(i + 1) - dataZ.at(i)) * sampleRate;
+ }
+ // cnpy needs a vector for the shape
+ cnpy::npz_save(filename_, std::to_string(event) + "X", electricFieldX.data(),
+ {electricFieldX.size()}, "a");
+ cnpy::npz_save(filename_, std::to_string(event) + "Y", electricFieldY.data(),
+ {electricFieldY.size()}, "a");
+ cnpy::npz_save(filename_, std::to_string(event) + "Z", electricFieldZ.data(),
+ {electricFieldZ.size()}, "a");
+ } else {
+ // cnpy needs a vector for the shape
+ // and write this event to the .npz archive
+ cnpy::npz_save(filename_, std::to_string(event) + "X", dataX.data(), {dataX.size()},
+ "a");
+ cnpy::npz_save(filename_, std::to_string(event) + "Y", dataY.data(), {dataY.size()},
+ "a");
+ cnpy::npz_save(filename_, std::to_string(event) + "Z", dataZ.data(), {dataZ.size()},
+ "a");
+ }
+ }
+
+ template <typename TAntennaImpl>
+ inline TAntennaImpl& Antenna<TAntennaImpl>::implementation() {
+ return static_cast<TAntennaImpl&>(*this);
+ }
+
+} // namespace corsika
diff --git a/corsika/detail/modules/radio/antennas/TimeDomainAntenna.inl b/corsika/detail/modules/radio/antennas/TimeDomainAntenna.inl
new file mode 100644
index 0000000000000000000000000000000000000000..920b303b9d63281e708ad59b278b8089a4a98cf8
--- /dev/null
+++ b/corsika/detail/modules/radio/antennas/TimeDomainAntenna.inl
@@ -0,0 +1,130 @@
+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+#pragma once
+
+#include <corsika/modules/radio/antennas/TimeDomainAntenna.hpp>
+
+namespace corsika {
+
+ inline TimeDomainAntenna::TimeDomainAntenna(std::string const& name,
+ Point const& location,
+ CoordinateSystemPtr coordinateSystem,
+ TimeType const& start_time,
+ TimeType const& duration,
+ InverseTimeType const& sample_rate,
+ TimeType const ground_hit_time)
+ : Antenna(name, location, coordinateSystem)
+ , start_time_(start_time)
+ , duration_(duration)
+ , sample_rate_(sample_rate)
+ , ground_hit_time_(ground_hit_time)
+ , num_bins_(static_cast<std::size_t>(duration * sample_rate + 1.5l))
+ , waveformEX_(num_bins_, 0)
+ , waveformEY_(num_bins_, 0)
+ , waveformEZ_(num_bins_, 0)
+ , time_axis_(createTimeAxis()){};
+
+ inline void TimeDomainAntenna::receive(const TimeType time,
+ const Vector<dimensionless_d>& receive_vector,
+ const ElectricFieldVector& 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>(
+ std::floor((time - start_time_) * sample_rate_ + 0.5l))};
+
+ // store the x,y,z electric field components.
+ auto const& Electric_field_components{efield.getComponents(coordinateSystem_)};
+ waveformEX_.at(timebin_) += (Electric_field_components.getX() * (1_m / 1_V));
+ waveformEY_.at(timebin_) += (Electric_field_components.getY() * (1_m / 1_V));
+ waveformEZ_.at(timebin_) += (Electric_field_components.getZ() * (1_m / 1_V));
+ // TODO: Check how they are stored in memory, row-wise or column-wise? Probably use
+ // a 3D object
+ }
+ }
+
+ inline void TimeDomainAntenna::receive(const TimeType time,
+ const Vector<dimensionless_d>& receive_vector,
+ const VectorPotential& vectorP) {
+
+ 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>(
+ std::floor((time - start_time_) * sample_rate_ + 0.5l))};
+
+ // store the x,y,z electric field components.
+ auto const& Vector_potential_components{vectorP.getComponents(coordinateSystem_)};
+ waveformEX_.at(timebin_) +=
+ (Vector_potential_components.getX() * (1_m / (1_V * 1_s)));
+ waveformEY_.at(timebin_) +=
+ (Vector_potential_components.getY() * (1_m / (1_V * 1_s)));
+ waveformEZ_.at(timebin_) +=
+ (Vector_potential_components.getZ() * (1_m / (1_V * 1_s)));
+ // TODO: Check how they are stored in memory, row-wise or column-wise? Probably use
+ // a 3D object
+ }
+ }
+
+ inline auto const& TimeDomainAntenna::getWaveformX() const { return waveformEX_; }
+
+ inline auto const& TimeDomainAntenna::getWaveformY() const { return waveformEY_; }
+
+ inline auto const& TimeDomainAntenna::getWaveformZ() const { return waveformEZ_; }
+
+ inline std::vector<long double> TimeDomainAntenna::createTimeAxis() const {
+
+ // create a 1-D xtensor to store time values so we can print them later.
+ std::vector<long double> times(num_bins_, 0);
+
+ // calculate the sample_period
+ auto sample_period{1 / sample_rate_};
+
+ // fill in every time-value
+ for (std::size_t i = 0; i < num_bins_; i++) {
+ // create the current time in nanoseconds
+ times.at(i) = static_cast<long double>(
+ ((start_time_ - ground_hit_time_) + i * sample_period) / 1_ns);
+ }
+
+ return times;
+ }
+
+ inline auto const& TimeDomainAntenna::getAxis() const { return time_axis_; }
+
+ inline InverseTimeType const& TimeDomainAntenna::getSampleRate() const {
+ return sample_rate_;
+ }
+
+ inline TimeType const& TimeDomainAntenna::getStartTime() const { return start_time_; }
+
+ inline void TimeDomainAntenna::reset() {
+ std::fill(waveformEX_.begin(), waveformEX_.end(), 0);
+ std::fill(waveformEY_.begin(), waveformEY_.end(), 0);
+ std::fill(waveformEZ_.begin(), waveformEZ_.end(), 0);
+ }
+
+ inline YAML::Node TimeDomainAntenna::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;
+ }
+
+} // namespace corsika
diff --git a/corsika/detail/modules/radio/detectors/AntennaCollection.inl b/corsika/detail/modules/radio/detectors/AntennaCollection.inl
new file mode 100644
index 0000000000000000000000000000000000000000..d985bb98027ccd9f89b202c8352355dccd3da35f
--- /dev/null
+++ b/corsika/detail/modules/radio/detectors/AntennaCollection.inl
@@ -0,0 +1,45 @@
+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+#pragma once
+
+#include <corsika/modules/radio/detectors/AntennaCollection.hpp>
+
+namespace corsika {
+
+ template <typename TAntennaImpl>
+ inline void AntennaCollection<TAntennaImpl>::addAntenna(TAntennaImpl const& antenna) {
+ antennas_.push_back(antenna);
+ }
+
+ template <typename TAntennaImpl>
+ inline TAntennaImpl& AntennaCollection<TAntennaImpl>::at(std::size_t const i) {
+ return antennas_.at(i);
+ }
+
+ template <typename TAntennaImpl>
+ inline TAntennaImpl const& AntennaCollection<TAntennaImpl>::at(
+ std::size_t const i) const {
+ return antennas_.at(i);
+ }
+
+ template <typename TAntennaImpl>
+ inline int AntennaCollection<TAntennaImpl>::size() const {
+ return antennas_.size();
+ }
+
+ template <typename TAntennaImpl>
+ inline std::vector<TAntennaImpl>& AntennaCollection<TAntennaImpl>::getAntennas() {
+ return antennas_;
+ }
+
+ template <typename TAntennaImpl>
+ inline void AntennaCollection<TAntennaImpl>::reset() {
+ std::for_each(antennas_.begin(), antennas_.end(), std::mem_fn(&TAntennaImpl::reset));
+ }
+
+} // namespace corsika
diff --git a/corsika/detail/modules/radio/propagators/RadioPropagator.inl b/corsika/detail/modules/radio/propagators/RadioPropagator.inl
new file mode 100644
index 0000000000000000000000000000000000000000..e5c65f97ef3518a5bf47467c8b3650956fcf3f65
--- /dev/null
+++ b/corsika/detail/modules/radio/propagators/RadioPropagator.inl
@@ -0,0 +1,18 @@
+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+#pragma once
+
+#include <corsika/modules/radio/propagators/RadioPropagator.hpp>
+
+namespace corsika {
+
+ template <typename TImpl, typename TEnvironment>
+ inline RadioPropagator<TImpl, TEnvironment>::RadioPropagator(TEnvironment const& env)
+ : env_(env) {}
+
+} // namespace corsika
diff --git a/corsika/detail/modules/radio/propagators/SignalPath.inl b/corsika/detail/modules/radio/propagators/SignalPath.inl
new file mode 100644
index 0000000000000000000000000000000000000000..7de9c8419a9e75bde2fd556d717e1135028fa0cb
--- /dev/null
+++ b/corsika/detail/modules/radio/propagators/SignalPath.inl
@@ -0,0 +1,29 @@
+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <corsika/modules/radio/propagators/SignalPath.hpp>
+
+namespace corsika {
+
+ inline SignalPath::SignalPath(
+ TimeType const propagation_time, double const average_refractive_index,
+ double const refractive_index_source, double const refractive_index_destination,
+ Vector<dimensionless_d> const& emit, Vector<dimensionless_d> const& receive,
+ LengthType const R_distance, std::deque<Point> const& points)
+ : Path(points)
+ , propagation_time_(propagation_time)
+ , average_refractive_index_(average_refractive_index)
+ , refractive_index_source_(refractive_index_source)
+ , refractive_index_destination_(refractive_index_destination)
+ , emit_(emit)
+ , receive_(receive)
+ , R_distance_(R_distance) {}
+
+} // namespace corsika
\ No newline at end of file
diff --git a/corsika/detail/modules/radio/propagators/SimplePropagator.inl b/corsika/detail/modules/radio/propagators/SimplePropagator.inl
new file mode 100644
index 0000000000000000000000000000000000000000..4101878983cb7f3060e0743392a08caca9689508
--- /dev/null
+++ b/corsika/detail/modules/radio/propagators/SimplePropagator.inl
@@ -0,0 +1,70 @@
+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+#pragma once
+
+#include <corsika/modules/radio/propagators/SimplePropagator.hpp>
+
+namespace corsika {
+
+ template <typename TEnvironment>
+ inline SimplePropagator<TEnvironment>::SimplePropagator(TEnvironment const& env)
+ : RadioPropagator<SimplePropagator, TEnvironment>(env){};
+
+ template <typename TEnvironment>
+ inline typename SimplePropagator<TEnvironment>::SignalPathCollection
+ SimplePropagator<TEnvironment>::propagate(Point const& source, Point const& destination,
+ LengthType const stepsize) const {
+
+ /**
+ * This is the simplest case of straight propagator
+ * where no integration takes place.
+ * This can be used for fast tests and checks of the radio module.
+ *
+ */
+
+ // these are used for the direction of emission and reception of signal at the antenna
+ auto const emit_{(destination - source).normalized()};
+ auto const receive_{-emit_};
+
+ // the geometrical distance from the point of emission to an observer
+ auto const distance_{(destination - source).getNorm()};
+
+ // get the universe for this environment
+ auto const* const universe{Base::env_.getUniverse().get()};
+
+ // the points that consist the signal path (source & destination).
+ std::deque<Point> points;
+
+ // store value of the refractive index at points.
+ std::vector<double> rindex;
+ rindex.reserve(2);
+
+ // get and store the refractive index of the first point 'source'.
+ auto const* const nodeSource{universe->getContainingNode(source)};
+ auto const ri_source{nodeSource->getModelProperties().getRefractiveIndex(source)};
+ rindex.push_back(ri_source);
+ points.push_back(source);
+
+ // add the refractive index of last point 'destination' and store it.
+ auto const* const node{universe->getContainingNode(destination)};
+ auto const ri_destination{node->getModelProperties().getRefractiveIndex(destination)};
+ rindex.push_back(ri_destination);
+ points.push_back(destination);
+
+ // compute the average refractive index.
+ auto const averageRefractiveIndex_ = (ri_source + ri_destination) / 2;
+
+ // compute the total time delay.
+ TimeType const time = averageRefractiveIndex_ * (distance_ / constants::c);
+
+ return {SignalPath(time, averageRefractiveIndex_, ri_source, ri_destination, emit_,
+ receive_, distance_, points)};
+
+ } // END: propagate()
+
+} // namespace corsika
diff --git a/corsika/detail/modules/radio/propagators/StraightPropagator.inl b/corsika/detail/modules/radio/propagators/StraightPropagator.inl
new file mode 100644
index 0000000000000000000000000000000000000000..9962fc133003d855ab682aa15b12c275e225c035
--- /dev/null
+++ b/corsika/detail/modules/radio/propagators/StraightPropagator.inl
@@ -0,0 +1,169 @@
+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+#pragma once
+
+#include <corsika/modules/radio/propagators/StraightPropagator.hpp>
+
+namespace corsika {
+
+ template <typename TEnvironment>
+ // TODO: maybe the constructor doesn't take any arguments for the environment (?)
+ inline StraightPropagator<TEnvironment>::StraightPropagator(TEnvironment const& env)
+ : RadioPropagator<StraightPropagator, TEnvironment>(env){};
+
+ template <typename TEnvironment>
+ inline typename StraightPropagator<TEnvironment>::SignalPathCollection
+ StraightPropagator<TEnvironment>::propagate(Point const& source,
+ Point const& destination,
+ LengthType const stepsize) const {
+
+ /*
+ * get the normalized (unit) vector from `source` to `destination'.
+ * this is also the `emit` and `receive` vectors in the SignalPath class.
+ * in this case emit and receive unit vectors should be the same
+ * so they are both called direction
+ */
+
+ // these are used for the direction of emission and reception of signal at the antenna
+ auto const emit_{(destination - source).normalized()};
+ auto const receive_{-emit_};
+
+ // the distance from the point of emission to an observer
+ auto const distance_{(destination - source).getNorm()};
+
+ try {
+ if (stepsize <= 0.5 * distance_) {
+
+ // "step" is the direction vector with length `stepsize`
+ auto const step{emit_ * stepsize};
+
+ // calculate the number of points (roughly) for the numerical integration
+ auto const n_points{(destination - source).getNorm() / stepsize};
+
+ // get the universe for this environment
+ auto const* const universe{Base::env_.getUniverse().get()};
+
+ // the points that we build along the way for the numerical integration
+ std::deque<Point> points;
+
+ // store value of the refractive index at points
+ std::vector<double> rindex;
+ rindex.reserve(n_points);
+
+ // get and store the refractive index of the first point 'source'
+ auto const* const nodeSource{universe->getContainingNode(source)};
+ auto const ri_source{nodeSource->getModelProperties().getRefractiveIndex(source)};
+ rindex.push_back(ri_source);
+ points.push_back(source);
+
+ // loop from `source` to `destination` to store values before Simpson's rule.
+ // this loop skips the last point 'destination' and "misses" the extra point
+ for (auto point = source + step; (point - destination).getNorm() > 0.6 * stepsize;
+ point = point + step) {
+
+ // get the environment node at this specific 'point'
+ auto const* const node{universe->getContainingNode(point)};
+
+ // get the associated refractivity at 'point'
+ auto const refractive_index{
+ node->getModelProperties().getRefractiveIndex(point)};
+ // auto const refractive_index{1.000327};
+ rindex.push_back(refractive_index);
+
+ // add this 'point' to our deque collection
+ points.push_back(point);
+ }
+
+ // Get the extra points that the for loop misses until the destination
+ auto const extrapoint_{points.back() + step};
+
+ // add the refractive index of last point 'destination' and store it
+ auto const* const node{universe->getContainingNode(destination)};
+ auto const ri_destination{
+ node->getModelProperties().getRefractiveIndex(destination)};
+ // auto const ri_destination{1.000327};
+ rindex.push_back(ri_destination);
+ points.push_back(destination);
+
+ auto N = rindex.size();
+ std::size_t index = 0;
+ double sum = rindex.at(index);
+ auto refra_ = rindex.at(index);
+ TimeType time{0_s};
+
+ if ((N - 1) % 2 == 0) {
+ // Apply the standard Simpson's rule
+ auto const h = ((destination - source).getNorm()) / (N - 1);
+
+ for (std::size_t index = 1; index < (N - 1); index += 2) {
+ sum += 4 * rindex.at(index);
+ refra_ += rindex.at(index);
+ }
+ for (std::size_t index = 2; index < (N - 1); index += 2) {
+ sum += 2 * rindex.at(index);
+ refra_ += rindex.at(index);
+ }
+ index = N - 1;
+ sum = sum + rindex.at(index);
+ refra_ += rindex.at(index);
+
+ // compute the total time delay.
+ time = sum * (h / (3 * constants::c));
+ } else {
+ // Apply Simpson's rule for one "extra" point and then subtract the difference
+ points.pop_back();
+ rindex.pop_back();
+ auto const* const node{universe->getContainingNode(extrapoint_)};
+ auto const ri_extrapoint{
+ node->getModelProperties().getRefractiveIndex(extrapoint_)};
+ rindex.push_back(ri_extrapoint);
+ points.push_back(extrapoint_);
+ auto const extrapoint2_{extrapoint_ + step};
+ auto const* const node2{universe->getContainingNode(extrapoint2_)};
+ auto const ri_extrapoint2{
+ node2->getModelProperties().getRefractiveIndex(extrapoint2_)};
+ rindex.push_back(ri_extrapoint2);
+ points.push_back(extrapoint2_);
+ N = rindex.size();
+ auto const h = ((extrapoint2_ - source).getNorm()) / (N - 1);
+ for (std::size_t index = 1; index < (N - 1); index += 2) {
+ sum += 4 * rindex.at(index);
+ refra_ += rindex.at(index);
+ }
+ for (std::size_t index = 2; index < (N - 1); index += 2) {
+ sum += 2 * rindex.at(index);
+ refra_ += rindex.at(index);
+ }
+ index = N - 1;
+ sum = sum + rindex.at(index);
+ refra_ += rindex.at(index);
+
+ // compute the total time delay including the correction
+ time =
+ sum * (h / (3 * constants::c)) -
+ (ri_extrapoint2 * ((extrapoint2_ - destination).getNorm()) / constants::c);
+ }
+
+ // uncomment the following if you want to skip the integration for fast tests
+ // TimeType time = ri_destination * (distance_ / constants::c);
+
+ // compute the average refractive index.
+ auto const averageRefractiveIndex_ = refra_ / N;
+
+ return {SignalPath(time, averageRefractiveIndex_, ri_source, ri_destination,
+ emit_, receive_, distance_, points)};
+ } else {
+ throw stepsize;
+ }
+ } catch (const LengthType& s) {
+ CORSIKA_LOG_ERROR("Please choose a smaller stepsize for the numerical integration");
+ }
+
+ } // END: propagate()
+
+} // namespace corsika
diff --git a/corsika/framework/core/PhysicalConstants.hpp b/corsika/framework/core/PhysicalConstants.hpp
index 92df20c2018ef382aac191517116d98363155f25..5fd079bd2985f9ece15619964da7a2483ee9576f 100644
--- a/corsika/framework/core/PhysicalConstants.hpp
+++ b/corsika/framework/core/PhysicalConstants.hpp
@@ -35,6 +35,10 @@ namespace corsika::constants {
// elementary charge
constexpr quantity<electric_charge_d> e{Rep(1.6021766208e-19L) * coulomb};
+ // vacuum permittivity
+ constexpr quantity<dimensions<-3, -1, 4, 2>> epsilonZero{Rep(8.8541878128e-12L) *
+ farad / meter};
+
// electronvolt
// constexpr quantity<hepenergy_d> eV{e / coulomb * joule};
diff --git a/corsika/framework/core/PhysicalGeometry.hpp b/corsika/framework/core/PhysicalGeometry.hpp
index 6a1227de457dc5febde7f4d730ca74893972fdd6..e0ea84b10d38899161b8dade986618311060e648 100644
--- a/corsika/framework/core/PhysicalGeometry.hpp
+++ b/corsika/framework/core/PhysicalGeometry.hpp
@@ -43,4 +43,14 @@ namespace corsika {
**/
using LengthVector = Vector<length_d>;
+ /**
+ * A 3D vector defined in a specific coordinate system with units ElectricFieldType
+ **/
+ typedef Vector<ElectricFieldType::dimension_type> ElectricFieldVector;
+
+ /**
+ * A 3D vector defined in a specific coordinate system with units VectorPotentialType
+ **/
+ typedef Vector<VectorPotentialType::dimension_type> VectorPotential;
+
} // namespace corsika
diff --git a/corsika/framework/core/PhysicalUnits.hpp b/corsika/framework/core/PhysicalUnits.hpp
index 82373acc9057817ccde077ad93adac8f3584360b..ba04c819969448fc586733c8e642cdc5a6b680d2 100644
--- a/corsika/framework/core/PhysicalUnits.hpp
+++ b/corsika/framework/core/PhysicalUnits.hpp
@@ -96,6 +96,10 @@ namespace corsika::units::si {
phys::units::quantity<phys::units::dimensions<2, -1, 0>, double>;
using MagneticFluxType =
phys::units::quantity<phys::units::magnetic_flux_density_d, double>;
+ using ElectricFieldType =
+ phys::units::quantity<phys::units::dimensions<1, 1, -3, -1>, double>;
+ using VectorPotentialType =
+ phys::units::quantity<phys::units::dimensions<1, 1, -2, -1>, double>;
template <typename DimFrom, typename DimTo>
auto constexpr conversion_factor_HEP_to_SI() {
diff --git a/corsika/modules/TimeCut.hpp b/corsika/modules/TimeCut.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..2c846ad0ca5dbb92266321469a46b79789ec79c7
--- /dev/null
+++ b/corsika/modules/TimeCut.hpp
@@ -0,0 +1,44 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/process/ContinuousProcess.hpp>
+#include <corsika/framework/core/Step.hpp>
+
+#include <corsika/setup/SetupStack.hpp>
+#include <corsika/setup/SetupTrajectory.hpp>
+
+namespace corsika {
+
+ /**
+ * Simple TimeCut process, removes particles older than specified cut time.
+ */
+ class TimeCut : public ContinuousProcess<TimeCut> {
+
+ public:
+ TimeCut(TimeType const time);
+
+ template <typename Particle>
+ ProcessReturn doContinuous(
+ Step<Particle> const& step,
+ const bool limitFlag = false); // this is not used for TimeCut
+ template <typename Particle, typename Track>
+ LengthType getMaxStepLength(Particle const&, Track const&) {
+ return meter * std::numeric_limits<double>::infinity();
+ }
+
+ private:
+ TimeType time_;
+ };
+
+} // namespace corsika
+
+#include <corsika/detail/modules/TimeCut.inl>
diff --git a/corsika/modules/radio/CoREAS.hpp b/corsika/modules/radio/CoREAS.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..0756a78af2b932093312fd3416780929c319379b
--- /dev/null
+++ b/corsika/modules/radio/CoREAS.hpp
@@ -0,0 +1,59 @@
+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+#pragma once
+
+#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>
+#include <corsika/modules/radio/propagators/SignalPath.hpp>
+#include <cmath>
+
+namespace corsika {
+
+ template <typename TRadioDetector, typename TPropagator>
+ class CoREAS final
+ : public RadioProcess<TRadioDetector, CoREAS<TRadioDetector, TPropagator>,
+ TPropagator> {
+
+ public:
+ // an identifier for which algorithm was used
+ static constexpr auto algorithm = "CoREAS";
+
+ /**
+ * 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.
+ *
+ * This must be provided by the TRadioImpl.
+ *
+ * @param particle The current particle.
+ * @param track The current track.
+ *
+ */
+ template <typename Particle>
+ ProcessReturn simulate(Step<Particle> const& step);
+
+ using Base =
+ RadioProcess<TRadioDetector, CoREAS<TRadioDetector, TPropagator>, TPropagator>;
+ using Base::antennas_;
+
+ }; // end of class CoREAS
+
+} // namespace corsika
+
+#include <corsika/detail/modules/radio/CoREAS.inl>
\ No newline at end of file
diff --git a/corsika/modules/radio/RadioProcess.hpp b/corsika/modules/radio/RadioProcess.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..9ae9087bde8e81b06c143f9fff320cacd7fe279f
--- /dev/null
+++ b/corsika/modules/radio/RadioProcess.hpp
@@ -0,0 +1,111 @@
+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+#pragma once
+
+#include <istream>
+#include <fstream>
+#include <iostream>
+#include <string>
+#include <corsika/output/BaseOutput.hpp>
+#include <corsika/framework/process/ContinuousProcess.hpp>
+#include <corsika/framework/core/Step.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.
+ * TAntennaCollection is the detector instance that stores antennas
+ * and is responsible for managing the output writing.
+ */
+ template <typename TAntennaCollection, typename TRadioImpl, typename TPropagator>
+ class RadioProcess : public ContinuousProcess<
+ RadioProcess<TAntennaCollection, TRadioImpl, TPropagator>>,
+ public BaseOutput {
+
+ /*
+ * 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.
+ */
+ TRadioImpl& implementation();
+
+ /**
+ * Get a const reference to the underlying implementation.
+ */
+ TRadioImpl const& implementation() const;
+
+ protected:
+ 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(TAntennaCollection& antennas, TArgs&&... args);
+
+ /**
+ * Perform the continuous process (radio emission).
+ *
+ * This handles filtering individual particle tracks
+ * before passing them to `Simulate`.`
+ *
+ * @param particle The current particle.
+ * @param track The current track.
+ */
+ template <typename Particle>
+ ProcessReturn doContinuous(Step<Particle> const& step, bool const);
+
+ /**
+ * 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.
+ */
+ template <typename Particle, typename Track>
+ LengthType getMaxStepLength(Particle const& vParticle, Track const& vTrack) const;
+
+ /**
+ * Called at the start of each library.
+ */
+ void startOfLibrary(boost::filesystem::path const& directory) final override;
+
+ /**
+ * Called at the end of each shower.
+ */
+ virtual void endOfShower(unsigned int const) final override;
+
+ /**
+ * Called at the end of each library.
+ *
+ */
+ void endOfLibrary() final override {}
+
+ /**
+ * Get the configuration of this output.
+ */
+ YAML::Node getConfig() const final;
+
+ }; // END: class RadioProcess
+
+} // namespace corsika
+
+#include <corsika/detail/modules/radio/RadioProcess.inl>
\ No newline at end of file
diff --git a/corsika/modules/radio/ZHS.hpp b/corsika/modules/radio/ZHS.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..5ec38895b8a8c75a9cb6974b996253a430de1f3e
--- /dev/null
+++ b/corsika/modules/radio/ZHS.hpp
@@ -0,0 +1,62 @@
+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+#pragma once
+
+#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>
+#include <corsika/modules/radio/propagators/SignalPath.hpp>
+#include <cmath>
+
+namespace corsika {
+
+ /**
+ * A concrete implementation of the ZHS algorithm.
+ */
+ template <typename TRadioDetector, typename TPropagator>
+ class ZHS final : public RadioProcess<TRadioDetector, ZHS<TRadioDetector, TPropagator>,
+ TPropagator> {
+
+ public:
+ // an identifier for which algorithm was used
+ static constexpr auto algorithm = "ZHS";
+
+ /**
+ * Construct a new ZHS instance.
+ *
+ * This forwards the detector and other arguments to
+ * the RadioProcess parent.
+ *
+ */
+ template <typename... TArgs>
+ ZHS(TRadioDetector& detector, TArgs&&... args)
+ : RadioProcess<TRadioDetector, ZHS, TPropagator>(detector, args...){};
+
+ /**
+ * Simulate the radio emission from a particle across a track.
+ *
+ * This must be provided by the TRadioImpl.
+ *
+ * @param particle The current particle.
+ * @param track The current track.
+ *
+ */
+ template <typename Particle>
+ ProcessReturn simulate(Step<Particle> const& step) const;
+
+ private:
+ using Base =
+ RadioProcess<TRadioDetector, ZHS<TRadioDetector, TPropagator>, TPropagator>;
+ using Base::antennas_;
+
+ }; // END: class ZHS
+
+} // namespace corsika
+
+#include <corsika/detail/modules/radio/ZHS.inl>
\ No newline at end of file
diff --git a/corsika/modules/radio/antennas/Antenna.hpp b/corsika/modules/radio/antennas/Antenna.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..b5aeb8fb99be1dd69e0ce57536a2e0e731921d09
--- /dev/null
+++ b/corsika/modules/radio/antennas/Antenna.hpp
@@ -0,0 +1,127 @@
+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+#pragma once
+
+#include <cnpy.hpp>
+#include <boost/filesystem.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/core/PhysicalGeometry.hpp>
+
+namespace corsika {
+
+ /**
+ * A common abstract interface for radio antennas.
+ *
+ * All concrete antenna implementations should be of
+ * type Antenna<T> where T is a concrete antenna implementation.
+ *
+ */
+ template <typename TAntennaImpl>
+ class Antenna {
+
+ protected:
+ std::string const name_; ///< The name/identifier of this antenna.
+ Point const location_; ///< The location of this antenna.
+ CoordinateSystemPtr const coordinateSystem_; ///< The coordinate system of the antenna
+ std::string filename_ = ""; ///< The filename for the output file for this antenna.
+
+ public:
+ using axistype = std::vector<long double>;
+
+ /**
+ * \brief Construct a base antenna instance.
+ *
+ * @param name A name for this antenna.
+ * @param location The location of this antenna.
+ *
+ */
+ Antenna(std::string const& name, Point const& location,
+ CoordinateSystemPtr const& coordinateSystem);
+
+ /**
+ * Receive a signal at this antenna.
+ *
+ * This is a general implementation call that must be specialized
+ * for the particular antenna implementation and usage.
+ *
+ */
+ template <typename... TVArgs>
+ void receive(TVArgs&&... args);
+
+ /**
+ * Get the location of this antenna.
+ */
+ Point const& getLocation() const;
+
+ /**
+ * Get the name of this name antenna.
+ *
+ * This is used in producing the output data file.
+ */
+ std::string const& getName() const;
+
+ /**
+ * Reset the antenna before starting a new simulation.
+ */
+ 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.
+ */
+ axistype getAxis() const;
+
+ /**
+ * Return a reference to the underlying waveform data for X polarization.
+ *
+ * This is used when writing the antenna information to disk
+ * and will be converted to a 32-bit float before writing.
+ */
+ std::vector<double> const& getWaveformX() const;
+
+ /**
+ * Return a reference to the underlying waveform data for Y polarization.
+ *
+ * This is used when writing the antenna information to disk
+ * and will be converted to a 32-bit float before writing.
+ */
+ std::vector<double> const& getWaveformY() const;
+
+ /**
+ * Return a reference to the underlying waveform data for Z polarization.
+ *
+ * This is used when writing the antenna information to disk
+ * and will be converted to a 32-bit float before writing.
+ */
+ std::vector<double> const& getWaveformZ() const;
+
+ /**
+ * Prepare for the start of the library.
+ */
+ void startOfLibrary(boost::filesystem::path const& directory,
+ std::string const radioImplementation);
+
+ /**
+ * Flush the data from this shower to disk.
+ */
+ void endOfShower(int const event, std::string const& radioImplementation,
+ double const sampleRate);
+
+ protected:
+ /**
+ * Get a reference to the underlying radio implementation.
+ */
+ TAntennaImpl& implementation();
+
+ }; // END: class Antenna final
+
+} // namespace corsika
+
+#include <corsika/detail/modules/radio/antennas/Antenna.inl>
diff --git a/corsika/modules/radio/antennas/TimeDomainAntenna.hpp b/corsika/modules/radio/antennas/TimeDomainAntenna.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..801c529b87a035a021ab1c3aeb52cd86adf9e14d
--- /dev/null
+++ b/corsika/modules/radio/antennas/TimeDomainAntenna.hpp
@@ -0,0 +1,138 @@
+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+#pragma once
+
+#include <corsika/modules/radio/antennas/Antenna.hpp>
+#include <vector>
+
+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> {
+
+ 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 const num_bins_; ///< The number of bins used.
+ std::vector<double> waveformEX_; ///< EX polarization.
+ std::vector<double> waveformEY_; ///< EY polarization.
+ std::vector<double> waveformEZ_; ///< EZ polarization.
+ TimeType const ground_hit_time_; ///< The time the primary particle hits the ground.
+ std::vector<long double> const
+ time_axis_; ///< The time axis corresponding to the electric field.
+
+ public:
+ // import the methods from the antenna
+
+ // label this as a time-domain antenna.
+ static constexpr bool is_time_domain{true};
+
+ 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 coordinateSystem The coordinate system 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 ground_hit_time The time the primary particle hits the ground on a
+ * straight vertical line.
+ *
+ */
+ TimeDomainAntenna(std::string const& name, Point const& location,
+ CoordinateSystemPtr coordinateSystem, TimeType const& start_time,
+ TimeType const& duration, InverseTimeType const& sample_rate,
+ TimeType const ground_hit_time);
+
+ /**
+ * 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);
+
+ void receive(TimeType const time, Vector<dimensionless_d> const& receive_vector,
+ VectorPotential const& vectorP);
+
+ /**
+ * Return the time-units of each waveform for X polarization
+ *
+ * This returns them in nanoseconds for ease of use.
+ */
+ auto const& getWaveformX() const;
+
+ /**
+ * Return the time-units of each waveform for Y polarization
+ *
+ * This returns them in nanoseconds for ease of use.
+ */
+ auto const& getWaveformY() const;
+
+ /**
+ * Return the time-units of each waveform for Z polarization
+ *
+ * This returns them in nanoseconds for ease of use.
+ */
+ auto const& getWaveformZ() const;
+
+ /**
+ * Creates time-units of each waveform.
+ *
+ * It creates them in nanoseconds for ease of use.
+ */
+ std::vector<long double> createTimeAxis() const;
+
+ /**
+ * Return the time-units of each waveform.
+ *
+ * This returns them in nanoseconds for ease of use.
+ */
+ auto const& getAxis() const;
+
+ /**
+ * Returns the sampling rate of the time domain antenna.
+ */
+ InverseTimeType const& getSampleRate() const;
+
+ /**
+ * Returns the start time of detection for the time domain antenna.
+ */
+ TimeType const& getStartTime() const;
+
+ /**
+ * Reset the antenna before starting a new simulation.
+ */
+ void reset();
+
+ /**
+ * Return a YAML configuration for this antenna.
+ */
+ YAML::Node getConfig() const;
+
+ }; // END: class TimeDomainAntenna
+
+} // namespace corsika
+
+#include <corsika/detail/modules/radio/antennas/TimeDomainAntenna.inl>
diff --git a/corsika/modules/radio/detectors/AntennaCollection.hpp b/corsika/modules/radio/detectors/AntennaCollection.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..368be18d934d2a4093a9871c6e51e4157523c058
--- /dev/null
+++ b/corsika/modules/radio/detectors/AntennaCollection.hpp
@@ -0,0 +1,69 @@
+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+#pragma once
+
+namespace corsika {
+
+ /**
+ * The base interface for radio detectors.
+ * At the moment it is a collection of antennas with the same implementation.
+ */
+
+ template <typename TAntennaImpl>
+ class AntennaCollection {
+
+ /**
+ * The collection of antennas used in this simulation.
+ */
+ std::vector<TAntennaImpl> antennas_;
+
+ public:
+ /**
+ * Add an antenna to this radio process.
+ *
+ * @param antenna The antenna to add
+ */
+ void addAntenna(TAntennaImpl const& antenna);
+
+ /**
+ * Get the specific antenna at that place in the collection
+ *
+ * @param index in the collection
+ */
+ TAntennaImpl& at(std::size_t const i);
+
+ TAntennaImpl const& at(std::size_t const i) const;
+
+ /**
+ * Get the number of antennas in the collection
+ */
+ int size() const;
+
+ /**
+ * Get a *non*-const reference to the collection of antennas.
+ *
+ * @returns An iterable mutable reference to the antennas.
+ */
+ std::vector<TAntennaImpl>& getAntennas();
+
+ /**
+ * Get a const reference to the collection of antennas.
+ *
+ * @returns An iterable mutable reference to the antennas.
+ */
+ std::vector<TAntennaImpl> const& getAntennas() const;
+
+ /**
+ * Reset all the antenna waveforms.
+ */
+ void reset();
+ }; // END: class RadioDetector
+
+} // namespace corsika
+
+#include <corsika/detail/modules/radio/detectors/AntennaCollection.inl>
diff --git a/corsika/modules/radio/propagators/RadioPropagator.hpp b/corsika/modules/radio/propagators/RadioPropagator.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..599ada07f8029beb776aeeab9ada233c9cb22f26
--- /dev/null
+++ b/corsika/modules/radio/propagators/RadioPropagator.hpp
@@ -0,0 +1,47 @@
+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+#pragma once
+
+#include <vector>
+
+#include <corsika/framework/geometry/Path.hpp>
+#include <corsika/modules/radio/propagators/SignalPath.hpp>
+
+namespace corsika {
+
+ /**
+ * Radio propagators are used to calculate the propagation
+ * paths from particles to antennas. Any class that wants
+ * to be used as a RadioPropagator must implement the
+ * following methods:
+ *
+ * SignalPathCollection Propagate(Point const& start,
+ * Point const& end,
+ * LengthType const stepsize);
+ */
+ template <typename TImpl, typename TEnvironment>
+ class RadioPropagator {
+
+ protected:
+ // Since we typically know roughly how many paths will
+ // be computed for a given propagator, we use a std::vector here.
+ using SignalPathCollection = std::vector<SignalPath> const;
+
+ TEnvironment const& env_; ///< The environment.
+
+ public:
+ /**
+ * Construct a new RadioPropagator instance.
+ */
+ RadioPropagator(TEnvironment const& env);
+
+ }; // class RadioPropagator
+
+} // namespace corsika
+
+#include <corsika/detail/modules/radio/propagators/RadioPropagator.inl>
\ No newline at end of file
diff --git a/corsika/modules/radio/propagators/SignalPath.hpp b/corsika/modules/radio/propagators/SignalPath.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..9a1239223838b5c56dd49384d9152bb42ee9032d
--- /dev/null
+++ b/corsika/modules/radio/propagators/SignalPath.hpp
@@ -0,0 +1,51 @@
+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/geometry/Path.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+
+namespace corsika {
+
+ /**
+ * Store the photon signal path between two points.
+ *
+ * This is basically a container class
+ */
+ struct SignalPath final : private Path {
+
+ // TODO: discuss if we need average refractivity or average refractive index
+ TimeType const propagation_time_; ///< The total propagation time.
+ double const average_refractive_index_; ///< The average refractive index.
+ double const refractive_index_source_; ///< The refractive index at the source.
+ double const
+ refractive_index_destination_; ///< The refractive index at the destination point.
+ Vector<dimensionless_d> const emit_; ///< The (unit-length) emission vector.
+ Vector<dimensionless_d> const receive_; ///< The (unit-length) receive vector.
+ std::deque<Point> const
+ points_; ///< A collection of points that make up the geometrical path.
+ LengthType const
+ R_distance_; ///< The distance from the point of emission to an observer. TODO:
+ ///< optical path, not geometrical! (probably)
+
+ /**
+ * Create a new SignalPath instance.
+ */
+ SignalPath(TimeType const propagation_time, double const average_refractive_index,
+ double const refractive_index_source,
+ double const refractive_index_destination,
+ Vector<dimensionless_d> const& emit,
+ Vector<dimensionless_d> const& receive, LengthType const R_distance,
+ std::deque<Point> const& points);
+
+ }; // END: class SignalPath final
+
+} // namespace corsika
+
+#include <corsika/detail/modules/radio/propagators/SignalPath.inl>
\ No newline at end of file
diff --git a/corsika/modules/radio/propagators/SimplePropagator.hpp b/corsika/modules/radio/propagators/SimplePropagator.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..33c402774178560248a7e5949bb69f1794ca0f7f
--- /dev/null
+++ b/corsika/modules/radio/propagators/SimplePropagator.hpp
@@ -0,0 +1,51 @@
+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+#pragma once
+
+#include <corsika/media/Environment.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+#include <corsika/framework/core/PhysicalConstants.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/modules/radio/propagators/RadioPropagator.hpp>
+
+namespace corsika {
+
+ /**
+ * This class implements a simple propagator that uses
+ * the straight-line (vector) between the particle
+ * location and the antenna as the trajectory.
+ *
+ */
+ template <typename TEnvironment>
+ class SimplePropagator final
+ : public RadioPropagator<SimplePropagator<TEnvironment>, TEnvironment> {
+
+ using Base = RadioPropagator<SimplePropagator<TEnvironment>, TEnvironment>;
+ using SignalPathCollection = typename Base::SignalPathCollection;
+
+ public:
+ /**
+ * Construct a new SimplePropagator with a given environment.
+ *
+ */
+ SimplePropagator(TEnvironment const& env);
+
+ /**
+ * Return the collection of paths from `source` to `destination`.
+ * Hence, the signal propagated from the
+ * emission point to the antenna location.
+ *
+ */
+ SignalPathCollection propagate(Point const& source, Point const& destination,
+ LengthType const stepsize) const;
+ }; // End: SimplePropagator
+
+} // namespace corsika
+
+#include <corsika/detail/modules/radio/propagators/SimplePropagator.inl>
\ No newline at end of file
diff --git a/corsika/modules/radio/propagators/StraightPropagator.hpp b/corsika/modules/radio/propagators/StraightPropagator.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..d6c156336efadd83e94434360e14caa8ed72bb07
--- /dev/null
+++ b/corsika/modules/radio/propagators/StraightPropagator.hpp
@@ -0,0 +1,52 @@
+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+#pragma once
+
+#include <corsika/media/Environment.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+#include <corsika/framework/core/PhysicalConstants.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/modules/radio/propagators/RadioPropagator.hpp>
+
+namespace corsika {
+
+ /**
+ * This class implements a basic propagator that uses
+ * the straight-line (vector) between the particle
+ * location and the antenna as the trajectory.
+ *
+ * This is what is used in ZHAireS and CoREAS in C7.
+ */
+ template <typename TEnvironment>
+ class StraightPropagator final
+ : public RadioPropagator<StraightPropagator<TEnvironment>, TEnvironment> {
+
+ using Base = RadioPropagator<StraightPropagator<TEnvironment>, TEnvironment>;
+ using SignalPathCollection = typename Base::SignalPathCollection;
+
+ public:
+ /**
+ * Construct a new StraightPropagator with a given environment.
+ *
+ */
+ StraightPropagator(TEnvironment const& env);
+
+ /**
+ * Return the collection of paths from `start` to `end`.
+ * or from 'source' which is the emission point to 'destination'
+ * which is the location of the antenna
+ */
+ SignalPathCollection propagate(Point const& source, Point const& destination,
+ LengthType const stepsize) const;
+
+ }; // End: StraightPropagator
+
+} // namespace corsika
+
+#include <corsika/detail/modules/radio/propagators/StraightPropagator.inl>
\ No newline at end of file
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index 7b89a91772b63dfc731b7ba24d33536158d356f9..7e1df4b56ae2ab93962396f77d941166c619c43b 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -68,4 +68,20 @@ add_executable (water water.cpp)
target_link_libraries (water CORSIKA8::CORSIKA8)
add_executable (environment environment.cpp)
+
target_link_libraries (environment CORSIKA8::CORSIKA8)
+add_executable (synchrotron_test_manual_tracking synchrotron_test_manual_tracking.cpp)
+target_link_libraries (synchrotron_test_manual_tracking CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (synchrotron_test_manual_tracking)
+
+add_executable (synchrotron_test_C8tracking synchrotron_test_C8tracking.cpp)
+target_link_libraries (synchrotron_test_C8tracking CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (synchrotron_test_C8tracking)
+
+add_executable (clover_leaf clover_leaf.cpp)
+target_link_libraries (clover_leaf CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (clover_leaf)
+
+add_executable (radio_em_shower radio_em_shower.cpp)
+target_link_libraries (radio_em_shower CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (radio_em_shower RUN_OPTIONS 10 1121673489)
diff --git a/examples/boundary_example.cpp b/examples/boundary_example.cpp
index 4ede912094f627f1e590996a183ba5ce4b98e5db..b9037bf019255e7d4232f3cbb4c256811a30766d 100644
--- a/examples/boundary_example.cpp
+++ b/examples/boundary_example.cpp
@@ -80,7 +80,7 @@ private:
//
int main() {
- logging::set_level(logging::level::info);
+ logging::set_level(logging::level::warn);
CORSIKA_LOG_INFO("boundary_example");
diff --git a/examples/cascade_example.cpp b/examples/cascade_example.cpp
index 8f846ac8130a3336ef6ec242d8398e7f72f7c214..c27511c18ef03c561838ee45fe21f13763ed002a 100644
--- a/examples/cascade_example.cpp
+++ b/examples/cascade_example.cpp
@@ -57,7 +57,7 @@ using namespace std;
//
int main() {
- logging::set_level(logging::level::info);
+ logging::set_level(logging::level::warn);
CORSIKA_LOG_INFO("cascade_example");
diff --git a/examples/cascade_proton_example.cpp b/examples/cascade_proton_example.cpp
index 62ce71079ddbeebf50e39115cfa7729ac42f9bb0..efb99770ad9e22dea7b94c71ed735a0ad1c680b4 100644
--- a/examples/cascade_proton_example.cpp
+++ b/examples/cascade_proton_example.cpp
@@ -60,7 +60,7 @@ using namespace std;
//
int main() {
- logging::set_level(logging::level::info);
+ logging::set_level(logging::level::warn);
std::cout << "cascade_proton_example" << std::endl;
diff --git a/examples/clover_leaf.cpp b/examples/clover_leaf.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ccb8c7e79aed238d9bfeabe5994dac0efe41cbff
--- /dev/null
+++ b/examples/clover_leaf.cpp
@@ -0,0 +1,228 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <corsika/framework/core/Cascade.hpp>
+#include <corsika/framework/geometry/Sphere.hpp>
+#include <corsika/framework/geometry/PhysicalGeometry.hpp>
+#include <corsika/framework/process/ProcessSequence.hpp>
+#include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/utility/CorsikaFenv.hpp>
+#include <corsika/framework/core/Logging.hpp>
+
+#include <corsika/output/OutputManager.hpp>
+
+#include <corsika/media/Environment.hpp>
+#include <corsika/media/HomogeneousMedium.hpp>
+#include <corsika/media/IMagneticFieldModel.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+#include <corsika/media/MediumPropertyModel.hpp>
+#include <corsika/media/UniformMagneticField.hpp>
+#include <corsika/media/UniformRefractiveIndex.hpp>
+#include <corsika/media/ShowerAxis.hpp>
+
+#include <corsika/setup/SetupStack.hpp>
+#include <corsika/setup/SetupTrajectory.hpp>
+
+#include <corsika/modules/radio/RadioProcess.hpp>
+#include <corsika/modules/radio/CoREAS.hpp>
+#include <corsika/modules/radio/ZHS.hpp>
+#include <corsika/modules/radio/antennas/Antenna.hpp>
+#include <corsika/modules/radio/antennas/TimeDomainAntenna.hpp>
+#include <corsika/modules/radio/detectors/AntennaCollection.hpp>
+#include <corsika/modules/radio/propagators/StraightPropagator.hpp>
+#include <corsika/modules/radio/propagators/SimplePropagator.hpp>
+#include <corsika/modules/radio/propagators/SignalPath.hpp>
+#include <corsika/modules/radio/propagators/RadioPropagator.hpp>
+#include <corsika/modules/TrackWriter.hpp>
+
+#include <corsika/modules/StackInspector.hpp>
+#include <corsika/modules/ParticleCut.hpp>
+#include <corsika/modules/TimeCut.hpp>
+//#include <corsika/modules/TrackWriter.hpp>
+
+/*
+ NOTE, WARNING, ATTENTION
+
+ The .../Random.hpppp implement the hooks of external modules to the C8 random
+ number generator. It has to occur excatly ONCE per linked
+ executable. If you include the header below multiple times and
+ link this together, it will fail.
+ */
+#include <corsika/modules/sibyll/Random.hpp>
+#include <corsika/modules/urqmd/Random.hpp>
+
+#include <iomanip>
+#include <iostream>
+#include <limits>
+#include <string>
+#include <typeinfo>
+
+using namespace corsika;
+using namespace std;
+
+//
+// A simple shower to get the electric field trace of an electron (& a positron)
+// in order to reveal the "clover leaf" pattern of energy fluence to the ground.
+//
+int main() {
+
+ logging::set_level(logging::level::warn);
+ corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
+
+ CORSIKA_LOG_INFO("Synchrotron radiation");
+
+ feenableexcept(FE_INVALID);
+ RNGManager<>::getInstance().registerRandomStream("cascade");
+ std::random_device rd;
+ auto seed = rd();
+ RNGManager<>::getInstance().setSeed(seed);
+
+ OutputManager output("1 electron - 1 positron");
+
+ // create a suitable environment
+ using IModelInterface =
+ IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
+ using AtmModel = UniformRefractiveIndex<
+ MediumPropertyModel<UniformMagneticField<HomogeneousMedium<IModelInterface>>>>;
+ using EnvType = Environment<AtmModel>;
+ EnvType env;
+ CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
+ Point const center{rootCS, 0_m, 0_m, 0_m};
+ // a refractive index for the vacuum
+ const double ri_{1.};
+ // the composition we use for the homogeneous medium
+ NuclearComposition const Composition({Code::Nitrogen}, {1.});
+ // create magnetic field vector
+ auto const Bmag{0.00005_T};
+ MagneticFieldVector B(rootCS, 0_T, Bmag, 0_T);
+ // create a Sphere for the medium
+ auto world = EnvType::createNode<Sphere>(center, 150_km);
+ // set the environment properties
+ world->setModelProperties<AtmModel>(ri_, Medium::AirDry1Atm, B,
+ 1_kg / (1_m * 1_m * 1_m), Composition);
+ // bind things together
+ env.getUniverse()->addChild(std::move(world));
+
+ Point injectionPos(rootCS, 0_m, 0_m, 5_km);
+
+ auto const centerX{center.getCoordinates().getX()};
+ auto const centerY{center.getCoordinates().getY()};
+ auto const centerZ{center.getCoordinates().getZ()};
+ auto const injectionPosX_{injectionPos.getCoordinates().getX()};
+ auto const injectionPosY_{injectionPos.getCoordinates().getY()};
+ auto const injectionPosZ_{injectionPos.getCoordinates().getZ()};
+ auto const triggerpoint_{Point(rootCS, injectionPosX_, injectionPosY_, injectionPosZ_)};
+ std::cout << "Trigger Point is: " << triggerpoint_ << std::endl;
+
+ // the antenna characteristics
+ const TimeType duration_{2e-6_s}; // 0.8e-4_s
+ const InverseTimeType sampleRate_{1e+11_Hz}; // 1e+9_Hz
+
+ // the detectors
+ AntennaCollection<TimeDomainAntenna> detectorCoREAS;
+ // AntennaCollection<TimeDomainAntenna> detectorZHS;
+
+ std::string name_center = "CoREAS_R=0_m--Phi=0degrees";
+ auto triggertime_center{((triggerpoint_ - center).getNorm() / constants::c) - 500_ns};
+ TimeDomainAntenna antenna_center(name_center, center, rootCS, triggertime_center,
+ duration_, sampleRate_, triggertime_center);
+ detectorCoREAS.addAntenna(antenna_center);
+
+ for (auto radius_1 = 25_m; radius_1 <= 1000_m; radius_1 += 25_m) {
+ for (auto phi_1 = 0; phi_1 <= 315; phi_1 += 45) {
+ auto phiRad_1 = phi_1 / 180. * M_PI;
+ auto rr_1 = static_cast<int>(radius_1 / 1_m);
+ auto const point_1{Point(rootCS, centerX + radius_1 * cos(phiRad_1),
+ centerY + radius_1 * sin(phiRad_1), centerZ)};
+ std::cout << "Antenna point: " << point_1 << std::endl;
+ auto triggertime_1{((triggerpoint_ - point_1).getNorm() / constants::c) - 500_ns};
+ std::string name_1 = "CoREAS_R=" + std::to_string(rr_1) +
+ "_m--Phi=" + std::to_string(phi_1) + "degrees";
+ TimeDomainAntenna antenna_1(name_1, point_1, rootCS, triggertime_1, duration_,
+ sampleRate_, triggertime_1);
+ detectorCoREAS.addAntenna(antenna_1);
+ }
+ }
+
+ // setup particle stack, and add primary particle
+ setup::Stack<EnvType> stack;
+ stack.clear();
+
+ // electron
+ const Code beamCode = Code::Electron;
+ auto const charge = get_charge(beamCode);
+ auto const mass = get_mass(beamCode);
+ auto const gyroradius = 100_m;
+ auto const pLabMag = convert_SI_to_HEP(charge * Bmag * gyroradius);
+ auto const omega_inv =
+ convert_HEP_to_SI<MassType::dimension_type>(mass) / (abs(charge) * Bmag);
+ MomentumVector const plab{rootCS, 0_MeV, 0_MeV, -10_MeV};
+ auto const Elab = sqrt(plab.getSquaredNorm() + static_pow<2>(mass));
+ auto gamma = Elab / mass;
+ TimeType const period = 2 * M_PI * omega_inv * gamma;
+
+ // positron
+ const Code beamCodeP = Code::Positron;
+ auto const chargeP = get_charge(beamCodeP);
+ auto const massP = get_mass(beamCodeP);
+ // auto const gyroradius = 100_m;
+ auto const pLabMagP = convert_SI_to_HEP(chargeP * Bmag * gyroradius);
+ auto const omega_invP =
+ convert_HEP_to_SI<MassType::dimension_type>(massP) / (abs(chargeP) * Bmag);
+ MomentumVector const plabP{rootCS, 0_MeV, 0_MeV, -10_MeV};
+ auto const ElabP = sqrt(plabP.getSquaredNorm() + static_pow<2>(massP));
+ auto gammaP = ElabP / massP;
+ TimeType const periodP = 2 * M_PI * omega_invP * gammaP;
+
+ stack.addParticle(std::make_tuple(beamCode,
+ calculate_kinetic_energy(plab.getNorm(), mass),
+ plab.normalized(), injectionPos, 0_ns));
+ stack.addParticle(std::make_tuple(beamCodeP,
+ calculate_kinetic_energy(plabP.getNorm(), massP),
+ plabP.normalized(), injectionPos, 0_ns));
+
+ // setup relevant processes
+ setup::Tracking tracking;
+
+ // put radio processes here
+ RadioProcess<decltype(detectorCoREAS),
+ CoREAS<decltype(detectorCoREAS), decltype(SimplePropagator(env))>,
+ decltype(SimplePropagator(env))>
+ coreas(detectorCoREAS, env);
+ output.add("CoREAS", coreas);
+
+ // RadioProcess<decltype(detectorZHS), ZHS<decltype(detectorZHS),
+ // decltype(SimplePropagator(env))>, decltype(SimplePropagator(env))>
+ // zhs(detectorZHS, env);
+ // output.add("ZHS", zhs);
+
+ TimeCut cut(period / 4);
+
+ // TrackWriter trackWriter;
+ // output.add("tracks", trackWriter); // register TrackWriter
+
+ // assemble all processes into an ordered process list
+ auto sequence = make_sequence(coreas, cut);
+
+ // define air shower object, run simulation
+ Cascade EAS(env, tracking, sequence, output, stack);
+ output.startOfShower();
+ EAS.run();
+ output.endOfShower();
+
+ CORSIKA_LOG_INFO("|p| electron = {} and E electron = {}", plab.getNorm(), Elab);
+ CORSIKA_LOG_INFO("period: {}", period);
+ CORSIKA_LOG_INFO("gamma: {}", gamma);
+
+ CORSIKA_LOG_INFO("|p| positron = {} and E positron = {}", plabP.getNorm(), ElabP);
+ CORSIKA_LOG_INFO("period: {}", periodP);
+ CORSIKA_LOG_INFO("gamma: {}", gammaP);
+
+ output.endOfLibrary();
+}
diff --git a/examples/em_shower.cpp b/examples/em_shower.cpp
index 1516bac2046e818d33e457ba61442db7ac0331cf..0c8e4c448cdf1d878e9a300bf5d8e0475af2ba56 100644
--- a/examples/em_shower.cpp
+++ b/examples/em_shower.cpp
@@ -81,7 +81,7 @@ using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;
int main(int argc, char** argv) {
- logging::set_level(logging::level::info);
+ logging::set_level(logging::level::warn);
if (!(argc == 2 || argc == 3)) {
std::cerr << "usage: em_shower <energy/GeV> [seed]" << std::endl
diff --git a/examples/geometry_example.cpp b/examples/geometry_example.cpp
index 596e6c1842011d56d8db27dc6e02fe81ab8bdbb9..059f6d849c5d6758b6fd5ffbdb99ddfdd71bd0e0 100644
--- a/examples/geometry_example.cpp
+++ b/examples/geometry_example.cpp
@@ -20,7 +20,7 @@ using namespace corsika;
int main() {
- logging::set_level(logging::level::info);
+ logging::set_level(logging::level::warn);
CORSIKA_LOG_INFO("geometry_example");
diff --git a/examples/helix_example.cpp b/examples/helix_example.cpp
index e980931b74ace074bf6d0459cf69e6d158f32b91..9ce484b67636bb68a693e5451e969ace248983df 100644
--- a/examples/helix_example.cpp
+++ b/examples/helix_example.cpp
@@ -20,7 +20,7 @@ using namespace corsika;
int main() {
- logging::set_level(logging::level::info);
+ logging::set_level(logging::level::warn);
CORSIKA_LOG_INFO("helix_example");
diff --git a/examples/particle_list_example.cpp b/examples/particle_list_example.cpp
index 241fda60d53fbeb8d5d9e7862c3a7df51b7a1235..36020bc1d6ac4fc9f9a10e2496441c15312585e2 100644
--- a/examples/particle_list_example.cpp
+++ b/examples/particle_list_example.cpp
@@ -32,7 +32,7 @@ using namespace std;
//
int main() {
- logging::set_level(logging::level::info);
+ logging::set_level(logging::level::warn);
corsika_logger->set_pattern("[%n:%^%-8l%$] %v");
logging::info(
diff --git a/examples/radio_em_shower.cpp b/examples/radio_em_shower.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f151b1e893d0444c47b54365f6145852f53dc0ff
--- /dev/null
+++ b/examples/radio_em_shower.cpp
@@ -0,0 +1,305 @@
+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <corsika/framework/process/ProcessSequence.hpp>
+#include <corsika/framework/process/SwitchProcessSequence.hpp>
+#include <corsika/framework/process/InteractionCounter.hpp>
+#include <corsika/framework/core/Cascade.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
+#include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/geometry/Sphere.hpp>
+#include <corsika/framework/geometry/Plane.hpp>
+#include <corsika/framework/geometry/Sphere.hpp>
+#include <corsika/framework/geometry/PhysicalGeometry.hpp>
+#include <corsika/framework/utility/CorsikaFenv.hpp>
+#include <corsika/framework/utility/SaveBoostHistogram.hpp>
+
+#include <corsika/output/OutputManager.hpp>
+#include <corsika/modules/writers/SubWriter.hpp>
+#include <corsika/modules/writers/EnergyLossWriter.hpp>
+#include <corsika/modules/writers/LongitudinalWriter.hpp>
+
+#include <corsika/media/Environment.hpp>
+#include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+#include <corsika/media/ShowerAxis.hpp>
+#include <corsika/media/MediumPropertyModel.hpp>
+#include <corsika/media/UniformMagneticField.hpp>
+#include <corsika/media/UniformRefractiveIndex.hpp>
+#include <corsika/media/CORSIKA7Atmospheres.hpp>
+
+#include <corsika/modules/LongitudinalProfile.hpp>
+#include <corsika/modules/ObservationPlane.hpp>
+#include <corsika/modules/ParticleCut.hpp>
+#include <corsika/modules/TrackWriter.hpp>
+#include <corsika/modules/Sibyll.hpp>
+#include <corsika/modules/PROPOSAL.hpp>
+
+#include <corsika/modules/radio/RadioProcess.hpp>
+#include <corsika/modules/radio/CoREAS.hpp>
+#include <corsika/modules/radio/ZHS.hpp>
+#include <corsika/modules/radio/antennas/Antenna.hpp>
+#include <corsika/modules/radio/antennas/TimeDomainAntenna.hpp>
+#include <corsika/modules/radio/detectors/AntennaCollection.hpp>
+#include <corsika/modules/radio/propagators/StraightPropagator.hpp>
+#include <corsika/modules/radio/propagators/SimplePropagator.hpp>
+#include <corsika/modules/radio/propagators/SignalPath.hpp>
+#include <corsika/modules/radio/propagators/RadioPropagator.hpp>
+
+#include <corsika/setup/SetupStack.hpp>
+#include <corsika/setup/SetupTrajectory.hpp>
+
+#include <iomanip>
+#include <iostream>
+#include <limits>
+#include <string>
+#include <typeinfo>
+
+/*
+ NOTE, WARNING, ATTENTION
+
+ The .../Random.hpppp implement the hooks of external modules to the C8 random
+ number generator. It has to occur exactly ONCE per linked
+ executable. If you include the header below multiple times and
+ link this together, it will fail.
+*/
+#include <corsika/modules/Random.hpp>
+
+using namespace corsika;
+using namespace std;
+
+void registerRandomStreams(int seed) {
+ RNGManager<>::getInstance().registerRandomStream("cascade");
+ RNGManager<>::getInstance().registerRandomStream("proposal");
+ RNGManager<>::getInstance().registerRandomStream("sibyll");
+ if (seed == 0) {
+ std::random_device rd;
+ seed = rd();
+ cout << "new random seed (auto) " << seed << endl;
+ }
+ RNGManager<>::getInstance().setSeed(seed);
+}
+
+template <typename TInterface>
+using MyExtraEnv =
+ UniformRefractiveIndex<MediumPropertyModel<UniformMagneticField<TInterface>>>;
+
+int main(int argc, char** argv) {
+
+ logging::set_level(logging::level::warn);
+
+ if (argc != 3) {
+ std::cerr
+ << "usage: radio_em_shower <energy/GeV> <seed> - put seed=0 to use random seed"
+ << std::endl;
+ return 1;
+ }
+
+ int seed{static_cast<int>(std::stof(std::string(argv[2])))};
+ std::cout << "Seed: " << seed << std::endl;
+ feenableexcept(FE_INVALID);
+ // initialize random number sequence(s)
+ registerRandomStreams(seed);
+
+ // setup environment, geometry
+ using EnvironmentInterface =
+ IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
+ using EnvType = Environment<EnvironmentInterface>;
+ EnvType env;
+ CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
+ Point const center{rootCS, 0_m, 0_m, 0_m};
+
+ create_5layer_atmosphere<EnvironmentInterface, MyExtraEnv>(
+ env, AtmosphereId::LinsleyUSStd, center, 1.000327, Medium::AirDry1Atm,
+ MagneticFieldVector{rootCS, 50_uT, 0_T, 0_T});
+
+ std::unordered_map<Code, HEPEnergyType> energy_resolution = {
+ {Code::Electron, 5_MeV},
+ {Code::Positron, 5_MeV},
+ {Code::Photon, 5_MeV},
+ };
+ for (auto [pcode, energy] : energy_resolution)
+ set_energy_production_threshold(pcode, energy);
+
+ // setup particle stack, and add primary particle
+ setup::Stack<EnvType> stack;
+ stack.clear();
+ const Code beamCode = Code::Electron;
+ auto const mass = get_mass(beamCode);
+ const HEPEnergyType E0 = 1_GeV * std::stof(std::string(argv[1]));
+ double theta = 0.;
+ auto const thetaRad = theta / 180. * M_PI;
+
+ HEPMomentumType P0 = calculate_momentum(E0, mass);
+ auto momentumComponents = [](double thetaRad, HEPMomentumType ptot) {
+ return std::make_tuple(ptot * sin(thetaRad), 0_eV, -ptot * cos(thetaRad));
+ };
+
+ auto const [px, py, pz] = momentumComponents(thetaRad, P0);
+ auto plab = MomentumVector(rootCS, {px, py, pz});
+ cout << "input particle: " << beamCode << endl;
+ cout << "input angles: theta=" << theta << endl;
+ cout << "input momentum: " << plab.getComponents() / 1_GeV
+ << ", norm = " << plab.getNorm() << endl;
+
+ auto const observationHeight = 1.4_km + constants::EarthRadius::Mean;
+ auto const injectionHeight = 112.75_km + constants::EarthRadius::Mean;
+ auto const t = -observationHeight * cos(thetaRad) +
+ sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) +
+ static_pow<2>(injectionHeight));
+ Point const showerCore{rootCS, 0_m, 0_m, observationHeight};
+ Point const injectionPos =
+ showerCore + DirectionVector{rootCS, {-sin(thetaRad), 0, cos(thetaRad)}} * t;
+
+ std::cout << "point of injection: " << injectionPos.getCoordinates() << std::endl;
+
+ stack.addParticle(std::make_tuple(
+ beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
+ plab.normalized(), injectionPos, 0_ns));
+
+ CORSIKA_LOG_INFO("shower axis length: {} ",
+ (showerCore - injectionPos).getNorm() * 1.02);
+
+ ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env,
+ false, 1000};
+
+ TimeType const groundHitTime{(showerCore - injectionPos).getNorm() / constants::c};
+
+ std::string outname_ = "radio_em_shower_outputs"; // + std::to_string(rr_);
+ OutputManager output(outname_);
+
+ // Radio antennas and relevant information
+ // the antenna time variables
+ const TimeType duration_{1e-6_s};
+ const InverseTimeType sampleRate_{1e+9_Hz};
+
+ // the detector (aka antenna collection) for CoREAS and ZHS
+ AntennaCollection<TimeDomainAntenna> detectorCoREAS;
+ AntennaCollection<TimeDomainAntenna> detectorZHS;
+
+ auto const showerCoreX_{showerCore.getCoordinates().getX()};
+ auto const showerCoreY_{showerCore.getCoordinates().getY()};
+ auto const injectionPosX_{injectionPos.getCoordinates().getX()};
+ auto const injectionPosY_{injectionPos.getCoordinates().getY()};
+ auto const injectionPosZ_{injectionPos.getCoordinates().getZ()};
+ auto const triggerpoint_{Point(rootCS, injectionPosX_, injectionPosY_, injectionPosZ_)};
+ std::cout << "Trigger Point is: " << triggerpoint_ << std::endl;
+
+ // // setup CoREAS antennas - use the for loop for star shape pattern
+ // for (auto radius_1 = 25_m; radius_1 <= 500_m; radius_1 += 25_m) {
+ // for (auto phi_1 = 0; phi_1 <= 315; phi_1 += 45) {
+ auto radius_1 = 200_m;
+ auto phi_1 = 45;
+ auto phiRad_1 = phi_1 / 180. * M_PI;
+ auto rr_1 = static_cast<int>(radius_1 / 1_m);
+ auto const point_1{Point(rootCS, showerCoreX_ + radius_1 * cos(phiRad_1),
+ showerCoreY_ + radius_1 * sin(phiRad_1),
+ constants::EarthRadius::Mean)};
+ std::cout << "Antenna point: " << point_1 << std::endl;
+ auto triggertime_1{(triggerpoint_ - point_1).getNorm() / constants::c};
+ std::string name_1 =
+ "CoREAS_R=" + std::to_string(rr_1) + "_m--Phi=" + std::to_string(phi_1) + "degrees";
+ TimeDomainAntenna antenna_1(name_1, point_1, rootCS, triggertime_1, duration_,
+ sampleRate_, triggertime_1);
+ detectorCoREAS.addAntenna(antenna_1);
+ // }
+ // }
+
+ // // setup ZHS antennas - use the for loop for star shape pattern
+ // for (auto radius_ = 25_m; radius_ <= 500_m; radius_ += 25_m) {
+ // for (auto phi_ = 0; phi_ <= 315; phi_ += 45) {
+ auto radius_ = 200_m;
+ auto phi_ = 45;
+ auto phiRad_ = phi_ / 180. * M_PI;
+ auto rr_ = static_cast<int>(radius_ / 1_m);
+ auto const point_{Point(rootCS, showerCoreX_ + radius_ * cos(phiRad_),
+ showerCoreY_ + radius_ * sin(phiRad_),
+ constants::EarthRadius::Mean)};
+ auto triggertime_{(triggerpoint_ - point_).getNorm() / constants::c};
+ std::string name_ =
+ "ZHS_R=" + std::to_string(rr_) + "_m--Phi=" + std::to_string(phi_) + "degrees";
+ TimeDomainAntenna antenna_(name_, point_, rootCS, triggertime_, duration_, sampleRate_,
+ triggertime_);
+ detectorZHS.addAntenna(antenna_);
+ // }
+ // }
+
+ // setup processes, decays and interactions
+
+ EnergyLossWriter dEdX{showerAxis, 10_g / square(1_cm), 200};
+ // register energy losses as output
+ output.add("dEdX", dEdX);
+
+ ParticleCut<SubWriter<decltype(dEdX)>> cut(5_MeV, 5_MeV, 100_GeV, 100_GeV, true, dEdX);
+
+ corsika::sibyll::Interaction sibyll{env};
+ HEPEnergyType heThresholdNN = 80_GeV;
+ corsika::proposal::Interaction emCascade(env, sibyll.getHadronInteractionModel(),
+ heThresholdNN);
+ corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>> emContinuous(env, dEdX);
+ // BetheBlochPDG<SubWriter<decltype(dEdX)>> emContinuous{dEdX};
+
+ // NOT possible right now, due to interface differenc in PROPOSAL
+ // InteractionCounter emCascadeCounted(emCascade);
+
+ TrackWriter tracks;
+ output.add("tracks", tracks);
+
+ // long. profile
+ LongitudinalWriter profile{showerAxis, 10_g / square(1_cm)};
+ output.add("profile", profile);
+ LongitudinalProfile<SubWriter<decltype(profile)>> longprof{profile};
+
+ // initiate CoREAS
+ RadioProcess<decltype(detectorCoREAS),
+ CoREAS<decltype(detectorCoREAS), decltype(SimplePropagator(env))>,
+ decltype(SimplePropagator(env))>
+ coreas(detectorCoREAS, env);
+
+ // register CoREAS with the output manager
+ output.add("CoREAS", coreas);
+
+ // initiate ZHS
+ RadioProcess<decltype(detectorZHS),
+ ZHS<decltype(detectorZHS), decltype(SimplePropagator(env))>,
+ decltype(SimplePropagator(env))>
+ zhs(detectorZHS, env);
+
+ // // register ZHS with the output manager
+ output.add("ZHS", zhs);
+
+ Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
+ ObservationPlane<setup::Tracking, ParticleWriterParquet> observationLevel{
+ obsPlane, DirectionVector(rootCS, {1., 0., 0.})};
+ output.add("particles", observationLevel);
+
+ // auto sequence = make_sequence(emCascade, emContinuous, longprof, cut, coreas, zhs);
+ auto sequence = make_sequence(emCascade, emContinuous, longprof, cut, coreas, zhs,
+ observationLevel, tracks);
+ // define air shower object, run simulation
+ setup::Tracking tracking;
+
+ output.startOfLibrary();
+ Cascade EAS(env, tracking, sequence, output, stack);
+
+ // to fix the point of first interaction, uncomment the following two lines:
+ // EAS.forceInteraction();
+
+ EAS.run();
+
+ HEPEnergyType const Efinal = dEdX.getEnergyLost() + observationLevel.getEnergyGround();
+
+ CORSIKA_LOG_INFO(
+ "total energy budget (GeV): {}, "
+ "relative difference (%): {}",
+ Efinal / 1_GeV, (Efinal / E0 - 1) * 100);
+
+ output.endOfLibrary();
+}
\ No newline at end of file
diff --git a/examples/stack_example.cpp b/examples/stack_example.cpp
index da97201afad3001f6a4b9e5cec6f2701df227091..8f7f6cdeb9a5e24e2547c309b69506836a4fc615 100644
--- a/examples/stack_example.cpp
+++ b/examples/stack_example.cpp
@@ -43,7 +43,7 @@ void read(VectorStack& s) {
int main() {
- logging::set_level(logging::level::info);
+ logging::set_level(logging::level::warn);
CORSIKA_LOG_INFO("stack_example");
VectorStack s;
diff --git a/examples/stopping_power.cpp b/examples/stopping_power.cpp
index 4aa5a0f5a504aa4a6baf5ebb2bffd14f87c01c26..8b39786e1380f7f25c0821635218f220e472094f 100644
--- a/examples/stopping_power.cpp
+++ b/examples/stopping_power.cpp
@@ -30,7 +30,7 @@ using namespace std;
//
int main() {
- logging::set_level(logging::level::info);
+ logging::set_level(logging::level::warn);
CORSIKA_LOG_INFO("stopping_power");
diff --git a/examples/synchrotron_test_C8tracking.cpp b/examples/synchrotron_test_C8tracking.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..cf1792188e4e463c33c2f5a9078b8121bdfb9ef3
--- /dev/null
+++ b/examples/synchrotron_test_C8tracking.cpp
@@ -0,0 +1,172 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <corsika/framework/core/Cascade.hpp>
+#include <corsika/framework/geometry/Sphere.hpp>
+#include <corsika/framework/geometry/PhysicalGeometry.hpp>
+#include <corsika/framework/process/ProcessSequence.hpp>
+#include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/utility/CorsikaFenv.hpp>
+#include <corsika/framework/core/Logging.hpp>
+
+#include <corsika/output/OutputManager.hpp>
+
+#include <corsika/media/Environment.hpp>
+#include <corsika/media/HomogeneousMedium.hpp>
+#include <corsika/media/IMagneticFieldModel.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+#include <corsika/media/MediumPropertyModel.hpp>
+#include <corsika/media/UniformMagneticField.hpp>
+#include <corsika/media/UniformRefractiveIndex.hpp>
+
+#include <corsika/setup/SetupStack.hpp>
+#include <corsika/setup/SetupTrajectory.hpp>
+
+#include <corsika/modules/radio/RadioProcess.hpp>
+#include <corsika/modules/radio/CoREAS.hpp>
+#include <corsika/modules/radio/ZHS.hpp>
+#include <corsika/modules/radio/antennas/TimeDomainAntenna.hpp>
+#include <corsika/modules/radio/detectors/AntennaCollection.hpp>
+#include <corsika/modules/radio/propagators/SimplePropagator.hpp>
+
+#include <corsika/modules/TimeCut.hpp>
+
+/*
+ NOTE, WARNING, ATTENTION
+
+ The .../Random.hpppp implement the hooks of external modules to the C8 random
+ number generator. It has to occur excatly ONCE per linked
+ executable. If you include the header below multiple times and
+ link this together, it will fail.
+ */
+#include <corsika/modules/sibyll/Random.hpp>
+#include <corsika/modules/urqmd/Random.hpp>
+
+#include <iomanip>
+#include <iostream>
+#include <limits>
+#include <string>
+#include <typeinfo>
+
+using namespace corsika;
+using namespace std;
+
+//
+// A simple shower to get the electric field trace of an electron using C8 tracking
+//
+int main() {
+
+ logging::set_level(logging::level::warn);
+ corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
+
+ CORSIKA_LOG_INFO("Synchrotron radiation");
+
+ feenableexcept(FE_INVALID);
+ RNGManager<>::getInstance().registerRandomStream("cascade");
+ std::random_device rd;
+ auto seed = rd();
+ RNGManager<>::getInstance().setSeed(seed);
+
+ OutputManager output("synchrotron_radiation_C8tracking-output");
+
+ // set up the environment
+ using EnvironmentInterface =
+ IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
+ using EnvType = Environment<EnvironmentInterface>;
+ EnvType env;
+ auto& universe = *(env.getUniverse());
+ CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
+
+ auto world = EnvType::createNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 150_km);
+
+ using MyHomogeneousModel = UniformRefractiveIndex<
+ MediumPropertyModel<UniformMagneticField<HomogeneousMedium<EnvironmentInterface>>>>;
+
+ auto const Bmag{0.0003809_T};
+ MagneticFieldVector B{rootCS, 0_T, 0_T, Bmag};
+
+ // the composition we use for the homogeneous medium
+ NuclearComposition const nitrogenComposition({Code::Nitrogen}, {1.});
+
+ world->setModelProperties<MyHomogeneousModel>(
+ 1, Medium::AirDry1Atm, B, 1_kg / (1_m * 1_m * 1_m), nitrogenComposition);
+
+ universe.addChild(std::move(world));
+
+ // the antenna locations
+ const auto point1{Point(rootCS, 30000_m, 0_m, 0_m)};
+
+ // the antenna time variables
+ const TimeType t1{0.994e-4_s};
+ const TimeType t2{1.07e-4_s - 0.994e-4_s};
+ const InverseTimeType t3{5e+11_Hz};
+
+ // the antennas
+ TimeDomainAntenna ant1("antenna CoREAS", point1, rootCS, t1, t2, t3, t1);
+ TimeDomainAntenna ant2("antenna ZHS", point1, rootCS, t1, t2, t3, t1);
+
+ // the detectors
+ AntennaCollection<TimeDomainAntenna> detectorCoREAS;
+ AntennaCollection<TimeDomainAntenna> detectorZHS;
+ detectorCoREAS.addAntenna(ant1);
+ detectorZHS.addAntenna(ant2);
+
+ // setup particle stack, and add primary particle
+ setup::Stack<EnvType> stack;
+ stack.clear();
+ const Code beamCode = Code::Electron;
+ auto const charge = get_charge(beamCode);
+ auto const mass = get_mass(beamCode);
+ auto const gyroradius = 100_m;
+ auto const pLabMag = convert_SI_to_HEP(charge * Bmag * gyroradius);
+ auto const omega_inv =
+ convert_HEP_to_SI<MassType::dimension_type>(mass) / (abs(charge) * Bmag);
+ MomentumVector const plab{rootCS, pLabMag, 0_MeV, 0_MeV};
+ auto const Elab = sqrt(plab.getSquaredNorm() + static_pow<2>(mass));
+ auto gamma = Elab / mass;
+ TimeType const period = 2 * M_PI * omega_inv * gamma;
+
+ Point injectionPos(rootCS, 0_m, 100_m, 0_m);
+ stack.addParticle(std::make_tuple(
+ beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
+ plab.normalized(), injectionPos, 0_ns));
+
+ // setup relevant processes
+ setup::Tracking tracking;
+
+ // put radio processes here
+ RadioProcess<decltype(detectorCoREAS),
+ CoREAS<decltype(detectorCoREAS), decltype(SimplePropagator(env))>,
+ decltype(SimplePropagator(env))>
+ coreas(detectorCoREAS, env);
+ output.add("CoREAS", coreas);
+
+ RadioProcess<decltype(detectorZHS),
+ ZHS<decltype(detectorZHS), decltype(SimplePropagator(env))>,
+ decltype(SimplePropagator(env))>
+ zhs(detectorZHS, env);
+ output.add("ZHS", zhs);
+
+ TimeCut cut(period);
+
+ // assemble all processes into an ordered process list
+ auto sequence = make_sequence(coreas, zhs, cut);
+
+ // define air shower object, run simulation
+ Cascade EAS(env, tracking, sequence, output, stack);
+ output.startOfShower();
+ EAS.run();
+ output.endOfShower();
+
+ CORSIKA_LOG_INFO("|p| = {} and E = {}", plab.getNorm(), Elab);
+ CORSIKA_LOG_INFO("period: {}", period);
+ CORSIKA_LOG_INFO("gamma: {}", gamma);
+
+ output.endOfLibrary();
+}
\ No newline at end of file
diff --git a/examples/synchrotron_test_manual_tracking.cpp b/examples/synchrotron_test_manual_tracking.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8921329f5d365d7f0e54003b7ed51e9417435315
--- /dev/null
+++ b/examples/synchrotron_test_manual_tracking.cpp
@@ -0,0 +1,176 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <corsika/framework/geometry/Sphere.hpp>
+#include <corsika/framework/geometry/PhysicalGeometry.hpp>
+#include <corsika/framework/process/ProcessSequence.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+#include <corsika/output/OutputManager.hpp>
+
+#include <corsika/media/Environment.hpp>
+#include <corsika/media/HomogeneousMedium.hpp>
+#include <corsika/media/IMagneticFieldModel.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+#include <corsika/media/MediumPropertyModel.hpp>
+#include <corsika/media/UniformMagneticField.hpp>
+#include <corsika/media/UniformRefractiveIndex.hpp>
+
+#include <corsika/setup/SetupStack.hpp>
+#include <corsika/setup/SetupTrajectory.hpp>
+
+#include <corsika/modules/radio/RadioProcess.hpp>
+#include <corsika/modules/radio/CoREAS.hpp>
+#include <corsika/modules/radio/ZHS.hpp>
+#include <corsika/modules/radio/antennas/TimeDomainAntenna.hpp>
+#include <corsika/modules/radio/detectors/AntennaCollection.hpp>
+#include <corsika/modules/radio/propagators/SimplePropagator.hpp>
+
+/*
+ NOTE, WARNING, ATTENTION
+
+ The .../Random.hpppp implement the hooks of external modules to the C8 random
+ number generator. It has to occur excatly ONCE per linked
+ executable. If you include the header below multiple times and
+ link this together, it will fail.
+ */
+#include <corsika/modules/sibyll/Random.hpp>
+#include <corsika/modules/urqmd/Random.hpp>
+
+#include <iomanip>
+#include <iostream>
+#include <limits>
+#include <string>
+
+using namespace corsika;
+using namespace std;
+
+//
+// A simple example to get the electric field trace of an electron using manual tracking
+//
+int main() {
+
+ // create a suitable environment
+ using IModelInterface =
+ IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
+ using AtmModel = UniformRefractiveIndex<
+ MediumPropertyModel<UniformMagneticField<HomogeneousMedium<IModelInterface>>>>;
+ using EnvType = Environment<AtmModel>;
+ EnvType env;
+ CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
+ Point const center{rootCS, 0_m, 0_m, 0_m};
+ // a refractive index for the vacuum
+ const double ri_{1};
+ // the constant density
+ const auto density{19.2_g / cube(1_cm)};
+ // the composition we use for the homogeneous medium
+ NuclearComposition const Composition({Code::Nitrogen}, {1.});
+ // create magnetic field vector
+ Vector B1(rootCS, 0_T, 0_T, 0.3809_T);
+ // create a Sphere for the medium
+ auto Medium =
+ EnvType::createNode<Sphere>(center, 1_km * std::numeric_limits<double>::infinity());
+ // set the environment properties
+ auto const props = Medium->setModelProperties<AtmModel>(ri_, Medium::AirDry1Atm, B1,
+ density, Composition);
+ // bind things together
+ env.getUniverse()->addChild(std::move(Medium));
+
+ // the antennas location
+ const auto point1{Point(rootCS, 30000_m, 0_m, 0_m)};
+
+ // create times for the antenna
+ const TimeType start{0.994e-4_s};
+ const TimeType duration{1.07e-4_s - 0.994e-4_s};
+ // 3 km antenna
+ // const TimeType start{0.994e-5_s};
+ // const TimeType duration{1.7e-5_s - 0.994e-5_s};
+ const InverseTimeType sampleRate_{5e+11_Hz};
+
+ std::cout << "number of points in time: " << duration * sampleRate_ << std::endl;
+
+ // create 2 antennas
+ TimeDomainAntenna ant1("CoREAS antenna", point1, rootCS, start, duration, sampleRate_,
+ start);
+ TimeDomainAntenna ant2("ZHS antenna", point1, rootCS, start, duration, sampleRate_,
+ start);
+
+ // construct a radio detector instance to store our antennas
+ AntennaCollection<TimeDomainAntenna> detectorCoREAS;
+ AntennaCollection<TimeDomainAntenna> detectorZHS;
+
+ // add the antennas to the detector
+ detectorCoREAS.addAntenna(ant1);
+ detectorZHS.addAntenna(ant2);
+
+ // create a new stack for each trial
+ setup::Stack<EnvType> stack;
+ stack.clear();
+
+ const Code particle{Code::Electron};
+ const HEPMassType pmass{get_mass(particle)};
+
+ // construct an energy // move in the for loop
+ const HEPEnergyType E0{11.4_MeV};
+
+ // construct the output manager
+ OutputManager outputs("synchrotron_radiation_manual_tracking-output");
+
+ // create a radio process instance using CoREAS
+ RadioProcess<
+ AntennaCollection<TimeDomainAntenna>,
+ CoREAS<AntennaCollection<TimeDomainAntenna>, decltype(SimplePropagator(env))>,
+ decltype(SimplePropagator(env))>
+ coreas(detectorCoREAS, env);
+ // register CoREAS to the output manager
+ outputs.add("CoREAS", coreas);
+
+ // create a radio process instance using ZHS
+ RadioProcess<AntennaCollection<TimeDomainAntenna>,
+ ZHS<AntennaCollection<TimeDomainAntenna>, decltype(SimplePropagator(env))>,
+ decltype(SimplePropagator(env))>
+ zhs(detectorZHS, env);
+ // register ZHS to the output manager
+ outputs.add("ZHS", zhs);
+
+ // trigger the start of the library and the first event
+ outputs.startOfLibrary();
+ outputs.startOfShower();
+
+ // the number of points that make up the circle
+ int const n_points{100000};
+ LengthType const radius{100_m};
+ TimeType timeCounter{0._s};
+
+ // loop over all the tracks twice (this produces 2 pulses)
+ for (size_t i = 0; i <= (n_points)*2; i++) {
+ Point const point_1(rootCS, {radius * cos(M_PI * 2 * i / n_points),
+ radius * sin(M_PI * 2 * i / n_points), 0_m});
+ Point const point_2(rootCS, {radius * cos(M_PI * 2 * (i + 1) / n_points),
+ radius * sin(M_PI * 2 * (i + 1) / n_points), 0_m});
+ TimeType t{(point_2 - point_1).getNorm() / (0.999 * constants::c)};
+ timeCounter = timeCounter + t;
+ VelocityVector v{(point_2 - point_1) / t};
+ auto beta{v / constants::c};
+ auto gamma{E0 / pmass};
+ auto plab{beta * pmass * gamma};
+ Line l{point_1, v};
+ StraightTrajectory track{l, t};
+ auto particle1{stack.addParticle(std::make_tuple(
+ particle, calculate_kinetic_energy(plab.getNorm(), get_mass(particle)),
+ plab.normalized(), point_1, timeCounter))};
+ Step step(particle1, track);
+ coreas.doContinuous(step, true);
+ zhs.doContinuous(step, true);
+ stack.clear();
+ }
+
+ // trigger the manager to write the data to disk
+ outputs.endOfShower();
+ outputs.endOfLibrary();
+}
\ No newline at end of file
diff --git a/examples/vertical_EAS.cpp b/examples/vertical_EAS.cpp
index 1a24ab355ae5180cebf04fc709f37e243e089e4b..e24a632fadef6f9f399d6e3a7604909e50a4156f 100644
--- a/examples/vertical_EAS.cpp
+++ b/examples/vertical_EAS.cpp
@@ -105,7 +105,7 @@ using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;
int main(int argc, char** argv) {
- logging::set_level(logging::level::info);
+ logging::set_level(logging::level::warn);
CORSIKA_LOG_INFO("vertical_EAS");
diff --git a/python/corsika/io/library.py b/python/corsika/io/library.py
index 236fc97ac972a11c0fb366876bd92cae20d6ab18..1c2ff38e025e9a728bde4b9c7a5e0cad5982deb1 100644
--- a/python/corsika/io/library.py
+++ b/python/corsika/io/library.py
@@ -10,8 +10,7 @@
import logging
import os
import os.path as op
-import re
-from typing import Any, Dict, Optional, List
+from typing import Any, Dict, List, Optional
import yaml
diff --git a/python/corsika/io/outputs/__init__.py b/python/corsika/io/outputs/__init__.py
index 63b76b533d41085aa77af57df5439c5dd3c68335..7678ab32b2fa79be594ca0b34e2846a166c933c4 100644
--- a/python/corsika/io/outputs/__init__.py
+++ b/python/corsika/io/outputs/__init__.py
@@ -14,6 +14,7 @@ from .bethe_bloch import BetheBlochPDG
from .particle_cut import ParticleCut
from .energy_loss import EnergyLoss
from .output import Output
+from .radio_process import RadioProcess
__all__ = [
"Output",
@@ -23,4 +24,5 @@ __all__ = [
"BetheBlochPDG",
"ParticleCut",
"EnergyLoss"
+ "RadioProcess",
]
diff --git a/python/corsika/io/outputs/output.py b/python/corsika/io/outputs/output.py
index fed4368aac5d1cb8be7c1cbe6a955f1f799df4df..5b59f1099261f10d7270622655d4829c41b32201 100644
--- a/python/corsika/io/outputs/output.py
+++ b/python/corsika/io/outputs/output.py
@@ -114,7 +114,7 @@ class Output(ABC):
Any:
The data in its default format.
"""
- return self.astype()
+ return self.astype() # type: ignore
@staticmethod
def load_config(path: str) -> Dict[str, Any]:
diff --git a/python/corsika/io/outputs/radio_process.py b/python/corsika/io/outputs/radio_process.py
new file mode 100644
index 0000000000000000000000000000000000000000..4936f92e8a25cfbce852c046ee0ec34b39d6eb89
--- /dev/null
+++ b/python/corsika/io/outputs/radio_process.py
@@ -0,0 +1,157 @@
+"""
+ Read data written by a RadioProcess
+
+ (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+
+ This software is distributed under the terms of the GNU General Public
+ Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ the license.
+"""
+import logging
+import os.path as op
+from typing import Any
+
+import numpy as np
+import xarray as xr
+
+from .output import Output
+
+
+class RadioProcess(Output):
+ """
+ Read particle data from an RadioProcess.
+
+ This *currently* can be used to read data written by ZHS
+ or CoREAS implementations of the RadioProcess.
+ """
+
+ def __init__(self, path: str):
+ """
+ Initialize this radio process reader (and load the data).
+
+ Since each antenna can have a different sample rate and duration,
+ we can't load them into a shared array. Therefore, we load each
+ of the antennas into an XArray dataset.
+
+ Parameters
+ ----------
+ path: str
+ The path to the directory containing this output.
+ """
+ super().__init__(path)
+
+ # try and load our data
+ try:
+ self.__data = self.load_data(path)
+ except Exception as e:
+ logging.getLogger("corsika").warn(
+ f"An error occured loading a RadioProcess: {e}"
+ )
+
+ def load_data(self, path: str) -> xr.Dataset:
+ """
+ Load the data associated with this radio process.
+
+ Parameters
+ ----------
+ path: str
+ The path to the directory containing this output.
+
+ """
+
+ # get the list of antenna names
+ antennas = list(self.config["antennas"].keys())
+
+ # if there are no antennas,
+ if len(antennas) == 0:
+ logging.warn(f"No antennas were found for {self.config['name']}")
+
+ # we build up the XArray Dataset in this dictionary
+ dataset = {}
+
+ # loop over each of the antennas
+ for iant, name in enumerate(antennas):
+
+ # load the data file associated with this antenna
+ try:
+ data = np.load(f"{op.join(path, name)}.npz")
+ except Exception as e:
+ raise RuntimeError(
+ (
+ f"Unable to open file for antenna {name}"
+ f"in {self.config['name']} as {e}"
+ )
+ )
+
+ # if we get here, we have successfully loaded the antennas data file
+
+ # extract the sample times (in ns)
+ times = data["Time"]
+
+ # calculate the number of showers for this antenna
+ nshowers = len(list(data.keys())) - 1
+
+ # check that we got some events
+ if nshowers == 0:
+ logging.warn(f"Antenna {name} contains data for 0 showers.")
+
+ # create the array to store the waveforms for all the events
+ waveforms = np.zeros((nshowers, *data["0"].shape), dtype=np.float32)
+
+ # fill in the 'waveforms' array
+ for iev in np.arange(nshowers):
+ waveforms[iev, ...] = data[str(iev)]
+
+ # create the data array
+ showers = xr.DataArray(
+ waveforms,
+ coords=(np.arange(nshowers), ["x", "y", "z"], times), # type: ignore
+ dims=["shower", "pol", "time"],
+ )
+
+ # save this data array
+ dataset[name] = showers
+
+ return xr.Dataset(dataset)
+
+ def is_good(self) -> bool:
+ """
+ Returns true if this output has been read successfully
+ and has the correct files/state/etc.
+
+ Returns
+ -------
+ bool:
+ True if this is a good output.
+ """
+ return self.__data is not None
+
+ def astype(self, dtype: str = "xarray", **kwargs: Any) -> Any:
+ """
+ Load the antenna data from this process.
+
+ Parameters
+ ----------
+ dtype: str
+ The data format to return the data in (i.e. numpy, pandas, etc.)
+
+ Returns
+ -------
+ Any:
+ The return type of this method is determined by `dtype`.
+ """
+ if dtype == "xarray":
+ return self.__data
+ else:
+ raise ValueError(
+ (
+ f"Unknown format '{dtype}' for RadioProcess. "
+ "We currently only support ['xarray']."
+ )
+ )
+
+ def __repr__(self) -> str:
+ """
+ Return a string representation of this class.
+ """
+ return f"RadioProcess('{self.config['name']}', '{self.config['algorithm']}')"
diff --git a/python/setup.cfg b/python/setup.cfg
index f58d04eeb061f000b2d890aa26f542f70f3aa0c0..2c8d668ffdf9a3a8135dc52fc6e4c79b66c06e72 100644
--- a/python/setup.cfg
+++ b/python/setup.cfg
@@ -1,6 +1,6 @@
[flake8]
# use a slightly longer line to be consistent with black
-max-line-length = 88
+max-line-length = 90
# E231 is missing whitespace after a comma
# this contradicts the black formatting rules
diff --git a/python/setup.py b/python/setup.py
index 4bc73b3260b617ac3c55b2273d399b546ffdd3f6..0e7b1babd8d98fc1228627dd3f70dcb7910af8fc 100644
--- a/python/setup.py
+++ b/python/setup.py
@@ -32,7 +32,7 @@ setup(
keywords=["cosmic ray", "physics", "air shower", "simulation"],
packages=find_packages(),
python_requires=">=3.6*, <4",
- install_requires=["numpy", "pyyaml", "pyarrow", "boost_histogram"],
+ install_requires=["numpy", "pyyaml", "pyarrow", "boost_histogram", "xarray"],
extras_require={
"test": [
"pytest",
diff --git a/python/tests/io/test_hist.py b/python/tests/io/test_hist.py
index 7cf5b38fb7798c81680cf14d108e3b31dd3dcaf1..3bd11003ca663337f322ed9fea4a7f7559e75d03 100644
--- a/python/tests/io/test_hist.py
+++ b/python/tests/io/test_hist.py
@@ -16,7 +16,7 @@ import corsika
from .. import build_directory
# the directory containing 'testSaveBoostHistogram'
-bindir = op.join(build_directory, "Framework", "Utilities")
+bindir = op.join(build_directory, "bin")
def generate_hist() -> str:
@@ -32,14 +32,14 @@ def generate_hist() -> str:
"""
# we construct the name of the bin
- bin = op.join(bindir, "testSaveBoostHistogram")
+ bin = op.join(bindir, "testFramework")
# check if a histogram already exists
if op.exists(op.join(bin, "hist.npz")):
return op.join(bin, "hist.npz"), False
else:
# run the program - this generates "hist.npz" in the CWD
- subprocess.call(bin)
+ subprocess.call([bin, "saveHistogram"])
return op.join(os.getcwd(), "hist.npz"), True
diff --git a/tests/modules/CMakeLists.txt b/tests/modules/CMakeLists.txt
index 2073a7d50f44a6647163cd40618fd73b8293d95f..78e3d5faff9e66ccd4a2d9720fc294717e1cd850 100644
--- a/tests/modules/CMakeLists.txt
+++ b/tests/modules/CMakeLists.txt
@@ -12,6 +12,7 @@ set (test_modules_sources
testParticleCut.cpp
testSibyll.cpp
testEpos.cpp
+ testRadio.cpp
)
CORSIKA_ADD_TEST (testModules SOURCES ${test_modules_sources})
diff --git a/tests/modules/testRadio.cpp b/tests/modules/testRadio.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4d019d9ecad783ee3e501d5b22a8130ca73026d2
--- /dev/null
+++ b/tests/modules/testRadio.cpp
@@ -0,0 +1,800 @@
+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+#include <catch2/catch.hpp>
+
+#include <corsika/modules/radio/ZHS.hpp>
+#include <corsika/modules/radio/CoREAS.hpp>
+#include <corsika/modules/radio/antennas/TimeDomainAntenna.hpp>
+#include <corsika/modules/radio/detectors/AntennaCollection.hpp>
+#include <corsika/modules/radio/propagators/StraightPropagator.hpp>
+#include <corsika/modules/radio/propagators/SimplePropagator.hpp>
+#include <corsika/modules/radio/propagators/SignalPath.hpp>
+#include <corsika/modules/radio/propagators/RadioPropagator.hpp>
+
+#include <vector>
+#include <istream>
+#include <fstream>
+#include <iostream>
+
+#include <corsika/media/Environment.hpp>
+#include <corsika/media/FlatExponential.hpp>
+#include <corsika/media/HomogeneousMedium.hpp>
+#include <corsika/media/IMagneticFieldModel.hpp>
+#include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+#include <corsika/media/MediumPropertyModel.hpp>
+#include <corsika/media/UniformMagneticField.hpp>
+#include <corsika/media/SlidingPlanarExponential.hpp>
+#include <corsika/media/IMediumModel.hpp>
+#include <corsika/media/IRefractiveIndexModel.hpp>
+#include <corsika/media/UniformRefractiveIndex.hpp>
+#include <corsika/media/ExponentialRefractiveIndex.hpp>
+#include <corsika/media/VolumeTreeNode.hpp>
+#include <corsika/media/CORSIKA7Atmospheres.hpp>
+
+#include <corsika/framework/geometry/CoordinateSystem.hpp>
+#include <corsika/framework/geometry/Line.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+#include <corsika/setup/SetupStack.hpp>
+#include <corsika/setup/SetupTrajectory.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/PhysicalConstants.hpp>
+
+#include <corsika/output/OutputManager.hpp>
+
+using namespace corsika;
+
+double constexpr absMargin = 1.0e-7;
+
+template <typename TInterface>
+using MyExtraEnv =
+ UniformRefractiveIndex<MediumPropertyModel<UniformMagneticField<TInterface>>>;
+
+TEST_CASE("Radio", "[processes]") {
+
+ logging::set_level(logging::level::debug);
+
+ SECTION("CoREAS process") {
+
+ // This serves as a compiler test for any changes in the CoREAS algorithm
+ // Environment
+ using EnvironmentInterface =
+ IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
+
+ // using EnvType = setup::Environment;
+ using EnvType = Environment<EnvironmentInterface>;
+ EnvType envCoREAS;
+ CoordinateSystemPtr const& rootCSCoREAS = envCoREAS.getCoordinateSystem();
+ Point const center{rootCSCoREAS, 0_m, 0_m, 0_m};
+
+ // 1.000327,
+ create_5layer_atmosphere<EnvironmentInterface, MyExtraEnv>(
+ envCoREAS, AtmosphereId::LinsleyUSStd, center, 1.000327, Medium::AirDry1Atm,
+ MagneticFieldVector{rootCSCoREAS, 0_T, 50_uT, 0_T});
+
+ // now create antennas and detectors
+ // the antennas location
+ const auto point1{Point(envCoREAS.getCoordinateSystem(), 100_m, 2_m, 3_m)};
+ const auto point2{Point(envCoREAS.getCoordinateSystem(), 4_m, 80_m, 6_m)};
+ const auto point3{Point(envCoREAS.getCoordinateSystem(), 7_m, 8_m, 9_m)};
+ const auto point4{Point(envCoREAS.getCoordinateSystem(), 5_m, 5_m, 10_m)};
+
+ // create times for the antenna
+ const TimeType t1{0_s};
+ const TimeType t2{10_s};
+ const InverseTimeType t3{1e+3_Hz};
+ const TimeType t4{11_s};
+
+ // check that I can create an antenna at (1, 2, 3)
+ TimeDomainAntenna ant1("antenna_name", point1, rootCSCoREAS, t1, t2, t3, t1);
+ TimeDomainAntenna ant2("antenna_name2", point2, rootCSCoREAS, t1, t2, t3, t1);
+
+ // construct a radio detector instance to store our antennas
+ AntennaCollection<TimeDomainAntenna> detector;
+
+ // add the antennas to the detector
+ detector.addAntenna(ant1);
+ detector.addAntenna(ant2);
+
+ // create a particle
+ const Code particle{Code::Electron};
+ // const Code particle{Code::Proton};
+ const auto pmass{get_mass(particle)};
+
+ VelocityVector v0(rootCSCoREAS, {5e+2_m / second, 5e+2_m / second, 5e+2_m / second});
+
+ Vector B0(rootCSCoREAS, 5_T, 5_T, 5_T);
+
+ Line const line(point3, v0);
+
+ auto const k{1_m * ((1_m) / ((1_s * 1_s) * 1_V))};
+
+ auto const t = 1e-12_s;
+ LeapFrogTrajectory base(point4, v0, B0, k, t);
+ // std::cout << "Leap Frog Trajectory is: " << base << std::endl;
+
+ // create a new stack for each trial
+ setup::Stack<EnvType> stack;
+
+ // construct an energy
+ const HEPEnergyType E0{1_TeV};
+
+ // compute the necessary momentumn
+ const HEPMomentumType P0{sqrt(E0 * E0 - pmass * pmass)};
+
+ // and create the momentum vector
+ const auto plab{MomentumVector(rootCSCoREAS, {0_GeV, 0_GeV, P0})};
+
+ // and create the location of the particle in this coordinate system
+ const Point pos(rootCSCoREAS, 50_m, 10_m, 80_m);
+
+ // add the particle to the stack
+ auto const particle1{stack.addParticle(std::make_tuple(
+ particle, calculate_kinetic_energy(plab.getNorm(), get_mass(particle)),
+ plab.normalized(), pos, 0_ns))};
+
+ auto const charge_{get_charge(particle1.getPID())};
+
+ // create a radio process instance using CoREAS
+ RadioProcess<decltype(detector),
+ CoREAS<decltype(detector), decltype(StraightPropagator(envCoREAS))>,
+ decltype(StraightPropagator(envCoREAS))>
+ coreas(detector, envCoREAS);
+
+ Step step(particle1, base);
+ // check doContinuous and simulate methods
+ coreas.doContinuous(step, true);
+ } // END: SECTION("CoREAS process")
+
+ SECTION("ZHS process") {
+
+ // This section serves as a compiler test for any changes in the ZHS algorithm
+ // Environment
+ using IModelInterface =
+ IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
+ using AtmModel = UniformRefractiveIndex<
+ MediumPropertyModel<UniformMagneticField<HomogeneousMedium<IModelInterface>>>>;
+ using EnvType = Environment<AtmModel>;
+ EnvType envZHS;
+ CoordinateSystemPtr const& rootCSZHS = envZHS.getCoordinateSystem();
+ // get the center point
+ Point const center{rootCSZHS, 0_m, 0_m, 0_m};
+ // a refractive index
+ const double ri_{1.000327};
+
+ // the constant density
+ const auto density{19.2_g / cube(1_cm)};
+
+ // the composition we use for the homogeneous medium
+ NuclearComposition const protonComposition({Code::Proton}, {1.});
+
+ // create magnetic field vector
+ Vector B1(rootCSZHS, 0_T, 0_T, 1_T);
+
+ auto Medium = EnvType::createNode<Sphere>(
+ center, 1_km * std::numeric_limits<double>::infinity());
+
+ auto const props = Medium->setModelProperties<AtmModel>(ri_, Medium::AirDry1Atm, B1,
+ density, protonComposition);
+ envZHS.getUniverse()->addChild(std::move(Medium));
+
+ // the antennas location
+ const auto point1{Point(envZHS.getCoordinateSystem(), 100_m, 2_m, 3_m)};
+ const auto point2{Point(envZHS.getCoordinateSystem(), 4_m, 80_m, 6_m)};
+ const auto point3{Point(envZHS.getCoordinateSystem(), 7_m, 8_m, 9_m)};
+ const auto point4{Point(envZHS.getCoordinateSystem(), 5_m, 5_m, 10_m)};
+
+ // create times for the antenna
+ const TimeType t1{0_s};
+ const TimeType t2{10_s};
+ const InverseTimeType t3{1e+3_Hz};
+ const TimeType t4{11_s};
+
+ // check that I can create an antenna at (1, 2, 3)
+ TimeDomainAntenna ant1("antenna_zhs", point1, rootCSZHS, t1, t2, t3, t1);
+ TimeDomainAntenna ant2("antenna_zhs2", point2, rootCSZHS, t1, t2, t3, t1);
+
+ // construct a radio detector instance to store our antennas
+ AntennaCollection<TimeDomainAntenna> detector;
+
+ // add the antennas to the detector
+ detector.addAntenna(ant1);
+ detector.addAntenna(ant2);
+
+ // create a particle
+ auto const particle{Code::Electron};
+ const auto pmass{get_mass(particle)};
+
+ VelocityVector v0(rootCSZHS, {5e+2_m / second, 5e+2_m / second, 5e+2_m / second});
+
+ Vector B0(rootCSZHS, 5_T, 5_T, 5_T);
+
+ Line const line(point3, v0);
+
+ auto const k{1_m * ((1_m) / ((1_s * 1_s) * 1_V))};
+
+ auto const t = 1e-12_s;
+ LeapFrogTrajectory base(point4, v0, B0, k, t);
+ // std::cout << "Leap Frog Trajectory is: " << base << std::endl;
+
+ // create a new stack for each trial
+ setup::Stack<EnvType> stack;
+
+ // construct an energy
+ const HEPEnergyType E0{1_TeV};
+
+ // compute the necessary momentum
+ const HEPMomentumType P0{sqrt(E0 * E0 - pmass * pmass)};
+
+ // and create the momentum vector
+ const auto plab{MomentumVector(rootCSZHS, {0_GeV, 0_GeV, P0})};
+
+ // and create the location of the particle in this coordinate system
+ const Point pos(rootCSZHS, 50_m, 10_m, 80_m);
+
+ // add the particle to the stack
+ auto const particle1{stack.addParticle(std::make_tuple(
+ particle, calculate_kinetic_energy(plab.getNorm(), get_mass(particle)),
+ plab.normalized(), pos, 0_ns))};
+
+ auto const charge_{get_charge(particle1.getPID())};
+
+ // create a radio process instance using ZHS
+ RadioProcess<
+ AntennaCollection<TimeDomainAntenna>,
+ ZHS<AntennaCollection<TimeDomainAntenna>, decltype(StraightPropagator(envZHS))>,
+ decltype(StraightPropagator(envZHS))>
+ zhs(detector, envZHS);
+
+ Step step(particle1, base);
+ // check doContinuous and simulate methods
+ zhs.doContinuous(step, true);
+
+ } // END: SECTION("ZHS process")
+
+ SECTION("Radio extreme cases") {
+
+ // Environment
+ using EnvironmentInterface =
+ IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
+
+ using EnvType = Environment<EnvironmentInterface>;
+ EnvType envRadio;
+ CoordinateSystemPtr const& rootCSRadio = envRadio.getCoordinateSystem();
+ Point const center{rootCSRadio, 0_m, 0_m, 0_m};
+
+ create_5layer_atmosphere<EnvironmentInterface, MyExtraEnv>(
+ envRadio, AtmosphereId::LinsleyUSStd, center, 1.415, Medium::AirDry1Atm,
+ MagneticFieldVector{rootCSRadio, 0_T, 50_uT, 0_T});
+
+ // now create antennas and detectors
+ // the antennas location
+ const auto point1{Point(envRadio.getCoordinateSystem(), 0_m, 0_m, 0_m)};
+
+ // track points
+ Point const point_1(rootCSRadio, {1_m, 1_m, 0_m});
+ Point const point_2(rootCSRadio, {2_km, 1_km, 0_m});
+ Point const point_4(rootCSRadio, {0_m, 1_m, 0_m});
+
+ // create times for the antenna
+ const TimeType start{0_s};
+ const TimeType duration{100_ns};
+ const InverseTimeType sample{1e+12_Hz};
+ const TimeType duration_dummy{2_s};
+ const InverseTimeType sample_dummy{1_Hz};
+
+ // check that I can create an antenna at (1, 2, 3)
+ TimeDomainAntenna ant1("antenna_name", point1, rootCSRadio, start, duration, sample,
+ start);
+ TimeDomainAntenna ant2("dummy", point1, rootCSRadio, start, duration_dummy,
+ sample_dummy, start);
+ // construct a radio detector instance to store our antennas
+ AntennaCollection<TimeDomainAntenna> detector;
+ AntennaCollection<TimeDomainAntenna> detector_dummy;
+ // add the antennas to the detector
+ detector.addAntenna(ant1);
+ detector_dummy.addAntenna(ant2);
+
+ // create a new stack for each trial
+ setup::Stack<EnvType> stack;
+ // create a particle
+ const Code particle{Code::Electron};
+ const Code particle2{Code::Proton};
+
+ const auto pmass{get_mass(particle)};
+ const auto pmass2{get_mass(particle2)};
+ // construct an energy
+ const HEPEnergyType E0{1_TeV};
+ // compute the necessary momentumn
+ const HEPMomentumType P0{sqrt(E0 * E0 - pmass * pmass)};
+ // and create the momentum vector
+ const auto plab{MomentumVector(rootCSRadio, {P0, 0_GeV, 0_GeV})};
+ // add the particle to the stack
+ auto const particle_stack{stack.addParticle(std::make_tuple(
+ particle, calculate_kinetic_energy(plab.getNorm(), get_mass(particle)),
+ plab.normalized(), point_1, 0_ns))};
+
+ // particle stack with proton
+ auto const particle_stack_proton{stack.addParticle(std::make_tuple(
+ particle2, calculate_kinetic_energy(plab.getNorm(), get_mass(particle2)),
+ plab.normalized(), point_1, 0_ns))};
+
+ // feed radio with a proton track to check that it skips that track.
+ TimeType tp{(point_2 - point_1).getNorm() / (0.999 * constants::c)};
+ VelocityVector vp{(point_2 - point_1) / tp};
+ Line lp{point_1, vp};
+ StraightTrajectory track_p{lp, tp};
+ Step step_proton(particle_stack_proton, track_p);
+
+ // feed radio with a track that triggers zhs like approx in coreas and fraunhofer
+ // limit check for zhs
+ TimeType th{(point_4 - point1).getNorm() / constants::c};
+ VelocityVector vh{(point_4 - point1) / th};
+ Line lh{point1, vh};
+ StraightTrajectory track_h{lh, th};
+ Step step_h(particle_stack, track_h);
+
+ // create radio process instances
+ RadioProcess<decltype(detector),
+ CoREAS<decltype(detector), decltype(SimplePropagator(envRadio))>,
+ decltype(SimplePropagator(envRadio))>
+ coreas(detector, envRadio);
+
+ RadioProcess<decltype(detector),
+ ZHS<decltype(detector), decltype(SimplePropagator(envRadio))>,
+ decltype(SimplePropagator(envRadio))>
+ zhs(detector, envRadio);
+
+ coreas.doContinuous(step_proton, true);
+ zhs.doContinuous(step_proton, true);
+ coreas.doContinuous(step_h, true);
+ zhs.doContinuous(step_h, true);
+
+ // create radio processes with "dummy" antenna to trigger extreme time-binning
+ RadioProcess<decltype(detector_dummy),
+ CoREAS<decltype(detector_dummy), decltype(SimplePropagator(envRadio))>,
+ decltype(SimplePropagator(envRadio))>
+ coreas_dummy(detector_dummy, envRadio);
+
+ RadioProcess<decltype(detector_dummy),
+ ZHS<decltype(detector_dummy), decltype(SimplePropagator(envRadio))>,
+ decltype(SimplePropagator(envRadio))>
+ zhs_dummy(detector_dummy, envRadio);
+
+ coreas_dummy.doContinuous(step_h, true);
+ zhs_dummy.doContinuous(step_h, true);
+
+ } // END: SECTION("Radio extreme cases")
+
+} // END: TEST_CASE("Radio", "[processes]")
+
+TEST_CASE("Antennas") {
+
+ SECTION("TimeDomainAntenna") {
+
+ // create an environment so we can get a coordinate system
+ using EnvType = Environment<IRefractiveIndexModel<IMediumModel>>;
+ EnvType env6;
+
+ using UniRIndex =
+ UniformRefractiveIndex<HomogeneousMedium<IRefractiveIndexModel<IMediumModel>>>;
+
+ // the antenna location
+ const auto point1{Point(env6.getCoordinateSystem(), 1_m, 2_m, 3_m)};
+ const auto point2{Point(env6.getCoordinateSystem(), 4_m, 5_m, 6_m)};
+
+ // get a coordinate system
+ const CoordinateSystemPtr rootCS6 = env6.getCoordinateSystem();
+
+ auto Medium6 = EnvType::createNode<Sphere>(
+ Point{rootCS6, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
+
+ auto const props6 = Medium6->setModelProperties<UniRIndex>(
+ 1, 1_kg / (1_m * 1_m * 1_m), NuclearComposition({Code::Nitrogen}, {1.}));
+
+ env6.getUniverse()->addChild(std::move(Medium6));
+
+ // create times for the antenna
+ const TimeType t1{10_s};
+ const TimeType t2{10_s};
+ const InverseTimeType t3{1 / 1_s};
+ const TimeType t4{11_s};
+
+ // check that I can create an antenna at (1, 2, 3)
+ TimeDomainAntenna ant1("antenna_name", point1, rootCS6, t1, t2, t3, t1);
+ TimeDomainAntenna ant2("antenna_name2", point2, rootCS6, t4, t2, t3, t4);
+
+ // assert that the antenna name is correct
+ REQUIRE(ant1.getName() == "antenna_name");
+ REQUIRE(ant2.getName() == "antenna_name2");
+
+ // and check that the antenna is at the right location
+ REQUIRE((ant1.getLocation() - point1).getNorm() < 1e-12 * 1_m);
+ REQUIRE((ant2.getLocation() - point2).getNorm() < 1e-12 * 1_m);
+
+ // construct a radio detector instance to store our antennas
+ AntennaCollection<TimeDomainAntenna> detector;
+
+ // add this antenna to the process
+ detector.addAntenna(ant1);
+ detector.addAntenna(ant2);
+ CHECK(detector.size() == 2);
+
+ // get a unit vector
+ Vector<dimensionless_d> v1(rootCS6, {0, 0, 1});
+ Vector<ElectricFieldType::dimension_type> v11(rootCS6,
+ {10_V / 1_m, 10_V / 1_m, 10_V / 1_m});
+
+ Vector<dimensionless_d> v2(rootCS6, {0, 1, 0});
+ Vector<ElectricFieldType::dimension_type> v22(rootCS6,
+ {20_V / 1_m, 20_V / 1_m, 20_V / 1_m});
+
+ // use receive methods
+ ant1.receive(15_s, v1, v11);
+ ant2.receive(16_s, v2, v22);
+
+ // use getDataX,Y,Z() and getAxis() methods
+ auto Ex = ant1.getWaveformX();
+ CHECK(Ex[5] - 10 == 0);
+ auto tx = ant1.getAxis();
+ CHECK(tx[5] - 5 * 1_s / 1_ns == Approx(0.0));
+ auto Ey = ant1.getWaveformY();
+ CHECK(Ey[5] - 10 == 0);
+ auto Ez = ant1.getWaveformZ();
+ CHECK(Ez[5] - 10 == 0);
+ auto ty = ant1.getAxis();
+ auto tz = ant1.getAxis();
+ CHECK(tx[5] - ty[5] == 0);
+ CHECK(ty[5] - tz[5] == 0);
+ auto Ex2 = ant2.getWaveformX();
+ CHECK(Ex2[5] - 20 == 0);
+ auto Ey2 = ant2.getWaveformY();
+ CHECK(Ey2[5] - 20 == 0);
+ auto Ez2 = ant2.getWaveformZ();
+ CHECK(Ez2[5] - 20 == 0);
+
+ // the following creates a star-shaped pattern of antennas in the ground
+ AntennaCollection<TimeDomainAntenna> detector__;
+ const auto point11{Point(env6.getCoordinateSystem(), 1000_m, 20_m, 30_m)};
+ const TimeType t2222{1e-6_s};
+ const InverseTimeType t3333{1e+9_Hz};
+
+ std::vector<std::string> antenna_names;
+ std::vector<Point> antenna_locations;
+ for (auto radius_ = 100_m; radius_ <= 200_m; radius_ += 100_m) {
+ for (auto phi_ = 0; phi_ <= 315; phi_ += 45) {
+ auto phiRad_ = phi_ / 180. * M_PI;
+ auto const point_{Point(env6.getCoordinateSystem(), radius_ * cos(phiRad_),
+ radius_ * sin(phiRad_), 0_m)};
+ antenna_locations.push_back(point_);
+ auto time__{(point11 - point_).getNorm() / constants::c};
+ const int rr_ = static_cast<int>(radius_ / 1_m);
+ std::string var_ = "antenna_R=" + std::to_string(rr_) +
+ "_m-Phi=" + std::to_string(phi_) + "degrees";
+ antenna_names.push_back(var_);
+ TimeDomainAntenna ant111(var_, point_, rootCS6, time__, t2222, t3333, time__);
+ detector__.addAntenna(ant111);
+ }
+ }
+
+ CHECK(detector__.size() == 16);
+ CHECK(detector__.getAntennas().size() == 16);
+ int i = 0;
+ // this prints out the antenna names and locations
+ for (auto const antenna : detector__.getAntennas()) {
+ CHECK(antenna.getName() == antenna_names[i]);
+ CHECK(distance(antenna.getLocation(), antenna_locations[i]) / 1_m == 0);
+ i++;
+ }
+
+ // Check the .at() method for radio detectors
+ for (size_t i = 0; i <= (detector__.size() - 1); i++) {
+ CHECK(detector__.at(i).getName() == antenna_names[i]);
+ CHECK(distance(detector__.at(i).getLocation(), antenna_locations[i]) / 1_m == 0);
+ }
+
+ } // END: SECTION("TimeDomainAntenna")
+
+} // END: TEST_CASE("Antennas")
+
+TEST_CASE("Propagators") {
+
+ SECTION("Simple Propagator w/ Uniform Refractive Index") {
+
+ // create a suitable environment
+ using IModelInterface =
+ IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
+ using AtmModel = UniformRefractiveIndex<
+ MediumPropertyModel<UniformMagneticField<HomogeneousMedium<IModelInterface>>>>;
+ using EnvType = Environment<AtmModel>;
+ EnvType env;
+ CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
+ // get the center point
+ Point const center{rootCS, 0_m, 0_m, 0_m};
+ // a refractive index for the vacuum
+ const double ri_{1};
+ // the constant density
+ const auto density{19.2_g / cube(1_cm)};
+ // the composition we use for the homogeneous medium
+ NuclearComposition const Composition({Code::Nitrogen}, {1.});
+ // create magnetic field vector
+ Vector B1(rootCS, 0_T, 0_T, 0.3809_T);
+ // create a Sphere for the medium
+ auto Medium = EnvType::createNode<Sphere>(
+ center, 1_km * std::numeric_limits<double>::infinity());
+ // set the environment properties
+ auto const props = Medium->setModelProperties<AtmModel>(ri_, Medium::AirDry1Atm, B1,
+ density, Composition);
+ // bind things together
+ env.getUniverse()->addChild(std::move(Medium));
+
+ // get some points
+ Point p0(rootCS, {0_m, 0_m, 0_m});
+ Point p10(rootCS, {0_m, 0_m, 10_m});
+
+ // get a unit vector
+ Vector<dimensionless_d> v1(rootCS, {0, 0, 1});
+ Vector<dimensionless_d> v2(rootCS, {0, 0, -1});
+
+ // get a geometrical path of points
+ Path P1({p0, p10});
+
+ // construct a Straight Propagator given the uniform refractive index environment
+ SimplePropagator SP(env);
+
+ // store the outcome of the Propagate method to paths_
+ auto const paths_ = SP.propagate(p0, p10, 1_m);
+
+ // perform checks to paths_ components
+ for (auto const& path : paths_) {
+ CHECK((path.propagation_time_ / 1_s) -
+ (((p10 - p0).getNorm() / constants::c) / 1_s) ==
+ Approx(0));
+ CHECK(path.average_refractive_index_ == Approx(1));
+ CHECK(path.refractive_index_source_ == Approx(1));
+ CHECK(path.refractive_index_destination_ == Approx(1));
+ CHECK(path.emit_.getComponents() == v1.getComponents());
+ CHECK(path.receive_.getComponents() == v2.getComponents());
+ CHECK(path.R_distance_ == 10_m);
+ CHECK(std::equal(P1.begin(), P1.end(), Path(path.points_).begin(),
+ [](Point a, Point b) { return (a - b).getNorm() / 1_m < 1e-5; }));
+ }
+
+ } // END: SECTION("Simple Propagator w/ Uniform Refractive Index")
+
+ // check that I can create working Straight Propagators in different environments
+ SECTION("Straight Propagator w/ Uniform Refractive Index") {
+
+ // create a suitable environment
+ using IModelInterface =
+ IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
+ using AtmModel = UniformRefractiveIndex<
+ MediumPropertyModel<UniformMagneticField<HomogeneousMedium<IModelInterface>>>>;
+ using EnvType = Environment<AtmModel>;
+ EnvType env;
+ CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
+ // get the center point
+ Point const center{rootCS, 0_m, 0_m, 0_m};
+ // a refractive index for the vacuum
+ const double ri_{1};
+ // the constant density
+ const auto density{19.2_g / cube(1_cm)};
+ // the composition we use for the homogeneous medium
+ NuclearComposition const Composition({Code::Nitrogen}, {1.});
+ // create magnetic field vector
+ Vector B1(rootCS, 0_T, 0_T, 0.3809_T);
+ // create a Sphere for the medium
+ auto Medium = EnvType::createNode<Sphere>(
+ center, 1_km * std::numeric_limits<double>::infinity());
+ // set the environment properties
+ auto const props = Medium->setModelProperties<AtmModel>(ri_, Medium::AirDry1Atm, B1,
+ density, Composition);
+ // bind things together
+ env.getUniverse()->addChild(std::move(Medium));
+
+ // get some points
+ Point p0(rootCS, {0_m, 0_m, 0_m});
+ Point p1(rootCS, {0_m, 0_m, 1_m});
+ Point p2(rootCS, {0_m, 0_m, 2_m});
+ Point p3(rootCS, {0_m, 0_m, 3_m});
+ Point p4(rootCS, {0_m, 0_m, 4_m});
+ Point p5(rootCS, {0_m, 0_m, 5_m});
+ Point p6(rootCS, {0_m, 0_m, 6_m});
+ Point p7(rootCS, {0_m, 0_m, 7_m});
+ Point p8(rootCS, {0_m, 0_m, 8_m});
+ Point p9(rootCS, {0_m, 0_m, 9_m});
+ Point p10(rootCS, {0_m, 0_m, 10_m});
+ Point p30(rootCS, {0_m, 0_m, 30000_m});
+
+ // get a unit vector
+ Vector<dimensionless_d> v1(rootCS, {0, 0, 1});
+ Vector<dimensionless_d> v2(rootCS, {0, 0, -1});
+
+ // get a geometrical path of points
+ Path P1({p0, p1, p2, p3, p4, p5, p6, p7, p8, p9, p10});
+
+ // construct a Straight Propagator given the uniform refractive index environment
+ StraightPropagator SP(env);
+
+ // store the outcome of the Propagate method to paths_
+ auto const paths_ = SP.propagate(p0, p10, 1_m);
+
+ // perform checks to paths_ components
+ for (auto const& path : paths_) {
+ CHECK((path.propagation_time_ / 1_s) -
+ (((p10 - p0).getNorm() / constants::c) / 1_s) ==
+ Approx(0).margin(absMargin));
+ CHECK(path.average_refractive_index_ == Approx(1));
+ CHECK(path.refractive_index_source_ == Approx(1));
+ CHECK(path.refractive_index_destination_ == Approx(1));
+ CHECK(path.emit_.getComponents() == v1.getComponents());
+ CHECK(path.receive_.getComponents() == v2.getComponents());
+ CHECK(path.R_distance_ == 10_m);
+ CHECK(std::equal(P1.begin(), P1.end(), Path(path.points_).begin(),
+ [](Point a, Point b) { return (a - b).getNorm() / 1_m < 1e-5; }));
+ }
+
+ // get another path to different points
+ auto const paths2_{SP.propagate(p0, p30, 909_m)};
+
+ for (auto const& path : paths2_) {
+ CHECK((path.propagation_time_ / 1_s) -
+ (((p30 - p0).getNorm() / constants::c) / 1_s) ==
+ Approx(0).margin(absMargin));
+ CHECK(path.average_refractive_index_ == Approx(1));
+ CHECK(path.refractive_index_source_ == Approx(1));
+ CHECK(path.refractive_index_destination_ == Approx(1));
+ CHECK(path.R_distance_ == 30000_m);
+ }
+
+ // get a third path using a weird stepsize
+ auto const paths3_{SP.propagate(p0, p30, 731.89_m)};
+
+ for (auto const& path : paths3_) {
+ CHECK((path.propagation_time_ / 1_s) -
+ (((p30 - p0).getNorm() / constants::c) / 1_s) ==
+ Approx(0).margin(absMargin));
+ CHECK(path.average_refractive_index_ == Approx(1));
+ CHECK(path.refractive_index_source_ == Approx(1));
+ CHECK(path.refractive_index_destination_ == Approx(1));
+ CHECK(path.R_distance_ == 30000_m);
+ }
+
+ CHECK(paths_.size() == 1);
+ CHECK(paths2_.size() == 1);
+ CHECK(paths3_.size() == 1);
+ } // END: SECTION("Straight Propagator w/ Uniform Refractive Index")
+
+ SECTION("Straight Propagator w/ Exponential Refractive Index") {
+
+ // create an environment with exponential refractive index (n_0 = 1 & lambda = 0)
+ using ExpoRIndex = ExponentialRefractiveIndex<
+ HomogeneousMedium<IRefractiveIndexModel<IMediumModel>>>;
+
+ using EnvType = Environment<IRefractiveIndexModel<IMediumModel>>;
+ EnvType env1;
+
+ // get another coordinate system
+ const CoordinateSystemPtr rootCS1 = env1.getCoordinateSystem();
+
+ // the center of the earth
+ Point const center1_{rootCS1, 0_m, 0_m, 0_m};
+ LengthType const radius_{0_m};
+
+ auto Medium1 = EnvType::createNode<Sphere>(
+ Point{rootCS1, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
+
+ auto const props1 = Medium1->setModelProperties<ExpoRIndex>(
+ 1, 0 / 1_m, center1_, radius_, 1_kg / (1_m * 1_m * 1_m),
+ NuclearComposition({Code::Nitrogen}, {1.}));
+
+ env1.getUniverse()->addChild(std::move(Medium1));
+
+ // get some points
+ Point pp0(rootCS1, {0_m, 0_m, 0_m});
+ Point pp1(rootCS1, {0_m, 0_m, 1_m});
+ Point pp2(rootCS1, {0_m, 0_m, 2_m});
+ Point pp3(rootCS1, {0_m, 0_m, 3_m});
+ Point pp4(rootCS1, {0_m, 0_m, 4_m});
+ Point pp5(rootCS1, {0_m, 0_m, 5_m});
+ Point pp6(rootCS1, {0_m, 0_m, 6_m});
+ Point pp7(rootCS1, {0_m, 0_m, 7_m});
+ Point pp8(rootCS1, {0_m, 0_m, 8_m});
+ Point pp9(rootCS1, {0_m, 0_m, 9_m});
+ Point pp10(rootCS1, {0_m, 0_m, 10_m});
+
+ // get a unit vector
+ Vector<dimensionless_d> vv1(rootCS1, {0, 0, 1});
+ Vector<dimensionless_d> vv2(rootCS1, {0, 0, -1});
+
+ // get a geometrical path of points
+ Path PP1({pp0, pp1, pp2, pp3, pp4, pp5, pp6, pp7, pp8, pp9, pp10});
+
+ // construct a Straight Propagator given the exponential refractive index environment
+ StraightPropagator SP1(env1);
+
+ // store the outcome of Propagate method to paths1_
+ auto const paths1_ = SP1.propagate(pp0, pp10, 1_m);
+
+ // perform checks to paths1_ components (this is just a sketch for now)
+ for (auto const& path : paths1_) {
+ CHECK((path.propagation_time_ / 1_s) -
+ (((pp10 - pp0).getNorm() / constants::c) / 1_s) ==
+ Approx(0).margin(absMargin));
+ CHECK(path.average_refractive_index_ == Approx(1));
+ CHECK(path.refractive_index_source_ == Approx(1));
+ CHECK(path.refractive_index_destination_ == Approx(1));
+ CHECK(path.emit_.getComponents() == vv1.getComponents());
+ CHECK(path.receive_.getComponents() == vv2.getComponents());
+ CHECK(path.R_distance_ == 10_m);
+ CHECK(std::equal(PP1.begin(), PP1.end(), Path(path.points_).begin(),
+ [](Point a, Point b) { return (a - b).getNorm() / 1_m < 1e-5; }));
+ }
+
+ CHECK(paths1_.size() == 1);
+
+ /*
+ * A second environment with another exponential refractive index
+ */
+
+ // create an environment with exponential refractive index (n_0 = 2 & lambda = 2)
+ using ExpoRIndex = ExponentialRefractiveIndex<
+ HomogeneousMedium<IRefractiveIndexModel<IMediumModel>>>;
+
+ using EnvType = Environment<IRefractiveIndexModel<IMediumModel>>;
+ EnvType env2;
+
+ // get another coordinate system
+ const CoordinateSystemPtr rootCS2 = env2.getCoordinateSystem();
+
+ // the center of the earth
+ Point const center2_{rootCS2, 0_m, 0_m, 0_m};
+
+ auto Medium2 = EnvType::createNode<Sphere>(
+ Point{rootCS2, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
+
+ auto const props2 = Medium2->setModelProperties<ExpoRIndex>(
+ 2, 2 / 1_m, center2_, radius_, 1_kg / (1_m * 1_m * 1_m),
+ NuclearComposition({Code::Nitrogen}, {1.}));
+
+ env2.getUniverse()->addChild(std::move(Medium2));
+
+ // get some points
+ Point ppp0(rootCS2, {0_m, 0_m, 0_m});
+ Point ppp10(rootCS2, {0_m, 0_m, 10_m});
+
+ // get a unit vector
+ Vector<dimensionless_d> vvv1(rootCS2, {0, 0, 1});
+ Vector<dimensionless_d> vvv2(rootCS2, {0, 0, -1});
+
+ // construct a Straight Propagator given the exponential refractive index environment
+ StraightPropagator SP2(env2);
+
+ // store the outcome of Propagate method to paths1_
+ auto const paths2_ = SP2.propagate(ppp0, ppp10, 1_m);
+
+ // perform checks to paths1_ components (this is just a sketch for now)
+ for (auto const& path : paths2_) {
+ CHECK((path.propagation_time_ / 1_s) -
+ ((3.177511688_m / (3 * constants::c)) / 1_s) ==
+ Approx(0).margin(absMargin));
+ CHECK(path.average_refractive_index_ == Approx(0.210275935));
+ CHECK(path.refractive_index_source_ == Approx(2));
+ // CHECK(path.refractive_index_destination_ == Approx(0.0000000041));
+ CHECK(path.emit_.getComponents() == vvv1.getComponents());
+ CHECK(path.receive_.getComponents() == vvv2.getComponents());
+ CHECK(path.R_distance_ == 10_m);
+ }
+
+ CHECK(paths2_.size() == 1);
+
+ } // END: SECTION("Straight Propagator w/ Exponential Refractive Index")
+
+} // END: TEST_CASE("Propagators")