IAP GITLAB

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

Introduce a coordinate system for the antennas

parent 699214a4
No related branches found
No related tags found
1 merge request!329Radio interface
/*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
* (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
......@@ -14,9 +14,10 @@
namespace corsika {
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, CoordinateSystemPtr const& coordinateSystem)
: name_(name)
, location_(location){};
, location_(location)
, coordinateSystem_(coordinateSystem) {};
template <typename TAntennaImpl>
inline Point const& Antenna<TAntennaImpl>::getLocation() const {
......@@ -68,15 +69,15 @@ namespace corsika {
template <typename TAntennaImpl>
inline void Antenna<TAntennaImpl>::endOfShower(const int event,
const std::string radioImplementation,
std::string const& radioImplementation,
const double sampleRate) {
// get the copy of the waveform data for this event
// we transpose it so that we can match dimensions with the
// time array that is already in the output file
std::vector<double> dataX = this->implementation().getDataX();
std::vector<double> dataY = this->implementation().getDataY();
std::vector<double> dataZ = this->implementation().getDataZ();
std::vector<double> const& dataX = this->implementation().getDataX();
std::vector<double> const& dataY = this->implementation().getDataY();
std::vector<double> const& dataZ = this->implementation().getDataZ();
if (radioImplementation == "ZHS") {
std::vector<double> electricFieldX(dataX.size() - 1,
......@@ -89,24 +90,18 @@ namespace corsika {
electricFieldZ.at(i) = -(dataZ.at(i + 1) - dataZ.at(i)) * sampleRate;
}
// cnpy needs a vector for the shape
std::vector<size_t> shapeX = {electricFieldX.size()};
std::vector<size_t> shapeY = {electricFieldY.size()};
std::vector<size_t> shapeZ = {electricFieldZ.size()};
cnpy::npz_save(filename_, std::to_string(event) + "X", electricFieldX.data(),
shapeX, "a");
{electricFieldX.size()}, "a");
cnpy::npz_save(filename_, std::to_string(event) + "Y", electricFieldY.data(),
shapeY, "a");
{electricFieldY.size()}, "a");
cnpy::npz_save(filename_, std::to_string(event) + "Z", electricFieldZ.data(),
shapeZ, "a");
{electricFieldZ.size()}, "a");
} else {
// cnpy needs a vector for the shape
std::vector<size_t> shapeX = {dataX.size()};
std::vector<size_t> shapeY = {dataY.size()};
std::vector<size_t> shapeZ = {dataZ.size()};
// 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");
cnpy::npz_save(filename_, std::to_string(event) + "X", dataX.data(), {dataX.size()}, "a");
cnpy::npz_save(filename_, std::to_string(event) + "Y", dataY.data(), {dataY.size()}, "a");
cnpy::npz_save(filename_, std::to_string(event) + "Z", dataZ.data(), {dataZ.size()}, "a");
}
}
......
/*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
* (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
......@@ -28,6 +28,7 @@ namespace corsika {
public:
std::string const name_; ///< The name/identifier of this antenna.
Point const location_; ///< The location of this antenna.
CoordinateSystemPtr const coordinateSystem_; ///< The coordinate system of the antenna
std::string filename_ = ""; ///< The filename for the output file for this antenna.
// this stores the polarization vector of an electric field
......@@ -41,7 +42,7 @@ namespace corsika {
* @param location The location of this antenna.
*
*/
Antenna(std::string const& name, Point const& location);
Antenna(std::string const& name, Point const& location, CoordinateSystemPtr const& coordinateSystem);
/**
* Receive a signal at this antenna.
......@@ -84,7 +85,7 @@ namespace corsika {
* This is used when writing the antenna information to disk
* and will be converted to a 32-bit float before writing.
*/
std::vector<double>& getDataX() const;
std::vector<double> const& getDataX() const;
/**
* Return a reference to the underlying data for Y polarization.
......@@ -92,7 +93,7 @@ namespace corsika {
* This is used when writing the antenna information to disk
* and will be converted to a 32-bit float before writing.
*/
std::vector<double>& getDataY() const;
std::vector<double> const& getDataY() const;
/**
* Return a reference to the underlying data for Z polarization.
......@@ -100,7 +101,7 @@ namespace corsika {
* This is used when writing the antenna information to disk
* and will be converted to a 32-bit float before writing.
*/
std::vector<double>& getDataZ() const;
std::vector<double> const& getDataZ() const;
/**
* Prepare for the start of the library.
......@@ -111,7 +112,7 @@ namespace corsika {
/**
* Flush the data from this shower to disk.
*/
void endOfShower(int const event, std::string const radioImplementation,
void endOfShower(int const event, std::string const& radioImplementation,
double const sampleRate);
protected:
......
......@@ -210,7 +210,7 @@ int main(int argc, char** argv) {
auto triggertime_1{(triggerpoint_ - point_1).getNorm() / constants::c};
std::string name_1 = "CoREAS_R=" + std::to_string(rr_1) +
"_m--Phi=" + std::to_string(phi_1) + "degrees";
TimeDomainAntenna antenna_1(name_1, point_1, triggertime_1, duration_, sampleRate_,
TimeDomainAntenna antenna_1(name_1, point_1, rootCS, triggertime_1, duration_, sampleRate_,
triggertime_1);
detectorCoREAS.addAntenna(antenna_1);
}
......@@ -230,7 +230,7 @@ int main(int argc, char** argv) {
auto triggertime_{(triggerpoint_ - point_).getNorm() / constants::c};
std::string name_ =
"ZHS_R=" + std::to_string(rr_) + "_m--Phi=" + std::to_string(phi_) + "degrees";
TimeDomainAntenna antenna_(name_, point_, triggertime_, duration_, sampleRate_,
TimeDomainAntenna antenna_(name_, point_, rootCS, triggertime_, duration_, sampleRate_,
triggertime_);
detectorZHS.addAntenna(antenna_);
}
......
......@@ -107,8 +107,8 @@ int main() {
const InverseTimeType t3{5e+11_Hz};
// the antennas
TimeDomainAntenna ant1("antenna CoREAS", point1, t1, t2, t3, t1);
TimeDomainAntenna ant2("antenna ZHS", point1, t1, t2, t3, t1);
TimeDomainAntenna ant1("antenna CoREAS", point1, rootCS, t1, t2, t3, t1);
TimeDomainAntenna ant2("antenna ZHS", point1, rootCS, t1, t2, t3, t1);
// the detectors
AntennaCollection<TimeDomainAntenna> detectorCoREAS;
......
......@@ -95,8 +95,8 @@ int main() {
std::cout << "number of points in time: " << duration * sampleRate_ << std::endl;
// create 2 antennas
TimeDomainAntenna ant1("CoREAS antenna", point1, start, duration, sampleRate_, start);
TimeDomainAntenna ant2("ZHS antenna", point1, start, duration, sampleRate_, start);
TimeDomainAntenna ant1("CoREAS antenna", point1, rootCS, start, duration, sampleRate_, start);
TimeDomainAntenna ant2("ZHS antenna", point1, rootCS, start, duration, sampleRate_, start);
// construct a radio detector instance to store our antennas
AntennaCollection<TimeDomainAntenna> detectorCoREAS;
......
/*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
* (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
......@@ -93,8 +93,8 @@ TEST_CASE("Radio", "[processes]") {
const TimeType t4{11_s};
// check that I can create an antenna at (1, 2, 3)
TimeDomainAntenna ant1("antenna_name", point1, t1, t2, t3, t1);
TimeDomainAntenna ant2("antenna_name2", point2, t1, t2, t3, t1);
TimeDomainAntenna ant1("antenna_name", point1, rootCSCoREAS, t1, t2, t3, t1);
TimeDomainAntenna ant2("antenna_name2", point2, rootCSCoREAS, t1, t2, t3, t1);
// construct a radio detector instance to store our antennas
AntennaCollection<TimeDomainAntenna> detector;
......@@ -197,8 +197,8 @@ TEST_CASE("Radio", "[processes]") {
const TimeType t4{11_s};
// check that I can create an antenna at (1, 2, 3)
TimeDomainAntenna ant1("antenna_zhs", point1, t1, t2, t3, t1);
TimeDomainAntenna ant2("antenna_zhs2", point2, t1, t2, t3, t1);
TimeDomainAntenna ant1("antenna_zhs", point1, rootCSZHS, t1, t2, t3, t1);
TimeDomainAntenna ant2("antenna_zhs2", point2, rootCSZHS, t1, t2, t3, t1);
// construct a radio detector instance to store our antennas
AntennaCollection<TimeDomainAntenna> detector;
......@@ -292,8 +292,8 @@ TEST_CASE("Antennas") {
const TimeType t4{11_s};
// check that I can create an antenna at (1, 2, 3)
TimeDomainAntenna ant1("antenna_name", point1, t1, t2, t3, t1);
TimeDomainAntenna ant2("antenna_name2", point2, t4, t2, t3, t4);
TimeDomainAntenna ant1("antenna_name", point1, rootCS6, t1, t2, t3, t1);
TimeDomainAntenna ant2("antenna_name2", point2, rootCS6, t4, t2, t3, t4);
// assert that the antenna name is correct
REQUIRE(ant1.getName() == "antenna_name");
......@@ -324,21 +324,24 @@ TEST_CASE("Antennas") {
ant1.receive(15_s, v1, v11);
ant2.receive(16_s, v2, v22);
// use getWaveform() methods
auto[tx, Ex] = ant1.getWaveformX();
// use getDataX,Y,Z() and getAxis() methods
auto Ex = ant1.getDataX();
CHECK(Ex[5] - 10 == 0);
auto tx = ant1.getAxis();
CHECK(tx[5] - 5 * 1_s / 1_ns == Approx(0.0));
auto[ty, Ey] = ant1.getWaveformY();
auto Ey = ant1.getDataY();
CHECK(Ey[5] - 10 == 0);
auto[tz, Ez] = ant1.getWaveformZ();
auto Ez = ant1.getDataZ();
CHECK(Ez[5] - 10 == 0);
auto ty = ant1.getAxis();
auto tz = ant1.getAxis();
CHECK(tx[5] - ty[5] == 0);
CHECK(ty[5] - tz[5] == 0);
auto[tx2, Ex2] = ant2.getWaveformX();
auto Ex2 = ant2.getDataX();
CHECK(Ex2[5] - 20 == 0);
auto[ty2, Ey2] = ant2.getWaveformY();
auto Ey2 = ant2.getDataY();
CHECK(Ey2[5] - 20 == 0);
auto[tz2, Ez2] = ant2.getWaveformZ();
auto Ez2 = ant2.getDataZ();
CHECK(Ez2[5] - 20 == 0);
// the following creates a star-shaped pattern of antennas in the ground
......@@ -360,7 +363,7 @@ TEST_CASE("Antennas") {
std::string var_ = "antenna_R=" + std::to_string(rr_) +
"_m-Phi=" + std::to_string(phi_) + "degrees";
antenna_names.push_back(var_);
TimeDomainAntenna ant111(var_, point_, time__, t2222, t3333, time__);
TimeDomainAntenna ant111(var_, point_, rootCS6, time__, t2222, t3333, time__);
detector__.addAntenna(ant111);
}
}
......
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