IAP GITLAB

Skip to content
Snippets Groups Projects
Commit dce89485 authored by Maximilian Reininghaus's avatar Maximilian Reininghaus :vulcan:
Browse files

Merge gitlab.ikp.kit.edu:AirShowerPhysics/corsika

parents 41ff19e9 744b984d
No related branches found
No related tags found
No related merge requests found
......@@ -31,9 +31,9 @@ int main()
auto const norm = diff.squaredNorm(); // squared length with the right dimension
// print the components of the vector as given in the different CS
std::cout << "p2-p1 components in root: " << diff.getComponents(root) << std::endl;
std::cout << "p2-p1 components in cs2: " << diff.getComponents(cs2) << std::endl; // by definition invariant under translations
std::cout << "p2-p1 components in cs3: " << diff.getComponents(cs3) << std::endl; // but not under rotations
std::cout << "p2-p1 components in root: " << diff.GetComponents(root) << std::endl;
std::cout << "p2-p1 components in cs2: " << diff.GetComponents(cs2) << std::endl; // by definition invariant under translations
std::cout << "p2-p1 components in cs3: " << diff.GetComponents(cs3) << std::endl; // but not under rotations
std::cout << "p2-p1 norm^2: " << norm << std::endl;
Sphere s(p1, 10_m); // define a sphere around a point with a radius
......@@ -54,10 +54,10 @@ int main()
// if a CS is not given as parameter for getComponents(), the components
// in the "home" CS are returned
std::cout << "v1: " << v1.getComponents() << std::endl;
std::cout << "v2: " <<v2.getComponents() << std::endl;
std::cout << "parallel projection of v1 onto v2: " << v3.getComponents() << std::endl;
std::cout << "normalized cross product of v1 x v2" << cross.getComponents() << std::endl;
std::cout << "v1: " << v1.GetComponents() << std::endl;
std::cout << "v2: " <<v2.GetComponents() << std::endl;
std::cout << "parallel projection of v1 onto v2: " << v3.GetComponents() << std::endl;
std::cout << "normalized cross product of v1 x v2" << cross.GetComponents() << std::endl;
return EXIT_SUCCESS;
}
......@@ -32,7 +32,7 @@ int main()
for (auto [t, i] = std::tuple{t0, 0}; t < t1; t += dt, ++i)
{
auto const r = h.getPosition(t).getCoordinates();
auto const r = h.GetPosition(t).GetCoordinates();
positions[i][0] = t / 1_s;
positions[i][1] = r[0] / 1_m;
......
......@@ -2,25 +2,44 @@
namespace cascade;
template<typename Sequence, typename Trajectory>
void
Cascade::Cascade()
{
}
template<typename Sequence, typename Trajectory>
void
Cascade::Process()
Cascade::Init()
{
Stack s;
if (!s.IsEmpty()) {
fStack.Init();
fProcesseList.Init();
}
s
template<typename Sequence, typename Trajectory>
void
Cascade::Run()
{
if (!fStack.IsEmpty()) {
if (!fStack.IsEmpty()) {
Particle& p = fStack.GetNextParticle();
Step(p);
}
// do cascade equations, which can put new particles on Stack,
// thus, the double loop
// DoCascadeEquations(); //
}
}
template<typename Trajectory>
template<typename Sequence, typename Trajectory>
void
Cascade::Step(auto& sequence, Particle& particle)
Cascade::Step(Particle& particle)
{
double nextStep = sequence.MinStepLength(particle);
Trajectory trajectory = sequence.Transport(particle, nextStep);
double nextStep = fProcesseList.MinStepLength(particle);
Trajectory trajectory = fProcesseList.Transport(particle, nextStep);
sequence.DoContinuous(particle, trajectory);
sequence.DoDiscrete(particle);
}
......
#ifndef _include_Cascade_h_
#define _include_Cascade_h_
namespace cascade {
template<typename Processes, typename Trajectory, typename Stack>
class Cascade {
public:
Cascade();
void Init();
void Run();
void Step(Particle& particle);
private:
Stack fStack;
Processes fProcesseList;
};
}
#endif
#include <Geometry/CoordinateSystem.h>
EigenTransform CoordinateSystem::getTransformation(CoordinateSystem const& c1, CoordinateSystem const& c2)
EigenTransform CoordinateSystem::GetTransformation(CoordinateSystem const& c1, CoordinateSystem const& c2)
{
CoordinateSystem const* a{&c1};
CoordinateSystem const* b{&c2};
......@@ -12,13 +12,13 @@ EigenTransform CoordinateSystem::getTransformation(CoordinateSystem const& c1, C
while (a != b && a != nullptr)
{
a = a->getReference();
a = a->GetReference();
}
if (a == b)
break;
b = b->getReference();
b = b->GetReference();
}
if (a == b && a != nullptr)
......@@ -38,16 +38,16 @@ EigenTransform CoordinateSystem::getTransformation(CoordinateSystem const& c1, C
while (p != commonBase)
{
t = p->getTransform() * t;
p = p->getReference();
t = p->GetTransform() * t;
p = p->GetReference();
}
p = &c2;
while (p != commonBase)
{
t = p->getTransform().inverse(Eigen::TransformTraits::Isometry) * t;
p = p->getReference();
t = p->GetTransform().inverse(Eigen::TransformTraits::Isometry) * t;
p = p->GetReference();
}
return t;
......
......@@ -22,7 +22,7 @@ class CoordinateSystem
}
public:
static EigenTransform getTransformation(CoordinateSystem const& c1, CoordinateSystem const& c2);
static EigenTransform GetTransformation(CoordinateSystem const& c1, CoordinateSystem const& c2);
CoordinateSystem() : // for creating the root CS
transf(EigenTransform::Identity())
......@@ -58,12 +58,12 @@ public:
return CoordinateSystem(*this, transf);
}
auto const* getReference() const
auto const* GetReference() const
{
return reference;
}
auto const& getTransform() const
auto const& GetTransform() const
{
return transf;
}
......
......@@ -27,12 +27,12 @@ public:
}
auto getPosition(Time t) const
auto GetPosition(Time t) const
{
return r0 + vPar * t + (vPerp * (cos(omegaC * t) - 1) + uPerp * sin(omegaC * t)) / omegaC;
}
auto getRadius() const
auto GetRadius() const
{
return radius;
}
......
......@@ -17,7 +17,7 @@ class LineTrajectory // TODO: inherit from Trajectory
{
}
auto getPosition(Time t) const
auto GetPosition(Time t) const
{
return r0 + v0 * t;
}
......
......@@ -23,12 +23,12 @@ public:
{
}
auto getCoordinates() const
auto GetCoordinates() const
{
return BaseVector<phys::units::length_d>::qVector;
}
auto getCoordinates(CoordinateSystem const& pCS) const
auto GetCoordinates(CoordinateSystem const& pCS) const
{
if (&pCS == BaseVector<phys::units::length_d>::cs)
{
......@@ -36,7 +36,7 @@ public:
}
else
{
return QuantityVector<phys::units::length_d>(CoordinateSystem::getTransformation(*BaseVector<phys::units::length_d>::cs, pCS) * BaseVector<phys::units::length_d>::qVector.eVector);
return QuantityVector<phys::units::length_d>(CoordinateSystem::GetTransformation(*BaseVector<phys::units::length_d>::cs, pCS) * BaseVector<phys::units::length_d>::qVector.eVector);
}
}
......@@ -46,13 +46,13 @@ public:
*/
void rebase(CoordinateSystem const& pCS)
{
BaseVector<phys::units::length_d>::qVector = getCoordinates(pCS);
BaseVector<phys::units::length_d>::qVector = GetCoordinates(pCS);
BaseVector<phys::units::length_d>::cs = &pCS;
}
Point operator+(Vector<phys::units::length_d> const& pVec) const
{
return Point(*BaseVector<phys::units::length_d>::cs, getCoordinates() + pVec.getComponents(*BaseVector<phys::units::length_d>::cs));
return Point(*BaseVector<phys::units::length_d>::cs, GetCoordinates() + pVec.GetComponents(*BaseVector<phys::units::length_d>::cs));
}
/*!
......@@ -61,7 +61,7 @@ public:
Vector<phys::units::length_d> operator-(Point const& pB) const
{
auto& cs = *BaseVector<phys::units::length_d>::cs;
return Vector<phys::units::length_d>(cs, getCoordinates() - pB.getCoordinates(cs));
return Vector<phys::units::length_d>(cs, GetCoordinates() - pB.GetCoordinates(cs));
}
};
......
......@@ -35,7 +35,7 @@ public:
* returns a QuantityVector with the components given in the "home"
* CoordinateSystem of the Vector
*/
auto getComponents() const
auto GetComponents() const
{
return BaseVector<dim>::qVector;
}
......@@ -44,7 +44,7 @@ public:
* returns a QuantityVector with the components given in an arbitrary
* CoordinateSystem
*/
auto getComponents(CoordinateSystem const& pCS) const
auto GetComponents(CoordinateSystem const& pCS) const
{
if (&pCS == BaseVector<dim>::cs)
{
......@@ -52,7 +52,7 @@ public:
}
else
{
return QuantityVector<dim>(CoordinateSystem::getTransformation(*BaseVector<dim>::cs, pCS).linear() * BaseVector<dim>::qVector.eVector);
return QuantityVector<dim>(CoordinateSystem::GetTransformation(*BaseVector<dim>::cs, pCS).linear() * BaseVector<dim>::qVector.eVector);
}
}
......@@ -62,7 +62,7 @@ public:
*/
void rebase(CoordinateSystem const& pCS)
{
BaseVector<dim>::qVector = getComponents(pCS);
BaseVector<dim>::qVector = GetComponents(pCS);
BaseVector<dim>::cs = &pCS;
}
......@@ -94,8 +94,8 @@ public:
template <typename dim2>
auto parallelProjectionOnto(Vector<dim2> const& pVec, CoordinateSystem const& pCS) const
{
auto const ourCompVec = getComponents(pCS);
auto const otherCompVec = pVec.getComponents(pCS);
auto const ourCompVec = GetComponents(pCS);
auto const otherCompVec = pVec.GetComponents(pCS);
auto const& a = ourCompVec.eVector;
auto const& b = otherCompVec.eVector;
......@@ -110,13 +110,13 @@ public:
auto operator+(Vector<dim> const& pVec) const
{
auto const components = getComponents(*BaseVector<dim>::cs) + pVec.getComponents(*BaseVector<dim>::cs);
auto const components = GetComponents(*BaseVector<dim>::cs) + pVec.GetComponents(*BaseVector<dim>::cs);
return Vector<dim>(*BaseVector<dim>::cs, components);
}
auto operator-(Vector<dim> const& pVec) const
{
auto const components = getComponents() - pVec.getComponents(*BaseVector<dim>::cs);
auto const components = GetComponents() - pVec.GetComponents(*BaseVector<dim>::cs);
return Vector<dim>(*BaseVector<dim>::cs, components);
}
......@@ -159,13 +159,13 @@ public:
auto& operator+=(Vector<dim> const& pVec)
{
BaseVector<dim>::qVector += pVec.getComponents(*BaseVector<dim>::cs);
BaseVector<dim>::qVector += pVec.GetComponents(*BaseVector<dim>::cs);
return *this;
}
auto& operator-=(Vector<dim> const& pVec)
{
BaseVector<dim>::qVector -= pVec.getComponents(*BaseVector<dim>::cs);
BaseVector<dim>::qVector -= pVec.GetComponents(*BaseVector<dim>::cs);
return *this;
}
......@@ -182,8 +182,8 @@ public:
template <typename dim2>
auto cross(Vector<dim2> pV) const
{
auto const c1 = getComponents().eVector;
auto const c2 = pV.getComponents(*BaseVector<dim>::cs).eVector;
auto const c1 = GetComponents().eVector;
auto const c2 = pV.GetComponents(*BaseVector<dim>::cs).eVector;
auto const bareResult = c1.cross(c2);
using ProdQuantity = phys::units::detail::Product<dim, dim2, double, double>;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment