IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 557ea1de authored by Nikos Karastathis's avatar Nikos Karastathis :ocean:
Browse files

clang format

parent 00914233
No related branches found
No related tags found
1 merge request!329Radio interface
Showing
with 1228 additions and 1148 deletions
...@@ -54,8 +54,8 @@ namespace corsika { ...@@ -54,8 +54,8 @@ namespace corsika {
} }
template <typename TTimeType, typename TSpaceVecType> template <typename TTimeType, typename TSpaceVecType>
inline FourVector<TTimeType, TSpaceVecType>& FourVector<TTimeType, TSpaceVecType>:: inline FourVector<TTimeType, TSpaceVecType>&
operator+=(FourVector const& b) { FourVector<TTimeType, TSpaceVecType>::operator+=(FourVector const& b) {
timeLike_ += b.timeLike_; timeLike_ += b.timeLike_;
spaceLike_ += b.spaceLike_; spaceLike_ += b.spaceLike_;
...@@ -63,39 +63,39 @@ namespace corsika { ...@@ -63,39 +63,39 @@ namespace corsika {
} }
template <typename TTimeType, typename TSpaceVecType> template <typename TTimeType, typename TSpaceVecType>
inline FourVector<TTimeType, TSpaceVecType>& FourVector<TTimeType, TSpaceVecType>:: inline FourVector<TTimeType, TSpaceVecType>&
operator-=(FourVector const& b) { FourVector<TTimeType, TSpaceVecType>::operator-=(FourVector const& b) {
timeLike_ -= b.timeLike_; timeLike_ -= b.timeLike_;
spaceLike_ -= b.spaceLike_; spaceLike_ -= b.spaceLike_;
return *this; return *this;
} }
template <typename TTimeType, typename TSpaceVecType> template <typename TTimeType, typename TSpaceVecType>
inline FourVector<TTimeType, TSpaceVecType>& FourVector<TTimeType, TSpaceVecType>:: inline FourVector<TTimeType, TSpaceVecType>&
operator*=(double const b) { FourVector<TTimeType, TSpaceVecType>::operator*=(double const b) {
timeLike_ *= b; timeLike_ *= b;
spaceLike_ *= b; spaceLike_ *= b;
return *this; return *this;
} }
template <typename TTimeType, typename TSpaceVecType> template <typename TTimeType, typename TSpaceVecType>
inline FourVector<TTimeType, TSpaceVecType>& FourVector<TTimeType, TSpaceVecType>:: inline FourVector<TTimeType, TSpaceVecType>&
operator/=(double const b) { FourVector<TTimeType, TSpaceVecType>::operator/=(double const b) {
timeLike_ /= b; timeLike_ /= b;
spaceLike_.getComponents() /= b; spaceLike_.getComponents() /= b;
return *this; return *this;
} }
template <typename TTimeType, typename TSpaceVecType> template <typename TTimeType, typename TSpaceVecType>
inline FourVector<TTimeType, TSpaceVecType>& FourVector<TTimeType, TSpaceVecType>:: inline FourVector<TTimeType, TSpaceVecType>&
operator/(double const b) { FourVector<TTimeType, TSpaceVecType>::operator/(double const b) {
*this /= b; *this /= b;
return *this; return *this;
} }
template <typename TTimeType, typename TSpaceVecType> template <typename TTimeType, typename TSpaceVecType>
inline typename FourVector<TTimeType, TSpaceVecType>::norm_type inline typename FourVector<TTimeType, TSpaceVecType>::norm_type
FourVector<TTimeType, TSpaceVecType>::operator*(FourVector const& b) { FourVector<TTimeType, TSpaceVecType>::operator*(FourVector const& b) {
if constexpr (std::is_same<time_type, decltype(std::declval<space_type>() / meter * if constexpr (std::is_same<time_type, decltype(std::declval<space_type>() / meter *
second)>::value) second)>::value)
return timeLike_ * b.timeLike_ * constants::cSquared - spaceLike_.norm(); return timeLike_ * b.timeLike_ * constants::cSquared - spaceLike_.norm();
......
...@@ -17,8 +17,8 @@ ...@@ -17,8 +17,8 @@
namespace corsika { namespace corsika {
template <typename TDimension> template <typename TDimension>
inline typename QuantityVector<TDimension>::quantity_type QuantityVector<TDimension>:: inline typename QuantityVector<TDimension>::quantity_type
operator[](size_t const index) const { QuantityVector<TDimension>::operator[](size_t const index) const {
return quantity_type(phys::units::detail::magnitude_tag, eigenVector_[index]); return quantity_type(phys::units::detail::magnitude_tag, eigenVector_[index]);
} }
......
...@@ -20,7 +20,7 @@ namespace corsika { ...@@ -20,7 +20,7 @@ namespace corsika {
template <class AConstIterator, class BConstIterator> template <class AConstIterator, class BConstIterator>
inline typename WeightProviderIterator<AConstIterator, BConstIterator>::value_type inline typename WeightProviderIterator<AConstIterator, BConstIterator>::value_type
WeightProviderIterator<AConstIterator, BConstIterator>::operator*() const { WeightProviderIterator<AConstIterator, BConstIterator>::operator*() const {
return ((*aIter_) * (*bIter_)).magnitude(); return ((*aIter_) * (*bIter_)).magnitude();
} }
......
...@@ -13,12 +13,11 @@ ...@@ -13,12 +13,11 @@
namespace corsika { namespace corsika {
inline TimeCut::TimeCut(const TimeType time) inline TimeCut::TimeCut(const TimeType time)
: time_(time) {} : time_(time) {}
template <typename Particle, typename Track> template <typename Particle, typename Track>
inline ProcessReturn TimeCut::doContinuous( inline ProcessReturn TimeCut::doContinuous(Particle& particle, Track const&,
Particle& particle, Track const&, bool const) {
bool const) {
CORSIKA_LOG_TRACE("TimeCut::doContinuous"); CORSIKA_LOG_TRACE("TimeCut::doContinuous");
if (particle.getTime() >= time_) { if (particle.getTime() >= time_) {
CORSIKA_LOG_TRACE("stopping continuous process"); CORSIKA_LOG_TRACE("stopping continuous process");
......
...@@ -56,8 +56,8 @@ namespace corsika::proposal { ...@@ -56,8 +56,8 @@ namespace corsika::proposal {
PROPOSAL::InterpolationSettings::TABLES_PATH = corsika_data("PROPOSAL").c_str(); PROPOSAL::InterpolationSettings::TABLES_PATH = corsika_data("PROPOSAL").c_str();
} }
inline size_t ProposalProcessBase::hash::operator()(const calc_key_t& p) const inline size_t ProposalProcessBase::hash::operator()(
noexcept { const calc_key_t& p) const noexcept {
return p.first ^ std::hash<Code>{}(p.second); return p.first ^ std::hash<Code>{}(p.second);
} }
......
This diff is collapsed.
...@@ -11,106 +11,116 @@ ...@@ -11,106 +11,116 @@
namespace corsika { namespace corsika {
template <typename TAntennaCollection, typename TRadioImpl, typename TPropagator> template <typename TAntennaCollection, typename TRadioImpl, typename TPropagator>
inline TRadioImpl& RadioProcess<TAntennaCollection, TRadioImpl, TPropagator>::implementation() inline TRadioImpl&
{ return static_cast<TRadioImpl&>(*this); } 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
{ template <typename TAntennaCollection, typename TRadioImpl, typename TPropagator>
return static_cast<TRadioImpl const&>(*this); 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, typename Track>
inline ProcessReturn RadioProcess<TAntennaCollection, TRadioImpl,
TPropagator>::doContinuous(const Particle& particle,
const Track& track,
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_{particle.getPID()};
if ((particleID_ == Code::Electron) || (particleID_ == Code::Positron)) {
CORSIKA_LOG_DEBUG("Particle for radio calculation: {} ", particleID_);
return this->implementation().simulate(particle, track);
} else {
CORSIKA_LOG_DEBUG("Particle {} is irrelevant for radio", particleID_);
return ProcessReturn::Ok;
} }
//}
template <typename TAntennaCollection, typename TRadioImpl, typename TPropagator> }
template <typename... TArgs>
inline RadioProcess<TAntennaCollection, TRadioImpl, TPropagator>::RadioProcess( template <typename TAntennaCollection, typename TRadioImpl, typename TPropagator>
TAntennaCollection &antennas, TArgs &&...args) template <typename Particle, typename Track>
: antennas_(antennas) inline LengthType
, propagator_(args...) {} RadioProcess<TAntennaCollection, TRadioImpl, TPropagator>::getMaxStepLength(
const Particle& vParticle, const Track& vTrack) const {
template <typename TAntennaCollection, typename TRadioImpl, typename TPropagator> return meter * std::numeric_limits<double>::infinity();
template <typename Particle, typename Track> }
inline ProcessReturn RadioProcess<TAntennaCollection, TRadioImpl, TPropagator>::doContinuous(
const Particle &particle, const Track &track, const bool) { template <typename TAntennaCollection, typename TRadioImpl, typename TPropagator>
// we want the following particles: inline void RadioProcess<TAntennaCollection, TRadioImpl, TPropagator>::startOfLibrary(
// Code::Electron & Code::Positron & Code::Gamma const boost::filesystem::path& directory) {
// we wrap Simulate() in doContinuous as the plan is to add particle level // loop over every antenna and set the initial path
// filtering or thinning for calculation of the radio emission. This is // this also writes the time-bins to disk.
// important for controlling the runtime of radio (by ignoring particles for (auto& antenna : antennas_.getAntennas()) {
// that aren't going to contribute i.e. heavy hadrons) antenna.startOfLibrary(directory, this->implementation().algorithm);
// if (valid(particle, track)) {
auto const particleID_ {particle.getPID()};
if ((particleID_ == Code::Electron) || (particleID_ == Code::Positron)) {
CORSIKA_LOG_DEBUG("Particle for radio calculation: {} ", particleID_);
return this->implementation().simulate(particle, track);
} 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> template <typename TAntennaCollection, typename TRadioImpl, typename TPropagator>
inline LengthType RadioProcess<TAntennaCollection, TRadioImpl, TPropagator>::getMaxStepLength( inline void RadioProcess<TAntennaCollection, TRadioImpl, TPropagator>::endOfShower(
const Particle &vParticle, const Track &vTrack) const { const unsigned int) {
return meter * std::numeric_limits<double>::infinity();
// 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.sample_rate_ * 1_s);
antenna.reset();
} }
template <typename TAntennaCollection, typename TRadioImpl, typename TPropagator> // increment our event counter
inline void RadioProcess<TAntennaCollection, TRadioImpl, TPropagator>::startOfLibrary( event_++;
const boost::filesystem::path &directory) { }
// loop over every antenna and set the initial path template <typename TAntennaCollection, typename TRadioImpl, typename TPropagator>
// this also writes the time-bins to disk. inline YAML::Node RadioProcess<TAntennaCollection, TRadioImpl, TPropagator>::getConfig()
for (auto& antenna : antennas_.getAntennas()) { antenna.startOfLibrary(directory,this->implementation().algorithm);} const {
}
// top-level YAML node
template <typename TAntennaCollection, typename TRadioImpl, typename TPropagator> YAML::Node config;
inline void RadioProcess<TAntennaCollection, TRadioImpl, TPropagator>::endOfShower(const unsigned int) {
// fill in some basics
// loop over every antenna and instruct them to config["type"] = "RadioProcess";
// flush data to disk, and then reset the antenna config["algorithm"] = this->implementation().algorithm;
// before the next event config["units"]["time"] = "ns";
for (auto& antenna : antennas_.getAntennas()) { config["units"]["frequency"] = "GHz";
antenna.endOfShower(event_, this->implementation().algorithm, antenna.sample_rate_*1_s); config["units"]["electric field"] = "V/m";
antenna.reset(); config["units"]["distance"] = "m";
}
for (auto& antenna : antennas_.getAntennas()) {
// increment our event counter // get the name/location of this antenna
event_++; 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);
} }
template <typename TAntennaCollection, typename TRadioImpl, typename TPropagator> return config;
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;
}
} // namespace corsika } // namespace corsika
\ No newline at end of file
This diff is collapsed.
...@@ -13,96 +13,106 @@ ...@@ -13,96 +13,106 @@
namespace corsika { namespace corsika {
template <typename TAntennaImpl> template <typename TAntennaImpl>
inline Antenna<TAntennaImpl>::Antenna(std::string const& name, Point const& location) inline Antenna<TAntennaImpl>::Antenna(std::string const& name, Point const& location)
: name_(name) : name_(name)
, location_(location){}; , location_(location){};
template <typename TAntennaImpl> template <typename TAntennaImpl>
inline Point const& Antenna<TAntennaImpl>::getLocation() const { return location_; } inline Point const& Antenna<TAntennaImpl>::getLocation() const {
return location_;
}
template <typename TAntennaImpl> template <typename TAntennaImpl>
inline std::string const& Antenna<TAntennaImpl>::getName() const { return name_; } inline std::string const& Antenna<TAntennaImpl>::getName() const {
return name_;
}
template <typename TAntennaImpl> template <typename TAntennaImpl>
inline void Antenna<TAntennaImpl>::startOfLibrary(const boost::filesystem::path &directory, inline void Antenna<TAntennaImpl>::startOfLibrary(
const std::string radioImplementation) { const boost::filesystem::path& directory, const std::string radioImplementation) {
// calculate and save our filename // calculate and save our filename
filename_ = (directory / this->getName()).string() + ".npz"; filename_ = (directory / this->getName()).string() + ".npz";
// get the axis labels for this antenna and write the first row. // get the axis labels for this antenna and write the first row.
std::vector<long double> axis = this->implementation().getAxis(); std::vector<long double> axis = this->implementation().getAxis();
// check for the axis name // check for the axis name
std::string label = "Unknown"; std::string label = "Unknown";
if constexpr (TAntennaImpl::is_time_domain) { if constexpr (TAntennaImpl::is_time_domain) {
label = "Time"; label = "Time";
} else if constexpr (TAntennaImpl::is_freq_domain) { } else if constexpr (TAntennaImpl::is_freq_domain) {
label = "Frequency"; label = "Frequency";
} }
if (radioImplementation == "ZHS" && TAntennaImpl::is_time_domain) { if (radioImplementation == "ZHS" && TAntennaImpl::is_time_domain) {
for (size_t i=0; i<axis.size()-1;i++) for (size_t i = 0; i < axis.size() - 1; i++) {
{ axis.at(i) = (axis.at(i + 1) + axis.at(i)) / 2.;
axis.at(i) = (axis.at(i+1)+axis.at(i))/2.; }
} // explicitly convert the arrays to the needed type for cnpy
// explicitly convert the arrays to the needed type for cnpy axis.pop_back();
axis.pop_back(); long double const* raw_data = axis.data();
long double const* raw_data = axis.data(); std::vector<size_t> N = {
std::vector<size_t> N = {axis.size()}; // cnpy needs a vector here -- this should be axis.size() - 1 -- axis.size()}; // cnpy needs a vector here -- this should be axis.size() - 1 --
// write the labels to the first row of the NumPy file // write the labels to the first row of the NumPy file
cnpy::npz_save(filename_, label, raw_data, N, "w"); cnpy::npz_save(filename_, label, raw_data, N, "w");
} else { } else {
// explicitly convert the arrays to the needed type for cnpy // explicitly convert the arrays to the needed type for cnpy
long double const* raw_data = axis.data(); long double const* raw_data = axis.data();
std::vector<size_t> N = {axis.size()}; // cnpy needs a vector here std::vector<size_t> N = {axis.size()}; // cnpy needs a vector here
// write the labels to the first row of the NumPy file // write the labels to the first row of the NumPy file
cnpy::npz_save(filename_, label, raw_data, N, "w"); cnpy::npz_save(filename_, label, raw_data, N, "w");
}
} }
}
template <typename TAntennaImpl> template <typename TAntennaImpl>
inline void Antenna<TAntennaImpl>::endOfShower(const int event, const std::string radioImplementation, inline void Antenna<TAntennaImpl>::endOfShower(const int event,
const double sampleRate) { const std::string radioImplementation,
const double sampleRate) {
// get the copy of the waveform data for this event // get the copy of the waveform data for this event
// we transpose it so that we can match dimensions with the // we transpose it so that we can match dimensions with the
// time array that is already in the output file // time array that is already in the output file
std::vector<double> dataX = this->implementation().getDataX(); std::vector<double> dataX = this->implementation().getDataX();
std::vector<double> dataY = this->implementation().getDataY(); std::vector<double> dataY = this->implementation().getDataY();
std::vector<double> dataZ = this->implementation().getDataZ(); std::vector<double> dataZ = this->implementation().getDataZ();
if (radioImplementation == "ZHS") { if (radioImplementation == "ZHS") {
std::vector<double> electricFieldX (dataX.size() - 1, 0); //num_bins_, std::vector<double>(3, 0) std::vector<double> electricFieldX(dataX.size() - 1,
std::vector<double> electricFieldY (dataY.size() - 1, 0); 0); // num_bins_, std::vector<double>(3, 0)
std::vector<double> electricFieldZ (dataZ.size() - 1, 0); std::vector<double> electricFieldY(dataY.size() - 1, 0);
for (size_t i = 0; i < electricFieldX.size(); i++) 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; electricFieldX.at(i) = -(dataX.at(i + 1) - dataX.at(i)) * sampleRate;
electricFieldY.at(i) = -(dataY.at(i+1)-dataY.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; electricFieldZ.at(i) = -(dataZ.at(i + 1) - dataZ.at(i)) * sampleRate;
} }
// cnpy needs a vector for the shape // cnpy needs a vector for the shape
std::vector<size_t> shapeX = {electricFieldX.size()}; std::vector<size_t> shapeX = {electricFieldX.size()};
std::vector<size_t> shapeY = {electricFieldY.size()}; std::vector<size_t> shapeY = {electricFieldY.size()};
std::vector<size_t> shapeZ = {electricFieldZ.size()}; std::vector<size_t> shapeZ = {electricFieldZ.size()};
cnpy::npz_save(filename_, std::to_string(event) + "X", electricFieldX.data(), shapeX, "a"); cnpy::npz_save(filename_, std::to_string(event) + "X", electricFieldX.data(),
cnpy::npz_save(filename_, std::to_string(event) + "Y", electricFieldY.data(), shapeY, "a"); shapeX, "a");
cnpy::npz_save(filename_, std::to_string(event) + "Z", electricFieldZ.data(), shapeZ, "a"); cnpy::npz_save(filename_, std::to_string(event) + "Y", electricFieldY.data(),
} else { shapeY, "a");
// cnpy needs a vector for the shape cnpy::npz_save(filename_, std::to_string(event) + "Z", electricFieldZ.data(),
std::vector<size_t> shapeX = {dataX.size()}; shapeZ, "a");
std::vector<size_t> shapeY = {dataY.size()}; } else {
std::vector<size_t> shapeZ = {dataZ.size()}; // cnpy needs a vector for the shape
// and write this event to the .npz archive std::vector<size_t> shapeX = {dataX.size()};
cnpy::npz_save(filename_, std::to_string(event) + "X", dataX.data(), shapeX, "a"); std::vector<size_t> shapeY = {dataY.size()};
cnpy::npz_save(filename_, std::to_string(event) + "Y", dataY.data(), shapeY, "a"); std::vector<size_t> shapeZ = {dataZ.size()};
cnpy::npz_save(filename_, std::to_string(event) + "Z", dataZ.data(), shapeZ, "a"); // and write this event to the .npz archive
} cnpy::npz_save(filename_, std::to_string(event) + "X", dataX.data(), shapeX, "a");
cnpy::npz_save(filename_, std::to_string(event) + "Y", dataY.data(), shapeY, "a");
cnpy::npz_save(filename_, std::to_string(event) + "Z", dataZ.data(), shapeZ, "a");
} }
}
template <typename TAntennaImpl> template <typename TAntennaImpl>
inline TAntennaImpl& Antenna<TAntennaImpl>::implementation() { return static_cast<TAntennaImpl&>(*this); } inline TAntennaImpl& Antenna<TAntennaImpl>::implementation() {
return static_cast<TAntennaImpl&>(*this);
}
} // namespace corsika } // namespace corsika
\ No newline at end of file
...@@ -13,108 +13,122 @@ ...@@ -13,108 +13,122 @@
namespace corsika { namespace corsika {
inline TimeDomainAntenna::TimeDomainAntenna(const std::string &name, const Point &location, inline TimeDomainAntenna::TimeDomainAntenna(const std::string& name,
const TimeType &start_time, const TimeType &duration, const Point& location,
const InverseTimeType &sample_rate, const TimeType &ground_hit_time) const TimeType& start_time,
: Antenna(name, location) const TimeType& duration,
, start_time_(start_time) const InverseTimeType& sample_rate,
, duration_(duration) const TimeType& ground_hit_time)
, sample_rate_(sample_rate) : Antenna(name, location)
, ground_hit_time_(ground_hit_time) , start_time_(start_time)
, num_bins_(static_cast<std::size_t>(duration * sample_rate + 1.5l)) , duration_(duration)
, waveformEX_(num_bins_, 0) , sample_rate_(sample_rate)
, waveformEY_(num_bins_, 0) , ground_hit_time_(ground_hit_time)
, waveformEZ_(num_bins_, 0) {}; , num_bins_(static_cast<std::size_t>(duration * sample_rate + 1.5l))
, waveformEX_(num_bins_, 0)
inline void TimeDomainAntenna::receive(const TimeType time, const Vector<dimensionless_d> &receive_vector, , waveformEY_(num_bins_, 0)
const ElectricFieldVector &efield) { , waveformEZ_(num_bins_, 0){};
if (time < start_time_ || time > (start_time_ + duration_)) { inline void TimeDomainAntenna::receive(const TimeType time,
return; const Vector<dimensionless_d>& receive_vector,
} else { const ElectricFieldVector& efield) {
// figure out the correct timebin to store the E-field value.
// NOTE: static cast is implicitly flooring if (time < start_time_ || time > (start_time_ + duration_)) {
auto timebin_{static_cast<std::size_t>(std::floor((time - start_time_) * sample_rate_ + 0.5l))}; return;
// CORSIKA_LOG_DEBUG("Timebin: {}", timebin_); } else {
// figure out the correct timebin to store the E-field value.
// ToDO: ask explicitly for a CS and use that specific on for writing the output // NOTE: static cast is implicitly flooring
auto timebin_{static_cast<std::size_t>(
// store the x,y,z electric field components. std::floor((time - start_time_) * sample_rate_ + 0.5l))};
waveformEX_.at(timebin_) += (efield.getComponents().getX() / (1_V / 1_m)); // CORSIKA_LOG_DEBUG("Timebin: {}", timebin_);
waveformEY_.at(timebin_) += (efield.getComponents().getY() / (1_V / 1_m));
waveformEZ_.at(timebin_) += (efield.getComponents().getZ() / (1_V / 1_m)); // ToDO: ask explicitly for a CS and use that specific on for writing the output
// TODO: Check how they are stored in memory, row-wise or column-wise?
} // store the x,y,z electric field components.
waveformEX_.at(timebin_) += (efield.getComponents().getX() / (1_V / 1_m));
waveformEY_.at(timebin_) += (efield.getComponents().getY() / (1_V / 1_m));
waveformEZ_.at(timebin_) += (efield.getComponents().getZ() / (1_V / 1_m));
// TODO: Check how they are stored in memory, row-wise or column-wise?
} }
}
inline void TimeDomainAntenna::receive(const TimeType time, const Vector<dimensionless_d> &receive_vector,
const VectorPotential &vectorP) { inline void TimeDomainAntenna::receive(const TimeType time,
const Vector<dimensionless_d>& receive_vector,
if (time < start_time_ || time > (start_time_ + duration_)) { const VectorPotential& vectorP) {
return;
} else { if (time < start_time_ || time > (start_time_ + duration_)) {
// figure out the correct timebin to store the E-field value. return;
// NOTE: static cast is implicitly flooring } else {
auto timebin_{static_cast<std::size_t>(std::floor((time - start_time_) * sample_rate_ + 0.5l))}; // figure out the correct timebin to store the E-field value.
// CORSIKA_LOG_DEBUG("Timebin: {}", timebin_); // NOTE: static cast is implicitly flooring
auto timebin_{static_cast<std::size_t>(
// ToDO: ask explicitly for a CS and use that specific on for writing the output std::floor((time - start_time_) * sample_rate_ + 0.5l))};
// CORSIKA_LOG_DEBUG("Timebin: {}", timebin_);
// store the x,y,z electric field components.
waveformEX_.at(timebin_) += (vectorP.getComponents().getX() / (1_V * 1_s / 1_m)); // ToDO: ask explicitly for a CS and use that specific on for writing the output
waveformEY_.at(timebin_) += (vectorP.getComponents().getY() / (1_V * 1_s / 1_m));
waveformEZ_.at(timebin_) += (vectorP.getComponents().getZ() / (1_V * 1_s / 1_m)); // store the x,y,z electric field components.
// TODO: Check how they are stored in memory, row-wise or column-wise? waveformEX_.at(timebin_) += (vectorP.getComponents().getX() / (1_V * 1_s / 1_m));
} waveformEY_.at(timebin_) += (vectorP.getComponents().getY() / (1_V * 1_s / 1_m));
waveformEZ_.at(timebin_) += (vectorP.getComponents().getZ() / (1_V * 1_s / 1_m));
// TODO: Check how they are stored in memory, row-wise or column-wise?
} }
}
inline auto& TimeDomainAntenna::getDataX() const { return waveformEX_; } inline auto& TimeDomainAntenna::getDataX() const { return waveformEX_; }
inline auto& TimeDomainAntenna::getDataY() const { return waveformEY_; }
inline auto& TimeDomainAntenna::getDataZ() const { return waveformEZ_; } inline auto& TimeDomainAntenna::getDataY() const { return waveformEY_; }
inline auto TimeDomainAntenna::getAxis() const { inline auto& TimeDomainAntenna::getDataZ() const { return waveformEZ_; }
// create a 1-D xtensor to store time values so we can print them later. inline auto TimeDomainAntenna::getAxis() const {
std::vector<long double> times(num_bins_, 0);
// calculate the sample_period // create a 1-D xtensor to store time values so we can print them later.
auto sample_period{1 / sample_rate_}; std::vector<long double> times(num_bins_, 0);
// fill in every time-value // calculate the sample_period
// TODO: Vectorize this auto sample_period{1 / sample_rate_};
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; // fill in every time-value
// TODO: Vectorize this
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);
} }
inline auto TimeDomainAntenna::getWaveformX() const { return std::make_pair(getAxis(), waveformEX_); } return times;
}
inline auto TimeDomainAntenna::getWaveformY() const { return std::make_pair(getAxis(), waveformEY_); } inline auto TimeDomainAntenna::getWaveformX() const {
return std::make_pair(getAxis(), waveformEX_);
}
inline auto TimeDomainAntenna::getWaveformZ() const { return std::make_pair(getAxis(), waveformEZ_); } inline auto TimeDomainAntenna::getWaveformY() const {
return std::make_pair(getAxis(), waveformEY_);
}
inline void TimeDomainAntenna::reset() { inline auto TimeDomainAntenna::getWaveformZ() const {
std::fill(waveformEX_.begin(), waveformEX_.end(), 0); return std::make_pair(getAxis(), waveformEZ_);
std::fill(waveformEY_.begin(), waveformEY_.end(), 0); }
std::fill(waveformEZ_.begin(), waveformEZ_.end(), 0);
}
inline YAML::Node TimeDomainAntenna::getConfig() const { 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);
}
// top-level config inline YAML::Node TimeDomainAntenna::getConfig() const {
YAML::Node config;
config["type"] = "TimeDomainAntenna"; // top-level config
config["start_time"] = start_time_ / 1_ns; YAML::Node config;
config["duration"] = duration_ / 1_ns;
config["sample_rate"] = sample_rate_ / 1_GHz;
return 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 } // namespace corsika
\ No newline at end of file
...@@ -13,32 +13,30 @@ ...@@ -13,32 +13,30 @@
namespace corsika { namespace corsika {
template <typename TAntennaImpl> template <typename TAntennaImpl>
inline AntennaCollection::AntennaCollection { inline AntennaCollection::AntennaCollection {
std::vector<TAntennaImpl> antennas_; std::vector<TAntennaImpl> antennas_;
} }
template <typename TAntennaImpl> template <typename TAntennaImpl>
inline void AntennaCollection::addAntenna(TAntennaImpl const& antenna) { inline void AntennaCollection::addAntenna(TAntennaImpl const& antenna) {
antennas_.push_back(antenna); antennas_.push_back(antenna);
} }
template <typename TAntennaImpl> template <typename TAntennaImpl>
inline TAntennaImpl AntennaCollection::at(std::size_t const i) { inline TAntennaImpl AntennaCollection::at(std::size_t const i) {
antennas_.at(i); antennas_.at(i);
} }
inline int AntennaCollection::size() { inline int AntennaCollection::size() { return antennas_.size(); }
return antennas_.size();
} template <typename TAntennaImpl>
inline std::vector<TAntennaImpl>& AntennaCollection::getAntennas() {
template <typename TAntennaImpl> return antennas_;
inline std::vector<TAntennaImpl>& AntennaCollection::getAntennas() { }
return antennas_;
} inline void AntennaCollection::reset() {
std::for_each(antennas_.begin(), antennas_.end(), std::mem_fn(&TAntennaImpl::reset));
inline void AntennaCollection::reset() { }
std::for_each(antennas_.begin(), antennas_.end(), std::mem_fn(&TAntennaImpl::reset));
}
} // namespace corsika } // namespace corsika
...@@ -11,8 +11,8 @@ ...@@ -11,8 +11,8 @@
namespace corsika { namespace corsika {
template <typename TImpl, typename TEnvironment> template <typename TImpl, typename TEnvironment>
inline RadioPropagator<TImpl, TEnvironment>::RadioPropagator(const TEnvironment &env) inline RadioPropagator<TImpl, TEnvironment>::RadioPropagator(const TEnvironment& env)
: env_(env) {} : env_(env) {}
} // namespace corsika } // namespace corsika
...@@ -12,17 +12,18 @@ ...@@ -12,17 +12,18 @@
namespace corsika { namespace corsika {
inline SignalPath::SignalPath(const TimeType propagation_time, const double average_refractive_index, inline SignalPath::SignalPath(
const double refractive_index_source, const double refractive_index_destination, const TimeType propagation_time, const double average_refractive_index,
const Vector <dimensionless_d> emit, const Vector <dimensionless_d> receive, const double refractive_index_source, const double refractive_index_destination,
const LengthType R_distance, const std::deque<Point> &points) const Vector<dimensionless_d> emit, const Vector<dimensionless_d> receive,
: Path(points) const LengthType R_distance, const std::deque<Point>& points)
, propagation_time_(propagation_time) : Path(points)
, average_refractive_index_(average_refractive_index) , propagation_time_(propagation_time)
, refractive_index_source_(refractive_index_source) , average_refractive_index_(average_refractive_index)
, refractive_index_destination_(refractive_index_destination) , refractive_index_source_(refractive_index_source)
, emit_(emit) , refractive_index_destination_(refractive_index_destination)
, receive_(receive) , emit_(emit)
, R_distance_(R_distance) {} , receive_(receive)
, R_distance_(R_distance) {}
} // namespace corsika } // namespace corsika
\ No newline at end of file
...@@ -11,61 +11,60 @@ ...@@ -11,61 +11,60 @@
namespace corsika { namespace corsika {
template <typename TEnvironment> template <typename TEnvironment>
inline SimplePropagator<TEnvironment>::SimplePropagator(const TEnvironment &env) inline SimplePropagator<TEnvironment>::SimplePropagator(const TEnvironment& env)
: RadioPropagator<SimplePropagator, TEnvironment>(env){}; : RadioPropagator<SimplePropagator, TEnvironment>(env){};
template <typename TEnvironment>
inline typename SimplePropagator<TEnvironment>::SignalPathCollection
SimplePropagator<TEnvironment>::propagate(const Point& source, const Point& destination,
const LengthType stepsize) const {
template <typename TEnvironment> /**
inline typename SimplePropagator<TEnvironment>::SignalPathCollection SimplePropagator<TEnvironment>::propagate(const Point &source, const Point &destination, * This is the simplest case of straight propagator
const LengthType stepsize) const { * 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
* This is the simplest case of straight propagator auto emit_{(destination - source).normalized()};
* where no integration takes place. auto receive_{-emit_};
* 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 // the geometrical distance from the point of emission to an observer
auto emit_{(destination - source).normalized()}; auto distance_{(destination - source).getNorm()};
auto receive_{ - emit_};
// the geometrical distance from the point of emission to an observer // get the universe for this environment
auto distance_ {(destination - source).getNorm()}; auto const* const universe{Base::env_.getUniverse().get()};
// get the universe for this environment // the points that consist the signal path (source & destination).
auto const* const universe{Base::env_.getUniverse().get()}; std::deque<Point> points;
// the points that consist the signal path (source & destination). // store value of the refractive index at points.
std::deque<Point> points; std::vector<double> rindex;
rindex.reserve(2);
// store value of the refractive index at points. // get and store the refractive index of the first point 'source'.
std::vector<double> rindex; auto const* nodeSource{universe->getContainingNode(source)};
rindex.reserve(2); auto const ri_source{nodeSource->getModelProperties().getRefractiveIndex(source)};
rindex.push_back(ri_source);
points.push_back(source);
// get and store the refractive index of the first point 'source'. // add the refractive index of last point 'destination' and store it.
auto const* nodeSource{universe->getContainingNode(source)}; auto const* node{universe->getContainingNode(destination)};
auto const ri_source{nodeSource->getModelProperties().getRefractiveIndex(source)}; auto const ri_destination{node->getModelProperties().getRefractiveIndex(destination)};
rindex.push_back(ri_source); rindex.push_back(ri_destination);
points.push_back(source); points.push_back(destination);
// add the refractive index of last point 'destination' and store it. // compute the average refractive index.
auto const* node{universe->getContainingNode(destination)}; auto averageRefractiveIndex_ = (ri_source + ri_destination) / 2;
auto const ri_destination{node->getModelProperties().getRefractiveIndex(destination)};
rindex.push_back(ri_destination);
points.push_back(destination);
// compute the average refractive index. // compute the total time delay.
auto averageRefractiveIndex_ = (ri_source + ri_destination) / 2; TimeType time = averageRefractiveIndex_ * (distance_ / constants::c);
// compute the total time delay. return {SignalPath(time, averageRefractiveIndex_, ri_source, ri_destination, emit_,
TimeType time = averageRefractiveIndex_ * (distance_ / constants::c); receive_, distance_, points)};
} // END: propagate()
return {SignalPath(time, averageRefractiveIndex_, ri_source, ri_destination,
emit_, receive_, distance_, points)};
} // END: propagate()
} // namespace corsika } // namespace corsika
...@@ -11,153 +11,159 @@ ...@@ -11,153 +11,159 @@
namespace corsika { namespace corsika {
template <typename TEnvironment> template <typename TEnvironment>
// TODO: maybe the constructor doesn't take any arguments for the environment (?) // TODO: maybe the constructor doesn't take any arguments for the environment (?)
inline StraightPropagator<TEnvironment>::StraightPropagator(const TEnvironment &env) inline StraightPropagator<TEnvironment>::StraightPropagator(const TEnvironment& env)
: RadioPropagator<StraightPropagator, TEnvironment>(env){}; : RadioPropagator<StraightPropagator, TEnvironment>(env){};
template <typename TEnvironment> template <typename TEnvironment>
inline typename StraightPropagator<TEnvironment>::SignalPathCollection StraightPropagator<TEnvironment>::propagate( inline typename StraightPropagator<TEnvironment>::SignalPathCollection
const Point &source, const Point &destination, const LengthType stepsize) const { StraightPropagator<TEnvironment>::propagate(const Point& source,
const Point& destination,
/** const LengthType 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 * get the normalized (unit) vector from `source` to `destination'.
* so they are both called direction * 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 emit_{(destination - source).normalized()};
auto receive_{ - emit_}; // these are used for the direction of emission and reception of signal at the antenna
auto emit_{(destination - source).normalized()};
// the distance from the point of emission to an observer auto receive_{-emit_};
auto distance_ {(destination - source).getNorm()};
// the distance from the point of emission to an observer
try { auto distance_{(destination - source).getNorm()};
if (stepsize <= 0.5 * distance_) {
try {
// "step" is the direction vector with length `stepsize` if (stepsize <= 0.5 * distance_) {
auto step{emit_ * stepsize};
// "step" is the direction vector with length `stepsize`
// calculate the number of points (roughly) for the numerical integration auto step{emit_ * stepsize};
auto n_points{(destination - source).getNorm() / stepsize};
// calculate the number of points (roughly) for the numerical integration
// get the universe for this environment auto n_points{(destination - source).getNorm() / stepsize};
auto const* const universe{Base::env_.getUniverse().get()};
// get the universe for this environment
// the points that we build along the way for the numerical integration auto const* const universe{Base::env_.getUniverse().get()};
std::deque<Point> points;
// the points that we build along the way for the numerical integration
// store value of the refractive index at points std::deque<Point> points;
std::vector<double> rindex;
rindex.reserve(n_points); // store value of the refractive index at points
std::vector<double> rindex;
// get and store the refractive index of the first point 'source' rindex.reserve(n_points);
auto const* nodeSource{universe->getContainingNode(source)};
auto const ri_source{ // get and store the refractive index of the first point 'source'
nodeSource->getModelProperties().getRefractiveIndex(source)}; auto const* nodeSource{universe->getContainingNode(source)};
rindex.push_back(ri_source); auto const ri_source{nodeSource->getModelProperties().getRefractiveIndex(source)};
points.push_back(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 // loop from `source` to `destination` to store values before Simpson's rule.
for (auto point = source + step; // this loop skips the last point 'destination' and "misses" the extra point
(point - destination).getNorm() > 0.6 * stepsize; point = point + step) { for (auto point = source + step; (point - destination).getNorm() > 0.6 * stepsize;
point = point + step) {
// get the environment node at this specific 'point'
auto const* node{universe->getContainingNode(point)}; // get the environment node at this specific 'point'
auto const* node{universe->getContainingNode(point)};
// get the associated refractivity at 'point'
auto const refractive_index{ // get the associated refractivity at 'point'
node->getModelProperties().getRefractiveIndex(point)}; auto const refractive_index{
// auto const refractive_index{1.000327}; node->getModelProperties().getRefractiveIndex(point)};
rindex.push_back(refractive_index); // auto const refractive_index{1.000327};
rindex.push_back(refractive_index);
// add this 'point' to our deque collection
points.push_back(point); // 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* 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 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* 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* 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 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 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() // 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* 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 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* 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* 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 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 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 } // namespace corsika
...@@ -36,7 +36,8 @@ namespace corsika::constants { ...@@ -36,7 +36,8 @@ namespace corsika::constants {
constexpr quantity<electric_charge_d> e{Rep(1.6021766208e-19L) * coulomb}; constexpr quantity<electric_charge_d> e{Rep(1.6021766208e-19L) * coulomb};
// vacuum permittivity // vacuum permittivity
constexpr quantity<dimensions<-3, -1, 4, 2>> epsilonZero{Rep(8.8541878128e-12L) * farad / meter}; constexpr quantity<dimensions<-3, -1, 4, 2>> epsilonZero{Rep(8.8541878128e-12L) *
farad / meter};
// electronvolt // electronvolt
// constexpr quantity<hepenergy_d> eV{e / coulomb * joule}; // constexpr quantity<hepenergy_d> eV{e / coulomb * joule};
......
...@@ -99,11 +99,10 @@ namespace corsika::units::si { ...@@ -99,11 +99,10 @@ namespace corsika::units::si {
using ElectricFieldType = using ElectricFieldType =
phys::units::quantity<phys::units::dimensions<1, 1, -3, -1>, double>; phys::units::quantity<phys::units::dimensions<1, 1, -3, -1>, double>;
using VectorPotentialType = using VectorPotentialType =
phys::units::quantity<phys::units::dimensions<1, 1, -2, -1>, double>; phys::units::quantity<phys::units::dimensions<1, 1, -2, -1>, double>;
using MagneticFieldType = using MagneticFieldType =
phys::units::quantity<phys::units::dimensions<-1, 0, 0, 1>, double>; phys::units::quantity<phys::units::dimensions<-1, 0, 0, 1>, double>;
template <typename DimFrom, typename DimTo> template <typename DimFrom, typename DimTo>
auto constexpr conversion_factor_HEP_to_SI() { auto constexpr conversion_factor_HEP_to_SI() {
static_assert(DimFrom::dim1 == 0 && DimFrom::dim2 == 0 && DimFrom::dim3 == 0 && static_assert(DimFrom::dim1 == 0 && DimFrom::dim2 == 0 && DimFrom::dim3 == 0 &&
......
...@@ -15,74 +15,74 @@ ...@@ -15,74 +15,74 @@
namespace corsika { namespace corsika {
/*!
* A Point represents a point in position space. It is defined by its
* coordinates with respect to some CoordinateSystem.
*/
class Point : public BaseVector<length_d> {
public:
Point(CoordinateSystemPtr const& pCS, QuantityVector<length_d> const& pQVector)
: BaseVector<length_d>(pCS, pQVector) {}
Point(CoordinateSystemPtr const& cs, LengthType x, LengthType y, LengthType z)
: BaseVector<length_d>(cs, {x, y, z}) {}
/** \todo TODO: this should be private or protected, we don NOT want to expose numbers
* without reference to outside:
*/
inline QuantityVector<length_d> const& getCoordinates() const;
inline QuantityVector<length_d>& getCoordinates();
/**
this always returns a QuantityVector as triple
\returns A value type QuantityVector, since it may have to create a temporary
object to transform to pCS.
**/
inline QuantityVector<length_d> getCoordinates(CoordinateSystemPtr const& pCS) const;
/**
* this always returns a QuantityVector as triple
*
* \returns A reference type QuantityVector&, but be aware, the underlying class data
* is actually transformed to pCS, if needed. Thus, there may be an implicit call to
* \ref rebase.
**/
inline QuantityVector<length_d>& getCoordinates(CoordinateSystemPtr const& pCS);
/**
* \name access coordinate components
* \{
*
* Note, if you access components in a different CoordinateSystem
* pCS than the stored data, internally a temporary object will be
* created and destroyed each call. This can be avoided by using
* \ref rebase first.
**/
inline LengthType getX(CoordinateSystemPtr const& pCS) const;
inline LengthType getY(CoordinateSystemPtr const& pCS) const;
inline LengthType getZ(CoordinateSystemPtr const& pCS) const;
/** \} **/
/*! /*!
* A Point represents a point in position space. It is defined by its * transforms the Point into another CoordinateSystem by changing its
* coordinates with respect to some CoordinateSystem. * coordinates interally
*/ */
class Point : public BaseVector<length_d> { inline void rebase(CoordinateSystemPtr const& pCS);
public: inline Point operator+(Vector<length_d> const& pVec) const;
Point(CoordinateSystemPtr const& pCS, QuantityVector<length_d> const& pQVector)
: BaseVector<length_d>(pCS, pQVector) {} /*!
* returns the distance Vector between two points
Point(CoordinateSystemPtr const& cs, LengthType x, LengthType y, LengthType z)
: BaseVector<length_d>(cs, {x, y, z}) {}
/** \todo TODO: this should be private or protected, we don NOT want to expose numbers
* without reference to outside:
*/
inline QuantityVector<length_d> const& getCoordinates() const;
inline QuantityVector<length_d>& getCoordinates();
/**
this always returns a QuantityVector as triple
\returns A value type QuantityVector, since it may have to create a temporary
object to transform to pCS.
**/
inline QuantityVector<length_d> getCoordinates(CoordinateSystemPtr const& pCS) const;
/**
* this always returns a QuantityVector as triple
*
* \returns A reference type QuantityVector&, but be aware, the underlying class data
* is actually transformed to pCS, if needed. Thus, there may be an implicit call to
* \ref rebase.
**/
inline QuantityVector<length_d>& getCoordinates(CoordinateSystemPtr const& pCS);
/**
* \name access coordinate components
* \{
*
* Note, if you access components in a different CoordinateSystem
* pCS than the stored data, internally a temporary object will be
* created and destroyed each call. This can be avoided by using
* \ref rebase first.
**/
inline LengthType getX(CoordinateSystemPtr const& pCS) const;
inline LengthType getY(CoordinateSystemPtr const& pCS) const;
inline LengthType getZ(CoordinateSystemPtr const& pCS) const;
/** \} **/
/*!
* transforms the Point into another CoordinateSystem by changing its
* coordinates interally
*/
inline void rebase(CoordinateSystemPtr const& pCS);
inline Point operator+(Vector<length_d> const& pVec) const;
/*!
* returns the distance Vector between two points
*/
inline Vector<length_d> operator-(Point const& pB) const;
};
/*
* calculates the distance between two points
*/ */
inline LengthType distance(Point const& p1, Point const& p2); inline Vector<length_d> operator-(Point const& pB) const;
};
/*
* calculates the distance between two points
*/
inline LengthType distance(Point const& p1, Point const& p2);
} // namespace corsika } // namespace corsika
......
...@@ -27,18 +27,15 @@ namespace corsika { ...@@ -27,18 +27,15 @@ namespace corsika {
template <typename Particle, typename Track> template <typename Particle, typename Track>
ProcessReturn doContinuous( ProcessReturn doContinuous(
Particle& vParticle, Particle& vParticle, Track const& vTrajectory,
Track const& vTrajectory,
const bool limitFlag = false); // this is not used for TimeCut const bool limitFlag = false); // this is not used for TimeCut
template <typename Particle, typename Track> template <typename Particle, typename Track>
LengthType getMaxStepLength(Particle const&, LengthType getMaxStepLength(Particle const&, Track const&) {
Track const&) {
return meter * std::numeric_limits<double>::infinity(); return meter * std::numeric_limits<double>::infinity();
} }
private: private:
TimeType time_; TimeType time_;
}; };
} // namespace corsika } // namespace corsika
......
...@@ -16,49 +16,50 @@ ...@@ -16,49 +16,50 @@
namespace corsika { namespace corsika {
template <typename TRadioDetector, typename TPropagator> template <typename TRadioDetector, typename TPropagator>
class CoREAS final : public RadioProcess<TRadioDetector, CoREAS<TRadioDetector, TPropagator>, TPropagator> { class CoREAS final
: public RadioProcess<TRadioDetector, CoREAS<TRadioDetector, TPropagator>,
TPropagator> {
public: public:
using ElectricFieldVector = Vector<ElectricFieldType::dimension_type>; using ElectricFieldVector = Vector<ElectricFieldType::dimension_type>;
// an identifier for which algorithm was used // an identifier for which algorithm was used
static constexpr auto algorithm = "CoREAS"; 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...){};
/** /**
* Construct a new CoREAS instance. * Simulate the radio emission from a particle across a track.
* *
* This forwards the detector and other arguments to * This must be provided by the TRadioImpl.
* the RadioProcess parent. *
* * @param particle The current particle.
*/ * @param track The current track.
template <typename... TArgs> *
CoREAS(TRadioDetector& detector, TArgs&&... args) */
: RadioProcess<TRadioDetector, CoREAS, TPropagator>(detector, args...) {}; template <typename Particle, typename Track>
ProcessReturn simulate(Particle const& particle, Track const& track);
private:
int tinycounter_{0};
int trackcounter_{0};
int zhscounter_{0};
/** using Base =
* Simulate the radio emission from a particle across a track. RadioProcess<TRadioDetector, CoREAS<TRadioDetector, TPropagator>, TPropagator>;
* using Base::antennas_;
* This must be provided by the TRadioImpl.
*
* @param particle The current particle.
* @param track The current track.
*
*/
template <typename Particle, typename Track>
ProcessReturn simulate(Particle const& particle, Track const& track);
private: }; // end of class CoREAS
int tinycounter_ {0};
int trackcounter_ {0};
int zhscounter_ {0};
using Base = RadioProcess<TRadioDetector, CoREAS<TRadioDetector, TPropagator>, TPropagator>;
using Base::antennas_;
}; // end of class CoREAS
} // namespace corsika } // namespace corsika
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment