/*
* (c) Copyright 2020 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
* the license.
*/
#include <catch2/catch.hpp>
#include <corsika/modules/radio/ZHS.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 <vector>
#include <xtensor/xtensor.hpp>
#include <xtensor/xbuilder.hpp>
#include <xtensor/xio.hpp>
#include <xtensor/xcsv.hpp>
#include <istream>
#include <fstream>
#include <iostream>
#include <corsika/media/Environment.hpp>
#include <corsika/media/FlatExponential.hpp>
#include <corsika/media/HomogeneousMedium.hpp>
#include <corsika/media/IMagneticFieldModel.hpp>
#include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp>
#include <corsika/media/NuclearComposition.hpp>
#include <corsika/media/MediumPropertyModel.hpp>
#include <corsika/media/UniformMagneticField.hpp>
#include <corsika/media/SlidingPlanarExponential.hpp>
#include <corsika/media/Environment.hpp>
#include <corsika/media/HomogeneousMedium.hpp>
#include <corsika/media/IMediumModel.hpp>
#include <corsika/media/IRefractiveIndexModel.hpp>
#include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp>
#include <corsika/media/UniformRefractiveIndex.hpp>
#include <corsika/media/ExponentialRefractiveIndex.hpp>
#include <corsika/media/VolumeTreeNode.hpp>
#include <corsika/framework/geometry/CoordinateSystem.hpp>
#include <corsika/framework/geometry/Line.hpp>
#include <corsika/framework/geometry/Point.hpp>
#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>
#include <corsika/media/UniformMagneticField.hpp>
using namespace corsika;
double constexpr absMargin = 1.0e-7;
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;
// 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 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));
// now create antennas and detectors
// the antennas location
const auto point1{Point(envCoREAS.getCoordinateSystem(), 100_m, 2_m, 3_m)};
const auto point2{Point(envCoREAS.getCoordinateSystem(), 4_m, 80_m, 6_m)};
const auto point3{Point(envCoREAS.getCoordinateSystem(), 7_m, 8_m, 9_m)};
const auto point4{Point(envCoREAS.getCoordinateSystem(), 5_m, 5_m, 10_m)};
// create times for the antenna
const TimeType t1{0_s};
const TimeType t2{100_s};
const InverseTimeType t3{1/1_s};
const TimeType t4{11_s};
// check that I can create an antenna at (1, 2, 3)
TimeDomainAntenna ant1("antenna_name", point1, t1, t2, t3);
TimeDomainAntenna ant2("antenna_name2", point2, 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);
// create a particle
auto const particle{Code::Electron};
// auto const particle{Code::Gamma};
const auto pmass{get_mass(particle)};
VelocityVector v0(rootCSCoREAS, {5e+2_m / second, 5e+2_m / second, 5e+2_m / second});
Vector B0(rootCSCoREAS, 5_T, 5_T, 5_T);
Line const line(point3, v0);
auto const k{1_m * ((1_m) / ((1_s * 1_s) * 1_V))};
auto const t = 10000_ms;
LeapFrogTrajectory base(point4, v0, B0, k, t);
// create a new stack for each trial
setup::Stack stack;
// construct an energy
const HEPEnergyType E0{1_TeV};
// compute the necessary momentumn
const HEPMomentumType P0{sqrt(E0 * E0 - pmass * pmass)};
// and create the momentum vector
const auto plab{MomentumVector(rootCSCoREAS, {0_GeV, 0_GeV, P0})};
// and create the location of the particle in this coordinate system
const Point pos(rootCSCoREAS, 50_m, 10_m, 80_m);
// add the particle to the stack
auto const particle1{stack.addParticle(std::make_tuple(particle, E0, plab, pos, 0_ns))};
auto const charge_ {get_charge(particle1.getPID())};
// std::cout << "charge: " << charge_ << std::endl;
// std::cout << "1 / c: " << 1. / constants::c << std::endl;
// set up a track object
// setup::Tracking tracking;
// auto startPoint_ {base.getPosition(0)};
// auto midPoint_ {base.getPosition(0.5)};
// auto endPoint_ {base.getPosition(1)};
// std::cout << "startPoint_: " << startPoint_ << std::endl;
// std::cout << "midPoint_: " << midPoint_ << std::endl;
// std::cout << "endPoint_: " << endPoint_ << std::endl;
// auto velo_ {base.getVelocity(0)};
// std::cout << "velocity: " << velo_ << std::endl;
// auto startTime_ {particle1.getTime() - base.getDuration()}; // time at the start point of the track hopefully.
// auto endTime_ {particle1.getTime()};
// std::cout << "startTime_: " << startTime_ << std::endl;
// std::cout << "endTime_: " << endTime_ << std::endl;
//
// auto beta_ {((endPoint_ - startPoint_) / (constants::c * (endTime_ - startTime_))).normalized()};
// std::cout << "beta_: " << beta_ << std::endl;
// Vector<dimensionless_d> v1(rootCSCoREAS, {0, 0, 1});
// std::cout << "v1: " << v1.getComponents() << std::endl;
//
// std::cout << "beta_.dot(v1): " << beta_.dot(v1) << std::endl;
//
// std::cout << "Pi: " << 1/M_PI << std::endl;
//
// std::cout << "speed of light: " << constants::c << std::endl;
//
// std::cout << "vacuum permitivity: " << constants::epsilonZero << std::endl;
// Create a CoREAS instance
// CoREAS<decltype(detector), decltype(StraightPropagator(envCoREAS))> coreas1(detector, envCoREAS);
// create a radio process instance using CoREAS
RadioProcess<decltype(detector), CoREAS<decltype(detector), decltype(StraightPropagator(envCoREAS))>, decltype(StraightPropagator(envCoREAS))>
coreas(detector, envCoREAS);
// check doContinuous and simulate methods
coreas.doContinuous(particle1, base);
// coreas1.simulate(particle1, base);
// check writeOutput method -> should produce 2 csv files for each antenna
coreas.writeOutput();
}
// SECTION("TimeDomainAntenna") {
//
// // create an environment so we can get a coordinate system
// using EnvType = Environment<IRefractiveIndexModel<IMediumModel>>;
// EnvType env6;
//
// using UniRIndex =
// UniformRefractiveIndex<HomogeneousMedium<IRefractiveIndexModel<IMediumModel>>>;
//
//
// // the antenna location
// const auto point1{Point(env6.getCoordinateSystem(), 1_m, 2_m, 3_m)};
// const auto point2{Point(env6.getCoordinateSystem(), 4_m, 5_m, 6_m)};
//
//
// // 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));
//
// // create times for the antenna
// const TimeType t1{10_s};
// const TimeType t2{10_s};
// const InverseTimeType t3{1/1_s};
// const TimeType t4{11_s};
//
// // check that I can create an antenna at (1, 2, 3)
// TimeDomainAntenna ant1("antenna_name", point1, t1, t2, t3);
// TimeDomainAntenna ant2("antenna_name2", point2, t4, t2, t3);
//
// // assert that the antenna name is correct
// REQUIRE(ant1.getName() == "antenna_name");
// REQUIRE(ant2.getName() == "antenna_name2");
//
// // and check that the antenna is at the right location
// REQUIRE((ant1.getLocation() - point1).getNorm() < 1e-12 * 1_m);
// REQUIRE((ant2.getLocation() - point2).getNorm() < 1e-12 * 1_m);
//
// // construct a radio detector instance to store our antennas
// AntennaCollection<TimeDomainAntenna> detector;
//
// // add this antenna to the process
// detector.addAntenna(ant1);
// detector.addAntenna(ant2);
// CHECK(detector.size() == 2);
//
// // 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};
//
// Vector<dimensionless_d> v2(rootCS6, {0, 1, 0});
// QuantityVector<ElectricFieldType::dimension_type> v22{20_V / 1_m, 20_V / 1_m, 20_V / 1_m};
//
// // use receive methods
// ant1.receive(15_s, v1, v11);
// ant2.receive(16_s, v2, v22);
//
// // use getWaveform() method
// auto [t111, E1] = ant1.getWaveform();
// CHECK(E1(5,0) - 10 == 0);
// auto [t222, E2] = ant2.getWaveform();
// CHECK(E2(5,0) -20 == 0);
//
// // use the receive method in a for loop. It works now!
// for (auto& xx : detector.getAntennas()) {
// xx.receive(15_s, v1, v11);
// }
//
// // t & E are correct! uncomment to see them
//// for (auto& xx : detector.getAntennas()) {
//// auto [t,E] = xx.getWaveform();
//// std::cout << t << std::endl;
//// std::cout << " ... "<< std::endl;
//// std::cout << E << std::endl;
//// std::cout << " ... "<< std::endl;
//// }
//
// // check output files. It works, uncomment to see.
//// int i = 1;
//// for (auto& antenna : detector.getAntennas()) {
////
//// auto [t,E] = antenna.getWaveform();
//// auto c0 = xt::hstack(xt::xtuple(t,E));
//// std::ofstream out_file ("antenna" + to_string(i) + "_output.csv");
//// xt::dump_csv(out_file, c0);
//// out_file.close();
//// ++i;
////
//// }
//
// // check reset method for antennas. Uncomment to see they are zero
//// ant1.reset();
//// ant2.reset();
////
//// std::cout << ant1.waveformE_ << std::endl;
//// std::cout << ant2.waveformE_ << std::endl;
////
//// std::cout << " ... "<< std::endl;
//// std::cout << " ... "<< std::endl;
//
// // check reset method for antenna collection. Uncomment to see they are zero
//// detector.reset();
//// for (auto& xx : detector.getAntennas()) {
//// std::cout << xx.waveformE_ << std::endl;
//// std::cout << " ... "<< std::endl;
//// }
//
//
// }
//
// // 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
//
//// std::cout << "path.total_time_: " << path.total_time_ << std::endl;
//// std::cout << "path.average_refractive_index_: " << path.average_refractive_index_ << std::endl;
//// std::cout << "path.emit_: " << path.emit_.getComponents() << std::endl;
//// std::cout << "path.R_distance_: " << path.R_distance_ << std::endl;
//
//
// }
//
// 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 );
//
// }
} // END: TEST_CASE("Radio", "[processes]")