/*
 * (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/SimplePropagator.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>

#include <corsika/output/OutputManager.hpp>

using namespace corsika;

double constexpr absMargin = 1.0e-7;

template <typename TInterface>
using MyExtraEnv =
UniformRefractiveIndex<MediumPropertyModel<UniformMagneticField<TInterface>>>;


TEST_CASE("Radio", "[processes]") {

logging::set_level(logging::level::debug);

  SECTION("CoREAS process") {

    // This serves as a compiler test for any changes in the CoREAS algorithm
    // Environment
    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
    // 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)};

//    auto css = point1.getCoordinateSystem();
//    std::cout << "CSS = " << css << "ENV = " << envCoREAS.getCoordinateSystem() << std::endl;


    // create times for the antenna
    const TimeType t1{0_s};
    const TimeType t2{10_s};
    const InverseTimeType t3{1e+3_Hz};
    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
    const Code particle{Code::Electron};
//    const Code particle{Code::Proton};
    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 = 1e-12_s;
    LeapFrogTrajectory base(point4, v0, B0, k, t);
//    std::cout << "Leap Frog Trajectory is: " << base << std::endl;

    // 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, plab, pos, 0_ns))};

    auto const charge_ {get_charge(particle1.getPID())};

    // 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, true);
    }


  SECTION("ZHS process") {

    // This section serves as a compiler test for any changes in the ZHS algorithm
    // Environment
    using IModelInterface = IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
    using AtmModel = UniformRefractiveIndex<MediumPropertyModel<UniformMagneticField<HomogeneousMedium
        <IModelInterface>>>>;
    using EnvType = Environment<AtmModel>;
    EnvType envZHS;
    CoordinateSystemPtr const& rootCSZHS = envZHS.getCoordinateSystem();
    // get the center point
    Point const center{rootCSZHS, 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(rootCSZHS, 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);
    envZHS.getUniverse()->addChild(std::move(Medium));

    // the antennas location
    const auto point1{Point(envZHS.getCoordinateSystem(), 100_m, 2_m, 3_m)};
    const auto point2{Point(envZHS.getCoordinateSystem(), 4_m, 80_m, 6_m)};
    const auto point3{Point(envZHS.getCoordinateSystem(), 7_m, 8_m, 9_m)};
    const auto point4{Point(envZHS.getCoordinateSystem(), 5_m, 5_m, 10_m)};

    // create times for the antenna
    const TimeType t1{0_s};
    const TimeType t2{10_s};
    const InverseTimeType t3{1e+3_Hz};
    const TimeType t4{11_s};

    // check that I can create an antenna at (1, 2, 3)
    TimeDomainAntenna ant1("antenna_zhs", point1, t1, t2, t3);
    TimeDomainAntenna ant2("antenna_zhs2", 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};
    const auto pmass{get_mass(particle)};

    VelocityVector v0(rootCSZHS, {5e+2_m / second, 5e+2_m / second, 5e+2_m / second});

    Vector B0(rootCSZHS, 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 = 1e-12_s;
    LeapFrogTrajectory base(point4, v0, B0, k, t);
//    std::cout << "Leap Frog Trajectory is: " << base << std::endl;

    // create a new stack for each trial
    setup::Stack stack;

    // construct an energy
    const HEPEnergyType E0{1_TeV};

    // compute the necessary momentum
    const HEPMomentumType P0{sqrt(E0 * E0 - pmass * pmass)};

    // and create the momentum vector
    const auto plab{MomentumVector(rootCSZHS, {0_GeV, 0_GeV, P0})};

    // and create the location of the particle in this coordinate system
    const Point pos(rootCSZHS, 50_m, 10_m, 80_m);

    // add the particle to the stack
    auto const particle1{stack.addParticle(std::make_tuple(particle, plab, pos, 0_ns))};

    auto const charge_ {get_charge(particle1.getPID())};

    // create a radio process instance using ZHS
    RadioProcess<decltype(detector), ZHS<decltype(detector), decltype(StraightPropagator(envZHS))>, decltype(StraightPropagator(envZHS))>
        zhs( detector, envZHS);

    // check doContinuous and simulate methods
    zhs.doContinuous(particle1, base, true);

  }

  SECTION("Synchrotron radiation") {

      // 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 the antenna collection
      // the antennas location
      const auto point1{Point(rootCS, 30000_m, 0_m, 0_m)};

      // create times for the antenna
      // 30 km antenna
      const TimeType start{0.994e-4_s};
      const TimeType duration{1.07e-4_s - 0.994e-4_s};
      // 3 km antenna
      //    const TimeType start{0.994e-5_s};
      //    const TimeType duration{1.7e-5_s - 0.994e-5_s};
      const InverseTimeType sampleRate_{5e+11_Hz};

      std::cout << "number of points in time: " << duration*sampleRate_ << std::endl;

      // create 4 cool antennas
      TimeDomainAntenna ant1("cool antenna", point1, start, duration, sampleRate_);

      // construct a radio detector instance to store our antennas
      AntennaCollection<TimeDomainAntenna> detector;

      // add the antennas to the detector
      detector.addAntenna(ant1);

      // create a new stack for each trial
      setup::Stack stack;
      stack.clear();

      const Code particle{Code::Electron};
//      const Code particle{Code::Proton};
      const HEPMassType pmass{get_mass(particle)};

      // construct an energy // move in the for loop
      const HEPEnergyType E0{11.4_MeV};

      // construct the output manager
      OutputManager outputs("radio_synchrotron_example");

      // create a radio process instance using CoREAS (to use ZHS simply change CoREAS with ZHS)
      RadioProcess<decltype(detector), CoREAS<decltype(detector), decltype(SimplePropagator(env))>, decltype(SimplePropagator(env))>
      coreas(detector, env);
      outputs.add("CoREAS", coreas); // register CoREAS with the output manager

      // trigger the start of the library and the first event
      outputs.startOfLibrary();
      outputs.startOfShower();

      // the number of points that make up the circle
      int const n_points {100000};
      LengthType const radius {100_m};
      TimeType timeCounter {0._s};

      // loop over all the tracks twice (this produces 2 pulses)
      for (size_t i = 0; i <= (n_points) * 2; 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)};
        timeCounter = timeCounter + t;
        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, plab, point_1, timeCounter))};
        coreas.doContinuous(particle1,track,true);
        stack.clear();
      }

      // trigger the manager to write the data to disk
      outputs.endOfShower();
      outputs.endOfLibrary();
}

  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});
    Vector<ElectricFieldType::dimension_type> v11(rootCS6, {10_V / 1_m, 10_V / 1_m, 10_V / 1_m});

    Vector<dimensionless_d> v2(rootCS6, {0, 1, 0});
    Vector<ElectricFieldType::dimension_type> v22(rootCS6, {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);

    // the following creates a star-shaped pattern of antennas in the ground
    AntennaCollection<TimeDomainAntenna> detector__;
    const auto point11{Point(env6.getCoordinateSystem(), 1000_m, 20_m, 30_m)};
    const TimeType t2222{1e-6_s};
    const InverseTimeType t3333{1e+9_Hz};

    for (auto radius_ = 100_m; radius_ <= 200_m; radius_ += 100_m) {
      for (auto phi_ = 0; phi_ <= 315; phi_ += 45) {
        auto phiRad_ = phi_ / 180. * M_PI;
        auto const point_ {Point(env6.getCoordinateSystem(), radius_ * cos(phiRad_), radius_ * sin(phiRad_), 0_m)};
        auto time__ {(point11 - point_).getNorm() / constants::c};
        const int rr_ = static_cast<int>(radius_ / 1_m);
        std::string var_ = "antenna_R=" + std::to_string(rr_) + "_m-Phi=" + std::to_string(phi_) + "degrees";
        TimeDomainAntenna ant111(var_, point_, time__, t2222, t3333);
        detector__.addAntenna(ant111);
      }
    }

    // this prints out the antenna names and locations
    for (auto const antenna : detector__.getAntennas()) {
      std::cout << antenna.getName() << " --++-- " << antenna.getLocation() << std::endl;
    }

  }

  SECTION("Simple Propagator w/ Uniform Refractive Index"){

      // 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));

      // get some points
      Point p0(rootCS, {0_m, 0_m, 0_m});
      Point p10(rootCS, {0_m, 0_m, 10_m});

      // get a unit vector
      Vector<dimensionless_d> v1(rootCS, {0, 0, 1});
      Vector<dimensionless_d> v2(rootCS, {0, 0, -1});

      // get a geometrical path of points
      Path P1({p0,p10});

      // construct a Straight Propagator given the uniform refractive index environment
      SimplePropagator 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.propagation_time_ / 1_s) - (((p10 - p0).getNorm() / constants::c) / 1_s) == Approx(0));
        CHECK(path.average_refractive_index_ == Approx(1));
        CHECK(path.refractive_index_source_ == Approx(1));
        CHECK(path.refractive_index_destination_ == Approx(1));
        CHECK(path.emit_.getComponents() == v1.getComponents());
        CHECK(path.receive_.getComponents() == v2.getComponents());
        CHECK(path.R_distance_ == 10_m);
        CHECK(std::equal(P1.begin(), P1.end(), Path(path.points_).begin(),[]
            (Point a, Point b) { return (a - b).getNorm() / 1_m < 1e-5;}));
      }

  }

    // check that I can create working Straight Propagators in different environments
  SECTION("Straight Propagator w/ Uniform Refractive Index") {

        // 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));

    // 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});
    Point p30(rootCS, {0_m, 0_m, 30000_m});

    // get a unit vector
    Vector<dimensionless_d> v1(rootCS, {0, 0, 1});
    Vector<dimensionless_d> v2(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.propagation_time_ / 1_s) - (((p10 - p0).getNorm() / constants::c) / 1_s) == Approx(0).margin(absMargin));
      CHECK(path.average_refractive_index_ == Approx(1));
      CHECK(path.refractive_index_source_ == Approx(1));
      CHECK(path.refractive_index_destination_ == Approx(1));
      CHECK(path.emit_.getComponents() == v1.getComponents());
      CHECK(path.receive_.getComponents() == v2.getComponents());
      CHECK(path.R_distance_ == 10_m);
      CHECK(std::equal(P1.begin(), P1.end(), Path(path.points_).begin(),[]
          (Point a, Point b) { return (a - b).getNorm() / 1_m < 1e-5;}));
    }

    // get another path to different points
    auto const paths2_ {SP.propagate(p0, p30, 909_m)};

    for (auto const& path : paths2_) {
      CHECK((path.propagation_time_ / 1_s) - (((p30 - p0).getNorm() / constants::c) / 1_s) == Approx(0).margin(absMargin));
      CHECK(path.average_refractive_index_ == Approx(1));
      CHECK(path.refractive_index_source_ == Approx(1));
      CHECK(path.refractive_index_destination_ == Approx(1));
      CHECK(path.R_distance_ == 30000_m);
    }

    // get a third path using a weird stepsize
    auto const paths3_ {SP.propagate(p0, p30, 731.89_m)};

    for (auto const& path : paths3_) {
      CHECK((path.propagation_time_ / 1_s) - (((p30 - p0).getNorm() / constants::c) / 1_s) == Approx(0).margin(absMargin));
      CHECK(path.average_refractive_index_ == Approx(1));
      CHECK(path.refractive_index_source_ == Approx(1));
      CHECK(path.refractive_index_destination_ == Approx(1));
      CHECK(path.R_distance_ == 30000_m);
    }

    CHECK(paths_.size() == 1);
    CHECK(paths2_.size() == 1);
    CHECK(paths3_.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});
      Vector<dimensionless_d> vv2(rootCS1, {0, 0, -1});

      // get a geometrical path of points
      Path PP1({pp0,pp1,pp2,pp3,pp4,pp5,pp6,pp7,pp8,pp9,pp10});

      // 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.propagation_time_ / 1_s) - (((pp10 - pp0).getNorm() / constants::c) / 1_s) == Approx(0).margin(absMargin));
        CHECK( path.average_refractive_index_ == Approx(1) );
        CHECK(path.refractive_index_source_ == Approx(1));
        CHECK(path.refractive_index_destination_ == Approx(1));
        CHECK( path.emit_.getComponents() == vv1.getComponents() );
        CHECK( path.receive_.getComponents() == vv2.getComponents() );
        CHECK( path.R_distance_ == 10_m );
        CHECK(std::equal(PP1.begin(), PP1.end(), Path(path.points_).begin(),[]
            (Point a, Point b) { return (a - b).getNorm() / 1_m < 1e-5;}));
      }

      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});
      Vector<dimensionless_d> vvv2(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.propagation_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.refractive_index_source_ == Approx(2));
        //      CHECK(path.refractive_index_destination_ == Approx(0.0000000041));
        CHECK( path.emit_.getComponents() == vvv1.getComponents() );
        CHECK( path.receive_.getComponents() == vvv2.getComponents() );
        CHECK( path.R_distance_ == 10_m );
      }

      CHECK( paths2_.size() == 1 );

    }

  } // END: TEST_CASE("Radio", "[processes]")