From 52add8c090b24efa6f86c83ccf0f1c3b24959500 Mon Sep 17 00:00:00 2001
From: Nikos Karastathis <n.karastathis@kit.edu>
Date: Wed, 28 Apr 2021 12:07:54 +0200
Subject: [PATCH] updates to sequence & another synchrotron implementation
---
corsika/modules/radio/CoREAS.hpp | 4 ++
corsika/modules/radio/ZHS.hpp | 27 +--------
examples/radio_shower2.cpp | 3 +-
tests/modules/testRadio.cpp | 101 ++++++++++++++++++++++++++++++-
4 files changed, 106 insertions(+), 29 deletions(-)
diff --git a/corsika/modules/radio/CoREAS.hpp b/corsika/modules/radio/CoREAS.hpp
index d6cf960b8..ef22fcd07 100755
--- a/corsika/modules/radio/CoREAS.hpp
+++ b/corsika/modules/radio/CoREAS.hpp
@@ -84,6 +84,8 @@ namespace corsika {
// loop over each antenna in the antenna collection (detector)
for (auto& antenna : detector_.getAntennas()) {
+ std::cout << "ANTENNA: " << antenna.getName() << std::endl;
+
// get the SignalPathCollection (path1) from the start "endpoint" to the antenna.
auto paths1{this->propagator_.propagate(startPoint_, antenna.getLocation(), 1_m)}; // TODO: Need to add the stepsize to .propagate()!!!!
@@ -370,6 +372,8 @@ namespace corsika {
} // End of loop over both paths to get signal info
} // End of looping over antennas
+ return ProcessReturn::Ok;
+
} // End of simulate method
}; // END: class CoREAS
diff --git a/corsika/modules/radio/ZHS.hpp b/corsika/modules/radio/ZHS.hpp
index ec75fdf02..a90cb0fba 100644
--- a/corsika/modules/radio/ZHS.hpp
+++ b/corsika/modules/radio/ZHS.hpp
@@ -185,31 +185,6 @@ namespace corsika {
return ProcessReturn::Ok;
} //end simulate
-
- /**
- * Return the maximum step length for this particle and track.
- *
- * This must be provided by the TRadioImpl.
- *
- * @param particle The current particle.
- * @param track The current track.
- *
- * @returns The maximum length of this track.
- */
- template <typename Particle, typename Track>
- LengthType MaxStepLength(Particle const& particle,
- Track const& track) const {
-
- // TODO : This is where we control the maximum step size
- // of a particle track in order to maintain the accuracy
- // of the particular formalism.
- //
- // This is part of the ZHS / CoReas formalisms and can
- // be related from the magnetic field / acceleration, charge,
- // etc. of the particle.
- return 1000000000_m;
- }
-
- }; // END: class RadioProcess
+ }; // END: class ZHS
} // namespace corsika
\ No newline at end of file
diff --git a/examples/radio_shower2.cpp b/examples/radio_shower2.cpp
index ab5ddbaab..4e07156d8 100644
--- a/examples/radio_shower2.cpp
+++ b/examples/radio_shower2.cpp
@@ -29,6 +29,7 @@
#include <corsika/modules/radio/RadioProcess.hpp>
#include <corsika/modules/radio/CoREAS.hpp>
+#include <corsika/modules/radio/ZHS.hpp>
#include <corsika/modules/radio/antennas/Antenna.hpp>
#include <corsika/modules/radio/antennas/TimeDomainAntenna.hpp>
#include <corsika/modules/radio/detectors/RadioDetector.hpp>
@@ -171,7 +172,7 @@ int main() {
TimeCut cut(period);
- TrackWriter trackWriter("tracks.dat");
+ TrackWriter trackWriter("synchrotron_tracks.dat");
// assemble all processes into an ordered process list
auto sequence = make_sequence(coreas, cut, trackWriter);
diff --git a/tests/modules/testRadio.cpp b/tests/modules/testRadio.cpp
index 081c50fcb..4ae4358e3 100644
--- a/tests/modules/testRadio.cpp
+++ b/tests/modules/testRadio.cpp
@@ -805,8 +805,6 @@ TEST_CASE("Radio", "[processes]") {
coreas.doContinuous(particle1,track,true);
}
- // ToDo: use just one electron, this way I have 400! (although it shouldn't affect the physics)
-
// get the last track
TimeType t {(points_[0] - points_[399]).getNorm() / (0.999 * constants::c)};
VelocityVector v { (points_[0] - points_[399]) / t };
@@ -823,6 +821,105 @@ TEST_CASE("Radio", "[processes]") {
}
+ SECTION("Synchrotron radiation 2") {
+
+ // create a suitable environment ///////////////////////////////////////////////////
+ using IModelInterface = IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
+ using AtmModel = UniformRefractiveIndex<MediumPropertyModel<UniformMagneticField<HomogeneousMedium
+ <IModelInterface>>>>;
+ using EnvType = Environment<AtmModel>;
+ EnvType env;
+ CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
+ // get the center point
+ Point const center{rootCS, 0_m, 0_m, 0_m};
+ // a refractive index for the vacuum
+ const double ri_{1};
+ // the constant density
+ const auto density{19.2_g / cube(1_cm)};
+ // the composition we use for the homogeneous medium
+ NuclearComposition const Composition(std::vector<Code>{Code::Nitrogen},
+ std::vector<float>{1.f});
+ // create magnetic field vector
+ Vector B1(rootCS, 0_T, 0_T, 0.3809_T);
+ // create a Sphere for the medium
+ auto Medium = EnvType::createNode<Sphere>(
+ center, 1_km * std::numeric_limits<double>::infinity());
+ // set the environment properties
+ auto const props = Medium->setModelProperties<AtmModel>(ri_, Medium::AirDry1Atm, B1, density, Composition);
+ // bind things together
+ env.getUniverse()->addChild(std::move(Medium));
+
+
+ // now create antennas and detectors/////////////////////////////////////////////
+ // the antennas location
+ const auto point1{Point(rootCS, 200_m, 0_m, 0_m)};
+// const auto point1{Point(rootCS, 30000_m, 0_m, 0_m)};
+// const auto point2{Point(rootCS, 5000_m, 100_m, 0_m)};
+// const auto point3{Point(rootCS, -100_m, -100_m, 0_m)};
+// const auto point4{Point(rootCS, -100_m, 100_m, 0_m)};
+
+ // create times for the antenna
+// const TimeType t1{0.998e-4_s};
+// const TimeType t2{1.0000e-4_s};
+// const InverseTimeType t3{1e+11_Hz};
+
+ const TimeType t1{0_s};
+ const TimeType t2{1e-6_s};
+ const InverseTimeType t3{1e+9_Hz};
+
+ // create 4 cool antennas
+ TimeDomainAntenna ant1("cool antenna", point1, t1, t2, t3);
+// TimeDomainAntenna ant2("cooler antenna", point2, t1, t2, t3);
+// TimeDomainAntenna ant3("coolest antenna", point3, t1, t2, t3);
+// TimeDomainAntenna ant4("No, I am the coolest antenna", point4, t1, t2, t3);
+
+ // construct a radio detector instance to store our antennas
+ AntennaCollection<TimeDomainAntenna> detector;
+
+ // add the antennas to the detector
+ detector.addAntenna(ant1);
+// detector.addAntenna(ant2);
+// detector.addAntenna(ant3);
+// detector.addAntenna(ant4);
+
+ //////////////////////////////////////////////////////////////////////////////////
+ // create a new stack for each trial
+ setup::Stack stack;
+ stack.clear();
+
+ const Code particle{Code::Electron};
+ const HEPMassType pmass{get_mass(particle)};
+
+ // construct an energy // move in the for loop
+ const HEPEnergyType E0{11.4_MeV};
+
+ // create a radio process instance using CoREAS or ZHS
+ RadioProcess<decltype(detector), CoREAS<decltype(detector), decltype(StraightPropagator(env))>, decltype(StraightPropagator(env))>
+ coreas(detector, env);
+
+ // loop over all the tracks except the last one
+ int const n_points {1000};
+ LengthType const radius {100_m};
+ for (size_t i = 0; i <= n_points; i++) {
+ Point const point_1(rootCS,{radius*cos(M_PI*2*i/n_points),radius*sin(M_PI*2*i/n_points), 0_m});
+ Point const point_2(rootCS,{radius*cos(M_PI*2*(i+1)/n_points),radius*sin(M_PI*2*(i+1)/n_points), 0_m});
+ TimeType t {(point_2 - point_1).getNorm() / (0.999 * constants::c)};
+ VelocityVector v { (point_2 - point_1) / t };
+ auto beta {v / constants::c};
+ auto gamma {E0/pmass};
+ auto plab {beta * pmass * gamma};
+ Line l {point_1,v};
+ StraightTrajectory track {l,t};
+ auto particle1{stack.addParticle(std::make_tuple(particle, E0, plab, point_1, t))}; //TODO: plab is inconsistent
+ coreas.doContinuous(particle1,track,true);
+ stack.clear();
+ }
+
+ // get the output
+ coreas.writeOutput();
+
+ }
+
SECTION("TimeDomainAntenna") {
--
GitLab