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

#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]") {

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


    // 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
    auto const particle{Code::Electron};
    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 = 1_s;
    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, 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 = 1_s;
    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 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 CoREAS
    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 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(StraightPropagator(env))>, decltype(StraightPropagator(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});
    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);

  }

//    // 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});
//
//    // 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));
//      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(paths_.size() == 1);
//  }
//
//  SECTION("Straight Propagator w/ Exponential Refractive Index") {
//
//    // 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 P1({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)  - ((34_m / (3 * 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(P1.begin(), P1.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_ - 0.0000000041 == 0);
//      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]")