IAP GITLAB

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

Antennas are tested and work. TODO: fix the foor loop for antennas.

parent 84f85357
No related branches found
No related tags found
1 merge request!329Radio interface
Hide whitespace changes
Inline Side-by-side
Showing
with 591 additions and 178 deletions
corsika/modules/radio/CoREAS.hpp
+ 294
− 120
View file @ 35158b33
......@@ -11,6 +11,7 @@
#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 {
......@@ -50,152 +51,325 @@ namespace corsika {
template <typename Particle, typename Track>
ProcessReturn simulate(Particle& particle, Track const& track) const {
// set threshold for application of ZHS-like approximation
const double approxThreshold_ {1.0e-3};
//get global simulation time for that track. (This is my best guess for now)
auto startTime_ {particle.getTime()
- track.getDuration()}; // time at start point of track.
auto endTime_ {particle.getTime()}; // time at end point of track.
// get the global simulation time for that track. (best guess for now)
auto startTime_ {particle.getTime() - track.getDuration()}; // time at the start point of the track hopefully.
auto endTime_ {particle.getTime()};
// beta is defined as velocity / speed of light
auto startBeta_ {track.getVelocity(0) / constants::c};
auto endBeta_ {track.getVelocity(1) / constants::c};
//TODO: check if they are different!!! They shouldn't!!!
// auto startBeta_ {(track.getPosition(0) / (particle.getTime()
// - track.getDuration()) ) / constants::c};
// auto endBeta_ {(track.getPosition(1) / particle.getTime() ) / constants::c};
// calculate gamma factor using beta (the proper way would be with energy over mass)
auto startGamma_ {1. / sqrt(1. - (startBeta_ * startBeta_))};
auto endGamma_ {1. / sqrt(1. - (endBeta_ * endBeta_))};
// gamma factor is calculated using beta
// auto startGamma_ {1. / sqrt(1. - (startBeta_ * startBeta_))};
// auto endGamma_ {1. / sqrt(1. - (endBeta_ * endBeta_))};
// get start and end position of the track
auto startPoint_ {track.getPosition(0)};
auto endPoint_ {track.getPosition(1)};
// beta is velocity / speed of light. Start & end should be the same!
auto beta_ {((endPoint_ - startPoint_) / (endTime_ - startTime_)).normalized()};
// get particle charge
auto const charge_ {get_charge(particle.getPID())};
// we loop over each antenna in the collection
for (auto& antenna : detector_.getAntennas()) {
// set threshold for application of ZHS-like approximation.
const double approxThreshold_ {1.0e-3};
ElectricFieldVector EV1_ {0_V / 0_m};
ElectricFieldVector EV2_ {0_V / 0_m};
ElectricFieldVector EV3_ {0_V / 0_m};
// we loop over each antenna in the collection.
for (auto& antenna : detector_.getAntennas()) {
auto startPointReceiveTime_ {0_ns};
auto endPointReceiveTime_ {0_ns};
std::vector<ElectricFieldVector> EVstart_;
std::vector<ElectricFieldVector> EVend_;
std::vector<TimeType> startTTimes_;
std::vector<TimeType> endTTimes_;
std::vector<double> preDoppler;
std::vector<double> postDoppler;
std::vector<QuantityVector<dimensionless_d>> ReceiveVectorsStart_;
std::vector<QuantityVector<dimensionless_d>> ReceiveVectorsEnd_;
// get the Path (path1) from the start "endpoint" to the antenna.
// This is a SignalPathCollection
auto paths1{this->propagator_.propagate(startPoint_, antenna.getLocation())};
// This is a Signal Path Collection
auto paths1 {this->propagator_.propagate(startPoint_, antenna.getLocation(), 1_nm)}; // TODO: Need to add the stepsize to .propagate()!!!!
// set postDoppler to approximate threshold in advance and check after loop
// if we need to perform ZHS-like approximation
auto preDoppler_ {approxThreshold_};
// now loop over the paths for startpoint that we got above
for (auto const& path : paths1) {
preDoppler_ = 1. - path.average_refractive_index_ *
startBeta_ * path.emit_;
// calculate preDoppler factor
double preDoppler_{1. - path.average_refractive_index_ *
beta_ * path.emit_};
// store it to the preDoppler std::vector for later comparisons
preDoppler.push_back(preDoppler_);
// calculate receive times for startpoint
auto startPointReceiveTime_ {path.total_time_ + startTime_};
// store it to startTTimes_ std::vector for later use
startTTimes_.push_back(startPointReceiveTime_);
if (preDoppler_ > approxThreshold_) {
startPointReceiveTime_ = path.total_time_ +
startTime_ ; // might do it on the fly
// store the receive unit vector
ReceiveVectorsStart_.push_back(path.receive_);
// CoREAS calculation -> get ElectricFieldVector1
EV1_= (charge_ / constants::c) *
path.receive_.cross(path.receive_.cross(startBeta_)) /
(path.R_distance_ * preDoppler_);
// calculate electric field vector for startpoint
auto EV1_= (charge_ / constants::c) *
path.receive_.cross(path.receive_.cross(beta_)) /
(path.R_distance_ * preDoppler_);
// pass it to the antenna
antenna.receive(startPointReceiveTime_, path.receive_, EV1_);
} else {
continue;
} // END preDoppler check
// store it to EVstart_ std::vector for later use
EVstart_.push_back(EV1_);
} // END: loop over paths for startpoint
} // End of looping over paths1
// get the Path (path2) from the end "endpoint" to the antenna.
// This is a SignalPathCollection
auto paths2{this->propagator_.propagate(endPoint_, antenna.getLocation())};
auto paths2 {this->propagator_.propagate(endPoint_, antenna.getLocation())};
// set postDoppler to approximate threshold in advance and check after loop
// if we need to perform ZHS-like approximation
auto postDoppler_ {approxThreshold_};
// now loop over the paths for endpoint that we got above
for (auto const& path : paths2) {
postDoppler_ = 1. - path.average_refractive_index_ *
endBeta_ * path.emit_;
if (preDoppler_ > approxThreshold_) {
auto endPointReceiveTime_ = path.total_time + endTime_ ; // might do it on the fly
//CoREAS calculation -> get ElectricFieldVector2
EV2_ = (charge_ / constants::c) *
path.receive_.cross(path.receive_.cross(endBeta_)) /
(path.R_distance_ * postDoppler_);
// pass it to the antenna
antenna.receive(endPointReceiveTime_, path.receive_, EV2_);
} else {
continue;
} // END postDoppler check
} // END: loop over paths for endpoint
// perform ZHS-like calculation close to Cherenkov angle
if (fabs(preDoppler_) <= approxThreshold_ || fabs(postDoppler_) <= approxThreshold_) {
// get global simulation time for the middle point of that track. (This is my best guess for now)
auto midTime_{particle.getTime() - (track.getDuration() / 2)};
// beta is defined as velocity / speed of light
auto midBeta_{track.getVelocity(0.5) / constants::c};
// auto midBeta_ {(track.getPosition(0.5) / (particle.getTime()
// - (track.getDuration() / 2) ) / constants::c};
// calculate gamma factor using beta (the proper way would be with energy over mass)
// auto midGamma_ {1. / sqrt(1. - (midBeta_ * midBeta_))};
// get start and end position of the track
auto midPoint_{track.getPosition(0.5)};
// get the Path (path3) from the middle "endpoint" to the antenna.
// This is a SignalPathCollection
auto paths3{this->propagator_.propagate(midPoint_, antenna.getLocation())};
// now loop over the paths for endpoint that we got above
for (auto const& path : paths3) {
midDoppler_ = 1. - path.average_refractive_index_ * midBeta_ * path.emit_;
auto const midPointReceiveTime_{path.total_time_ +
midTime_}; // might do it on the fly
// CoREAS calculation -> get ElectricFieldVector3 for "midPoint"
EV3_ = (charge_ / constants::c) *
path.receive_.cross(path.receive_.cross(midBeta_)) /
(path.R_distance_ * midDoppler_);
// pass it to the antenna but first check if "start" or "end" are double counted!!!
if (fabs(preDoppler_) > approxThreshold_ &&
fabs(postDoppler_) <= approxThreshold_) {
antenna.receive(startPointReceiveTime_, path.receive_, -EV1_);
antenna.receive(midPointReceiveTime_, path.receive_, EV3_);
} else if (fabs(preDoppler_) <= approxThreshold_ &&
fabs(postDoppler_) > approxThreshold_) {
antenna.receive(endPointReceiveTime_, path.receive_, -EV2_);
antenna.receive(midPointReceiveTime_, path.receive_, EV3_);
} else {
antenna.receive(midPointReceiveTime_, path.receive_, EV3_);
} // end deleting double-counted values
}
} // end of ZHS-like approximation
} // END: loop over antennas
}
double postDoppler_{1. - path.average_refractive_index_ *
beta_ * path.emit_}; // maybe this is path.receive_ (?)
// store it to the postDoppler std::vector for later comparisons
postDoppler.push_back(postDoppler_);
// calculate receive times for endpoint
auto endPointReceiveTime_ {path.total_time_ + endTime_};
// store it to endTTimes_ std::vector for later use
endTTimes_.push_back(endPointReceiveTime_);
// store the receive unit vector
ReceiveVectorsEnd_.push_back(path.receive_);
// calculate electric field vector for endpoint
auto EV2_= (charge_ / constants::c) *
path.receive_.cross(path.receive_.cross(beta_)) /
(path.R_distance_ * postDoppler_);
// store it to EVstart_ std::vector for later use
EVend_.push_back(EV2_);
} // End of looping over paths2
// start doing comparisons for preDoppler and postDoppler
// first check that start and end paths have the same number of paths
if (EVstart_.size() == EVend_.size()) {
// use this to access different elements of std::vectors
std::size_t index = 0;
for (auto& preDoppler__ : preDoppler) {
// redistribute contributions over time scale defined by the observation time resolution
// this is to make sure that "start" and "end" won't end up in the same bin (xtensor)!!
if ((preDoppler__ < 1.e-9) || (postDoppler.at(index) < 1.e-9)) {
auto gridResolution_ {antenna.duration_};
auto deltaT_ { endTTimes_.at(index) - startTTimes_.at(index) };
if (fabs(deltaT_) < gridResolution_) {
EVstart_.at(index) = EVstart_.at(index) * fabs(deltaT_ / gridResolution_);
EVend_.at(index) = EVend_.at(index) * fabs(deltaT_ / gridResolution_);
const long startBin = static_cast<long>(floor(startTTimes_.at(index)/gridResolution_+0.5l));
const long endBin = static_cast<long>(floor(endTTimes_.at(index)/gridResolution_+0.5l));
const double startBinFraction = (startTTimes_.at(index)/gridResolution_)-floor(startTTimes_.at(index)/gridResolution_);
const double endBinFraction = (endTTimes_.at(index)/gridResolution_)-floor(endTTimes_.at(index)/gridResolution_);
// only do timing modification if contributions would land in same bin
if (startBin == endBin) {
// if startE arrives before endE
if (deltaT_ >= 0) {
if ((startBinFraction >= 0.5) && (endBinFraction >= 0.5)) // both points left of bin center
{
startTTimes_.at(index) = startTTimes_.at(index) - gridResolution_; // shift EV1_ to previous gridpoint
}
else if ((startBinFraction < 0.5) && (endBinFraction < 0.5)) // both points right of bin center
{
endTTimes_.at(index) = endTTimes_.at(index) + gridResolution_; // shift EV2_ to next gridpoint
}
else // points on both sides of bin center
{
const double leftDist = 1.0-startBinFraction;
const double rightDist = endBinFraction;
// check if asymmetry to right or left
if (rightDist >= leftDist)
{
endTTimes_.at(index) = endTTimes_.at(index) + gridResolution_; // shift EV2_ to next gridpoint
}
else
{
startTTimes_.at(index) = startTTimes_.at(index) - 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
{
endTTimes_.at(index) = endTTimes_.at(index) - gridResolution_; // shift EV2_ to previous gridpoint
}
else if ((startBinFraction < 0.5) && (endBinFraction < 0.5)) // both points right of bin center
{
startTTimes_.at(index) = startTTimes_.at(index) + gridResolution_; // shift EV1_ to next gridpoint
}
else // points on both sides of bin center
{
const double leftDist = 1.0-endBinFraction;
const double rightDist = startBinFraction;
// check if asymmetry to right or left
if (rightDist >= leftDist)
{
startTTimes_.at(index) = startTTimes_.at(index) + gridResolution_; // shift EV1_ to next gridpoint
}
else
{
endTTimes_.at(index) = endTTimes_.at(index) - gridResolution_; // shift EV2_ to previous gridpoint
}
}
} // End of else statement
} // End of if for startbin == endbin
} // End of if deltaT < gridresolution
} // End of checking for very small doppler factors
// perform ZHS-like calculation close to Cherenkov angle
if (fabs(preDoppler__) <= approxThreshold_ || fabs(postDoppler.at(index)) <= approxThreshold_) {
// get global simulation time for the middle point of that track. (This is my best guess for now)
auto midTime_{particle.getTime() - (track.getDuration() / 2)};
// get "mid" position of the track (that may not work properly)
auto midPoint_{track.getPosition(0.5)};
// get the Path (path3) from the middle "endpoint" to the antenna.
// This is a SignalPathCollection
auto paths3{this->propagator_.propagate(midPoint_, antenna.getLocation())};
// std::size_t j_index {0}; // this will be useful for multiple paths (aka curved propagators)
// now loop over the paths for endpoint that we got above
for (auto const& path : paths3) {
// EVstart_.erase(EVstart_.begin() + index + j_index); // this should work for curved + curved propagators
// EVend_.erase(EVend_.begin() + index + j_index); // for now just use one index and not j_index since at the moment you are working with StraightPropagator
auto const midPointReceiveTime_{path.total_time_ + midTime_};
auto midDoppler_{1. - path.average_refractive_index_ * beta_ * path.emit_};
// change the values of the receive unit vectors of start and end
ReceiveVectorsStart_.at(index) = path.receive_;
ReceiveVectorsEnd_.at(index) = path.receive_;
// CoREAS calculation -> get ElectricFieldVector3 for "midPoint"
ElectricFieldVector EVmid_ = (charge_ / constants::c) *
path.receive_.cross(path.receive_.cross(beta_)) /
(path.R_distance_ * midDoppler_);
// EVstart_.insert(EVstart_.begin() + index + j_index, EVmid_); // this should work for curved + curved propagators
// EVend_.insert(EVend_.begin() + index + j_index, - EVmid_); // for now just use one index and not j_index since at the moment you are working with StraightPropagator
EVstart_.at(index) = EVmid_;
EVend_.at(index) = - EVmid_;
auto deltaT_{midPoint_.getNorm() / (constants::c * beta_ * fabs(midDoppler_))};
if (startTTimes_.at(index) < endTTimes_.at(index)) // EVstart_ arrives earlier
{
startTTimes_.at(index) = midPointReceiveTime_ - 0.5 * deltaT_;
endTTimes_.at(index) = midPointReceiveTime_ + 0.5 * deltaT_;
}
else // EVend_ arrives earlier
{
startTTimes_.at(index) = midPointReceiveTime_ + 0.5 * deltaT_;
endTTimes_.at(index) = midPointReceiveTime_ - 0.5 * deltaT_;
}
const long double gridResolution_{antenna.duration_};
deltaT_ = endTTimes_.at(index) - startTTimes_.at(index);
// redistribute contributions over time scale defined by the observation time resolution
if (fabs(deltaT_) < gridResolution_) {
EVstart_.at(index) = EVstart_.at(index) * fabs(deltaT_ / gridResolution_);
EVend_.at(index) = EVend_.at(index) * fabs(deltaT_ / gridResolution_);
const long startBin = static_cast<long>(floor(startTTimes_.at(index)/gridResolution_+0.5l));
const long endBin = static_cast<long>(floor(endTTimes_.at(index)/gridResolution_+0.5l));
const double startBinFraction = (startTTimes_.at(index)/gridResolution_)-floor(startTTimes_.at(index)/gridResolution_);
const double endBinFraction = (endTTimes_.at(index)/gridResolution_)-floor(endTTimes_.at(index)/gridResolution_);
// only do timing modification if contributions would land in same bin
if (startBin == endBin) {
// if startE arrives before endE
if (deltaT_ >= 0) {
if ((startBinFraction >= 0.5) && (endBinFraction >= 0.5)) // both points left of bin center
{
startTTimes_.at(index) = startTTimes_.at(index) - gridResolution_; // shift EV1_ to previous gridpoint
}
else if ((startBinFraction < 0.5) && (endBinFraction < 0.5)) // both points right of bin center
{
endTTimes_.at(index) = endTTimes_.at(index) + gridResolution_; // shift EV2_ to next gridpoint
}
else // points on both sides of bin center
{
const double leftDist = 1.0-startBinFraction;
const double rightDist = endBinFraction;
// check if asymmetry to right or left
if (rightDist >= leftDist)
{
endTTimes_.at(index) = endTTimes_.at(index) + gridResolution_; // shift EV2_ to next gridpoint
}
else
{
startTTimes_.at(index) = startTTimes_.at(index) - 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
{
endTTimes_.at(index) = endTTimes_.at(index) - gridResolution_; // shift EV2_ to previous gridpoint
}
else if ((startBinFraction < 0.5) && (endBinFraction < 0.5)) // both points right of bin center
{
startTTimes_.at(index) = startTTimes_.at(index) + gridResolution_; // shift EV1_ to next gridpoint
}
else // points on both sides of bin center
{
const double leftDist = 1.0-endBinFraction;
const double rightDist = startBinFraction;
// check if asymmetry to right or left
if (rightDist >= leftDist)
{
startTTimes_.at(index) = startTTimes_.at(index) + gridResolution_; // shift EV1_ to next gridpoint
}
else
{
endTTimes_.at(index) = endTTimes_.at(index) - gridResolution_; // shift EV2_ to previous gridpoint
}
}
} // End of else statement
} // End of if for startbin == endbin
} // End of if deltaT < gridresolution
} // End of looping over paths3
} // end of ZHS-like approximation
// Feed start and end to the antenna
antenna.receive(startTTimes_.at(index), ReceiveVectorsStart_.at(index), EVstart_.at(index));
antenna.receive(endTTimes_.at(index), ReceiveVectorsEnd_.at(index), EVend_.at(index));
// update index
index = index + 1;
} // End of for loop for preDoppler factor (this includes checking for postDoppler factors)
} // End of checking of vector sizes
} // End of looping over the antennas.
} // End of simulate method.
/**
* Return the maximum step length for this particle and track.
......@@ -218,7 +392,7 @@ namespace corsika {
// This is part of the ZHS / CoReas formalisms and can
// be related from the magnetic field / acceleration, charge,
// etc. of the particle.
return 1000000000000;
return 1000000000000_m;
}
}; // END: class RadioProcess
......
corsika/modules/radio/RadioProcess.hpp
+ 3
− 2
View file @ 35158b33
......@@ -111,13 +111,14 @@ namespace corsika {
* we wait for the true output formatting to be ready.
**/
bool writeOutput() const {
// this for loop still has some issues
int i = 1;
for (auto& antenna : detector_.getAntennas()) {
auto [t,E] = antenna.getWaveform();
std::ofstream out_file ("antenna" + to_string(i) + "_output.csv");
xt::dump_csv(out_file, t, E);
xt::dump_csv(out_file, t);
xt::dump_csv(out_file, E);
out_file.close();
++i;
......
corsika/modules/radio/ZHS.hpp
+ 14
− 11
View file @ 35158b33
......@@ -11,6 +11,7 @@
#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 <bits/stdc++.h>
namespace corsika {
......@@ -72,24 +73,26 @@ namespace corsika {
// get the Path from the track to the antenna
// This is a SignalPathCollection
auto paths{this->propagator_.propagate(startPoint_, antenna.getLocation())};
auto R_ {(startPoint_ - antenna.getLocation()).getNorm()};
auto paths{this->propagator_.propagate(startPoint_, antenna.getLocation(), 1_nm)};
// now loop over the paths that we got above
// Note: for the StraightPropagator, there will only be a single
// path but other propagators may return more than one.
for (auto const& path : paths) {
auto t1 = path.total_time_;
QuantityVector<ElectricFieldType::dimension_type> v11{10_V / 1_m, 10_V / 1_m, 10_V / 1_m};
// calculate the ZHS formalism for this particle-antenna
ElectricFieldVector EV_ = ((- constants) * trackVelocity_.dot(path.emit)) *
((midTime_ + path.total_time_ - (1 - path.average_refractivity_ * beta_ *
acos(track.getDirection(0).dot(path.emit_))) * startTime_)
- (midTime_ + path.total_time_ - (1 - path.average_refractivity_ * beta_ *
acos(track.getDirection(0).dot(path.emit_))) * endTime_))
/ (1 - path.average_refractivity_ * beta_ * acos(track.getDirection(0).dot(path.emit_)));
// pass it to the antenna
antenna.receive(midTime_ + path.total_time_, path.receive_, EV_);
// ElectricFieldVector EV_ = ((- constants) * trackVelocity_.dot(path.emit)) *
// ((midTime_ + path.total_time_ - (1 - path.average_refractivity_ * beta_ *
// acos(track.getDirection(0).dot(path.emit_))) * startTime_)
// - (midTime_ + path.total_time_ - (1 - path.average_refractivity_ * beta_ *
// acos(track.getDirection(0).dot(path.emit_))) * endTime_))
// / (1 - path.average_refractivity_ * beta_ * acos(track.getDirection(0).dot(path.emit_)));
//
// // pass it to the antenna
// antenna.receive(midTime_ + path.total_time_, path.receive_, EV_);
} // END: loop over paths
......
corsika/modules/radio/antennas/Antenna.hpp
+ 1
− 1
View file @ 35158b33
......@@ -34,7 +34,7 @@ namespace corsika {
// QuantityVector<MagneticFieldType::dimension_type>;
// a dimensionless vector used for the incident direction
using Vector = QuantityVector<dimensionless_d>;
// using Vector = QuantityVector<dimensionless_d>;
/**
* \brief Construct a base antenna instance.
......
corsika/modules/radio/antennas/TimeDomainAntenna.hpp
+ 20
− 20
View file @ 35158b33
......@@ -18,26 +18,25 @@ namespace corsika {
/**
* An implementation of a time-domain antenna that has a customized
* start time, sampling period, and waveform duration.
* start time, sampling rate, and waveform duration.
*
*/
class TimeDomainAntenna final : public Antenna<TimeDomainAntenna> {
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 num_bins_; ///< The number of bins used.
xt::xtensor<double,3> waveformE_; ///< The waveform stored by this antenna.
std::pair<xt::xtensor<double, 1>,
xt::xtensor<double,3>> waveform_; ///< useful for .getWaveform()
xt::xtensor<double,2> waveformE_; ///< The waveform stored by this antenna.
std::pair<xt::xtensor<double, 2>,
xt::xtensor<double,2>> waveform_; ///< useful for .getWaveform()
protected:
// expose the CRTP interfaces constructor
public:
// import the methods from the antenna
using Antenna<TimeDomainAntenna>::Antenna;
using Antenna<TimeDomainAntenna>::receive;
using Antenna<TimeDomainAntenna>::getName;
using Antenna<TimeDomainAntenna>::getLocation;
......@@ -62,7 +61,7 @@ namespace corsika {
, duration_(duration)
, sample_rate_(sample_rate)
, num_bins_ (static_cast<int>(duration * sample_rate))
, waveformE_ (xt::zeros<double>({3, num_bins_}))
, waveformE_ (xt::zeros<double>({num_bins_, 3}))
{};
/**
......@@ -76,20 +75,20 @@ namespace corsika {
* @param field The incident electric field vector.
*
*/
void receive(TimeType const time, Vector const& receive_vector,
void receive(TimeType const time, Vector<dimensionless_d> const& receive_vector,
ElectricFieldVector const& efield) {
if (time < start_time_ || time > start_time_ + duration_) {
if (time < start_time_ || time >= start_time_ + duration_) {
return;
} else {
// figure out the correct timebin to store the E-field value.
auto timebin_ {static_cast<std::size_t>((time - start_time_) * sample_rate_)};
// store the x,y,z electric field components.
waveformE_.at(0, timebin_) += (efield.getX().magnitude());
waveformE_.at(1, timebin_) += (efield.getY().magnitude());
waveformE_.at(2, timebin_) += (efield.getZ().magnitude());
waveformE_.at(timebin_, 0) += efield.getX().magnitude();
waveformE_.at(timebin_, 1) += efield.getY().magnitude();
waveformE_.at(timebin_, 2) += efield.getZ().magnitude();
//TODO: Check how they are stored in memory, row-wise or column-wise?
}
}
......@@ -100,17 +99,17 @@ namespace corsika {
*
* @returns A pair of the sample times, and the field
*/
std::pair<xt::xtensor<double, 1>, xt::xtensor<double,3>> getWaveform() const {
std::pair<xt::xtensor<double, 2>, xt::xtensor<double,2>> getWaveform() const {
// TODO: divide by Δt for CoREAS ONLY!
// create a 1-D xtensor to store time values so we can print them later.
xt::xtensor<double, 1> times_ {xt::zeros<double>({num_bins_, 1})};
xt::xtensor<double, 2> times_ (xt::zeros<double>({num_bins_, 1}));
for (auto i = 0; i <= num_bins_; i++) {
times_.at(i) = static_cast<double>(start_time_ / 1_ns) +
static_cast<double>(i * sample_rate_ * 1_ns);
for (int i = 0; i < num_bins_; i++) {
times_.at(i,0) = static_cast<double>(start_time_ / 1_s + i * sample_rate_ * 1_s);
}
return waveform_;
return std::make_pair(times_, waveformE_);
};
/**
......@@ -118,6 +117,7 @@ namespace corsika {
*/
void reset() {
waveformE_ = xt::zeros_like(waveformE_);
// times_ = xt::zeros_like(times_);
};
}; // END: class Antenna final
......
corsika/modules/radio/detectors/RadioDetector.hpp
+ 1
− 1
View file @ 35158b33
......@@ -37,7 +37,7 @@ namespace corsika {
*
* @returns An iterable mutable reference to the antennas.
*/
std::vector<TAntennaImpl> const& getAntennas() { return antennas_; }
std::vector<TAntennaImpl> const& getAntennas() { return antennas_; } // maybe the const& here is an issue
/**
* Reset all the antenna waveforms.
......
tests/modules/testRadio.cpp
+ 258
− 23
View file @ 35158b33
......@@ -8,12 +8,22 @@
#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/TimeDomainDetector.hpp>
#include <corsika/modules/radio/detectors/RadioDetector.hpp>
#include <corsika/modules/radio/propagators/StraightPropagator.hpp>
#include <corsika/modules/radio/propagators/SignalPath.hpp>
#include <corsika/modules/radio/propagators/RadioPropagator.hpp>
#include <vector>
#include <xtensor/xtensor.hpp>
#include <xtensor/xbuilder.hpp>
#include <xtensor/xio.hpp>
#include <xtensor/xcsv.hpp>
#include <istream>
#include <fstream>
#include <iostream>
#include <corsika/media/Environment.hpp>
#include <corsika/media/HomogeneousMedium.hpp>
......@@ -39,45 +49,158 @@ double constexpr absMargin = 1.0e-7;
TEST_CASE("Radio", "[processes]") {
// logging::set_level(logging::level::info);
// corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
logging::set_level(logging::level::info);
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
// check that I can create and reset a TimeDomain process
SECTION("TimeDomainDetector") {
// construct a time domain detector
TimeDomainDetector<TimeDomainAntenna> detector;
// // and construct some time domain radio process
// TimeDomain<ZHS<>, TimeDomainDetector<TimeDomainAntenna>> process(detector);
// TimeDomain<ZHS<CPU>, TimeDomainDetector<TimeDomainAntenna>> process2(detector);
// TimeDomain<ZHS<CPU>, decltype(detector)> process3(detector);
}
// SECTION("TimeDomainDetector") {
//
// // construct a time domain detector
// AntennaCollection<TimeDomainAntenna> detector;
//
// // // and construct some time domain radio process
// // TimeDomain<ZHS<>, TimeDomainDetector<TimeDomainAntenna>> process(detector);
// // TimeDomain<ZHS<CPU>, TimeDomainDetector<TimeDomainAntenna>> process2(detector);
// // TimeDomain<ZHS<CPU>, decltype(detector)> process3(detector);
// }
// check that I can create time domain antennas
SECTION("TimeDomainDetector") {
SECTION("TimeDomainAntenna") {
// create an environment so we can get a coordinate system
Environment<IMediumModel> env;
// the antenna location
const auto point{Point(env.getCoordinateSystem(), 1_m, 2_m, 3_m)};
const auto point1{Point(env.getCoordinateSystem(), 1_m, 2_m, 3_m)};
const auto point2{Point(env.getCoordinateSystem(), 4_m, 5_m, 6_m)};
// create times for the antenna
const TimeType t1{10_s};
const TimeType t2{10_s};
const InverseTimeType t3{1/1_s};
// check that I can create an antenna at (1, 2, 3)
const auto ant{TimeDomainAntenna("antenna_name", point)};
TimeDomainAntenna ant1("antenna_name", point1, t1, t2, t3);
TimeDomainAntenna ant2("antenna_name2", point2, 11_s, t2, t3);
// assert that the antenna name is correct
REQUIRE(ant.getName() == "antenna_name");
REQUIRE(ant1.getName() == "antenna_name");
// and check that the antenna is at the right location
REQUIRE((ant.getLocation() - point).getNorm() < 1e-12 * 1_m);
REQUIRE((ant1.getLocation() - point1).getNorm() < 1e-12 * 1_m);
// construct a radio detector instance to store our antennas
TimeDomainDetector<TimeDomainAntenna> detector;
AntennaCollection<TimeDomainAntenna> detector;
// add this antenna to the process
detector.addAntenna(ant);
detector.addAntenna(ant1);
detector.addAntenna(ant2);
// create an environment with uniform refractive index of 1
using UniRIndex =
UniformRefractiveIndex<HomogeneousMedium<IRefractiveIndexModel<IMediumModel>>>;
using EnvType = Environment<IRefractiveIndexModel<IMediumModel>>;
EnvType env6;
// 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(
std::vector<Code>{Code::Nitrogen},
std::vector<float>{1.f}));
env6.getUniverse()->addChild(std::move(Medium6));
// get a unit vector
Vector<dimensionless_d> v1(rootCS6, {0, 0, 1});
QuantityVector<ElectricFieldType::dimension_type> v11{10_V / 1_m, 10_V / 1_m, 10_V / 1_m};
Vector<dimensionless_d> v2(rootCS6, {0, 1, 0});
QuantityVector<ElectricFieldType::dimension_type> v22{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 getWaveform() method
auto [t11, E1] = ant1.getWaveform();
CHECK(E1(5,0) - 10 == 0);
auto [t222, E2] = ant2.getWaveform();
CHECK(E2(5,0) -20 == 0);
// the rest is just random ideas
//////////////////////////////////////////////////////////////////////////////////////
// for (auto& kostas : detector.getAntennas()) {
// std::cout << kostas.getName() << std::endl;
// kostas.receive(15_s, v1, v11);
// std::cout << kostas.waveformE_;
// auto [t, E] = kostas.getWaveform();
// std::cout << E(5,0) << std::endl;
// kostas.receive(16_s, v2, v22);
// std::cout << E(5,0) << std::endl;
// }
// auto t33{static_cast<int>(t1 / t2)};
// std::cout << t33 << std::endl;
// xt::xtensor<double,2> waveformE_ = xt::zeros<double>({2, 3});
// xt::xtensor<double,2> res = xt::view(waveformE_);
// std::cout << ant1.waveformE_ << std::endl;
// std::cout << " " << std::endl;
// std::cout << ant2.waveformE_ << std::endl;
// std::cout << " " << std::endl;
// std::cout << ant1.times_ << std::endl;
// std::cout << " " << std::endl;
// std::cout << ant2.times_ << std::endl;
// auto [t, E] = ant2.getWaveform();
//
//
//
// std::ofstream out_file("antenna_output.csv");
// xt::dump_csv(out_file, t);
// xt::dump_csv(out_file, E);
// out_file.close();
// for (auto& antenna : detector.getAntennas()) {
// auto [t, E] = antenna.getWaveform();
// }
// int i = 1;
// auto x = detector.getAntennas();
// for (auto& antenna : x) {
//
// auto [t, E] = antenna.getWaveform();
// std::ofstream out_file("antenna" + to_string(i) + "_output.csv");
// std::cout << E(5,0) << std::endl;
// xt::dump_csv(out_file, t);
// xt::dump_csv(out_file, E);
// out_file.close();
// ++i;
// }
// detector.writeOutput();
// ant1.reset();
// ant2.reset();
//
// std::cout << ant1.waveformE_ << std::endl;
// std::cout << ant2.waveformE_ << std::endl;
//////////////////////////////////////////////////////////////////////////////////////
}
// check that I can create working Straight Propagators in different environments
......@@ -133,9 +256,10 @@ TEST_CASE("Radio", "[processes]") {
for (auto const& path : paths_) {
CHECK((path.total_time_ / 1_s) - ((34_m / (3 * constants::c)) / 1_s) ==
Approx(0).margin(absMargin));
CHECK(path.average_refractivity_ == Approx(1));
CHECK(path.average_refractive_index_ == Approx(1));
CHECK(path.emit_.getComponents() == v1.getComponents());
CHECK(path.receive_.getComponents() == v1.getComponents());
CHECK(path.R_distance_ == 10_m);
// CHECK(std::equal(P1.begin(), P1.end(), path.points_.begin(),[]
// (Point a, Point b) { return (a - b).norm() / 1_m < 1e-5;}));
//TODO:THINK ABOUT THE POINTS IN THE SIGNALPATH.H
......@@ -198,9 +322,10 @@ TEST_CASE("Radio", "[processes]") {
for (auto const& path :paths1_) {
CHECK( (path.total_time_ / 1_s) - ((34_m / (3 * constants::c)) / 1_s)
== Approx(0).margin(absMargin) );
CHECK( path.average_refractivity_ == Approx(1) );
CHECK( path.average_refractive_index_ == Approx(1) );
CHECK( path.emit_.getComponents() == vv1.getComponents() );
CHECK( path.receive_.getComponents() == vv1.getComponents() );
CHECK( path.R_distance_ == 10_m );
}
CHECK( paths1_.size() == 1 );
......@@ -249,9 +374,10 @@ TEST_CASE("Radio", "[processes]") {
for (auto const& path :paths2_) {
CHECK( (path.total_time_ / 1_s) - ((3.177511688_m / (3 * constants::c)) / 1_s)
== Approx(0).margin(absMargin) );
CHECK( path.average_refractivity_ == Approx(0.210275935) );
CHECK( path.average_refractive_index_ == Approx(0.210275935) );
CHECK( path.emit_.getComponents() == vvv1.getComponents() );
CHECK( path.receive_.getComponents() == vvv1.getComponents() );
CHECK( path.R_distance_ == 10_m );
}
CHECK( paths2_.size() == 1 );
......@@ -262,6 +388,115 @@ TEST_CASE("Radio", "[processes]") {
//
// // TODO: construct the environment for the propagator
//////////////////////////////////////////////////////////////////////////////////////
// // useful information
// // create an environment so we can get a coordinate system
// environment::Environment<environment::IMediumModel> env;
//
// // the antenna location
// const auto point{geometry::Point(env.GetCoordinateSystem(), 1_m, 2_m, 3_m)};
//
// // check that I can create an antenna at (1, 2, 3)
// const auto ant{TimeDomainAntenna("antenna_name", point)};
//
// // assert that the antenna name is correct
// REQUIRE(ant.GetName() == "antenna_name");
//
// // and check that the antenna is at the right location
// REQUIRE((ant.GetLocation() - point).norm() < 1e-12 * 1_m);
//
// // construct a radio detector instance to store our antennas
// TimeDomainDetector<TimeDomainAntenna> detector;
//
// // add this antenna to the process
// detector.AddAntenna(ant);
//
// // create a TimeDomain process
// auto zhs{TimeDomain<ZHS<CPU>, decltype(detector)>(detector)};
//
// // get the antenna back from the process and check the properties
// REQUIRE(detector.GetAntennas().cbegin()->GetName() == "antenna_name");
//
// // and check the location is the same
// REQUIRE((detector.GetAntennas().cbegin()->GetLocation() - point).norm() <
// 1e-12 * 1_m);
//
// // and ANOTHER antenna
// const auto ant2{TimeDomainAntenna("antenna_name", point)};
// detector.AddAntenna(ant2);
// // and check that the number of antennas visible to the process has increased
// REQUIRE(zhs.GetDetector().GetAntennas().size() == 2);
//////////////////////////////////////////////////////////////////////////////////////////
// using VelocityVec = Vector<corsika::units::si::SpeedType::dimension_type>;
//
// // setup environment, geometry
// environment::Environment<environment::IMediumModel> env;
//
// // and get the coordinate systm
// geometry::CoordinateSystem const& cs{env.GetCoordinateSystem()};
//
// // a test particle
// const auto pcode{particles::Code::Electron};
// const auto mass{particles::Electron::GetMass()};
//
// // create a new stack for each trial
// setup::Stack stack;
//
// // construct an energy
// const HEPEnergyType E0{1_TeV};
//
// // compute the necessary momentumn
// const HEPMomentumType P0{sqrt(E0 * E0 - mass * mass)};
//
// // and create the momentum vector
// const auto plab{corsika::stack::MomentumVector(cs, {0_GeV, 0_GeV, P0})};
//
// // and create the location of the particle in this coordinate system
// const geometry::Point pos(cs, 5_m, 1_m, 8_m);
//
// // add the particle to the stack
// auto particle{stack.AddParticle(std::make_tuple(pcode, E0, plab, pos, 0_ns))};
//
// // get a view into the stack
// corsika::stack::SecondaryView stackview(particle);
//
// // create a trajectory
// const auto& track{
// Trajectory(Line(pos, VelocityVec(cs, 0_m / 1_s, 0_m / 1_s, 0.1_m / 1_ns)), 1_ns)};
//
// // and create a view vector
// const auto& view{Vector(cs, 0_m, 0_m, 1_m)};
//
// // construct a radio detector instance to store our antennas
// TimeDomainDetector<TimeDomainAntenna> detector;
//
// // the antenna location
// const auto point{geometry::Point(cs, 1_m, 2_m, 3_m)};
//
// // check that I can create an antenna at (1, 2, 3)
// const auto ant{TimeDomainAntenna("antenna_name", point)};
//
// // add this antenna to the process
// detector.AddAntenna(ant);
//
// // create a TimeDomain process
// auto zhs{TimeDomain<ZHS<CPU>, decltype(detector)>(detector)};
//
// // get the projectile
// auto projectile{stackview.GetProjectile()};
//
// // call ZHS over the track
// zhs.DoContinuous(projectile, track);
// zhs.Simulate(projectile, track);
//
// // try and call the ZHS implementation directly for a given view direction
// ZHS<CPU>::Emit(projectile, track, view);
//////////////////////////////////////////////////////////////////////////////////////////
// create an environment with uniform refractive index of 1
// using UniRIndex =
// UniformRefractiveIndex<HomogeneousMedium<IRefractiveIndexModel<IMediumModel>>>;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

data privacy - impressum