IAP GITLAB

Commit b98680ef authored by Ralf Ulrich's avatar Ralf Ulrich

migrated Helix

parent 55748265
......@@ -17,23 +17,24 @@
namespace corsika {
Point Helix::GetPosition(TimeType t) const {
return r0 + vPar * t +
(vPerp * (cos(omegaC * t) - 1) + uPerp * sin(omegaC * t)) / omegaC;
}
LengthType Helix::getRadius() const { return radius_; }
Point Helix::PositionFromArclength(LengthType l) const {
return GetPosition(TimeFromArclength(l));
Point Helix::getPosition(TimeType t) const {
return r0_ + vPar_ * t +
(vPerp_ * (std::cos(omegaC_ * t) - 1) + uPerp_ * std::sin(omegaC_ * t)) /
omegaC_;
}
LengthType Helix::GetRadius() const { return radius; }
Point Helix::getPositionFromArclength(LengthType l) const {
return getPosition(getTimeFromArclength(l));
}
LengthType Helix::ArcLength(TimeType t1, TimeType t2) const {
return (vPar + vPerp).norm() * (t2 - t1);
LengthType Helix::getArcLength(TimeType t1, TimeType t2) const {
return (vPar_ + vPerp_).norm() * (t2 - t1);
}
TimeType Helix::TimeFromArclength(LengthType l) const {
return l / (vPar + vPerp).norm();
TimeType Helix::getTimeFromArclength(LengthType l) const {
return l / (vPar_ + vPerp_).norm();
}
} // namespace corsika
......@@ -14,8 +14,9 @@ n/*
#include <corsika/framework/geometry/Vector.hpp>
namespace corsika {
/*!
* \class Helix
* Defines a helical path
*
* A Helix is defined by the cyclotron frequency \f$ \omega_c \f$, the initial
* Point r0 and
......@@ -30,34 +31,36 @@ namespace corsika {
class Helix {
using VelocityVec = Vector<SpeedType::dimension_type> ;
Point const r0;
FrequencyType const omegaC;
VelocityVec const vPar;
VelocityVec const vPerp, uPerp;
LengthType const radius;
///! \todo move VelocityVec into PhysicalUnits
using VelocityVec = Vector<SpeedType::dimension_type>;
public:
Helix(Point const& pR0, FrequencyType pOmegaC, VelocityVec const& pvPar,
VelocityVec const& pvPerp)
: r0(pR0)
, omegaC(pOmegaC)
, vPar(pvPar)
, vPerp(pvPerp)
, uPerp(vPerp.cross(vPar.normalized()))
, radius(pvPar.norm() / abs(pOmegaC)) {}
: r0_(pR0)
, omegaC_(pOmegaC)
, vPar_(pvPar)
, vPerp_(pvPerp)
, uPerp_(vPerp_.cross(vPar_.normalized()))
, radius_(pvPar.norm() / abs(pOmegaC)) {}
inline Point GetPosition(TimeType t) const;
inline LengthType getRadius() const;
inline Point PositionFromArclength(LengthType l) const;
inline Point getPosition(TimeType t) const;
inline Point getPositionFromArclength(LengthType l) const;
inline LengthType GetRadius() const;
inline LengthType getArcLength(TimeType t1, TimeType t2) const;
inline LengthType ArcLength(TimeType t1, TimeType t2) const;
inline TimeType getTimeFromArclength(LengthType l) const;
inline TimeType TimeFromArclength(LengthType l) const;
private:
Point r0_; ///! origin of helix, but this is in the center of the
///! "cylinder" on which the helix rotates
FrequencyType omegaC_; ///! speed of angular rotation
VelocityVec vPar_; ///! speed along direction of "cylinder"
VelocityVec vPerp_, uPerp_;
LengthType radius_;
};
} // namespace corsika
......
......@@ -6,6 +6,7 @@ set (test_framework_sources
#testCOMBoost.cpp #FIXME: reenable this
#testCorsikaFenv.cpp # does not work because of use of exceptions in catch2
testFourVector.cpp
testHelix.cpp
testFunctionTimer.cpp
testGeometry.cpp
testLogging.cpp
......
......@@ -11,7 +11,6 @@
#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>
......@@ -19,7 +18,6 @@
#include <corsika/framework/geometry/Trajectory.hpp>
using namespace corsika;
using namespace corsika::units::si;
double constexpr absMargin = 1.0e-8;
......@@ -233,31 +231,4 @@ TEST_CASE("Trajectories") {
.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));
}
}
/*
* (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/framework/geometry/Helix.hpp>
#include <corsika/framework/core/PhysicalUnits.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>
using namespace corsika;
double constexpr absMargin = 1.0e-8;
TEST_CASE("Helix class") {
CoordinateSystem& rootCS =
RootCoordinateSystem::getInstance().GetRootCoordinateSystem();
Point r0(rootCS, {0_m, 0_m, 0_m});
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.getPositionFromArclength(helix.getArcLength(0_s, 7_s)))
.norm()
.magnitude() == Approx(0).margin(absMargin));
/*
// we have to consider this, if we need it
auto const t = 1234_s;
Trajectory<Helix> const base(helix, t);
CHECK(helix.getPosition(t).GetCoordinates() == base.GetPosition(1.).GetCoordinates());
CHECK(base.ArcLength(0_s, 1_s) / 1_m == Approx(5));
*/
}
}
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