From 253a12254eaa667c088d063a0ab53e2c256cd3fa Mon Sep 17 00:00:00 2001
From: Nikos Karastathis <n.karastathis@kit.edu>
Date: Fri, 5 Mar 2021 21:08:16 +0100
Subject: [PATCH] first attempt of radio shower example
---
corsika/modules/radio/RadioProcess.hpp | 4 +-
examples/radio_shower.cpp | 194 ++++++++-
tests/modules/testRadio.cpp | 529 +++++++++----------------
3 files changed, 378 insertions(+), 349 deletions(-)
diff --git a/corsika/modules/radio/RadioProcess.hpp b/corsika/modules/radio/RadioProcess.hpp
index 3bdbf21ec..93e3dd75b 100644
--- a/corsika/modules/radio/RadioProcess.hpp
+++ b/corsika/modules/radio/RadioProcess.hpp
@@ -83,7 +83,9 @@ namespace corsika {
// important for controlling the runtime of radio (by ignoring particles
// that aren't going to contribute i.e. heavy hadrons)
//if (valid(particle, track)) {
- if (particle.is_em) { return this->implementation().simulate(particle, track); }
+ if (particle == Code::Electron || particle == Code::Positron) {
+ return this->implementation().simulate(particle, track);
+ }
//}
}
diff --git a/examples/radio_shower.cpp b/examples/radio_shower.cpp
index 93b9e79dd..35e162254 100644
--- a/examples/radio_shower.cpp
+++ b/examples/radio_shower.cpp
@@ -49,6 +49,15 @@
#include <corsika/modules/UrQMD.hpp>
#include <corsika/modules/PROPOSAL.hpp>
+#include <corsika/modules/radio/RadioProcess.hpp>
+#include <corsika/modules/radio/CoREAS.hpp>
+#include <corsika/modules/radio/antennas/TimeDomainAntenna.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 <corsika/setup/SetupStack.hpp>
#include <corsika/setup/SetupTrajectory.hpp>
@@ -110,15 +119,198 @@ int main(int argc, char** argv) {
// initialize random number sequence(s)
registerRandomStreams(seed);
- // setup environment (idea 2)
+ // setup environment
using EnvType = setup::Environment;
EnvType env;
CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
Point const center{rootCS, 0_m, 0_m, 0_m};
+ // the antenna location
+ const auto point1{Point(env.getCoordinateSystem(), 50_m, 50_m, 50_m)};
+ const auto point2{Point(env.getCoordinateSystem(), 25_m, 25_m, 25_m)};
+ // the antennas
+ TimeDomainAntenna ant1("antenna1", point1, 0_s, 100_s, 1/1e-8_s);
+ TimeDomainAntenna ant2("antenna2", point2, 0_s, 100_s, 1/1e-8_s);
+ // the detector
+ std::vector<TimeDomainAntenna> detector;
+ detector.push_back(ant1);
+ detector.push_back(ant2);
auto builder = make_layered_spherical_atmosphere_builder<
setup::EnvironmentInterface, MyExtraEnv>::create(center,
constants::EarthRadius::Mean, 1.,
Medium::AirDry1Atm,
MagneticFieldVector{rootCS, 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(env);
+
+ // setup particle stack, and add primary particle
+ setup::Stack stack;
+ stack.clear();
+ const Code beamCode = Code::Nucleus;
+ unsigned short const A = std::stoi(std::string(argv[1]));
+ unsigned short Z = std::stoi(std::string(argv[2]));
+ auto const mass = get_nucleus_mass(A, Z);
+ const HEPEnergyType E0 = 1_GeV * std::stof(std::string(argv[3]));
+ double theta = 0.;
+ auto const thetaRad = theta / 180. * M_PI;
+
+ auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
+ return sqrt((Elab - m) * (Elab + m));
+ };
+ HEPMomentumType P0 = elab2plab(E0, mass);
+ auto momentumComponents = [](double thetaRad, HEPMomentumType ptot) {
+ return std::make_tuple(ptot * sin(thetaRad), 0_eV, -ptot * cos(thetaRad));
+ };
+
+ auto const [px, py, pz] = momentumComponents(thetaRad, P0);
+ auto plab = MomentumVector(rootCS, {px, py, pz});
+ cout << "input particle: " << beamCode << endl;
+ cout << "input angles: theta=" << theta << endl;
+ cout << "input momentum: " << plab.getComponents() / 1_GeV
+ << ", norm = " << plab.getNorm() << endl;
+
+ auto const observationHeight = 0_km + builder.getEarthRadius();
+ auto const injectionHeight = 112.75_km + builder.getEarthRadius();
+ auto const t = -observationHeight * cos(thetaRad) +
+ sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) +
+ static_pow<2>(injectionHeight));
+ Point const showerCore{rootCS, 0_m, 0_m, observationHeight};
+ Point const injectionPos =
+ showerCore + DirectionVector{rootCS, {-sin(thetaRad), 0, cos(thetaRad)}} * t;
+
+ std::cout << "point of injection: " << injectionPos.getCoordinates() << std::endl;
+
+ if (A != 1) {
+ stack.addParticle(std::make_tuple(beamCode, E0, plab, injectionPos, 0_ns, A, Z));
+
+ } else {
+ if (Z == 1) {
+ stack.addParticle(std::make_tuple(Code::Proton, E0, plab, injectionPos, 0_ns));
+ } else if (Z == 0) {
+ stack.addParticle(std::make_tuple(Code::Neutron, E0, plab, injectionPos, 0_ns));
+ } else {
+ std::cerr << "illegal parameters" << std::endl;
+ return EXIT_FAILURE;
+ }
+ }
+
+ // we make the axis much longer than the inj-core distance since the
+ // profile will go beyond the core, depending on zenith angle
+ std::cout << "shower axis length: " << (showerCore - injectionPos).getNorm() * 1.5
+ << std::endl;
+
+ ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.5, env};
+
+ // setup processes, decays and interactions
+
+ corsika::sibyll::Interaction sibyll;
+ InteractionCounter sibyllCounted(sibyll);
+
+ corsika::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
+ InteractionCounter sibyllNucCounted(sibyllNuc);
+
+ corsika::pythia8::Decay decayPythia;
+
+ // use sibyll decay routine for decays of particles unknown to pythia
+ corsika::sibyll::Decay decaySibyll{{
+ Code::N1440Plus,
+ Code::N1440MinusBar,
+ Code::N1440_0,
+ Code::N1440_0Bar,
+ Code::N1710Plus,
+ Code::N1710MinusBar,
+ Code::N1710_0,
+ Code::N1710_0Bar,
+
+ Code::Pi1300Plus,
+ Code::Pi1300Minus,
+ Code::Pi1300_0,
+
+ Code::KStar0_1430_0,
+ Code::KStar0_1430_0Bar,
+ Code::KStar0_1430_Plus,
+ Code::KStar0_1430_MinusBar,
+ }};
+
+ decaySibyll.printDecayConfig();
+
+ ParticleCut cut{60_GeV, 60_GeV, 60_GeV, 60_GeV, true};
+ corsika::proposal::Interaction emCascade(env);
+ corsika::proposal::ContinuousProcess emContinuous(env);
+ InteractionCounter emCascadeCounted(emCascade);
+ // put radio here
+ CoREAS<TimeDomainAntenna, StraightPropagator(env)> coreas;
+
+ OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
+ TrackWriter trackWriter("tracks.dat");
+
+ LongitudinalProfile longprof{showerAxis};
+
+ Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
+ ObservationPlane observationLevel(obsPlane, DirectionVector(rootCS, {1., 0., 0.}),
+ "particles.dat");
+
+ corsika::urqmd::UrQMD urqmd;
+ InteractionCounter urqmdCounted{urqmd};
+ StackInspector<setup::Stack> stackInspect(1000, false, E0);
+
+ // assemble all processes into an ordered process list
+ struct EnergySwitch {
+ HEPEnergyType cutE_;
+ EnergySwitch(HEPEnergyType cutE)
+ : cutE_(cutE) {}
+ SwitchResult operator()(const Particle& p) {
+ if (p.getEnergy() < cutE_)
+ return SwitchResult::First;
+ else
+ return SwitchResult::Second;
+ }
+ };
+ auto hadronSequence = make_select(
+ urqmdCounted, make_sequence(sibyllNucCounted, sibyllCounted), EnergySwitch(55_GeV));
+ auto decaySequence = make_sequence(decayPythia, decaySibyll);
+ auto sequence =
+ make_sequence(stackInspect, hadronSequence, reset_particle_mass, decaySequence,
+ emContinuous, cut, coreas, trackWriter, observationLevel, longprof);
+
+ // define air shower object, run simulation
+ setup::Tracking tracking;
+ Cascade EAS(env, tracking, sequence, stack);
+
+ // to fix the point of first interaction, uncomment the following two lines:
+ // EAS.forceInteraction();
+
+ EAS.run();
+
+ cut.showResults();
+ emContinuous.showResults();
+ observationLevel.showResults();
+ const HEPEnergyType Efinal = cut.getCutEnergy() + cut.getInvEnergy() +
+ cut.getEmEnergy() + emContinuous.getEnergyLost() +
+ observationLevel.getEnergyGround();
+ cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
+ << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
+ observationLevel.reset();
+ cut.reset();
+ emContinuous.reset();
+
+ // get radio pulse
+ coreas.writeOutput();
+
+ auto const hists = sibyllCounted.getHistogram() + sibyllNucCounted.getHistogram() +
+ urqmdCounted.getHistogram();
+
+ save_hist(hists.labHist(), "inthist_lab_verticalEAS.npz", true);
+ save_hist(hists.CMSHist(), "inthist_cms_verticalEAS.npz", true);
+ longprof.save("longprof_verticalEAS.txt");
+
}
diff --git a/tests/modules/testRadio.cpp b/tests/modules/testRadio.cpp
index d761aafd8..bbc54aaa5 100644
--- a/tests/modules/testRadio.cpp
+++ b/tests/modules/testRadio.cpp
@@ -49,6 +49,7 @@
#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
#include <corsika/framework/geometry/Vector.hpp>
#include <corsika/setup/SetupStack.hpp>
+#include <corsika/setup/SetupEnvironment.hpp>
#include <corsika/setup/SetupTrajectory.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/core/PhysicalConstants.hpp>
@@ -60,48 +61,14 @@ using namespace corsika;
double constexpr absMargin = 1.0e-7;
-template <typename T>
-using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;
-
-template <typename T>
-using UniRIndex =
-UniformRefractiveIndex<HomogeneousMedium<IRefractiveIndexModel<MediumPropertyModel<UniformMagneticField<T>>>>>;
-
-
TEST_CASE("Radio", "[processes]") {
- // create an environment with uniform refractive index of 1
-
- using EnvType = Environment<IRefractiveIndexModel<IMediumModel>>;
- EnvType env9;
-
-
- using MyHomogeneousModel = MediumPropertyModel<
- UniformMagneticField<HomogeneousMedium<UniformRefractiveIndex<IRefractiveIndexModel<IMediumModel>>>>>;
-
- auto& universe = *(env9.getUniverse());
- CoordinateSystemPtr const& rootCS9 = env9.getCoordinateSystem();
-
- auto world = EnvType::createNode<Sphere>(Point{rootCS9, 0_m, 0_m, 0_m}, 150_km);
-
- world->setModelProperties<MyHomogeneousModel>(
- Medium::AirDry1Atm, MagneticFieldVector(rootCS9, 0_T, 0_T, 1_T),
- 1_kg / (1_m * 1_m * 1_m),
- NuclearComposition(std::vector<Code>{Code::Hydrogen},
- std::vector<float>{(float)1.}), 1);
-
- universe.addChild(std::move(world));
-
-// world->setModelProperties<UniRIndex>(
-// 1);
-//
-// universe.addChild(std::move(world));
-
SECTION("TimeDomainAntenna") {
// create an environment so we can get a coordinate system
- Environment<IMediumModel> env;
+ using EnvType = setup::Environment;
+ EnvType env;
// the antenna location
const auto point1{Point(env.getCoordinateSystem(), 1_m, 2_m, 3_m)};
@@ -132,138 +99,6 @@ TEST_CASE("Radio", "[processes]") {
// detector.addAntenna(ant1);
//// detector.addAntenna(ant2);
-
-// EnvType = Environment<IRefractiveIndexModel<IMediumModel>>;
-// using UniRIndex =
-// UniformRefractiveIndex<HomogeneousMedium<IRefractiveIndexModel<IMediumModel>>>;
-
-// using UniRIndex =
-// UniformRefractiveIndex<HomogeneousMedium<IRefractiveIndexModel<IMediumModel>>>;
-//
-// using EnvType = Environment<IRefractiveIndexModel<IMediumModel>>;
-// EnvType env;
-
-// 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));
-
-// using EnvType = setup::Environment;
-// EnvType env6;
-//
-// CoordinateSystemPtr const& rootCS = env6.getCoordinateSystem();
-//
-// Point const center{rootCS, 0_m, 0_m, 0_m};
-//
-// auto builder = make_layered_spherical_atmosphere_builder<
-// setup::EnvironmentInterface, UniRIndex>::create(center,
-// constants::EarthRadius::Mean,
-// Medium::AirDry1Atm,
-// 1,
-// MagneticFieldVector{rootCS, 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(env6);
-
- // setup particle stack, and add primary particle
-// setup::Stack stack;
-// stack.clear();
-// const Code beamCode = Code::Nucleus;
-// unsigned short const A = std::stoi(std::string(argv[1]));
-// unsigned short Z = std::stoi(std::string(argv[2]));
-// auto const mass = get_nucleus_mass(A, Z);
-// const HEPEnergyType E0 = 1_GeV * std::stof(std::string(argv[3]));
-// double theta = 0.;
-// auto const thetaRad = theta / 180. * M_PI;
-//
-// auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
-// return sqrt((Elab - m) * (Elab + m));
-// };
-// HEPMomentumType P0 = elab2plab(E0, mass);
-// auto momentumComponents = [](double thetaRad, HEPMomentumType ptot) {
-// return std::make_tuple(ptot * sin(thetaRad), 0_eV, -ptot * cos(thetaRad));
-// };
-//
-// auto const [px, py, pz] = momentumComponents(thetaRad, P0);
-// auto plab = MomentumVector(rootCS, {px, py, pz});
-// std::cout << "input particle: " << beamCode << std::endl;
-// std::cout << "input angles: theta=" << theta << std::endl;
-// std::cout << "input momentum: " << plab.getComponents() / 1_GeV
-// << ", norm = " << plab.getNorm() << std::endl;
-//
-// auto const observationHeight = 0_km + builder.getEarthRadius();
-// auto const injectionHeight = 112.75_km + builder.getEarthRadius();
-// auto const t = -observationHeight * cos(thetaRad) +
-// sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) +
-// static_pow<2>(injectionHeight));
-// Point const showerCore{rootCS, 0_m, 0_m, observationHeight};
-// Point const injectionPos =
-// showerCore + DirectionVector{rootCS, {-sin(thetaRad), 0, cos(thetaRad)}} * t;
-//
-// std::cout << "point of injection: " << injectionPos.getCoordinates() << std::endl;
-
-
-
-
-
-
-
-
-// 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};
@@ -361,185 +196,185 @@ TEST_CASE("Radio", "[processes]") {
}
// check that I can create working Straight Propagators in different environments
- SECTION("Straight Propagator w/ Uniform Refractive Index") {
-
- // create an environment with uniform refractive index of 1
- using UniRIndex =
- UniformRefractiveIndex<HomogeneousMedium<IRefractiveIndexModel<IMediumModel>>>;
-
- using EnvType = Environment<IRefractiveIndexModel<IMediumModel>>;
- EnvType env;
-
- // get a coordinate system
- const CoordinateSystemPtr rootCS = env.getCoordinateSystem();
-
- auto Medium = EnvType::createNode<Sphere>(
- Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
-
- auto const props = Medium->setModelProperties<UniRIndex>(
- 1, 1_kg / (1_m * 1_m * 1_m),
- NuclearComposition(
- std::vector<Code>{Code::Nitrogen},
- std::vector<float>{1.f}));
-
- env.getUniverse()->addChild(std::move(Medium));
-
- // get some points
- Point p0(rootCS, {0_m, 0_m, 0_m});
- // Point p1(rootCS, {0_m, 0_m, 1_m});
- // Point p2(rootCS, {0_m, 0_m, 2_m});
- // Point p3(rootCS, {0_m, 0_m, 3_m});
- // Point p4(rootCS, {0_m, 0_m, 4_m});
- // Point p5(rootCS, {0_m, 0_m, 5_m});
- // Point p6(rootCS, {0_m, 0_m, 6_m});
- // Point p7(rootCS, {0_m, 0_m, 7_m});
- // Point p8(rootCS, {0_m, 0_m, 8_m});
- // Point p9(rootCS, {0_m, 0_m, 9_m});
- Point p10(rootCS, {0_m, 0_m, 10_m});
-
- // get a unit vector
- Vector<dimensionless_d> v1(rootCS, {0, 0, 1});
-
- // // get a geometrical path of points
- // Path P1({p0,p1,p2,p3,p4,p5,p6,p7,p8,p9,p10});
-
- // construct a Straight Propagator given the uniform refractive index environment
- StraightPropagator SP(env);
-
- // store the outcome of the Propagate method to paths_
- auto const paths_ = SP.propagate(p0, p10, 1_m);
-
- // perform checks to paths_ components
- for (auto const& path : paths_) {
- CHECK((path.total_time_ / 1_s) - ((34_m / (3 * constants::c)) / 1_s) ==
- Approx(0).margin(absMargin));
- 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
- }
-
- CHECK(paths_.size() == 1);
- }
-
- SECTION("Straight Propagator w/ Exponential Refractive Index") {
-
-// logging::set_level(logging::level::info);
-// corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
-
- // create an environment with exponential refractive index (n_0 = 1 & lambda = 0)
- using ExpoRIndex = ExponentialRefractiveIndex<HomogeneousMedium
- <IRefractiveIndexModel<IMediumModel>>>;
-
- using EnvType = Environment<IRefractiveIndexModel<IMediumModel>>;
- EnvType env1;
-
- //get another coordinate system
- const CoordinateSystemPtr rootCS1 = env1.getCoordinateSystem();
-
- auto Medium1 = EnvType::createNode<Sphere>(
- Point{rootCS1, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
-
- auto const props1 =
- Medium1
- ->setModelProperties<ExpoRIndex>( 1, 0 / 1_m,
- 1_kg / (1_m * 1_m * 1_m),
- NuclearComposition(
- std::vector<Code>{Code::Nitrogen},
- std::vector<float>{1.f}));
-
- env1.getUniverse()->addChild(std::move(Medium1));
-
- // get some points
- Point pp0(rootCS1, {0_m, 0_m, 0_m});
-// Point pp1(rootCS1, {0_m, 0_m, 1_m});
-// Point pp2(rootCS1, {0_m, 0_m, 2_m});
-// Point pp3(rootCS1, {0_m, 0_m, 3_m});
-// Point pp4(rootCS1, {0_m, 0_m, 4_m});
-// Point pp5(rootCS1, {0_m, 0_m, 5_m});
-// Point pp6(rootCS1, {0_m, 0_m, 6_m});
-// Point pp7(rootCS1, {0_m, 0_m, 7_m});
-// Point pp8(rootCS1, {0_m, 0_m, 8_m});
-// Point pp9(rootCS1, {0_m, 0_m, 9_m});
- Point pp10(rootCS1, {0_m, 0_m, 10_m});
-
- // get a unit vector
- Vector<dimensionless_d> vv1(rootCS1, {0, 0, 1});
-
- // construct a Straight Propagator given the exponential refractive index environment
- StraightPropagator SP1(env1);
-
- // store the outcome of Propagate method to paths1_
- auto const paths1_ = SP1.propagate(pp0, pp10, 1_m);
-
- // perform checks to paths1_ components (this is just a sketch for now)
- for (auto const& path :paths1_) {
- CHECK( (path.total_time_ / 1_s) - ((34_m / (3 * constants::c)) / 1_s)
- == Approx(0).margin(absMargin) );
- 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 );
-
- /*
- * A second environment with another exponential refractive index
- */
-
- // create an environment with exponential refractive index (n_0 = 2 & lambda = 2)
- using ExpoRIndex = ExponentialRefractiveIndex<HomogeneousMedium
- <IRefractiveIndexModel<IMediumModel>>>;
-
- using EnvType = Environment<IRefractiveIndexModel<IMediumModel>>;
- EnvType env2;
-
- //get another coordinate system
- const CoordinateSystemPtr rootCS2 = env2.getCoordinateSystem();
-
- auto Medium2 = EnvType::createNode<Sphere>(
- Point{rootCS2, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
-
- auto const props2 =
- Medium2
- ->setModelProperties<ExpoRIndex>( 2, 2 / 1_m,
- 1_kg / (1_m * 1_m * 1_m),
- NuclearComposition(
- std::vector<Code>{Code::Nitrogen},
- std::vector<float>{1.f}));
-
- env2.getUniverse()->addChild(std::move(Medium2));
-
- // get some points
- Point ppp0(rootCS2, {0_m, 0_m, 0_m});
- Point ppp10(rootCS2, {0_m, 0_m, 10_m});
-
- // get a unit vector
- Vector<dimensionless_d> vvv1(rootCS2, {0, 0, 1});
-
- // construct a Straight Propagator given the exponential refractive index environment
- StraightPropagator SP2(env2);
-
- // store the outcome of Propagate method to paths1_
- auto const paths2_ = SP2.propagate(ppp0, ppp10, 1_m);
-
- // perform checks to paths1_ components (this is just a sketch for now)
- for (auto const& path :paths2_) {
- CHECK( (path.total_time_ / 1_s) - ((3.177511688_m / (3 * constants::c)) / 1_s)
- == Approx(0).margin(absMargin) );
- CHECK( path.average_refractive_index_ == Approx(0.210275935) );
- CHECK( path.emit_.getComponents() == vvv1.getComponents() );
- CHECK( path.receive_.getComponents() == vvv1.getComponents() );
- CHECK( path.R_distance_ == 10_m );
- }
-
- CHECK( paths2_.size() == 1 );
-
- }
+// SECTION("Straight Propagator w/ Uniform Refractive Index") {
+//
+// // create an environment with uniform refractive index of 1
+// using UniRIndex =
+// UniformRefractiveIndex<HomogeneousMedium<IRefractiveIndexModel<IMediumModel>>>;
+//
+// using EnvType = Environment<IRefractiveIndexModel<IMediumModel>>;
+// EnvType env;
+//
+// // get a coordinate system
+// const CoordinateSystemPtr rootCS = env.getCoordinateSystem();
+//
+// auto Medium = EnvType::createNode<Sphere>(
+// Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
+//
+// auto const props = Medium->setModelProperties<UniRIndex>(
+// 1, 1_kg / (1_m * 1_m * 1_m),
+// NuclearComposition(
+// std::vector<Code>{Code::Nitrogen},
+// std::vector<float>{1.f}));
+//
+// env.getUniverse()->addChild(std::move(Medium));
+//
+// // get some points
+// Point p0(rootCS, {0_m, 0_m, 0_m});
+// // Point p1(rootCS, {0_m, 0_m, 1_m});
+// // Point p2(rootCS, {0_m, 0_m, 2_m});
+// // Point p3(rootCS, {0_m, 0_m, 3_m});
+// // Point p4(rootCS, {0_m, 0_m, 4_m});
+// // Point p5(rootCS, {0_m, 0_m, 5_m});
+// // Point p6(rootCS, {0_m, 0_m, 6_m});
+// // Point p7(rootCS, {0_m, 0_m, 7_m});
+// // Point p8(rootCS, {0_m, 0_m, 8_m});
+// // Point p9(rootCS, {0_m, 0_m, 9_m});
+// Point p10(rootCS, {0_m, 0_m, 10_m});
+//
+// // get a unit vector
+// Vector<dimensionless_d> v1(rootCS, {0, 0, 1});
+//
+// // // get a geometrical path of points
+// // Path P1({p0,p1,p2,p3,p4,p5,p6,p7,p8,p9,p10});
+//
+// // construct a Straight Propagator given the uniform refractive index environment
+// StraightPropagator SP(env);
+//
+// // store the outcome of the Propagate method to paths_
+// auto const paths_ = SP.propagate(p0, p10, 1_m);
+//
+// // perform checks to paths_ components
+// for (auto const& path : paths_) {
+// CHECK((path.total_time_ / 1_s) - ((34_m / (3 * constants::c)) / 1_s) ==
+// Approx(0).margin(absMargin));
+// 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
+// }
+//
+// CHECK(paths_.size() == 1);
+// }
+//
+// SECTION("Straight Propagator w/ Exponential Refractive Index") {
+//
+//// logging::set_level(logging::level::info);
+//// corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
+//
+// // create an environment with exponential refractive index (n_0 = 1 & lambda = 0)
+// using ExpoRIndex = ExponentialRefractiveIndex<HomogeneousMedium
+// <IRefractiveIndexModel<IMediumModel>>>;
+//
+// using EnvType = Environment<IRefractiveIndexModel<IMediumModel>>;
+// EnvType env1;
+//
+// //get another coordinate system
+// const CoordinateSystemPtr rootCS1 = env1.getCoordinateSystem();
+//
+// auto Medium1 = EnvType::createNode<Sphere>(
+// Point{rootCS1, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
+//
+// auto const props1 =
+// Medium1
+// ->setModelProperties<ExpoRIndex>( 1, 0 / 1_m,
+// 1_kg / (1_m * 1_m * 1_m),
+// NuclearComposition(
+// std::vector<Code>{Code::Nitrogen},
+// std::vector<float>{1.f}));
+//
+// env1.getUniverse()->addChild(std::move(Medium1));
+//
+// // get some points
+// Point pp0(rootCS1, {0_m, 0_m, 0_m});
+//// Point pp1(rootCS1, {0_m, 0_m, 1_m});
+//// Point pp2(rootCS1, {0_m, 0_m, 2_m});
+//// Point pp3(rootCS1, {0_m, 0_m, 3_m});
+//// Point pp4(rootCS1, {0_m, 0_m, 4_m});
+//// Point pp5(rootCS1, {0_m, 0_m, 5_m});
+//// Point pp6(rootCS1, {0_m, 0_m, 6_m});
+//// Point pp7(rootCS1, {0_m, 0_m, 7_m});
+//// Point pp8(rootCS1, {0_m, 0_m, 8_m});
+//// Point pp9(rootCS1, {0_m, 0_m, 9_m});
+// Point pp10(rootCS1, {0_m, 0_m, 10_m});
+//
+// // get a unit vector
+// Vector<dimensionless_d> vv1(rootCS1, {0, 0, 1});
+//
+// // construct a Straight Propagator given the exponential refractive index environment
+// StraightPropagator SP1(env1);
+//
+// // store the outcome of Propagate method to paths1_
+// auto const paths1_ = SP1.propagate(pp0, pp10, 1_m);
+//
+// // perform checks to paths1_ components (this is just a sketch for now)
+// for (auto const& path :paths1_) {
+// CHECK( (path.total_time_ / 1_s) - ((34_m / (3 * constants::c)) / 1_s)
+// == Approx(0).margin(absMargin) );
+// 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 );
+//
+// /*
+// * A second environment with another exponential refractive index
+// */
+//
+// // create an environment with exponential refractive index (n_0 = 2 & lambda = 2)
+// using ExpoRIndex = ExponentialRefractiveIndex<HomogeneousMedium
+// <IRefractiveIndexModel<IMediumModel>>>;
+//
+// using EnvType = Environment<IRefractiveIndexModel<IMediumModel>>;
+// EnvType env2;
+//
+// //get another coordinate system
+// const CoordinateSystemPtr rootCS2 = env2.getCoordinateSystem();
+//
+// auto Medium2 = EnvType::createNode<Sphere>(
+// Point{rootCS2, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
+//
+// auto const props2 =
+// Medium2
+// ->setModelProperties<ExpoRIndex>( 2, 2 / 1_m,
+// 1_kg / (1_m * 1_m * 1_m),
+// NuclearComposition(
+// std::vector<Code>{Code::Nitrogen},
+// std::vector<float>{1.f}));
+//
+// env2.getUniverse()->addChild(std::move(Medium2));
+//
+// // get some points
+// Point ppp0(rootCS2, {0_m, 0_m, 0_m});
+// Point ppp10(rootCS2, {0_m, 0_m, 10_m});
+//
+// // get a unit vector
+// Vector<dimensionless_d> vvv1(rootCS2, {0, 0, 1});
+//
+// // construct a Straight Propagator given the exponential refractive index environment
+// StraightPropagator SP2(env2);
+//
+// // store the outcome of Propagate method to paths1_
+// auto const paths2_ = SP2.propagate(ppp0, ppp10, 1_m);
+//
+// // perform checks to paths1_ components (this is just a sketch for now)
+// for (auto const& path :paths2_) {
+// CHECK( (path.total_time_ / 1_s) - ((3.177511688_m / (3 * constants::c)) / 1_s)
+// == Approx(0).margin(absMargin) );
+// CHECK( path.average_refractive_index_ == Approx(0.210275935) );
+// CHECK( path.emit_.getComponents() == vvv1.getComponents() );
+// CHECK( path.receive_.getComponents() == vvv1.getComponents() );
+// CHECK( path.R_distance_ == 10_m );
+// }
+//
+// CHECK( paths2_.size() == 1 );
+//
+// }
// SECTION("ZHS process") {
// // first step is to construct an environment for the propagation (uni index)
--
GitLab