-
ralfulrich authored
rename framweork/sequence int framework/process remove all dynamic build files from main corsika directory, move into src re-added main/shower fixed wrong file locations Update corsika.hpp. Re-defining the version macros. Update corsika.hpp Update CONTRIBUTING.md Update MAINTAINERS.md Deleted FIXME.md Deleted AUTHORS Update COLLABORATION_AGREEMENT.md Update .clang-format Update .gitlab-ci.yml added packages: proposal, conex, cnpy, spdlog, lcov, corsika-data prevent in-source builds fotran language flags, sections build types fixed stack_example rename directories corsika8 interface target fixed few compilation, dependency issues in new system clang error simpler cmake files, using functions, re-introduce run_examples target clang format copyright notices fixed CI script, clang-format simplify ctest output better interfaces for urqmd and qgsjII, conex updated conex, data cmake integration build system cmake updates also updated corsika.hpp added earth radius moved static_pow from corsika::units::si::detail to corsika::units updated units::si namespaces throughout the project No python jobs removed a lot of stray using namespaces and corsika::units::si in entire codebase [refactory-2020] stack implementations: updated [refactory-2020] stack implementations: updated [refactory-2020] stack implementations [refactory-2020] stack implementations: !290 headers [refactory-2020] stack implementations: !290 headers
ralfulrich authoredrename framweork/sequence int framework/process remove all dynamic build files from main corsika directory, move into src re-added main/shower fixed wrong file locations Update corsika.hpp. Re-defining the version macros. Update corsika.hpp Update CONTRIBUTING.md Update MAINTAINERS.md Deleted FIXME.md Deleted AUTHORS Update COLLABORATION_AGREEMENT.md Update .clang-format Update .gitlab-ci.yml added packages: proposal, conex, cnpy, spdlog, lcov, corsika-data prevent in-source builds fotran language flags, sections build types fixed stack_example rename directories corsika8 interface target fixed few compilation, dependency issues in new system clang error simpler cmake files, using functions, re-introduce run_examples target clang format copyright notices fixed CI script, clang-format simplify ctest output better interfaces for urqmd and qgsjII, conex updated conex, data cmake integration build system cmake updates also updated corsika.hpp added earth radius moved static_pow from corsika::units::si::detail to corsika::units updated units::si namespaces throughout the project No python jobs removed a lot of stray using namespaces and corsika::units::si in entire codebase [refactory-2020] stack implementations: updated [refactory-2020] stack implementations: updated [refactory-2020] stack implementations [refactory-2020] stack implementations: !290 headers [refactory-2020] stack implementations: !290 headers
testGeometry.cpp 9.83 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 <cmath>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/CoordinateSystem.hpp>
#include <corsika/framework/geometry/Helix.hpp>
#include <corsika/framework/geometry/Line.hpp>
#include <corsika/framework/geometry/Point.hpp>
#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
#include <corsika/framework/geometry/Trajectory.hpp>
using namespace corsika;
using namespace corsika::units::si;
double constexpr absMargin = 1.0e-8;
TEST_CASE("transformations between CoordinateSystems") {
CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
REQUIRE(getTransformation(rootCS, rootCS).isApprox(EigenTransform::Identity()));
QuantityVector<length_d> const coordinates{0_m, 0_m, 0_m};
Point p1(rootCS, coordinates);
QuantityVector<magnetic_flux_density_d> components{1. * tesla, 0. * tesla, 0. * tesla};
Vector<magnetic_flux_density_d> v1(rootCS, components);
REQUIRE((p1.GetCoordinates() - coordinates).norm().magnitude() ==
Approx(0).margin(absMargin));
REQUIRE((p1.GetCoordinates(rootCS) - coordinates).norm().magnitude() ==
Approx(0).margin(absMargin));
/*
SECTION("unconnected CoordinateSystems") {
CoordinateSystem rootCS2 = CoordinateSystem::CreateRootCS();
REQUIRE_THROWS(CoordinateSystem::GetTransformation(rootCS, rootCS2));
}*/
SECTION("translations") {
QuantityVector<length_d> const translationVector{0_m, 4_m, 0_m};
CoordinateSystem translatedCS = rootCS.translate(translationVector);
REQUIRE(translatedCS.GetReference() == &rootCS);
REQUIRE((p1.GetCoordinates(translatedCS) + translationVector).norm().magnitude() ==
Approx(0).margin(absMargin));
// Vectors are not subject to translations
REQUIRE(
(v1.GetComponents(rootCS) - v1.GetComponents(translatedCS)).norm().magnitude() ==
Approx(0).margin(absMargin));
Point p2(translatedCS, {0_m, 0_m, 0_m});
REQUIRE(((p2 - p1).GetComponents() - translationVector).norm().magnitude() ==
Approx(0).margin(absMargin));
}
SECTION("multiple translations") {
QuantityVector<length_d> const tv1{0_m, 5_m, 0_m}; CoordinateSystem cs2 = rootCS.translate(tv1);
QuantityVector<length_d> const tv2{3_m, 0_m, 0_m};
CoordinateSystem cs3 = rootCS.translate(tv2);
QuantityVector<length_d> const tv3{0_m, 0_m, 2_m};
CoordinateSystem cs4 = cs3.translate(tv3);
REQUIRE(cs4.GetReference()->GetReference() == &rootCS);
REQUIRE(getTransformation(cs3, cs2).isApprox(
rootCS.translate({3_m, -5_m, 0_m}).GetTransform()));
REQUIRE(getTransformation(cs2, cs3).isApprox(
rootCS.translate({-3_m, +5_m, 0_m}).GetTransform()));
}
SECTION("rotations") {
QuantityVector<length_d> const axis{0_m, 0_m, 1_km};
double const angle = 90. / 180. * M_PI;
CoordinateSystem rotatedCS = rootCS.rotate(axis, angle);
REQUIRE(rotatedCS.GetReference() == &rootCS);
REQUIRE(v1.GetComponents(rotatedCS)[1].magnitude() ==
Approx((-1. * tesla).magnitude()));
// vector norm invariant under rotation
REQUIRE(v1.GetComponents(rotatedCS).norm().magnitude() ==
Approx(v1.GetComponents(rootCS).norm().magnitude()));
}
SECTION("multiple rotations") {
QuantityVector<length_d> const zAxis{0_m, 0_m, 1_km};
QuantityVector<length_d> const yAxis{0_m, 7_nm, 0_m};
QuantityVector<length_d> const xAxis{2_m, 0_nm, 0_m};
QuantityVector<magnetic_flux_density_d> components{1. * tesla, 2. * tesla,
3. * tesla};
Vector<magnetic_flux_density_d> v1(rootCS, components);
double const angle = 90. / 180. * M_PI;
CoordinateSystem rotated1 = rootCS.rotate(zAxis, angle);
CoordinateSystem rotated2 = rotated1.rotate(yAxis, angle);
CoordinateSystem rotated3 = rotated2.rotate(zAxis, -angle);
CoordinateSystem combined = rootCS.rotate(xAxis, -angle);
auto comp1 = v1.GetComponents(rotated3);
auto comp3 = v1.GetComponents(combined);
REQUIRE((comp1 - comp3).norm().magnitude() == Approx(0).margin(absMargin));
}
SECTION("RotateToZ positive") {
Vector const v{rootCS, 0_m, 1_m, 1_m};
auto const csPrime = rootCS.RotateToZ(v);
Vector const zPrime{csPrime, 0_m, 0_m, 5_m};
Vector const xPrime{csPrime, 5_m, 0_m, 0_m};
Vector const yPrime{csPrime, 0_m, 5_m, 0_m};
CHECK(xPrime.dot(v).magnitude() == Approx(0).margin(absMargin));
CHECK(yPrime.dot(v).magnitude() == Approx(0).margin(absMargin));
CHECK((zPrime.dot(v) / 1_m).magnitude() == Approx(5 * sqrt(2)));
CHECK(zPrime.GetComponents(rootCS)[1].magnitude() ==
Approx(zPrime.GetComponents(rootCS)[2].magnitude()));
CHECK(zPrime.GetComponents(rootCS)[0].magnitude() == Approx(0));
CHECK(xPrime.GetComponents(rootCS).eVector.dot(
yPrime.GetComponents(rootCS).eVector) == Approx(0)); CHECK(zPrime.GetComponents(rootCS).eVector.dot(
xPrime.GetComponents(rootCS).eVector) == Approx(0));
CHECK(yPrime.GetComponents(rootCS).eVector.dot(
zPrime.GetComponents(rootCS).eVector) == Approx(0));
CHECK(yPrime.GetComponents(rootCS).eVector.dot(
yPrime.GetComponents(rootCS).eVector) == Approx((5_m * 5_m).magnitude()));
CHECK(xPrime.GetComponents(rootCS).eVector.dot(
xPrime.GetComponents(rootCS).eVector) == Approx((5_m * 5_m).magnitude()));
CHECK(zPrime.GetComponents(rootCS).eVector.dot(
zPrime.GetComponents(rootCS).eVector) == Approx((5_m * 5_m).magnitude()));
}
SECTION("RotateToZ negative") {
Vector const v{rootCS, 0_m, 0_m, -1_m};
auto const csPrime = rootCS.RotateToZ(v);
Vector const zPrime{csPrime, 0_m, 0_m, 5_m};
Vector const xPrime{csPrime, 5_m, 0_m, 0_m};
Vector const yPrime{csPrime, 0_m, 5_m, 0_m};
CHECK(zPrime.dot(v).magnitude() > 0);
CHECK(xPrime.GetComponents(rootCS).eVector.dot(v.GetComponents().eVector) ==
Approx(0));
CHECK(yPrime.GetComponents(rootCS).eVector.dot(v.GetComponents().eVector) ==
Approx(0));
CHECK(xPrime.GetComponents(rootCS).eVector.dot(
yPrime.GetComponents(rootCS).eVector) == Approx(0));
CHECK(zPrime.GetComponents(rootCS).eVector.dot(
xPrime.GetComponents(rootCS).eVector) == Approx(0));
CHECK(yPrime.GetComponents(rootCS).eVector.dot(
zPrime.GetComponents(rootCS).eVector) == Approx(0));
CHECK(yPrime.GetComponents(rootCS).eVector.dot(
yPrime.GetComponents(rootCS).eVector) == Approx((5_m * 5_m).magnitude()));
CHECK(xPrime.GetComponents(rootCS).eVector.dot(
xPrime.GetComponents(rootCS).eVector) == Approx((5_m * 5_m).magnitude()));
CHECK(zPrime.GetComponents(rootCS).eVector.dot(
zPrime.GetComponents(rootCS).eVector) == Approx((5_m * 5_m).magnitude()));
}
}
TEST_CASE("Sphere") {
CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
Point center(rootCS, {0_m, 3_m, 4_m});
Sphere sphere(center, 5_m);
SECTION("GetCenter") {
CHECK((sphere.GetCenter().GetCoordinates(rootCS) -
QuantityVector<length_d>(0_m, 3_m, 4_m))
.norm()
.magnitude() == Approx(0).margin(absMargin));
CHECK(sphere.GetRadius() / 5_m == Approx(1));
}
SECTION("Contains") {
REQUIRE_FALSE(sphere.Contains(Point(rootCS, {100_m, 0_m, 0_m})));
REQUIRE(sphere.Contains(Point(rootCS, {2_m, 3_m, 4_m})));
}
}
TEST_CASE("Trajectories") {
CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
Point r0(rootCS, {0_m, 0_m, 0_m});
SECTION("Line") {
Vector<SpeedType::dimension_type> v0(rootCS,
{3_m / second, 0_m / second, 0_m / second});
Line const line(r0, v0);
CHECK(
(line.GetPosition(2_s).GetCoordinates() - QuantityVector<length_d>(6_m, 0_m, 0_m))
.norm()
.magnitude() == Approx(0).margin(absMargin));
CHECK((line.PositionFromArclength(4_m).GetCoordinates() -
QuantityVector<length_d>(4_m, 0_m, 0_m))
.norm()
.magnitude() == Approx(0).margin(absMargin));
CHECK((line.GetPosition(7_s) - line.PositionFromArclength(line.ArcLength(0_s, 7_s)))
.norm()
.magnitude() == Approx(0).margin(absMargin));
auto const t = 1_s;
LineTrajectory base(line, t);
CHECK(line.GetPosition(t).GetCoordinates() == base.GetPosition(1.).GetCoordinates());
CHECK((base.GetVelocity(0).normalized().GetComponents(rootCS) -
QuantityVector<dimensionless_d>{1, 0, 0})
.eVector.norm() == Approx(0).margin(absMargin));
}
SECTION("Helix") {
Vector<SpeedType::dimension_type> const vPar(
rootCS, {0_m / second, 0_m / second, 4_m / second});
Vector<SpeedType::dimension_type> const vPerp(
rootCS, {3_m / second, 0_m / second, 0_m / second});
auto const T = 1_s;
auto const omegaC = 2 * M_PI / T;
Helix const helix(r0, omegaC, vPar, vPerp);
CHECK((helix.GetPosition(1_s).GetCoordinates() -
QuantityVector<length_d>(0_m, 0_m, 4_m))
.norm()
.magnitude() == Approx(0).margin(absMargin));
CHECK((helix.GetPosition(0.25_s).GetCoordinates() -
QuantityVector<length_d>(-3_m / (2 * M_PI), -3_m / (2 * M_PI), 1_m))
.norm()
.magnitude() == Approx(0).margin(absMargin));
CHECK(
(helix.GetPosition(7_s) - helix.PositionFromArclength(helix.ArcLength(0_s, 7_s)))
.norm()
.magnitude() == Approx(0).margin(absMargin));
}
}