testRadio.cpp

/*
 * (c) Copyright 2022 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/RadioProcess.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/AntennaCollection.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 <boost/filesystem.hpp>
#include <vector>
#include <istream>
#include <filesystem>
#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/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 and
    // check the radio process output

    // Environment
    using EnvironmentInterface =
        IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;

    using EnvType = Environment<EnvironmentInterface>;
    //    using EnvType = setup::Environment;
    EnvType envCoREAS;
    CoordinateSystemPtr const& rootCS = envCoREAS.getCoordinateSystem();
    Point const center{rootCS, 0_m, 0_m, 0_m};

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

    // create the detector
    const auto ant1Loc{Point(rootCS, 100_m, 2_m, 3_m)};
    const auto ant2Loc{Point(rootCS, 4_m, 80_m, 6_m)};
    const TimeType t1{0_s};
    const TimeType t2{10_s};
    const InverseTimeType t3{1e+3_Hz};
    TimeDomainAntenna ant1("antenna_name", ant1Loc, rootCS, t1, t2, t3, t1);
    TimeDomainAntenna ant2("antenna_name2", ant2Loc, rootCS, t1, t2, t3, t1);
    AntennaCollection<TimeDomainAntenna> detector;
    detector.addAntenna(ant1);
    detector.addAntenna(ant2);

    const auto trackStart{Point(rootCS, 7_m, 8_m, 9_m)};
    const auto trackEnd{Point(rootCS, 5_m, 5_m, 10_m)};

    // create an electron
    const Code electron{Code::Electron};
    const auto pmass{get_mass(electron)};

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

    Vector B0(rootCS, 5_T, 5_T, 5_T);

    Line const line(trackStart, v0);

    auto const k{1_m * ((1_m) / ((1_s * 1_s) * 1_V))};

    auto const t = 1e-12_s;
    LeapFrogTrajectory base(trackEnd, v0, B0, k, t);

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

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

    // compute the necessary momentumn
    const HEPMomentumType P0{sqrt(E0 * E0 - pmass * pmass)};
    const auto plab{MomentumVector(rootCS, {0_GeV, 0_GeV, P0})};

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

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

    // create a radio process instance using CoREAS
    RadioProcess<decltype(detector),
                 CoREAS<decltype(detector), decltype(StraightPropagator(envCoREAS))>,
                 decltype(StraightPropagator(envCoREAS))>
        coreas(detector, envCoREAS);

    Step step(particle1, base);
    // check doContinuous and simulate methods
    auto const result = coreas.doContinuous(step, true);
    REQUIRE(ProcessReturn::Ok == result);

    for (auto const& ant : detector.getAntennas()) {
      // make sure something was put into the antenna
      auto totalX = ant.getWaveformX()[0];
      auto totalY = ant.getWaveformY()[0];
      auto totalZ = ant.getWaveformZ()[0];
      for (size_t i = 0; i < ant.getWaveformX().size(); i++) {
        totalX += ant.getWaveformX()[i];
        totalY += ant.getWaveformY()[i];
        totalZ += ant.getWaveformZ()[i];
      }
      REQUIRE((totalX + totalY + totalZ) > (totalX * 0));
    }

    //////////////////////////////////////
    // reset everything for a new particle
    //////////////////////////////////////
    ant1.reset();
    ant2.reset();
    stack.purge();

    // add the particle to the stack that is VERY late
    auto const particle2{stack.addParticle(std::make_tuple(
        electron, calculate_kinetic_energy(plab.getNorm(), get_mass(electron)),
        plab.normalized(), pos, t1 + t2 * 100000))};
    Step step2(particle2, base);
    auto const result2 = coreas.doContinuous(step2, true);
    REQUIRE(ProcessReturn::Ok == result2);
    for (auto const& ant : detector.getAntennas()) {
      // make sure something was put into the antenna
      auto total = ant.getWaveformX()[0];
      for (size_t i = 0; i < ant.getWaveformX().size(); i++) {
        total += ant.getWaveformX()[i] * ant.getWaveformX()[i];
        total += ant.getWaveformY()[i] * ant.getWaveformY()[i];
        total += ant.getWaveformZ()[i] * ant.getWaveformZ()[i];
      }
      REQUIRE(total < (1e-12 * ant.getWaveformX().size()));
    }

    // coreas output check
    std::string const implemencoreas{"CoREAS"};

    boost::filesystem::path const tempPathC{boost::filesystem::temp_directory_path() /
                                            ("test_corsika_radio_" + implemencoreas)};
    if (boost::filesystem::exists(tempPathC)) {
      boost::filesystem::remove_all(tempPathC);
    }

    boost::filesystem::create_directory(tempPathC);
    coreas.startOfLibrary(tempPathC);
    auto const outputFileC = tempPathC / ("antennas.parquet");
    CHECK(boost::filesystem::exists(outputFileC));
    // run end of shower and make sure that something extra was added
    auto const fileSizeC = boost::filesystem::file_size(outputFileC);
    coreas.endOfShower(0);
    CHECK(boost::filesystem::file_size(outputFileC) > fileSizeC);
    coreas.endOfLibrary();