testRadio.cpp

/*
 * (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 <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/IMediumModel.hpp>
#include <corsika/media/IRefractiveIndexModel.hpp>
#include <corsika/media/UniformRefractiveIndex.hpp>
#include <corsika/media/ExponentialRefractiveIndex.hpp>
#include <corsika/media/VolumeTreeNode.hpp>
#include <corsika/media/CORSIKA7Atmospheres.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/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 = setup::Environment;
    using EnvType = Environment<EnvironmentInterface>;
    EnvType envCoREAS;
    CoordinateSystemPtr const& rootCSCoREAS = envCoREAS.getCoordinateSystem();
    Point const center{rootCSCoREAS, 0_m, 0_m, 0_m};

//        1.000327,
        create_5layer_atmosphere<EnvironmentInterface, MyExtraEnv>(envCoREAS, AtmosphereId::LinsleyUSStd, center,
                                                                   1.000327, Medium::AirDry1Atm,
                                                  MagneticFieldVector{rootCSCoREAS, 0_T,
                                                                      50_uT, 0_T});



        // 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, t1);
    TimeDomainAntenna ant2("antenna_name2", point2, t1, t2, t3, t1);

    // 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,
                                                           calculate_kinetic_energy(plab.getNorm(), get_mass(particle)),
                                                           plab.normalized(), 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);
    } // END: SECTION("CoREAS process")


  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({Code::Proton}, {1.});

    // 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, t1);
    TimeDomainAntenna ant2("antenna_zhs2", point2, t1, t2, t3, t1);

    // 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, calculate_kinetic_energy(plab.getNorm(), get_mass(particle)),
                                                           plab.normalized(), pos, 0_ns))};

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

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

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

  } // END: SECTION("ZHS process")

  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({Code::Nitrogen}, {1.}));

    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, t1);
    TimeDomainAntenna ant2("antenna_name2", point2, t4, t2, t3, t4);

    // 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() methods
    auto [tx, Ex] = ant1.getWaveformX();
    CHECK(Ex[5] - 10 == 0);
    CHECK(tx[5] - 5 * 1_s / 1_ns == Approx(0.0));
    auto [ty, Ey] = ant1.getWaveformY();
    CHECK(Ey[5] - 10 == 0);
    auto [tz, Ez] = ant1.getWaveformZ();
    CHECK(Ez[5] - 10 == 0);
    CHECK(tx[5] - ty[5] == 0);
    CHECK(ty[5] - tz[5] == 0);
    auto [tx2, Ex2] = ant2.getWaveformX();
    CHECK(Ex2[5] - 20 == 0);
    auto [ty2, Ey2] = ant2.getWaveformY();
    CHECK(Ey2[5] - 20 == 0);
    auto [tz2, Ez2] = ant2.getWaveformZ();
    CHECK(Ez2[5] - 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, time__);
        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;
    }

  } // END: SECTION("TimeDomainAntenna")

  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({Code::Nitrogen}, {1.});
      // 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;}));
      }

  } // END: SECTION("Simple Propagator w/ Uniform Refractive Index")

    // 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({Code::Nitrogen}, {1.});
    // 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);
  } // END: SECTION("Straight Propagator w/ Uniform Refractive Index")

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

      // the center of the earth
      Point const center1_{rootCS1, 0_m, 0_m, 0_m};
      LengthType const radius_{0_m};

      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, center1_, radius_,
                                                1_kg / (1_m * 1_m * 1_m),
                                                NuclearComposition({Code::Nitrogen}, {1.}));

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

      // the center of the earth
      Point const center2_{rootCS2, 0_m, 0_m, 0_m};

      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, center2_, radius_,
                                                1_kg / (1_m * 1_m * 1_m),
                                                NuclearComposition({Code::Nitrogen}, {1.}));

      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: SECTION("Straight Propagator w/ Exponential Refractive Index")

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