IAP GITLAB

Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • AirShowerPhysics/corsika
  • rulrich/corsika
  • AAAlvesJr/corsika
  • Andre/corsika
  • arrabito/corsika
  • Nikos/corsika
  • olheiser73/corsika
  • AirShowerPhysics/papers/corsika
  • pranav/corsika
9 results
Show changes
Showing
with 0 additions and 1793 deletions
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* 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.
*/
namespace cascade;
void Cascade::Step(auto& sequence, Particle& particle) {
double nextStep = sequence.MinStepLength(particle);
Trajectory trajectory = sequence.Transport(particle, nextStep);
sequence.DoContinuous(particle, trajectory);
sequence.DoDiscrete(particle);
}
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* 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 <limits>
#include <corsika/environment/Environment.h>
#include <corsika/cascade/Cascade.h>
#include <corsika/process/ProcessSequence.h>
#include <corsika/process/stack_inspector/StackInspector.h>
#include <corsika/process/tracking_line/TrackingLine.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/geometry/Vector.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
using corsika::setup::Trajectory;
#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
// cpp file
#include <catch2/catch.hpp>
using namespace corsika;
using namespace corsika::process;
using namespace corsika::units;
using namespace corsika::units::si;
using namespace corsika::geometry;
#include <iostream>
using namespace std;
environment::Environment MakeDummyEnv() {
environment::Environment env; // dummy environment
auto& universe = *(env.GetUniverse());
auto theMedium = environment::Environment::CreateNode<Sphere>(
Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
1_km * std::numeric_limits<double>::infinity());
using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
theMedium->SetModelProperties<MyHomogeneousModel>(
1_g / (1_m * 1_m * 1_m),
environment::NuclearComposition(
std::vector<particles::Code>{particles::Code::Proton}, std::vector<float>{1.}));
universe.AddChild(std::move(theMedium));
return env;
}
class ProcessSplit : public process::ContinuousProcess<ProcessSplit> {
int fCount = 0;
int fCalls = 0;
HEPEnergyType fEcrit;
public:
ProcessSplit(HEPEnergyType e)
: fEcrit(e) {}
template <typename Particle, typename T>
LengthType MaxStepLength(Particle&, T&) const {
return 1_m;
}
template <typename Particle, typename T, typename Stack>
EProcessReturn DoContinuous(Particle& p, T&, Stack&) {
fCalls++;
HEPEnergyType E = p.GetEnergy();
if (E < fEcrit) {
p.Delete();
fCount++;
} else {
p.SetEnergy(E / 2);
p.AddSecondary(std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point,
units::si::TimeType>{p.GetPID(), E / 2, p.GetMomentum(),
p.GetPosition(), p.GetTime()});
}
return EProcessReturn::eOk;
}
void Init() {
fCount = 0;
fCalls = 0;
}
int GetCount() const { return fCount; }
int GetCalls() const { return fCalls; }
private:
};
TEST_CASE("Cascade", "[Cascade]") {
random::RNGManager& rmng = random::RNGManager::GetInstance();
rmng.RegisterRandomStream("cascade");
auto env = MakeDummyEnv();
tracking_line::TrackingLine<setup::Stack, setup::Trajectory> tracking(env);
stack_inspector::StackInspector<setup::Stack> p0(true);
const HEPEnergyType Ecrit = 85_MeV;
ProcessSplit p1(Ecrit);
auto sequence = p0 << p1;
setup::Stack stack;
cascade::Cascade EAS(env, tracking, sequence, stack);
CoordinateSystem const& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
stack.Clear();
HEPEnergyType E0 = 100_GeV;
stack.AddParticle(
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
particles::Code::Electron, E0,
corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, -1_GeV}),
Point(rootCS, {0_m, 0_m, 10_km}), 0_ns});
EAS.Init();
EAS.Run();
CHECK(p1.GetCount() == 2048);
CHECK(p1.GetCalls() == 4095);
}
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* 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.
*/
#ifndef _include_BASETRAJECTORY_H
#define _include_BASETRAJECTORY_H
#include <corsika/geometry/Point.h>
#include <corsika/geometry/Vector.h>
#include <corsika/units/PhysicalUnits.h>
#include <string>
namespace corsika::geometry {
/*!
* Interface / base class for trajectories.
*/
class BaseTrajectory {
BaseTrajectory() = delete;
public:
BaseTrajectory(corsika::units::si::TimeType start, corsika::units::si::TimeType end)
: fTStart(start)
, fTEnd(end) {}
//!< for \f$ t = 0 \f$, the starting Point shall be returned.
virtual Point GetPosition(corsika::units::si::TimeType) const = 0;
//!< the Point is return from u=0 (start) to u=1 (end)
virtual Point GetPosition(double u) const = 0;
/*!
* returns the length between two points of the trajectory
* parameterized by \arg t1 and \arg t2. Requires \arg t2 > \arg t1.
*/
virtual corsika::units::si::TimeType TimeFromArclength(
corsika::units::si::LengthType) const = 0;
virtual LengthType ArcLength(corsika::units::si::TimeType t1,
corsika::units::si::TimeType t2) const = 0;
virtual corsika::units::si::TimeType GetDuration(
corsika::units::si::TimeType t1, corsika::units::si::TimeType t2) const {
return t2 - t1;
}
virtual Point GetEndpoint() const { return GetPosition(fTEnd); }
virtual Point GetStartpoint() const { return GetPosition(fTStart); }
protected:
corsika::units::si::TimeType const fTStart, fTEnd;
};
} // namespace corsika::geometry
#endif
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* 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.
*/
#ifndef _include_BASEVECTOR_H_
#define _include_BASEVECTOR_H_
#include <corsika/geometry/CoordinateSystem.h>
#include <corsika/geometry/QuantityVector.h>
namespace corsika::geometry {
/*!
* Common base class for Vector and Point. Currently it does basically nothing.
*/
template <typename dim>
class BaseVector {
protected:
QuantityVector<dim> qVector;
CoordinateSystem const* cs;
public:
BaseVector(CoordinateSystem const& pCS, QuantityVector<dim> pQVector)
: qVector(pQVector)
, cs(&pCS) {}
auto const& GetCoordinateSystem() const { return *cs; }
};
} // namespace corsika::geometry
#endif
set (
GEOMETRY_SOURCES
CoordinateSystem.cc
)
set (
GEOMETRY_HEADERS
Vector.h
Point.h
Line.h
Sphere.h
Volume.h
CoordinateSystem.h
RootCoordinateSystem.h
Helix.h
BaseVector.h
QuantityVector.h
Trajectory.h
FourVector.h
)
set (
GEOMETRY_NAMESPACE
corsika/geometry
)
add_library (CORSIKAgeometry STATIC ${GEOMETRY_SOURCES})
CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAgeometry ${GEOMETRY_NAMESPACE} ${GEOMETRY_HEADERS})
set_target_properties (
CORSIKAgeometry
PROPERTIES
VERSION ${PROJECT_VERSION}
SOVERSION 1
PUBLIC_HEADER "${GEOMETRY_HEADERS}"
)
# target dependencies on other libraries (also the header onlys)
target_link_libraries (
CORSIKAgeometry
CORSIKAunits
)
target_include_directories (
CORSIKAgeometry
PUBLIC ${EIGEN3_INCLUDE_DIR}
INTERFACE ${EIGEN3_INCLUDE_DIR}
$<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<INSTALL_INTERFACE:include/include>
)
install (
TARGETS CORSIKAgeometry
LIBRARY DESTINATION lib
ARCHIVE DESTINATION lib
PUBLIC_HEADER DESTINATION include/${GEOMETRY_NAMESPACE}
)
# --------------------
# code unit testing
add_executable (testGeometry testGeometry.cc)
target_link_libraries (
testGeometry
CORSIKAgeometry
CORSIKAunits
CORSIKAthirdparty # for catch2
)
CORSIKA_ADD_TEST(testGeometry)
add_executable (testFourVector testFourVector.cc)
target_link_libraries (
testFourVector
CORSIKAgeometry
CORSIKAunits
CORSIKAthirdparty # for catch2
)
CORSIKA_ADD_TEST(testFourVector)
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* 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/geometry/CoordinateSystem.h>
#include <stdexcept>
using namespace corsika::geometry;
/**
* returns the transformation matrix necessary to transform primitives with coordinates
* in \a pFrom to \a pTo, e.g.
* \f$ \vec{v}^{\text{(to)}} = \mathcal{M} \vec{v}^{\text{(from)}} \f$
* (\f$ \vec{v}^{(.)} \f$ denotes the coordinates/components of the component in
* the indicated CoordinateSystem).
*/
EigenTransform CoordinateSystem::GetTransformation(CoordinateSystem const& pFrom,
CoordinateSystem const& pTo) {
CoordinateSystem const* a{&pFrom};
CoordinateSystem const* b{&pTo};
CoordinateSystem const* commonBase{nullptr};
while (a != b && b != nullptr) {
a = &pFrom;
while (a != b && a != nullptr) { a = a->GetReference(); }
if (a == b) break;
b = b->GetReference();
}
if (a == b && a != nullptr) {
commonBase = a;
} else {
throw std::runtime_error("no connection between coordinate systems found!");
}
EigenTransform t = EigenTransform::Identity();
auto* p = &pFrom;
while (p != commonBase) {
t = p->GetTransform() * t;
p = p->GetReference();
}
p = &pTo;
while (p != commonBase) {
t = t * p->GetTransform().inverse(Eigen::TransformTraits::Isometry);
p = p->GetReference();
}
return t;
}
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* 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.
*/
#ifndef _include_COORDINATESYSTEM_H_
#define _include_COORDINATESYSTEM_H_
#include <corsika/geometry/QuantityVector.h>
#include <corsika/units/PhysicalUnits.h>
#include <Eigen/Dense>
#include <stdexcept>
typedef Eigen::Transform<double, 3, Eigen::Affine> EigenTransform;
typedef Eigen::Translation<double, 3> EigenTranslation;
namespace corsika::geometry {
class RootCoordinateSystem;
using corsika::units::si::length_d;
class CoordinateSystem {
CoordinateSystem const* reference = nullptr;
EigenTransform transf;
CoordinateSystem(CoordinateSystem const& reference, EigenTransform const& transf)
: reference(&reference)
, transf(transf) {}
CoordinateSystem()
: // for creating the root CS
transf(EigenTransform::Identity()) {}
protected:
static auto CreateCS() { return CoordinateSystem(); }
friend corsika::geometry::RootCoordinateSystem; /// this is the only class that can
/// create ONE unique root CS
public:
static EigenTransform GetTransformation(CoordinateSystem const& c1,
CoordinateSystem const& c2);
auto& operator=(const CoordinateSystem& pCS) {
reference = pCS.reference;
transf = pCS.transf;
return *this;
}
auto translate(QuantityVector<length_d> vector) const {
EigenTransform const translation{EigenTranslation(vector.eVector)};
return CoordinateSystem(*this, translation);
}
auto rotate(QuantityVector<phys::units::length_d> axis, double angle) const {
if (axis.eVector.isZero()) {
throw std::runtime_error("null-vector given as axis parameter");
}
EigenTransform const rotation{Eigen::AngleAxisd(angle, axis.eVector.normalized())};
return CoordinateSystem(*this, rotation);
}
auto translateAndRotate(QuantityVector<phys::units::length_d> translation,
QuantityVector<phys::units::length_d> axis, double angle) {
if (axis.eVector.isZero()) {
throw std::runtime_error("null-vector given as axis parameter");
}
EigenTransform const transf{Eigen::AngleAxisd(angle, axis.eVector.normalized()) *
EigenTranslation(translation.eVector)};
return CoordinateSystem(*this, transf);
}
auto const* GetReference() const { return reference; }
auto const& GetTransform() const { return transf; }
};
} // namespace corsika::geometry
#endif
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* 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.
*/
#ifndef _include_corsika_framework_geometry_fourvector_h_
#define _include_corsika_framework_geometry_fourvector_h_
#include <corsika/geometry/Vector.h>
#include <corsika/units/PhysicalUnits.h>
#include <iostream>
#include <type_traits>
namespace corsika::geometry {
/**
FourVector supports "full" units, e.g. E in [GeV/c] and p in [GeV],
or also t in [s] and r in [m], etc.
However, for HEP applications it is also possible to use E and p
both in [GeV].
The FourVector can return NormSqr and Norm, whereas Norm is
sqrt(abs(NormSqr)). The physical units are always calculated and
returned properly.
FourVector can also return if it is TimeLike, SpaceLike or PhotonLike.
When a FourVector is initialized with a lvalue reference, this is
also used for the internal storage, which should lead to complete
disappearance of the FourVector class during optimization.
*/
template <typename TimeType, typename SpaceVecType>
class FourVector {
public:
using SpaceType = typename std::decay<SpaceVecType>::type::Quantity;
//! check the types and the physical units here:
static_assert(
std::is_same<typename std::decay<TimeType>::type, SpaceType>::value ||
std::is_same<typename std::decay<TimeType>::type,
decltype(std::declval<SpaceType>() / corsika::units::si::meter *
corsika::units::si::second)>::value,
"Units of time-like and space-like coordinates must either be idential "
"(e.g. GeV) or [E/c]=[p]");
public:
FourVector(const TimeType& eT, const SpaceVecType& eS)
: fTimeLike(eT)
, fSpaceLike(eS) {}
TimeType GetTimeLikeComponent() const { return fTimeLike; }
SpaceVecType& GetSpaceLikeComponents() { return fSpaceLike; }
const SpaceVecType& GetSpaceLikeComponents() const { return fSpaceLike; }
auto GetNormSqr() const { return GetTimeSquared() - fSpaceLike.squaredNorm(); }
SpaceType GetNorm() const { return sqrt(abs(GetNormSqr())); }
bool IsTimelike() const {
return GetTimeSquared() < fSpaceLike.squaredNorm();
} //! Norm2 < 0
bool IsSpacelike() const {
return GetTimeSquared() > fSpaceLike.squaredNorm();
} //! Norm2 > 0
/* this is not numerically stable
bool IsPhotonlike() const {
return GetTimeSquared() == fSpaceLike.squaredNorm();
} //! Norm2 == 0
*/
FourVector& operator+=(const FourVector& b) {
fTimeLike += b.fTimeLike;
fSpaceLike += b.fSpaceLike;
return *this;
}
FourVector& operator-=(const FourVector& b) {
fTimeLike -= b.fTimeLike;
fSpaceLike -= b.fSpaceLike;
return *this;
}
FourVector& operator*=(const double b) {
fTimeLike *= b;
fSpaceLike *= b;
return *this;
}
FourVector& operator/=(const double b) {
fTimeLike /= b;
fSpaceLike.GetComponents() /= b; // TODO: WHY IS THIS??????
return *this;
}
FourVector& operator/(const double b) {
*this /= b;
return *this;
}
/**
Note that the product between two 4-vectors assumes that you use
the same "c" convention for both. Only the LHS vector is checked
for this. You cannot mix different conventions due to
unit-checking.
*/
SpaceType operator*(const FourVector& b) {
if constexpr (std::is_same<typename std::decay<TimeType>::type,
decltype(std::declval<SpaceType>() /
corsika::units::si::meter *
corsika::units::si::second)>::value)
return fTimeLike * b.fTimeLike *
(corsika::units::constants::c * corsika::units::constants::c) -
fSpaceLike.norm();
else
return fTimeLike * fTimeLike - fSpaceLike.norm();
}
private:
/**
This function is automatically compiled to use of ignore the
extra factor of "c" for the time-like quantity
*/
auto GetTimeSquared() const {
if constexpr (std::is_same<typename std::decay<TimeType>::type,
decltype(std::declval<SpaceType>() /
corsika::units::si::meter *
corsika::units::si::second)>::value)
return fTimeLike * fTimeLike *
(corsika::units::constants::c * corsika::units::constants::c);
else
return fTimeLike * fTimeLike;
}
protected:
//! the data members
TimeType fTimeLike;
SpaceVecType fSpaceLike;
//! the friends: math operators
template <typename T, typename U>
friend FourVector<typename std::decay<T>::type, typename std::decay<U>::type>
operator+(const FourVector<T, U>&, const FourVector<T, U>&);
template <typename T, typename U>
friend FourVector<typename std::decay<T>::type, typename std::decay<U>::type>
operator-(const FourVector<T, U>&, const FourVector<T, U>&);
template <typename T, typename U>
friend FourVector<typename std::decay<T>::type, typename std::decay<U>::type>
operator*(const FourVector<T, U>&, const double);
template <typename T, typename U>
friend FourVector<typename std::decay<T>::type, typename std::decay<U>::type>
operator/(const FourVector<T, U>&, const double);
};
/**
The math operator+
*/
template <typename TimeType, typename SpaceVecType>
inline FourVector<typename std::decay<TimeType>::type,
typename std::decay<SpaceVecType>::type>
operator+(const FourVector<TimeType, SpaceVecType>& a,
const FourVector<TimeType, SpaceVecType>& b) {
return FourVector<typename std::decay<TimeType>::type,
typename std::decay<SpaceVecType>::type>(
a.fTimeLike + b.fTimeLike, a.fSpaceLike + b.fSpaceLike);
}
/**
The math operator-
*/
template <typename TimeType, typename SpaceVecType>
inline FourVector<typename std::decay<TimeType>::type,
typename std::decay<SpaceVecType>::type>
operator-(const FourVector<TimeType, SpaceVecType>& a,
const FourVector<TimeType, SpaceVecType>& b) {
return FourVector<typename std::decay<TimeType>::type,
typename std::decay<SpaceVecType>::type>(
a.fTimeLike - b.fTimeLike, a.fSpaceLike - b.fSpaceLike);
}
/**
The math operator*
*/
template <typename TimeType, typename SpaceVecType>
inline FourVector<typename std::decay<TimeType>::type,
typename std::decay<SpaceVecType>::type>
operator*(const FourVector<TimeType, SpaceVecType>& a, const double b) {
return FourVector<typename std::decay<TimeType>::type,
typename std::decay<SpaceVecType>::type>(a.fTimeLike * b,
a.fSpaceLike * b);
}
/**
The math operator/
*/
template <typename TimeType, typename SpaceVecType>
inline FourVector<typename std::decay<TimeType>::type,
typename std::decay<SpaceVecType>::type>
operator/(const FourVector<TimeType, SpaceVecType>& a, const double b) {
return FourVector<typename std::decay<TimeType>::type,
typename std::decay<SpaceVecType>::type>(a.fTimeLike / b,
a.fSpaceLike / b);
}
} // namespace corsika::geometry
#endif
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* 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.
*/
#ifndef _include_HELIX_H_
#define _include_HELIX_H_
#include <corsika/geometry/Point.h>
#include <corsika/geometry/Vector.h>
#include <corsika/units/PhysicalUnits.h>
#include <cmath>
namespace corsika::geometry {
/*!
* A Helix is defined by the cyclotron frequency \f$ \omega_c \f$, the initial
* Point r0 and
* the velocity vectors \f$ \vec{v}_{\parallel} \f$ and \f$ \vec{v}_{\perp} \f$
* denoting the projections of the initial velocity \f$ \vec{v}_0 \f$ parallel
* and perpendicular to the axis \f$ \vec{B} \f$, respectively, i.e.
* \f{align*}{
\vec{v}_{\parallel} &= \frac{\vec{v}_0 \cdot \vec{B}}{\vec{B}^2} \vec{B} \\
\vec{v}_{\perp} &= \vec{v}_0 - \vec{v}_{\parallel}
\f}
*/
class Helix {
using VelocityVec = Vector<corsika::units::si::SpeedType::dimension_type>;
Point const r0;
corsika::units::si::FrequencyType const omegaC;
VelocityVec const vPar;
VelocityVec const vPerp, uPerp;
corsika::units::si::LengthType const radius;
public:
Helix(Point const& pR0, corsika::units::si::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)) {}
Point GetPosition(corsika::units::si::TimeType t) const {
return r0 + vPar * t +
(vPerp * (cos(omegaC * t) - 1) + uPerp * sin(omegaC * t)) / omegaC;
}
Point PositionFromArclength(corsika::units::si::LengthType l) const {
return GetPosition(TimeFromArclength(l));
}
auto GetRadius() const { return radius; }
corsika::units::si::LengthType ArcLength(corsika::units::si::TimeType t1,
corsika::units::si::TimeType t2) const {
return (vPar + vPerp).norm() * (t2 - t1);
}
corsika::units::si::TimeType TimeFromArclength(
corsika::units::si::LengthType l) const {
return l / (vPar + vPerp).norm();
}
};
} // namespace corsika::geometry
#endif
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* 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.
*/
#ifndef _include_LINETRAJECTORY_H
#define _include_LINETRAJECTORY_H
#include <corsika/geometry/Point.h>
#include <corsika/geometry/Vector.h>
#include <corsika/units/PhysicalUnits.h>
namespace corsika::geometry {
class Line {
using VelocityVec = Vector<corsika::units::si::SpeedType::dimension_type>;
Point const r0;
VelocityVec const v0;
public:
Line(Point const& pR0, VelocityVec const& pV0)
: r0(pR0)
, v0(pV0) {}
Point GetPosition(corsika::units::si::TimeType t) const { return r0 + v0 * t; }
Point PositionFromArclength(corsika::units::si::LengthType l) const {
return r0 + v0.normalized() * l;
}
LengthType ArcLength(corsika::units::si::TimeType t1,
corsika::units::si::TimeType t2) const {
return v0.norm() * (t2 - t1);
}
corsika::units::si::TimeType TimeFromArclength(
corsika::units::si::LengthType t) const {
return t / v0.norm();
}
auto GetR0() const { return r0; }
auto GetV0() const { return v0; }
};
} // namespace corsika::geometry
#endif
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* 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.
*/
#ifndef _include_POINT_H_
#define _include_POINT_H_
#include <corsika/geometry/BaseVector.h>
#include <corsika/geometry/QuantityVector.h>
#include <corsika/geometry/Vector.h>
#include <corsika/units/PhysicalUnits.h>
namespace corsika::geometry {
using corsika::units::si::length_d;
using corsika::units::si::LengthType;
/*!
* A Point represents a point in position space. It is defined by its
* coordinates with respect to some CoordinateSystem.
*/
class Point : public BaseVector<length_d> {
public:
Point(CoordinateSystem const& pCS, QuantityVector<length_d> pQVector)
: BaseVector<length_d>(pCS, pQVector) {}
Point(CoordinateSystem const& cs, LengthType x, LengthType y, LengthType z)
: BaseVector<length_d>(cs, {x, y, z}) {}
// TODO: this should be private or protected, we don NOT want to expose numbers
// without reference to outside:
auto GetCoordinates() const { return BaseVector<length_d>::qVector; }
/// this always returns a QuantityVector as triple
auto GetCoordinates(CoordinateSystem const& pCS) const {
if (&pCS == BaseVector<length_d>::cs) {
return BaseVector<length_d>::qVector;
} else {
return QuantityVector<length_d>(
CoordinateSystem::GetTransformation(*BaseVector<length_d>::cs, pCS) *
BaseVector<length_d>::qVector.eVector);
}
}
/*!
* transforms the Point into another CoordinateSystem by changing its
* coordinates interally
*/
void rebase(CoordinateSystem const& pCS) {
BaseVector<length_d>::qVector = GetCoordinates(pCS);
BaseVector<length_d>::cs = &pCS;
}
Point operator+(Vector<length_d> const& pVec) const {
return Point(*BaseVector<length_d>::cs,
GetCoordinates() + pVec.GetComponents(*BaseVector<length_d>::cs));
}
/*!
* returns the distance Vector between two points
*/
Vector<length_d> operator-(Point const& pB) const {
auto& cs = *BaseVector<length_d>::cs;
return Vector<length_d>(cs, GetCoordinates() - pB.GetCoordinates(cs));
}
};
} // namespace corsika::geometry
#endif
This diff is collapsed.
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* 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.
*/
#ifndef _include_corsika_geometry_rootcoordinatesystem_h_
#define _include_corsika_geometry_rootcoordinatesystem_h_
#include <corsika/utl/Singleton.h>
#include <corsika/geometry/CoordinateSystem.h>
/*!
* This is the only way to get a root-coordinate system, and it is a
* singleton. All other CoordinateSystems must be relative to the
* RootCoordinateSystem
*/
namespace corsika::geometry {
class RootCoordinateSystem : public corsika::utl::Singleton<RootCoordinateSystem> {
friend class corsika::utl::Singleton<RootCoordinateSystem>;
protected:
RootCoordinateSystem() {}
public:
corsika::geometry::CoordinateSystem& GetRootCoordinateSystem() { return fRootCS; }
const corsika::geometry::CoordinateSystem& GetRootCoordinateSystem() const {
return fRootCS;
}
private:
corsika::geometry::CoordinateSystem fRootCS; // THIS IS IT
};
} // namespace corsika::geometry
#endif
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* 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.
*/
#ifndef _include_SPHERE_H_
#define _include_SPHERE_H_
#include <corsika/geometry/Point.h>
#include <corsika/geometry/Volume.h>
#include <corsika/units/PhysicalUnits.h>
namespace corsika::geometry {
class Sphere : public Volume {
Point const fCenter;
LengthType const fRadius;
public:
Sphere(Point const& pCenter, LengthType const pRadius)
: fCenter(pCenter)
, fRadius(pRadius) {}
//! returns true if the Point p is within the sphere
bool Contains(Point const& p) const override {
return fRadius * fRadius > (fCenter - p).squaredNorm();
}
auto& GetCenter() const { return fCenter; }
auto GetRadius() const { return fRadius; }
};
} // namespace corsika::geometry
#endif
This diff is collapsed.
This diff is collapsed.
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* 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.
*/
#ifndef _include_VOLUME_H_
#define _include_VOLUME_H_
#include <corsika/geometry/Point.h>
namespace corsika::geometry {
class Volume {
public:
//! returns true if the Point p is within the volume
virtual bool Contains(Point const& p) const = 0;
virtual ~Volume() = default;
};
} // namespace corsika::geometry
#endif
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.