-
Nikos Karastathis authoredNikos Karastathis authored
testRadio.cpp 25.18 KiB
/*
* (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()) {
CORSIKA_LOG_DEBUG("Antenna name: {} ", antenna.getName());
CORSIKA_LOG_DEBUG("Antenna location: {} ", antenna.getLocation());
}
} // 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]")