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
Hide whitespace changes
Inline Side-by-side
Showing
with 0 additions and 1921 deletions
Framework/Cascade/CMakeLists.txt deleted 100644 → 0
View file @ dd195476
# namespace of library -> location of header files
set (
CORSIKAcascade_NAMESPACE
corsika/cascade
)
# header files of this library
set (
CORSIKAcascade_HEADERS
Cascade.h
sibyll2.3c.h
SibStack.h
)
#set (
# CORSIKAcascade_SOURCES
# TrackingStep.cc
# Cascade.cc
# )
set (
CORSIKAsibyll_NAMESPACE
corsika/cascade
)
set (
CORSIKAsibyll_SOURCES
sibyll2.3c.f
gasdev.f
)
add_library (CORSIKAsibyll STATIC ${CORSIKAsibyll_SOURCES})
#add_library (CORSIKAcascade STATIC ${CORSIKAcascade_SOURCES})
add_library (CORSIKAcascade INTERFACE)
CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAcascade ${CORSIKAcascade_NAMESPACE} ${CORSIKAcascade_HEADERS})
#target_link_libraries (
# CORSIKAcascade
# CORSIKAparticles
# CORSIKAunits
# CORSIKAthirdparty # for catch2
# )
# include directive for upstream code
target_include_directories (
CORSIKAcascade
INTERFACE
$<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<INSTALL_INTERFACE:include/>
)
# install library
install (
FILES ${CORSIKAcascade_HEADERS}
DESTINATION include/${CORSIKAcascade_NAMESPACE}
)
# ----------------
# code unit testing
add_executable (
testCascade
testCascade.cc
)
target_link_libraries (
testCascade
# CORSIKAutls
CORSIKArandom
CORSIKAsibyll
CORSIKAcascade
ProcessStackInspector
CORSIKAstackinterface
CORSIKAprocesses
CORSIKAparticles
CORSIKAgeometry
CORSIKAprocesssequence
CORSIKAunits
CORSIKAthirdparty # for catch2
)
add_test (
NAME testCascade
COMMAND testCascade
)
Framework/Cascade/Cascade.h deleted 100644 → 0
View file @ dd195476
/**
* (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_cascade_Cascade_h_
#define _include_corsika_cascade_Cascade_h_
#include <corsika/process/ProcessReturn.h>
#include <corsika/units/PhysicalUnits.h>
#include <type_traits>
#include <corsika/setup/SetupTrajectory.h>
using namespace corsika::units::si;
namespace corsika::cascade {
template <typename Tracking, typename ProcessList, typename Stack>
class Cascade {
typedef typename Stack::ParticleType Particle;
Cascade() = delete;
public:
Cascade(Tracking& tr, ProcessList& pl, Stack& stack)
: fTracking(tr)
, fProcesseList(pl)
, fStack(stack) {}
void Init() {
fTracking.Init();
fProcesseList.Init();
fStack.Init();
}
void Run() {
while (!fStack.IsEmpty()) {
while (!fStack.IsEmpty()) {
Particle& pNext = *fStack.GetNextParticle();
Step(pNext);
}
// do cascade equations, which can put new particles on Stack,
// thus, the double loop
// DoCascadeEquations(); //
}
}
void Step(Particle& particle) {
corsika::setup::Trajectory step = fTracking.GetTrack(particle);
fProcesseList.MinStepLength(particle, step);
/// here the particle is actually moved along the trajectory to new position:
// std::visit(corsika::setup::ParticleUpdate<Particle>{particle}, step);
particle.SetPosition(step.GetPosition(1));
corsika::process::EProcessReturn status =
fProcesseList.DoContinuous(particle, step, fStack);
if (status == corsika::process::EProcessReturn::eParticleAbsorbed) {
fStack.Delete(particle); // TODO: check if this is really needed
} else {
fProcesseList.DoDiscrete(particle, fStack);
}
}
private:
Tracking& fTracking;
ProcessList& fProcesseList;
Stack& fStack;
};
} // namespace corsika::cascade
#endif
Framework/Cascade/SibStack.h deleted 100644 → 0
View file @ dd195476
#ifndef _include_sibstack_h_
#define _include_sibstack_h_
#include <string>
#include <vector>
#include <corsika/cascade/sibyll2.3c.h>
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/stack/Stack.h>
using namespace std;
using namespace corsika::stack;
class SibStackData {
public:
void Init();
void Clear() { s_plist_.np = 0; }
int GetSize() const { return s_plist_.np; }
int GetCapacity() const { return 8000; }
void SetId(const int i, const int v) { s_plist_.llist[i] = v; }
void SetEnergy(const int i, const double v) { s_plist_.p[3][i] = v; }
int GetId(const int i) const { return s_plist_.llist[i]; }
double GetEnergy(const int i) const { return s_plist_.p[3][i]; }
void Copy(const int i1, const int i2) {
s_plist_.llist[i1] = s_plist_.llist[i2];
s_plist_.p[3][i1] = s_plist_.p[3][i2];
}
protected:
void IncrementSize() { s_plist_.np++; }
void DecrementSize() {
if (s_plist_.np > 0) { s_plist_.np--; }
}
};
template <typename StackIteratorInterface>
class ParticleInterface : public ParticleBase<StackIteratorInterface> {
using ParticleBase<StackIteratorInterface>::GetStackData;
using ParticleBase<StackIteratorInterface>::GetIndex;
public:
void SetEnergy(const double v) { GetStackData().SetEnergy(GetIndex(), v); }
double GetEnergy() const { return GetStackData().GetEnergy(GetIndex()); }
void SetPID(const int v) { GetStackData().SetId(GetIndex(), v); }
corsika::process::sibyll::SibyllCode GetPID() const {
return static_cast<corsika::process::sibyll::SibyllCode>(
GetStackData().GetId(GetIndex()));
}
};
typedef Stack<SibStackData, ParticleInterface> SibStack;
#endif
Framework/Cascade/Step.cc deleted 100644 → 0
View file @ dd195476
/**
* (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);
}
Framework/Cascade/rndm_dbl.f deleted 100644 → 0
View file @ dd195476
C***********************************************************************
C
C interface to PHOJET double precision random number generator
C for SIBYLL \FR'14
C
C***********************************************************************
DOUBLE PRECISION FUNCTION S_RNDM(IDUMMY)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DUMMY = dble(IDUMMY)
S_RNDM= PHO_RNDM(DUMMY)
END
C***********************************************************************
C
C initialization routine for double precision random number generator
C calls PHO_RNDIN \FR'14
C
C***********************************************************************
SUBROUTINE RND_INI
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
COMMON /RNDMGAS/ ISET
ISET = 0
CALL PHO_RNDIN(12,34,56,78)
END
DOUBLE PRECISION FUNCTION GASDEV(Idum)
C***********************************************************************
C Gaussian deviation
C***********************************************************************
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
IMPLICIT INTEGER(I-N)
COMMON /RNDMGAS/ ISET
SAVE
DATA ISET/0/
gasdev=idum
IF (ISET.EQ.0) THEN
1 V1=2.D0*S_RNDM(0)-1.D0
V2=2.D0*S_RNDM(1)-1.D0
R=V1**2+V2**2
IF(R.GE.1.D0)GO TO 1
FAC=SQRT(-2.D0*LOG(R)/R)
GSET=V1*FAC
GASDEV=V2*FAC
ISET=1
ELSE
GASDEV=GSET
ISET=0
ENDIF
RETURN
END
C***********************************************************************
DOUBLE PRECISION FUNCTION PHO_RNDM(DUMMY)
C***********************************************************************
C
C random number generator
C
C initialization by call to PHO_RNDIN needed!
C
C the algorithm is taken from
C G.Marsaglia, A.Zaman: 'Toward a unversal random number generator'
C Florida State Univ. preprint FSU-SCRI-87-70
C
C implementation by K. Hahn (Dec. 88), changed to include possibility
C of saving / reading generator registers to / from file (R.E. 10/98)
C
C generator should not depend on the hardware (if a real has
C at least 24 significant bits in internal representation),
C the period is about 2**144,
C
C internal registers:
C U(97),C,CD,CM,I,J - seed values as initialized in PHO_RNDIN
C
C
C***********************************************************************
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
SAVE
COMMON /PORAND/ U(97),C,CD,CM,I,J
100 CONTINUE
RNDMI = DUMMY
RNDMI = U(I)-U(J)
IF ( RNDMI.LT.0.D0 ) RNDMI = RNDMI+1.D0
U(I) = RNDMI
I = I-1
IF ( I.EQ.0 ) I = 97
J = J-1
IF ( J.EQ.0 ) J = 97
C = C-CD
IF ( C.LT.0.D0 ) C = C+CM
RNDMI = RNDMI-C
IF ( RNDMI.LT.0.D0 ) RNDMI = RNDMI+1.D0
IF((ABS(RNDMI).LT.0.D0).OR.(ABS(RNDMI-1.D0).LT.1.D-10)) GOTO 100
PHO_RNDM = RNDMI
END
CDECK ID>, PHO_RNDIN
SUBROUTINE PHO_RNDIN(NA1,NA2,NA3,NB1)
C***********************************************************************
C
C initialization of PHO_RNDM, has to be called before using PHO_RNDM
C
C input:
C NA1,NA2,NA3,NB1 - values for initializing the generator
C NA? must be in 1..178 and not all 1;
C 12,34,56 are the standard values
C NB1 must be in 1..168;
C 78 is the standard value
C
C***********************************************************************
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
SAVE
COMMON /PORAND/ U(97),C,CD,CM,I,J
MA1 = NA1
MA2 = NA2
MA3 = NA3
MB1 = NB1
I = 97
J = 33
DO 20 II2 = 1,97
S = 0.D0
T = 0.5D0
DO 10 II1 = 1,24
MAT = MOD(MOD(MA1*MA2,179)*MA3,179)
MA1 = MA2
MA2 = MA3
MA3 = MAT
MB1 = MOD(53*MB1+1,169)
IF ( MOD(MB1*MAT,64).GE.32 ) S = S+T
T = 0.5D0*T
10 CONTINUE
U(II2) = S
20 CONTINUE
C = 362436.D0/16777216.D0
CD = 7654321.D0/16777216.D0
CM = 16777213.D0/16777216.D0
END
CDECK ID>, PHO_RNDSI
SUBROUTINE PHO_RNDSI(UIN,CIN,CDIN,CMIN,IIN,JIN)
C***********************************************************************
C
C updates internal random number generator registers using
C registers given as arguments
C
C***********************************************************************
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
SAVE
DIMENSION UIN(97)
COMMON /PORAND/ U(97),C,CD,CM,I,J
DO 10 KKK = 1,97
U(KKK) = UIN(KKK)
10 CONTINUE
C = CIN
CD = CDIN
CM = CMIN
I = IIN
J = JIN
END
CDECK ID>, PHO_RNDSO
SUBROUTINE PHO_RNDSO(UOUT,COUT,CDOUT,CMOUT,IOUT,JOUT)
C***********************************************************************
C
C copies internal registers from randon number generator
C to arguments
C
C***********************************************************************
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
SAVE
DIMENSION UOUT(97)
COMMON /PORAND/ U(97),C,CD,CM,I,J
DO 10 KKK = 1,97
UOUT(KKK) = U(KKK)
10 CONTINUE
COUT = C
CDOUT = CD
CMOUT = CM
IOUT = I
JOUT = J
END
CDECK ID>, PHO_RNDTE
SUBROUTINE PHO_RNDTE(IO)
C***********************************************************************
C
C test of random number generator PHO_RNDM
C
C input:
C IO defines output
C 0 output only if an error is detected
C 1 output independend on an error
C
C uses PHO_RNDSI and PHO_RNDSO to bring the random number generator
C to same status as it had before the test run
C
C***********************************************************************
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
SAVE
C input/output channels
INTEGER LI,LO
COMMON /POINOU/ LI,LO
DIMENSION UU(97)
DIMENSION U(6),X(6),D(6)
DATA U / 6533892.D0 , 14220222.D0 , 7275067.D0 ,
& 6172232.D0 , 8354498.D0 , 10633180.D0 /
CALL PHO_RNDSO(UU,CC,CCD,CCM,II,JJ)
CALL PHO_RNDIN(12,34,56,78)
DO 10 II1 = 1,20000
XX = PHO_RNDM(SD)
10 CONTINUE
SD = 0.D0
DO 20 II2 = 1,6
X(II2) = 4096.D0*(4096.D0*PHO_RNDM(XX))
D(II2) = X(II2)-U(II2)
SD = SD+ABS(D(II2))
20 CONTINUE
CALL PHO_RNDSI(UU,CC,CCD,CCM,II,JJ)
IF ((IO.EQ.1).OR.(ABS(SD).GT.0.D-10)) THEN
WRITE(LO,50) (U(I),X(I),D(I),I=1,6)
ENDIF
50 FORMAT(/,' PHO_RNDTE: test of the random number generator:',/,
& ' expected value calculated value difference',/,
& 6(F17.1,F20.1,F15.3,/),
& ' generator has the same status as before calling PHO_RNDTE',/)
END
CDECK ID>, PHO_RNDST
SUBROUTINE PHO_RNDST(MODE,FILENA)
C***********************************************************************
C
C read / write random number generator status from / to file
C
C input: MODE 1 read registers from file
C 2 dump registers to file
C
C FILENA file name
C
C***********************************************************************
IMPLICIT NONE
SAVE
INTEGER MODE
CHARACTER*(*) FILENA
C input/output channels
INTEGER LI,LO
COMMON /POINOU/ LI,LO
DOUBLE PRECISION UU,CC,CCD,CCM
DIMENSION UU(97)
INTEGER I,II,JJ
CHARACTER*80 CH_DUMMY
IF(MODE.EQ.1) THEN
WRITE(LO,'(/,1X,2A,A,/)') 'PHO_RNDST: ',
& 'reading random number registers from file ',FILENA
OPEN(12,FILE=FILENA,ERR=1010,STATUS='OLD')
READ(12,*,ERR=1010) CH_DUMMY
DO I=1,97
READ(12,*,ERR=1010) UU(I)
ENDDO
READ(12,*,ERR=1010) CC
READ(12,*,ERR=1010) CCD
READ(12,*,ERR=1010) CCM
READ(12,*,ERR=1010) II,JJ
CLOSE(12)
CALL PHO_RNDSI(UU,CC,CCD,CCM,II,JJ)
ELSE IF(MODE.EQ.2) THEN
WRITE(LO,'(/,1X,2A,A,/)') 'PHO_RNDST: ',
& 'dumping random number registers to file ',FILENA
OPEN(12,FILE=FILENA,ERR=1010,STATUS='UNKNOWN')
CALL PHO_RNDSO(UU,CC,CCD,CCM,II,JJ)
WRITE(12,'(1X,A)',ERR=1020) 'random number status registers:'
DO I=1,97
WRITE(12,'(1X,1P,E28.20)',ERR=1020) UU(I)
ENDDO
WRITE(12,'(1X,1P,E28.20)',ERR=1020) CC
WRITE(12,'(1X,1P,E28.20)',ERR=1020) CCD
WRITE(12,'(1X,1P,E28.20)',ERR=1020) CCM
WRITE(12,'(1X,2I4)',ERR=1020) II,JJ
CLOSE(12)
ELSE
WRITE(LO,'(/,1X,2A,I6,/)') 'PHO_RNDST: ',
& 'called with invalid mode, nothing done (mode)',MODE
ENDIF
RETURN
1010 CONTINUE
WRITE(LO,'(1X,2A,A,/)') 'PHO_RNDST: ',
& 'cannot open or read file ',FILENA
RETURN
1020 CONTINUE
WRITE(LO,'(1X,A,A,/)') 'PHO_RNDST: ',
& 'cannot open or write file ',FILENA
RETURN
END
C----------------------------------------
C standard generator
C----------------------------------------
REAL FUNCTION S_RNDM_std(IDUMMY)
C...Generator from the LUND montecarlo
C...Purpose: to generate random numbers uniformly distributed between
C...0 and 1, excluding the endpoints.
COMMON/LUDATR/MRLU(6),RRLU(100)
SAVE /LUDATR/
EQUIVALENCE (MRLU1,MRLU(1)),(MRLU2,MRLU(2)),(MRLU3,MRLU(3)),
&(MRLU4,MRLU(4)),(MRLU5,MRLU(5)),(MRLU6,MRLU(6)),
&(RRLU98,RRLU(98)),(RRLU99,RRLU(99)),(RRLU00,RRLU(100))
C... Initialize generation from given seed.
S_RNDM_std = real(idummy)
IF(MRLU2.EQ.0) THEN
IF (MRLU1 .EQ. 0) MRLU1 = 19780503 ! initial seed
IJ=MOD(MRLU1/30082,31329)
KL=MOD(MRLU1,30082)
I=MOD(IJ/177,177)+2
J=MOD(IJ,177)+2
K=MOD(KL/169,178)+1
L=MOD(KL,169)
DO 110 II=1,97
S=0.
T=0.5
DO 100 JJ=1,24
M=MOD(MOD(I*J,179)*K,179)
I=J
J=K
K=M
L=MOD(53*L+1,169)
IF(MOD(L*M,64).GE.32) S=S+T
T=0.5*T
100 CONTINUE
RRLU(II)=S
110 CONTINUE
TWOM24=1.
DO 120 I24=1,24
TWOM24=0.5*TWOM24
120 CONTINUE
RRLU98=362436.*TWOM24
RRLU99=7654321.*TWOM24
RRLU00=16777213.*TWOM24
MRLU2=1
MRLU3=0
MRLU4=97
MRLU5=33
ENDIF
C...Generate next random number.
130 RUNI=RRLU(MRLU4)-RRLU(MRLU5)
IF(RUNI.LT.0.) RUNI=RUNI+1.
RRLU(MRLU4)=RUNI
MRLU4=MRLU4-1
IF(MRLU4.EQ.0) MRLU4=97
MRLU5=MRLU5-1
IF(MRLU5.EQ.0) MRLU5=97
RRLU98=RRLU98-RRLU99
IF(RRLU98.LT.0.) RRLU98=RRLU98+RRLU00
RUNI=RUNI-RRLU98
IF(RUNI.LT.0.) RUNI=RUNI+1.
IF(RUNI.LE.0.OR.RUNI.GE.1.) GOTO 130
C...Update counters. Random number to output.
MRLU3=MRLU3+1
IF(MRLU3.EQ.1000000000) THEN
MRLU2=MRLU2+1
MRLU3=0
ENDIF
S_RNDM_std=RUNI
END
Framework/Cascade/testCascade.cc deleted 100644 → 0
View file @ dd195476
/**
* (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/cascade/Cascade.h>
#include <corsika/process/ProcessSequence.h>
#include <corsika/process/stack_inspector/StackInspector.h>
#include <corsika/process/tracking_line/TrackingLine.h>
#include <corsika/stack/super_stupid/SuperStupidStack.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/geometry/Vector.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::geometry;
#include <iostream>
using namespace std;
static int fCount = 0;
class ProcessSplit : public corsika::process::BaseProcess<ProcessSplit> {
public:
ProcessSplit() {}
template <typename Particle>
void MinStepLength(Particle&, setup::Trajectory&) const {}
template <typename Particle, typename Stack>
EProcessReturn DoContinuous(Particle&, setup::Trajectory&, Stack&) const {
return EProcessReturn::eOk;
}
template <typename Particle, typename Stack>
void DoDiscrete(Particle& p, Stack& s) const {
EnergyType E = p.GetEnergy();
if (E < 85_MeV) {
p.Delete();
fCount++;
} else {
p.SetEnergy(E / 2);
auto pnew = s.NewParticle();
// s.Copy(p, pnew);
pnew.SetEnergy(E / 2);
pnew.SetPosition(p.GetPosition());
pnew.SetMomentum(p.GetMomentum());
}
}
void Init() { fCount = 0; }
int GetCount() const { return fCount; }
private:
};
TEST_CASE("Cascade", "[Cascade]") {
tracking_line::TrackingLine<setup::Stack> tracking;
stack_inspector::StackInspector<setup::Stack> p0(true);
ProcessSplit p1;
const auto sequence = p0 + p1;
setup::Stack stack;
corsika::cascade::Cascade EAS(tracking, sequence, stack);
CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
stack.Clear();
auto particle = stack.NewParticle();
EnergyType E0 = 100_GeV;
particle.SetEnergy(E0);
particle.SetPosition(Point(rootCS, {0_m, 0_m, 10_km}));
particle.SetMomentum(corsika::stack::super_stupid::MomentumVector(
rootCS, {0 * newton * second, 0 * newton * second, -1 * newton * second}));
EAS.Init();
EAS.Run();
SECTION("sectionTwo") {
for (int i = 0; i < 0; ++i) {
stack.Clear();
auto particle = stack.NewParticle();
EnergyType E0 = 100_GeV * pow(10, i);
particle.SetEnergy(E0);
EAS.Init();
EAS.Run();
// cout << "Result: E0=" << E0 / 1_GeV << "GeV, count=" << p1.GetCount() << endl;
}
}
}
Framework/Geometry/BaseTrajectory.h deleted 100644 → 0
View file @ dd195476
/**
* (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 LengthType GetDistance(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
Framework/Geometry/BaseVector.h deleted 100644 → 0
View file @ dd195476
/**
* (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) {}
};
} // namespace corsika::geometry
#endif
Framework/Geometry/CMakeLists.txt deleted 100644 → 0
View file @ dd195476
set (
GEOMETRY_SOURCES
CoordinateSystem.cc
)
set (
GEOMETRY_HEADERS
Vector.h
Point.h
Line.h
Sphere.h
CoordinateSystem.h
RootCoordinateSystem.h
Helix.h
BaseVector.h
QuantityVector.h
Trajectory.h
# BaseTrajectory.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
)
add_test (NAME testGeometry COMMAND testGeometry)
Framework/Geometry/CoordinateSystem.cc deleted 100644 → 0
View file @ dd195476
/**
* (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>
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::string("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;
}
Framework/Geometry/CoordinateSystem.h deleted 100644 → 0
View file @ dd195476
/**
* (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>
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
/// creat 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::string("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::string("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
Framework/Geometry/Helix.h deleted 100644 → 0
View file @ dd195476
/**
* (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;
}
auto GetRadius() const { return radius; }
LengthType GetDistanceBetween(corsika::units::si::TimeType t1,
corsika::units::si::TimeType t2) const {
return (vPar + vPerp).norm() * (t2 - t1);
}
};
} // namespace corsika::geometry
#endif
Framework/Geometry/Line.h deleted 100644 → 0
View file @ dd195476
/**
* (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; }
LengthType GetDistanceBetween(corsika::units::si::TimeType t1,
corsika::units::si::TimeType t2) const {
// assert(t2 >= t1);
return v0.norm() * (t2 - t1);
}
};
} // namespace corsika::geometry
#endif
Framework/Geometry/Point.h deleted 100644 → 0
View file @ dd195476
/**
* (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
Framework/Geometry/QuantityVector.h deleted 100644 → 0
View file @ dd195476
/**
* (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_QUANTITYVECTOR_H_
#define _include_QUANTITYVECTOR_H_
#include <corsika/units/PhysicalUnits.h>
#include <Eigen/Dense>
#include <iostream>
#include <utility>
namespace corsika::geometry {
/*!
* A QuantityVector is a three-component container based on Eigen::Vector3d
* with a phys::units::si::dimension. Arithmethic operators are defined that
* propagate the dimensions by dimensional analysis.
*/
template <typename dim>
class QuantityVector {
protected:
// todo: check if we need to move "quantity" into namespace corsika::units
using Quantity = phys::units::quantity<dim, double>; //< the phys::units::quantity
// corresponding to the dimension
public:
Eigen::Vector3d eVector; //!< the actual container where the raw numbers are stored
typedef dim dimension; //!< should be a phys::units::dimension
QuantityVector(Quantity a, Quantity b, Quantity c)
: eVector{a.magnitude(), b.magnitude(), c.magnitude()} {}
QuantityVector(Eigen::Vector3d pBareVector)
: eVector(pBareVector) {}
auto operator[](size_t index) const {
return Quantity(phys::units::detail::magnitude_tag, eVector[index]);
}
auto norm() const {
return Quantity(phys::units::detail::magnitude_tag, eVector.norm());
}
auto squaredNorm() const {
using QuantitySquared =
decltype(std::declval<Quantity>() * std::declval<Quantity>());
return QuantitySquared(phys::units::detail::magnitude_tag, eVector.squaredNorm());
}
auto operator+(QuantityVector<dim> const& pQVec) const {
return QuantityVector<dim>(eVector + pQVec.eVector);
}
auto operator-(QuantityVector<dim> const& pQVec) const {
return QuantityVector<dim>(eVector - pQVec.eVector);
}
template <typename ScalarDim>
auto operator*(phys::units::quantity<ScalarDim, double> const p) const {
using ResQuantity = phys::units::detail::Product<ScalarDim, dim, double, double>;
if constexpr (std::is_same<ResQuantity, double>::value) // result dimensionless, not
// a "Quantity" anymore
{
return QuantityVector<phys::units::dimensionless_d>(eVector * p.magnitude());
} else {
return QuantityVector<typename ResQuantity::dimension_type>(eVector *
p.magnitude());
}
}
template <typename ScalarDim>
auto operator/(phys::units::quantity<ScalarDim, double> const p) const {
return (*this) * (1 / p);
}
auto operator*(double const p) const { return QuantityVector<dim>(eVector * p); }
auto operator/(double const p) const { return QuantityVector<dim>(eVector / p); }
auto& operator/=(double const p) {
eVector /= p;
return *this;
}
auto& operator*=(double const p) {
eVector *= p;
return *this;
}
auto& operator+=(QuantityVector<dim> const& pQVec) {
eVector += pQVec.eVector;
return *this;
}
auto& operator-=(QuantityVector<dim> const& pQVec) {
eVector -= pQVec.eVector;
return *this;
}
auto& operator-() const { return QuantityVector<dim>(-eVector); }
auto normalized() const { return (*this) * (1 / norm()); }
auto operator==(QuantityVector<dim> const& p) const { return eVector == p.eVector; }
};
} // namespace corsika::geometry
template <typename dim>
auto& operator<<(std::ostream& os, corsika::geometry::QuantityVector<dim> qv) {
using Quantity = phys::units::quantity<dim, double>;
os << '(' << qv.eVector(0) << ' ' << qv.eVector(1) << ' ' << qv.eVector(2) << ") "
<< phys::units::to_unit_symbol<dim, double>(
Quantity(phys::units::detail::magnitude_tag, 1));
return os;
}
#endif
Framework/Geometry/RootCoordinateSystem.h deleted 100644 → 0
View file @ dd195476
#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& GetRootCS() { return fRootCS; }
const corsika::geometry::CoordinateSystem& GetRootCS() const { return fRootCS; }
private:
corsika::geometry::CoordinateSystem fRootCS; // THIS IS IT
};
} // namespace corsika::geometry
#endif
Framework/Geometry/Sphere.h deleted 100644 → 0
View file @ dd195476
/**
* (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/units/PhysicalUnits.h>
namespace corsika::geometry {
class Sphere {
Point center;
LengthType const radius;
public:
Sphere(Point const& pCenter, LengthType const pRadius)
: center(pCenter)
, radius(pRadius) {}
//! returns true if the Point \a p is within the sphere
auto Contains(Point const& p) const {
return radius * radius > (center - p).squaredNorm();
}
};
} // namespace corsika::geometry
#endif
Framework/Geometry/Trajectory.h deleted 100644 → 0
View file @ dd195476
/**
* (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_TRAJECTORY_H
#define _include_TRAJECTORY_H
#include <corsika/units/PhysicalUnits.h>
using corsika::units::si::LengthType;
using corsika::units::si::TimeType;
namespace corsika::geometry {
template <typename T>
class Trajectory : public T {
corsika::units::si::TimeType fTimeLength;
public:
using T::GetDistanceBetween;
using T::GetPosition;
Trajectory(T const& theT, corsika::units::si::TimeType timeLength)
: T(theT)
, fTimeLength(timeLength) {}
/*Point GetPosition(corsika::units::si::TimeType t) const {
return fTraj.GetPosition(t + fTStart);
}*/
Point GetPosition(const double u) const { return T::GetPosition(fTimeLength * u); }
TimeType GetDuration() const { return fTimeLength; }
LengthType GetDistance(const corsika::units::si::TimeType t) const {
assert(t > fTimeLength);
assert(t >= 0 * corsika::units::si::second);
return T::DistanceBetween(0, t);
}
};
} // namespace corsika::geometry
#endif
Framework/Geometry/Vector.h deleted 100644 → 0
View file @ dd195476
/**
* (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_VECTOR_H_
#define _include_VECTOR_H_
#include <corsika/geometry/BaseVector.h>
#include <corsika/geometry/QuantityVector.h>
#include <corsika/units/PhysicalUnits.h>
/*!
* A Vector represents a 3-vector in Euclidean space. It is defined by components
* given in a specific CoordinateSystem. It has a physical dimension ("unit")
* as part of its type, so you cannot mix up e.g. electric with magnetic fields
* (but you could calculate their cross-product to get an energy flux vector).
*
* When transforming coordinate systems, a Vector is subject to the rotational
* part only and invariant under translations.
*/
namespace corsika::geometry {
template <typename dim>
class Vector : public BaseVector<dim> {
using Quantity = phys::units::quantity<dim, double>;
public:
Vector(CoordinateSystem const& pCS, QuantityVector<dim> pQVector)
: BaseVector<dim>(pCS, pQVector) {}
Vector(CoordinateSystem const& cs, Quantity x, Quantity y, Quantity z)
: BaseVector<dim>(cs, QuantityVector<dim>(x, y, z)) {}
/*!
* returns a QuantityVector with the components given in the "home"
* CoordinateSystem of the Vector
*/
auto GetComponents() const { return BaseVector<dim>::qVector; }
/*!
* returns a QuantityVector with the components given in an arbitrary
* CoordinateSystem
*/
auto GetComponents(CoordinateSystem const& pCS) const {
if (&pCS == BaseVector<dim>::cs) {
return BaseVector<dim>::qVector;
} else {
return QuantityVector<dim>(
CoordinateSystem::GetTransformation(*BaseVector<dim>::cs, pCS).linear() *
BaseVector<dim>::qVector.eVector);
}
}
/*!
* transforms the Vector into another CoordinateSystem by changing
* its components internally
*/
void rebase(CoordinateSystem const& pCS) {
BaseVector<dim>::qVector = GetComponents(pCS);
BaseVector<dim>::cs = &pCS;
}
/*!
* returns the norm/length of the Vector. Before using this method,
* think about whether squaredNorm() might be cheaper for your computation.
*/
auto norm() const { return BaseVector<dim>::qVector.norm(); }
/*!
* returns the squared norm of the Vector. Before using this method,
* think about whether norm() might be cheaper for your computation.
*/
auto squaredNorm() const { return BaseVector<dim>::qVector.squaredNorm(); }
/*!
* returns a Vector \f$ \vec{v}_{\parallel} \f$ which is the parallel projection
* of this vector \f$ \vec{v}_1 \f$ along another Vector \f$ \vec{v}_2 \f$ given by
* \f[
* \vec{v}_{\parallel} = \frac{\vec{v}_1 \cdot \vec{v}_2}{\vec{v}_2^2} \vec{v}_2
* \f]
*/
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& a = ourCompVec.eVector;
auto const& b = otherCompVec.eVector;
return Vector<dim>(pCS, QuantityVector<dim>(b * ((a.dot(b)) / b.squaredNorm())));
}
template <typename dim2>
auto parallelProjectionOnto(Vector<dim2> const& pVec) const {
return parallelProjectionOnto<dim2>(pVec, *BaseVector<dim>::cs);
}
auto operator+(Vector<dim> const& pVec) const {
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);
return Vector<dim>(*BaseVector<dim>::cs, components);
}
auto& operator*=(double const p) {
BaseVector<dim>::qVector *= p;
return *this;
}
template <typename ScalarDim>
auto operator*(phys::units::quantity<ScalarDim, double> const p) const {
using ProdQuantity = phys::units::detail::Product<dim, ScalarDim, double, double>;
if constexpr (std::is_same<ProdQuantity, double>::value) // result dimensionless,
// not a "Quantity" anymore
{
return Vector<phys::units::dimensionless_d>(*BaseVector<dim>::cs,
BaseVector<dim>::qVector * p);
} else {
return Vector<typename ProdQuantity::dimension_type>(
*BaseVector<dim>::cs, BaseVector<dim>::qVector * p);
}
}
template <typename ScalarDim>
auto operator/(phys::units::quantity<ScalarDim, double> const p) const {
return (*this) * (1 / p);
}
auto operator*(double const p) const {
return Vector<dim>(*BaseVector<dim>::cs, BaseVector<dim>::qVector * p);
}
auto operator/(double const p) const {
return Vector<dim>(*BaseVector<dim>::cs, BaseVector<dim>::qVector / p);
}
auto& operator+=(Vector<dim> const& pVec) {
BaseVector<dim>::qVector += pVec.GetComponents(*BaseVector<dim>::cs);
return *this;
}
auto& operator-=(Vector<dim> const& pVec) {
BaseVector<dim>::qVector -= pVec.GetComponents(*BaseVector<dim>::cs);
return *this;
}
auto& operator-() const {
return Vector<dim>(*BaseVector<dim>::cs, -BaseVector<dim>::qVector);
}
auto normalized() const { return (*this) * (1 / norm()); }
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 bareResult = c1.cross(c2);
using ProdQuantity = phys::units::detail::Product<dim, dim2, double, double>;
if constexpr (std::is_same<ProdQuantity, double>::value) // result dimensionless,
// not a "Quantity" anymore
{
return Vector<phys::units::dimensionless_d>(*BaseVector<dim>::cs, bareResult);
} else {
return Vector<typename ProdQuantity::dimension_type>(*BaseVector<dim>::cs,
bareResult);
}
}
};
} // namespace corsika::geometry
#endif
Framework/Geometry/testGeometry.cc deleted 100644 → 0
View file @ dd195476
/**
* (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.
*/
#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
// cpp file
#include <catch2/catch.hpp>
#include <corsika/geometry/CoordinateSystem.h>
#include <corsika/geometry/Helix.h>
#include <corsika/geometry/Line.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/geometry/Sphere.h>
#include <corsika/geometry/Trajectory.h>
#include <corsika/units/PhysicalUnits.h>
#include <cmath>
using namespace corsika::geometry;
using namespace corsika::units::si;
double constexpr absMargin = 1.0e-8;
TEST_CASE("transformations between CoordinateSystems") {
CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
REQUIRE(CoordinateSystem::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(CoordinateSystem::GetTransformation(cs3, cs2).isApprox(
rootCS.translate({3_m, -5_m, 0_m}).GetTransform()));
REQUIRE(CoordinateSystem::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));
}
}
TEST_CASE("Sphere") {
CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
Point center(rootCS, {0_m, 3_m, 4_m});
Sphere sphere(center, 5_m);
SECTION("isInside") {
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().GetRootCS();
Point r0(rootCS, {0_m, 0_m, 0_m});
SECTION("Line") {
Vector<SpeedType::dimension_type> v0(rootCS,
{1_m / second, 0_m / second, 0_m / second});
Line const line(r0, v0);
CHECK(
(line.GetPosition(2_s).GetCoordinates() - QuantityVector<length_d>(2_m, 0_m, 0_m))
.norm()
.magnitude() == Approx(0).margin(absMargin));
Trajectory<Line> base(line, 1_s);
CHECK(line.GetPosition(2_s).GetCoordinates() ==
base.GetPosition(2_s).GetCoordinates());
CHECK(base.GetDistanceBetween(1_s, 2_s) / 1_m == Approx(1));
}
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 omegaC = 2 * M_PI / 1_s;
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));
Trajectory<Helix> const base(helix, 1_s);
CHECK(helix.GetPosition(1234_s).GetCoordinates() ==
base.GetPosition(1234_s).GetCoordinates());
CHECK(base.GetDistanceBetween(1_s, 2_s) / 1_m == Approx(5));
}
}

data privacy - impressum