IAP GITLAB

Commit 676c72e6 authored by Remy Prechelt's avatar Remy Prechelt Committed by Ralf Ulrich

Port over refractive index models.

parent 36a1be21
/*
* (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.
*/
#pragma once
#include <corsika/media/IRefractiveIndexModel.hpp>
namespace corsika {
template <typename T>
template <typename... Args>
UniformRefractiveIndex<T>::UniformRefractiveIndex(double const n, Args&&... args)
: T(std::forward<Args>(args)...)
, n_(n) {}
template <typename T>
double UniformRefractiveIndex<T>::getRefractiveIndex(Point const&) const {
return n_;
}
template <typename T>
void UniformRefractiveIndex<T>::setRefractiveIndex(double const& n) {
n_ = n;
}
} // namespace corsika
/*
* (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.
*/
#pragma once
#include <corsika/framework/geometry/Point.hpp>
namespace corsika {
/**
* An interface for refractive index models.
*
* This is the base interface for refractive index mixins.
*
*/
template <typename TModel>
class IRefractiveIndexModel : public TModel {
public:
/**
* Evaluate the refractive index at a given location.
*
* @param point The location to evaluate at.
* @returns The refractive index at this point.
*/
virtual double getRefractiveIndex(Point const&) const = 0;
/**
* A virtual default destructor.
*/
virtual ~IRefractiveIndexModel() = default;
}; // END: class IRefractiveIndexModel
} // namespace corsika
/*
* (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.
*/
#pragma once
#include <corsika/media/IRefractiveIndexModel.hpp>
namespace corsika {
/**
* A uniform refractive index.
*
* This class returns the same refractive index
* for all evaluated locations.
*
*/
template <typename T>
class UniformRefractiveIndex final : public T {
double n_; ///< The constant refractive index that we use.
public:
/**
* Construct a UniformRefractiveIndex.
*
* This is initialized with a fixed refractive index
* and returns this refractive index at all locations.
*
* @param field The refractive index to return everywhere.
*/
template <typename... Args>
UniformRefractiveIndex(double const n, Args&&... args);
/**
* Evaluate the refractive index at a given location.
*
* @param point The location to evaluate at.
* @returns The refractive index at this point.
*/
double getRefractiveIndex(Point const&) const override;
/**
* Set the refractive index returned by this instance.
*
* @param point The location to evaluate at.
* @returns The refractive index at this location.
*/
void setRefractiveIndex(double const& n);
}; // END: class RefractiveIndex
} // namespace corsika
#include <corsika/detail/media/UniformRefractiveIndex.inl>
/*
* (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 <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/Line.hpp>
#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
#include <corsika/framework/geometry/Vector.hpp>
#include <corsika/media/FlatExponential.hpp>
#include <corsika/media/HomogeneousMedium.hpp>
#include <corsika/media/IMediumModel.hpp>
#include <corsika/media/InhomogeneousMedium.hpp>
#include <corsika/media/NuclearComposition.hpp>
#include <corsika/media/UniformRefractiveIndex.hpp>
#include <catch2/catch.hpp>
using namespace corsika;
using namespace corsika::units::si;
TEST_CASE("UniformRefractiveIndex w/ Homogeneous") {
CoordinateSystem const& gCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
Point const gOrigin(gCS, {0_m, 0_m, 0_m});
// setup our interface types
using IModelInterface = IRefractiveIndexModel<IMediumModel>;
using AtmModel = UniformRefractiveIndex<HomogeneousMedium<IModelInterface>>;
// the constant density
const auto density{19.2_g / cube(1_cm)};
// the composition we use for the homogenous medium
NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
std::vector<float>{1.f});
// the refrative index that we use
const double n{1.000327};
// create the atmospheric model
AtmModel medium(n, density, protonComposition);
// and require that it is constant
CHECK(n == medium.getRefractiveIndex(Point(gCS, -10_m, 4_m, 35_km)));
CHECK(n == medium.getRefractiveIndex(Point(gCS, +210_m, 0_m, 7_km)));
CHECK(n == medium.getRefractiveIndex(Point(gCS, 0_m, 0_m, 0_km)));
CHECK(n == medium.getRefractiveIndex(Point(gCS, 100_km, 400_km, 350_km)));
// a new refractive index
const double n2{2.3472123};
// update the refractive index of this atmospheric model
medium.setRefractiveIndex(n2);
// check that the returned refractive index is correct
CHECK(n2 == medium.getRefractiveIndex(Point(gCS, -10_m, 4_m, 35_km)));
CHECK(n2 == medium.getRefractiveIndex(Point(gCS, +210_m, 0_m, 7_km)));
CHECK(n2 == medium.getRefractiveIndex(Point(gCS, 0_m, 0_m, 0_km)));
CHECK(n2 == medium.getRefractiveIndex(Point(gCS, 100_km, 400_km, 350_km)));
// define our axis vector
Vector const axis(gCS, QuantityVector<dimensionless_d>(0, 0, 1));
// check the density and nuclear composition
CHECK(density == medium.getMassDensity(Point(gCS, 0_m, 0_m, 0_m)));
medium.getNuclearComposition();
// create a line of length 1 m
Line const line(gOrigin, Vector<SpeedType::dimension_type>(
gCS, {1_m / second, 0_m / second, 0_m / second}));
// the end time of our line
auto const tEnd = 1_s;
// and the associated trajectory
Trajectory<Line> const trajectory(line, tEnd);
// and check the integrated grammage
CHECK((medium.integratedGrammage(trajectory, 3_m) / (density * 3_m)) == Approx(1));
CHECK((medium.arclengthFromGrammage(trajectory, density * 5_m) / 5_m) == Approx(1));
}
Markdown is supported
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