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 1554 deletions
Environment/HomogeneousMedium.h deleted 100644 → 0
View file @ a4738ea3
/**
* (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_HomogeneousMedium_h_
#define _include_HomogeneousMedium_h_
#include <corsika/environment/NuclearComposition.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/geometry/Trajectory.h>
#include <corsika/geometry/Line.h>
#include <corsika/geometry/Point.h>
/**
* a homogeneous medium
*/
namespace corsika::environment {
template <class T>
class HomogeneousMedium : T {
corsika::units::si::MassDensityType const fDensity;
NuclearComposition const fNuclComp;
public:
HomogeneousMedium(corsika::units::si::MassDensityType pDensity,
NuclearComposition pNuclComp)
: fDensity(pDensity)
, fNuclComp(pNuclComp){};
corsika::units::si::MassDensityType GetMassDensity([
[maybe_unused]] corsika::geometry::Point const& p) const override {
return fDensity;
}
NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
corsika::units::si::GrammageType IntegratedGrammage(
corsika::geometry::Trajectory<corsika::geometry::Line> const& pTraj,
corsika::units::si::TimeType pTo) const override {
using namespace corsika::units::si;
return pTraj.ArcLength(0_s, pTo) * fDensity;
}
corsika::units::si::TimeType FromGrammage(
corsika::geometry::Trajectory<corsika::geometry::Line> const& pTraj,
corsika::units::si::GrammageType pGrammage) const override {
return pTraj.TimeFromArclength(pGrammage / fDensity);
}
};
} // namespace corsika::environment
#endif
Environment/IMediumModel.h deleted 100644 → 0
View file @ a4738ea3
#ifndef _include_IMediumModel_h
#define _include_IMediumModel_h
#include <corsika/environment/NuclearComposition.h>
#include <corsika/geometry/Trajectory.h>
#include <corsika/geometry/Point.h>
#include <corsika/units/PhysicalUnits.h>
namespace corsika::environment {
class IMediumModel {
public:
virtual ~IMediumModel() = default;
virtual corsika::units::si::MassDensityType GetMassDensity(
corsika::geometry::Point const&) const = 0;
// todo: think about the mixin inheritance of the trajectory vs the BaseTrajectory approach
// for now, only lines are supported
virtual corsika::units::si::GrammageType IntegratedGrammage(
corsika::geometry::Trajectory<corsika::geometry::Line> const&, corsika::units::si::TimeType) const = 0;
virtual corsika::units::si::TimeType FromGrammage(
corsika::geometry::Trajectory<corsika::geometry::Line> const&,
corsika::units::si::GrammageType) const = 0;
virtual NuclearComposition const& GetNuclearComposition() const = 0;
};
} // namespace corsika::environment
#endif
Environment/NuclearComposition.h deleted 100644 → 0
View file @ a4738ea3
#ifndef _include_NuclearComposition_h
#define _include_NuclearComposition_h
#include <corsika/particles/ParticleProperties.h>
#include <numeric>
#include <vector>
namespace corsika::environment {
class NuclearComposition {
std::vector<float> const fNumberFractions; //!< relative fractions of number density
std::vector<corsika::particles::Code> const
fComponents; //!< particle codes of consitutents
public:
NuclearComposition(std::vector<corsika::particles::Code> pComponents,
std::vector<float> pFractions)
: fNumberFractions(pFractions)
, fComponents(pComponents) {
auto const sumFractions =
std::accumulate(pFractions.cbegin(), pFractions.cend(), 0.f);
if (!(0.999f < sumFractions && sumFractions < 1.001f)) {
throw std::string("element fractions do not add up to 1");
}
}
auto size() const { return fNumberFractions.size(); }
auto const& GetFractions() const { return fNumberFractions; }
auto const& GetComponents() const { return fComponents; }
};
} // namespace corsika::environment
#endif
Environment/VolumeTreeNode.h deleted 100644 → 0
View file @ a4738ea3
/**
* (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_VolumeTreeNode_H
#define _include_VolumeTreeNode_H
#include <corsika/geometry/Volume.h>
#include <memory>
#include <vector>
namespace corsika::environment {
class Empty {}; //<! intended for usage as default template argument
template <typename IModelProperties = Empty>
class VolumeTreeNode {
public:
using VTNUPtr = std::unique_ptr<VolumeTreeNode<IModelProperties>>;
using IMPSharedPtr = std::shared_ptr<IModelProperties>;
using VolUPtr = std::unique_ptr<corsika::geometry::Volume>;
VolumeTreeNode(VolUPtr pVolume = nullptr)
: fGeoVolume(std::move(pVolume)) {}
bool Contains(corsika::geometry::Point const& p) const {
return fGeoVolume->Contains(p);
}
VolumeTreeNode<IModelProperties> const* Excludes(
corsika::geometry::Point const& p) const {
auto exclContainsIter =
std::find_if(fExcludedNodes.cbegin(), fExcludedNodes.cend(),
[&](auto const& s) { return bool(s->Contains(p)); });
return exclContainsIter != fExcludedNodes.cend() ? *exclContainsIter : nullptr;
}
/** returns a pointer to the sub-VolumeTreeNode which is "responsible" for the given
* \class Point \arg p, or nullptr iff \arg p is not contained in this volume.
*/
VolumeTreeNode<IModelProperties> const* GetContainingNode(
corsika::geometry::Point const& p) const {
if (!Contains(p)) { return nullptr; }
if (auto const childContainsIter =
std::find_if(fChildNodes.cbegin(), fChildNodes.cend(),
[&](auto const& s) { return bool(s->Contains(p)); });
childContainsIter == fChildNodes.cend()) // not contained in any of the children
{
if (auto const exclContainsIter = Excludes(p)) // contained in any excluded nodes
{
return exclContainsIter->GetContainingNode(p);
} else {
return this;
}
} else {
return (*childContainsIter)->GetContainingNode(p);
}
}
void addChild(VTNUPtr pChild) {
pChild->fParentNode = this;
fChildNodes.push_back(std::move(pChild));
// It is a bad idea to return an iterator to the inserted element
// because it might get invalidated when the vector needs to grow
// later and the caller won't notice.
}
void excludeOverlapWith(VTNUPtr const& pNode) {
fExcludedNodes.push_back(pNode.get());
}
auto GetParent() const { return fParentNode; };
auto const& GetVolume() const { return *fGeoVolume; }
auto const& GetModelProperties() const { return *fModelProperties; }
template <typename ModelProperties, typename... Args>
auto SetModelProperties(Args&&... args) {
static_assert(std::is_base_of_v<IModelProperties, ModelProperties>,
"unusable type provided");
fModelProperties = std::make_unique<ModelProperties>(std::forward<Args>(args)...);
}
void SetModelProperties(IMPSharedPtr ptr) { fModelProperties = ptr; }
template <class MediumType, typename... Args>
static auto CreateMedium(Args&&... args) {
static_assert(std::is_base_of_v<IMediumModel, MediumType>,
"unusable type provided, needs to be derived from \"IMediumModel\"");
return std::make_shared<MediumType>(std::forward<Args>(args)...);
}
// factory methods for creation of nodes
template <class VolumeType, typename... Args>
static auto CreateNode(Args&&... args) {
static_assert(std::is_base_of_v<corsika::geometry::Volume, VolumeType>,
"unusable type provided, needs to be derived from "
"\"corsika::geometry::Volume\"");
return std::make_unique<VolumeTreeNode<IModelProperties>>(
std::make_unique<VolumeType>(std::forward<Args>(args)...));
}
private:
std::vector<VTNUPtr> fChildNodes;
std::vector<VolumeTreeNode<IModelProperties> const*> fExcludedNodes;
VolumeTreeNode<IModelProperties> const* fParentNode = nullptr;
VolUPtr fGeoVolume;
std::shared_ptr<IModelProperties> fModelProperties;
};
} // namespace corsika::environment
#endif
Environment/testEnvironment.cc deleted 100644 → 0
View file @ a4738ea3
#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
// cpp file
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/IMediumModel.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/environment/VolumeTreeNode.h>
#include <corsika/particles/ParticleProperties.h>
#include <catch2/catch.hpp>
using namespace corsika::geometry;
using namespace corsika::environment;
using namespace corsika::units::si;
TEST_CASE("HomogeneousMedium") {
NuclearComposition const protonComposition(
std::vector<corsika::particles::Code>{corsika::particles::Code::Proton},
std::vector<float>{1.f});
HomogeneousMedium<IMediumModel> const medium(19.2_g / cube(1_cm), protonComposition);
}
Framework/CMakeLists.txt deleted 100644 → 0
View file @ a4738ea3
add_subdirectory (Utilities)
add_subdirectory (Units)
add_subdirectory (Geometry)
add_subdirectory (Particles)
add_subdirectory (Logging)
add_subdirectory (StackInterface)
add_subdirectory (ProcessSequence)
add_subdirectory (Cascade)
add_subdirectory (Random)
Framework/Cascade/CMakeLists.txt deleted 100644 → 0
View file @ a4738ea3
# 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.cc deleted 100644 → 0
View file @ a4738ea3
/**
* (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>
using namespace corsika::cascade;
template <typename ProcessList, typename Stack> //, typename Trajectory>
void Cascade::Init() {
fStack.Init();
fProcesseList.Init();
}
template <typename ProcessList, typename Stack> //, 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 ProcessList, typename Stack> //, typename Trajectory>
void Cascade::Step(Particle& particle) {
double nextStep = fProcesseList.MinStepLength(particle);
corsika::geometry::LineTrajectory trajectory =
fProcesseList.Transport(particle, nextStep);
sequence.DoContinuous(particle, trajectory);
// whats going on here? Everywhere else DoDiscrete is passed a Stack reference as well
sequence.DoDiscrete(particle);
}
Framework/Cascade/Cascade.h deleted 100644 → 0
View file @ a4738ea3
/**
* (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) {
// static_assert(std::is_member_function_pointer<decltype(&ProcessList::DoDiscrete)>::value,
//"ProcessList has not function DoDiscrete.");
// static_assert(std::is_member_function_pointer<decltype(&ProcessList::DoContinuous)>::value,
// "ProcessList has not function DoContinuous.");
}
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);
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 @ a4738ea3
#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 @ a4738ea3
/**
* (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 @ a4738ea3
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 @ a4738ea3
/**
* (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/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 @ a4738ea3
/**
* (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
Framework/Geometry/BaseVector.h deleted 100644 → 0
View file @ a4738ea3
/**
* (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 @ a4738ea3
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
# 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 @ a4738ea3
/**
* (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 @ a4738ea3
/**
* (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 @ a4738ea3
/**
* (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; }
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 t) const {
return t / (vPar + vPerp).norm();
}
};
} // namespace corsika::geometry
#endif
Framework/Geometry/Line.h deleted 100644 → 0
View file @ a4738ea3
/**
* (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 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();
}
};
} // namespace corsika::geometry
#endif

data privacy - impressum