diff --git a/corsika/modules/radio/propagators/StraightPropagator.hpp b/corsika/modules/radio/propagators/StraightPropagator.hpp index 6d3aceeaa1d0f2ce0683808a883274bfe743b64e..7bd009ca766658dc0f8f1f1e3e49ea7932a58b77 100644 --- a/corsika/modules/radio/propagators/StraightPropagator.hpp +++ b/corsika/modules/radio/propagators/StraightPropagator.hpp @@ -87,6 +87,7 @@ namespace corsika { // get the associated refractivity at 'point' auto const refractive_index{node->getModelProperties().getRefractiveIndex(point)}; +// auto const refractive_index{1.000327}; rindex.push_back(refractive_index); // add this 'point' to our deque collection @@ -96,6 +97,7 @@ namespace corsika { //add the refractive index of last point 'destination' and store it auto const* node{universe->getContainingNode(destination)}; auto const refractive_index{node->getModelProperties().getRefractiveIndex(destination)}; +// auto const refractive_index{1.000327}; rindex.push_back(refractive_index); points.push_back(destination); diff --git a/examples/radio_shower.cpp b/examples/radio_shower.cpp index 48d8b2f71cb26fcc0aa42113974a5779aab1a38a..af34300fa9cdd90fa0b83045de2e4fc973a72604 100644 --- a/examples/radio_shower.cpp +++ b/examples/radio_shower.cpp @@ -109,7 +109,7 @@ int main(int argc, char** argv) { corsika_logger->set_pattern("[%n:%^%-8l%$] %s:%#: %v"); logging::set_level(logging::level::info); - CORSIKA_LOG_INFO("vertical_EAS"); + CORSIKA_LOG_INFO("Vertical Radio Shower"); if (argc < 4) { std::cerr << "usage: vertical_EAS <A> <Z> <energy/GeV> [seed]" << std::endl; diff --git a/examples/radio_shower2.cpp b/examples/radio_shower2.cpp index 269b3dcbcacf5db524cf921ebc9b43566de49760..2135097ceb1826361076a0ea75838ab6ce2a9eef 100644 --- a/examples/radio_shower2.cpp +++ b/examples/radio_shower2.cpp @@ -21,6 +21,7 @@ #include <corsika/media/ShowerAxis.hpp> #include <corsika/media/MediumPropertyModel.hpp> #include <corsika/media/UniformMagneticField.hpp> +#include <corsika/media/UniformRefractiveIndex.hpp> #include <corsika/setup/SetupEnvironment.hpp> #include <corsika/setup/SetupStack.hpp> @@ -103,25 +104,24 @@ int main() { auto world = EnvType::createNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 150_km); - using MyHomogeneousModel = MediumPropertyModel< - UniformMagneticField<HomogeneousMedium<setup::EnvironmentInterface>>>; + using MyHomogeneousModel = UniformRefractiveIndex<MediumPropertyModel< + UniformMagneticField<HomogeneousMedium<setup::EnvironmentInterface>>>>; - world->setModelProperties<MyHomogeneousModel>( - Medium::AirDry1Atm, MagneticFieldVector(rootCS, 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.})); + world->setModelProperties<MyHomogeneousModel>(1.000327, + Medium::AirDry1Atm, MagneticFieldVector(rootCS, 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.})); universe.addChild(std::move(world)); + // setup particle stack, and add primary particle setup::Stack stack; stack.clear(); const Code beamCode = Code::Electron; const HEPMassType mass = Electron::mass; const HEPEnergyType E0 = 1000_GeV; - double theta = 0.; - double phi = 0.; Point injectionPos(rootCS, 0_m, 0_m, 0_m); { @@ -129,22 +129,16 @@ int main() { return sqrt(Elab * Elab - m * m); }; HEPMomentumType P0 = elab2plab(E0, mass); - auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) { - return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi), - -ptot * cos(theta)); - }; - auto const [px, py, pz] = - momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0); - auto plab = MomentumVector(rootCS, {px, py, pz}); + + auto plab = MomentumVector(rootCS, {0, P0, 0}); cout << "input particle: " << beamCode << endl; - cout << "input angles: theta=" << theta << " phi=" << phi << endl; cout << "input momentum: " << plab.getComponents() / 1_GeV << endl; stack.addParticle(std::make_tuple(beamCode, E0, plab, injectionPos, 0_ns)); } // setup processes, decays and interactions setup::Tracking tracking; - StackInspector<setup::Stack> stackInspect(1, true, E0); +// StackInspector<setup::Stack> stackInspect(1, true, E0); // put radio process here RadioProcess<decltype(detector), CoREAS<decltype(detector), @@ -153,10 +147,10 @@ int main() { TrackWriter trackWriter("tracks.dat"); - ShowerAxis const showerAxis{injectionPos, Vector{rootCS, 0_m, 0_m, -100_km}, env}; +// ShowerAxis const showerAxis{injectionPos, Vector{rootCS, 0_m, 0_m, -100_km}, env}; // assemble all processes into an ordered process list - auto sequence = make_sequence(coreas, trackWriter, stackInspect); + auto sequence = make_sequence(coreas, trackWriter); // define air shower object, run simulation Cascade EAS(env, tracking, sequence, stack);