IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 55dc75f4 authored by Nikos Karastathis's avatar Nikos Karastathis :ocean:
Browse files

built compiling environments & some notes for further testing

parent 02e98678
No related branches found
No related tags found
1 merge request!329Radio interface
Hide whitespace changes
Inline Side-by-side
corsika/modules/radio/CoREAS.hpp
+ 16
− 16
View file @ 55dc75f4
......@@ -192,15 +192,15 @@ namespace corsika {
auto gridResolution_ {antenna.duration_};
auto deltaT_ { endTTimes_.at(index) - startTTimes_.at(index) };
if (fabs(deltaT_ / 1_s) < gridResolution_ / 1_s) {
if (std::fabs(deltaT_ / 1_s) < gridResolution_ / 1_s) {
EVstart_.at(index) = EVstart_.at(index) * fabs(deltaT_ / gridResolution_);
EVend_.at(index) = EVend_.at(index) * fabs(deltaT_ / gridResolution_);
EVstart_.at(index) = EVstart_.at(index) * std::fabs(deltaT_ / gridResolution_);
EVend_.at(index) = EVend_.at(index) * std::fabs(deltaT_ / gridResolution_);
const long startBin = static_cast<long>(floor(startTTimes_.at(index)/gridResolution_+0.5l));
const long endBin = static_cast<long>(floor(endTTimes_.at(index)/gridResolution_+0.5l));
const double startBinFraction = (startTTimes_.at(index)/gridResolution_)-floor(startTTimes_.at(index)/gridResolution_);
const double endBinFraction = (endTTimes_.at(index)/gridResolution_)-floor(endTTimes_.at(index)/gridResolution_);
const long startBin = static_cast<long>(std::floor(startTTimes_.at(index)/gridResolution_+0.5l));
const long endBin = static_cast<long>(std::floor(endTTimes_.at(index)/gridResolution_+0.5l));
const double startBinFraction = (startTTimes_.at(index)/gridResolution_)-std::floor(startTTimes_.at(index)/gridResolution_);
const double endBinFraction = (endTTimes_.at(index)/gridResolution_)-std::floor(endTTimes_.at(index)/gridResolution_);
// only do timing modification if contributions would land in same bin
if (startBin == endBin) {
......@@ -260,7 +260,7 @@ namespace corsika {
} // End of checking for very small doppler factors
// perform ZHS-like calculation close to Cherenkov angle
if (fabs(preDoppler__) <= approxThreshold_ || fabs(postDoppler.at(index)) <= approxThreshold_) {
if (std::fabs(preDoppler__) <= approxThreshold_ || std::fabs(postDoppler.at(index)) <= approxThreshold_) {
// get global simulation time for the middle point of that track. (This is my best guess for now)
auto midTime_{particle.getTime() - (track.getDuration() / 2)};
......@@ -296,7 +296,7 @@ namespace corsika {
EVstart_.at(index) = EVmid_;
EVend_.at(index) = - EVmid_;
auto deltaT_{(endPoint_ - startPoint_).getNorm() / (constants::c * beta_.getNorm() * fabs(midDoppler_))}; // TODO: Caution with this!
auto deltaT_{(endPoint_ - startPoint_).getNorm() / (constants::c * beta_.getNorm() * std::fabs(midDoppler_))}; // TODO: Caution with this!
if (startTTimes_.at(index) < endTTimes_.at(index)) // EVstart_ arrives earlier
{
......@@ -313,15 +313,15 @@ namespace corsika {
deltaT_ = endTTimes_.at(index) - startTTimes_.at(index);
// redistribute contributions over time scale defined by the observation time resolution
if (fabs(deltaT_ / 1_s) < gridResolution_) {
if (std::fabs(deltaT_ / 1_s) < gridResolution_) {
EVstart_.at(index) = EVstart_.at(index) * fabs((deltaT_ / 1_s) / gridResolution_);
EVend_.at(index) = EVend_.at(index) * fabs((deltaT_ / 1_s) / gridResolution_);
EVstart_.at(index) = EVstart_.at(index) * std::fabs((deltaT_ / 1_s) / gridResolution_);
EVend_.at(index) = EVend_.at(index) * std::fabs((deltaT_ / 1_s) / gridResolution_);
const long startBin = static_cast<long>(floor((startTTimes_.at(index) / 1_s)/gridResolution_+0.5l));
const long endBin = static_cast<long>(floor((endTTimes_.at(index) / 1_s) /gridResolution_+0.5l));
const double startBinFraction = ((startTTimes_.at(index) / 1_s)/gridResolution_)-floor((startTTimes_.at(index) / 1_s)/gridResolution_);
const double endBinFraction = ((endTTimes_.at(index) / 1_s)/gridResolution_)-floor((endTTimes_.at(index) / 1_s)/gridResolution_);
const long startBin = static_cast<long>(std::floor((startTTimes_.at(index) / 1_s)/gridResolution_+0.5l));
const long endBin = static_cast<long>(std::floor((endTTimes_.at(index) / 1_s) /gridResolution_+0.5l));
const double startBinFraction = ((startTTimes_.at(index) / 1_s)/gridResolution_)-std::floor((startTTimes_.at(index) / 1_s)/gridResolution_);
const double endBinFraction = ((endTTimes_.at(index) / 1_s)/gridResolution_)-std::floor((endTTimes_.at(index) / 1_s)/gridResolution_);
// only do timing modification if contributions would land in same bin
if (startBin == endBin) {
......
tests/modules/testRadio.cpp
+ 82
− 16
View file @ 55dc75f4
......@@ -61,30 +61,96 @@ 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") {
// first step is to construct an environment for the propagation (uniform index 1)
using UniRIndex =
UniformRefractiveIndex<HomogeneousMedium<IRefractiveIndexModel<IMediumModel>>>;
using EnvType = Environment<IRefractiveIndexModel<IMediumModel>>;
EnvType envCoREAS;
// TODO: construct sychnotron radiation example with one electron
// // Environment 1 (works)
// // first step is to construct an environment for the propagation (uniform index 1)
// using UniRIndex =
// UniformRefractiveIndex<HomogeneousMedium<IRefractiveIndexModel<IMediumModel>>>;
//
// using EnvType = Environment<IRefractiveIndexModel<IMediumModel>>;
// EnvType envCoREAS;
//
// // get a coordinate system
// const CoordinateSystemPtr rootCSCoREAS = envCoREAS.getCoordinateSystem();
//
// auto MediumCoREAS = EnvType::createNode<Sphere>(
// Point{rootCSCoREAS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
//
// auto const propsCoREAS = MediumCoREAS->setModelProperties<UniRIndex>(
// 1.000327, 1_kg / (1_m * 1_m * 1_m),
// NuclearComposition(
// std::vector<Code>{Code::Nitrogen},
// std::vector<float>{1.f}));
//
// envCoREAS.getUniverse()->addChild(std::move(MediumCoREAS));
// get a coordinate system
const CoordinateSystemPtr rootCSCoREAS = envCoREAS.getCoordinateSystem();
auto MediumCoREAS = EnvType::createNode<Sphere>(
Point{rootCSCoREAS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
//////////////////////////////////////////////////////////////////////////////////////
// // Environment 2 (works)
// using IModelInterface = IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
// using AtmModel = UniformRefractiveIndex<MediumPropertyModel<UniformMagneticField<HomogeneousMedium
// <IModelInterface>>>>;
// using EnvType = Environment<AtmModel>;
// EnvType envCoREAS;
// CoordinateSystemPtr const& rootCSCoREAS = envCoREAS.getCoordinateSystem();
// // get the center point
// Point const center{rootCSCoREAS, 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(rootCSCoREAS, 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);
// envCoREAS.getUniverse()->addChild(std::move(Medium));
auto const propsZHS = MediumCoREAS->setModelProperties<UniRIndex>(
1, 1_kg / (1_m * 1_m * 1_m),
NuclearComposition(
std::vector<Code>{Code::Nitrogen},
std::vector<float>{1.f}));
envCoREAS.getUniverse()->addChild(std::move(MediumCoREAS));
//////////////////////////////////////////////////////////////////////////////////////
// Environment 3 (works)
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
......@@ -96,7 +162,7 @@ TEST_CASE("Radio", "[processes]") {
// create times for the antenna
const TimeType t1{0_s};
const TimeType t1{0_s}; // TODO: initialization of times to antennas! particle hits the observation level should be zero
const TimeType t2{100_s};
const InverseTimeType t3{1/1_s};
const TimeType t4{11_s};
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

data privacy - impressum