From 55dc75f4c8ec3c648edbd951de39139e356d68b4 Mon Sep 17 00:00:00 2001
From: Nikos Karastathis <n.karastathis@kit.edu>
Date: Mon, 15 Mar 2021 19:41:20 +0100
Subject: [PATCH] built compiling environments & some notes for further testing
---
corsika/modules/radio/CoREAS.hpp | 32 +++++------
tests/modules/testRadio.cpp | 98 ++++++++++++++++++++++++++------
2 files changed, 98 insertions(+), 32 deletions(-)
diff --git a/corsika/modules/radio/CoREAS.hpp b/corsika/modules/radio/CoREAS.hpp
index 627ec2063..7600232a6 100755
--- a/corsika/modules/radio/CoREAS.hpp
+++ b/corsika/modules/radio/CoREAS.hpp
@@ -192,15 +192,15 @@ namespace corsika {
auto gridResolution_ {antenna.duration_};
auto deltaT_ { endTTimes_.at(index) - startTTimes_.at(index) };
- if (fabs(deltaT_ / 1_s) < gridResolution_ / 1_s) {
+ if (std::fabs(deltaT_ / 1_s) < gridResolution_ / 1_s) {
- EVstart_.at(index) = EVstart_.at(index) * fabs(deltaT_ / gridResolution_);
- EVend_.at(index) = EVend_.at(index) * fabs(deltaT_ / gridResolution_);
+ EVstart_.at(index) = EVstart_.at(index) * std::fabs(deltaT_ / gridResolution_);
+ EVend_.at(index) = EVend_.at(index) * std::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_);
+ const long startBin = static_cast<long>(std::floor(startTTimes_.at(index)/gridResolution_+0.5l));
+ const long endBin = static_cast<long>(std::floor(endTTimes_.at(index)/gridResolution_+0.5l));
+ const double startBinFraction = (startTTimes_.at(index)/gridResolution_)-std::floor(startTTimes_.at(index)/gridResolution_);
+ const double endBinFraction = (endTTimes_.at(index)/gridResolution_)-std::floor(endTTimes_.at(index)/gridResolution_);
// only do timing modification if contributions would land in same bin
if (startBin == endBin) {
@@ -260,7 +260,7 @@ namespace corsika {
} // 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_) {
+ if (std::fabs(preDoppler__) <= approxThreshold_ || std::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)};
@@ -296,7 +296,7 @@ namespace corsika {
EVstart_.at(index) = EVmid_;
EVend_.at(index) = - EVmid_;
- auto deltaT_{(endPoint_ - startPoint_).getNorm() / (constants::c * beta_.getNorm() * fabs(midDoppler_))}; // TODO: Caution with this!
+ auto deltaT_{(endPoint_ - startPoint_).getNorm() / (constants::c * beta_.getNorm() * std::fabs(midDoppler_))}; // TODO: Caution with this!
if (startTTimes_.at(index) < endTTimes_.at(index)) // EVstart_ arrives earlier
{
@@ -313,15 +313,15 @@ namespace corsika {
deltaT_ = endTTimes_.at(index) - startTTimes_.at(index);
// redistribute contributions over time scale defined by the observation time resolution
- if (fabs(deltaT_ / 1_s) < gridResolution_) {
+ if (std::fabs(deltaT_ / 1_s) < gridResolution_) {
- EVstart_.at(index) = EVstart_.at(index) * fabs((deltaT_ / 1_s) / gridResolution_);
- EVend_.at(index) = EVend_.at(index) * fabs((deltaT_ / 1_s) / gridResolution_);
+ EVstart_.at(index) = EVstart_.at(index) * std::fabs((deltaT_ / 1_s) / gridResolution_);
+ EVend_.at(index) = EVend_.at(index) * std::fabs((deltaT_ / 1_s) / gridResolution_);
- const long startBin = static_cast<long>(floor((startTTimes_.at(index) / 1_s)/gridResolution_+0.5l));
- const long endBin = static_cast<long>(floor((endTTimes_.at(index) / 1_s) /gridResolution_+0.5l));
- const double startBinFraction = ((startTTimes_.at(index) / 1_s)/gridResolution_)-floor((startTTimes_.at(index) / 1_s)/gridResolution_);
- const double endBinFraction = ((endTTimes_.at(index) / 1_s)/gridResolution_)-floor((endTTimes_.at(index) / 1_s)/gridResolution_);
+ const long startBin = static_cast<long>(std::floor((startTTimes_.at(index) / 1_s)/gridResolution_+0.5l));
+ const long endBin = static_cast<long>(std::floor((endTTimes_.at(index) / 1_s) /gridResolution_+0.5l));
+ const double startBinFraction = ((startTTimes_.at(index) / 1_s)/gridResolution_)-std::floor((startTTimes_.at(index) / 1_s)/gridResolution_);
+ const double endBinFraction = ((endTTimes_.at(index) / 1_s)/gridResolution_)-std::floor((endTTimes_.at(index) / 1_s)/gridResolution_);
// only do timing modification if contributions would land in same bin
if (startBin == endBin) {
diff --git a/tests/modules/testRadio.cpp b/tests/modules/testRadio.cpp
index 19bfa885c..0cec7d9e1 100644
--- a/tests/modules/testRadio.cpp
+++ b/tests/modules/testRadio.cpp
@@ -61,30 +61,96 @@ using namespace corsika;
double constexpr absMargin = 1.0e-7;
+template <typename TInterface>
+using MyExtraEnv =
+UniformRefractiveIndex<MediumPropertyModel<UniformMagneticField<TInterface>>>;
+
TEST_CASE("Radio", "[processes]") {
SECTION("CoREAS process") {
- // first step is to construct an environment for the propagation (uniform index 1)
- using UniRIndex =
- UniformRefractiveIndex<HomogeneousMedium<IRefractiveIndexModel<IMediumModel>>>;
- using EnvType = Environment<IRefractiveIndexModel<IMediumModel>>;
- EnvType envCoREAS;
+ // TODO: construct sychnotron radiation example with one electron
+
+// // Environment 1 (works)
+// // first step is to construct an environment for the propagation (uniform index 1)
+// using UniRIndex =
+// UniformRefractiveIndex<HomogeneousMedium<IRefractiveIndexModel<IMediumModel>>>;
+//
+// using EnvType = Environment<IRefractiveIndexModel<IMediumModel>>;
+// EnvType envCoREAS;
+//
+// // get a coordinate system
+// const CoordinateSystemPtr rootCSCoREAS = envCoREAS.getCoordinateSystem();
+//
+// auto MediumCoREAS = EnvType::createNode<Sphere>(
+// Point{rootCSCoREAS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
+//
+// auto const propsCoREAS = MediumCoREAS->setModelProperties<UniRIndex>(
+// 1.000327, 1_kg / (1_m * 1_m * 1_m),
+// NuclearComposition(
+// std::vector<Code>{Code::Nitrogen},
+// std::vector<float>{1.f}));
+//
+// envCoREAS.getUniverse()->addChild(std::move(MediumCoREAS));
- // get a coordinate system
- const CoordinateSystemPtr rootCSCoREAS = envCoREAS.getCoordinateSystem();
- auto MediumCoREAS = EnvType::createNode<Sphere>(
- Point{rootCSCoREAS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
+ //////////////////////////////////////////////////////////////////////////////////////
+// // Environment 2 (works)
+// using IModelInterface = IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
+// using AtmModel = UniformRefractiveIndex<MediumPropertyModel<UniformMagneticField<HomogeneousMedium
+// <IModelInterface>>>>;
+// using EnvType = Environment<AtmModel>;
+// EnvType envCoREAS;
+// CoordinateSystemPtr const& rootCSCoREAS = envCoREAS.getCoordinateSystem();
+// // get the center point
+// Point const center{rootCSCoREAS, 0_m, 0_m, 0_m};
+// // a refractive index
+// const double ri_{1.000327};
+//
+// // the constant density
+// const auto density{19.2_g / cube(1_cm)};
+//
+// // the composition we use for the homogeneous medium
+// NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
+// std::vector<float>{1.f});
+//
+// // create magnetic field vector
+// Vector B1(rootCSCoREAS, 0_T, 0_T, 1_T);
+//
+// auto Medium = EnvType::createNode<Sphere>(
+// center, 1_km * std::numeric_limits<double>::infinity());
+//
+// auto const props = Medium->setModelProperties<AtmModel>(ri_, Medium::AirDry1Atm, B1, density, protonComposition);
+// envCoREAS.getUniverse()->addChild(std::move(Medium));
- auto const propsZHS = MediumCoREAS->setModelProperties<UniRIndex>(
- 1, 1_kg / (1_m * 1_m * 1_m),
- NuclearComposition(
- std::vector<Code>{Code::Nitrogen},
- std::vector<float>{1.f}));
- envCoREAS.getUniverse()->addChild(std::move(MediumCoREAS));
+ //////////////////////////////////////////////////////////////////////////////////////
+ // Environment 3 (works)
+ using EnvironmentInterface =
+ IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
+ using EnvType = Environment<EnvironmentInterface>;
+ EnvType envCoREAS;
+ CoordinateSystemPtr const& rootCSCoREAS = envCoREAS.getCoordinateSystem();
+ Point const center{rootCSCoREAS, 0_m, 0_m, 0_m};
+ auto builder = make_layered_spherical_atmosphere_builder<
+ EnvironmentInterface, MyExtraEnv>::create(center,
+ constants::EarthRadius::Mean, 1.000327,
+ Medium::AirDry1Atm,
+ MagneticFieldVector{rootCSCoREAS, 0_T,
+ 50_uT, 0_T});
+
+ builder.setNuclearComposition(
+ {{Code::Nitrogen, Code::Oxygen},
+ {0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
+
+// builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km);
+// builder.addExponentialLayer(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km);
+// builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km);
+// builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km);
+ builder.addLinearLayer(1e9_cm, 112.8_km);
+ builder.assemble(envCoREAS);
+//////////////////////////////////////////////////////////////////////////////////////////
// now create antennas and detectors
@@ -96,7 +162,7 @@ TEST_CASE("Radio", "[processes]") {
// create times for the antenna
- const TimeType t1{0_s};
+ const TimeType t1{0_s}; // TODO: initialization of times to antennas! particle hits the observation level should be zero
const TimeType t2{100_s};
const InverseTimeType t3{1/1_s};
const TimeType t4{11_s};
--
GitLab