IAP GITLAB

Skip to content
Snippets Groups Projects

Compare revisions

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

Source

Select target project
No results found

Target

Select target project
  • AirShowerPhysics/corsika
  • rulrich/corsika
  • AAAlvesJr/corsika
  • Andre/corsika
  • arrabito/corsika
  • Nikos/corsika
  • olheiser73/corsika
  • AirShowerPhysics/papers/corsika
  • pranav/corsika
9 results
Show changes
Showing
with 596 additions and 650 deletions
This diff is collapsed.
...@@ -3,9 +3,8 @@ ...@@ -3,9 +3,8 @@
# #
# See file AUTHORS for a list of contributors. # See file AUTHORS for a list of contributors.
# #
# This software is distributed under the terms of the GNU General Public # This software is distributed under the terms of the 3-clause BSD license.
# Licence version 3 (GPL Version 3). See file LICENSE for a full version of # See file LICENSE for a full version of the license.
# the license.
# #
set (CORSIKA8_VERSION @c8_version@) set (CORSIKA8_VERSION @c8_version@)
...@@ -31,10 +30,16 @@ set (COMPILE_OPTIONS @COMPILE_OPTIONS@) ...@@ -31,10 +30,16 @@ set (COMPILE_OPTIONS @COMPILE_OPTIONS@)
set (CMAKE_VERBOSE_MAKEFILE @CMAKE_VERBOSE_MAKEFILE@) set (CMAKE_VERBOSE_MAKEFILE @CMAKE_VERBOSE_MAKEFILE@)
#+++++++++++++++++++++++++++++ #+++++++++++++++++++++++++++++
# Setup external dependencies via conan # external dependencies
# same list as top-level CML.txt, except for Catch2 (not needed here)
# #
include (${CMAKE_CURRENT_LIST_DIR}/conanbuildinfo.cmake) find_package(Boost COMPONENTS filesystem REQUIRED)
conan_basic_setup (TARGETS) find_package(CLI11 REQUIRED)
find_package(Eigen3 REQUIRED)
find_package(spdlog REQUIRED)
find_package(yaml-cpp REQUIRED)
find_package(Arrow REQUIRED)
find_package(PROPOSAL REQUIRED)
#+++++++++++++++++++++++++++++ #+++++++++++++++++++++++++++++
# Import Pythia8 # Import Pythia8
...@@ -51,6 +56,35 @@ set_target_properties ( ...@@ -51,6 +56,35 @@ set_target_properties (
set (Pythia8_FOUND @Pythia8_FOUND@) set (Pythia8_FOUND @Pythia8_FOUND@)
message (STATUS "Pythia8 at: @Pythia8_PREFIX@") message (STATUS "Pythia8 at: @Pythia8_PREFIX@")
#+++++++++++++++++++++++++++++
# Import TAUOLA
# since we always import TAUOLA (ExternalProject_Add) we have to
# import here, too.
#
add_library (C8::ext::tauola::CxxInterface STATIC IMPORTED GLOBAL)
add_library (C8::ext::tauola::Fortran STATIC IMPORTED GLOBAL)
add_library(C8::ext::tauola INTERFACE IMPORTED)
set_property(TARGET C8::ext::tauola
PROPERTY
INTERFACE_LINK_LIBRARIES
C8::ext::tauola::CxxInterface
C8::ext::tauola::Fortran)
set_target_properties (
C8::ext::tauola::CxxInterface PROPERTIES
IMPORTED_LOCATION @TAUOLA_LIBDIR@/libTauolaCxxInterface.a
IMPORTED_LINK_INTERFACE_LIBRARIES dl
INTERFACE_INCLUDE_DIRECTORIES @TAUOLA_INCDIR@
)
set_target_properties (
C8::ext::tauola::Fortran PROPERTIES
IMPORTED_LOCATION @TAUOLA_LIBDIR@/libTauolaFortran.a
IMPORTED_LINK_INTERFACE_LIBRARIES dl
INTERFACE_INCLUDE_DIRECTORIES @TAUOLA_INCDIR@
)
set (TAUOLA_FOUND @TAUOLA_FOUND@)
message (STATUS "TAUOLA at: @TAUOLA_PREFIX@")
#++++++++++++++++++++++++++++++ #++++++++++++++++++++++++++++++
# import CORSIKA8 # import CORSIKA8
......
...@@ -3,9 +3,8 @@ ...@@ -3,9 +3,8 @@
# #
# See file AUTHORS for a list of contributors. # See file AUTHORS for a list of contributors.
# #
# This software is distributed under the terms of the GNU General Public # This software is distributed under the terms of the 3-clause BSD license.
# Licence version 3 (GPL Version 3). See file LICENSE for a full version of # See file LICENSE for a full version of the license.
# the license.
# #
...@@ -24,7 +23,7 @@ set (CMAKE_CXX_FLAGS_COVERAGE "-g --coverage") ...@@ -24,7 +23,7 @@ set (CMAKE_CXX_FLAGS_COVERAGE "-g --coverage")
set (CMAKE_EXE_LINKER_FLAGS_COVERAGE "--coverage") set (CMAKE_EXE_LINKER_FLAGS_COVERAGE "--coverage")
set (CMAKE_SHARED_LINKER_FLAGS_COVERAGE "--coverage") set (CMAKE_SHARED_LINKER_FLAGS_COVERAGE "--coverage")
# set a flag to inform code that we are in debug mode # set a flag to inform code that we are in debug mode
set (CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -DDEBUG") set (CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -D_C8_DEBUG_")
#+++++++++++++++++++++++++++++ #+++++++++++++++++++++++++++++
......
...@@ -3,9 +3,8 @@ ...@@ -3,9 +3,8 @@
# #
# See file AUTHORS for a list of contributors. # See file AUTHORS for a list of contributors.
# #
# This software is distributed under the terms of the GNU General Public # This software is distributed under the terms of the 3-clause BSD license.
# Licence version 3 (GPL Version 3). See file LICENSE for a full version of # See file LICENSE for a full version of the license.
# the license.
# #
################################################# #################################################
...@@ -55,9 +54,10 @@ function (CORSIKA_ADD_TEST) ...@@ -55,9 +54,10 @@ function (CORSIKA_ADD_TEST)
else () else ()
set(sanitize ${C8_ADD_TEST_SANITIZE}) set(sanitize ${C8_ADD_TEST_SANITIZE})
endif () endif ()
find_package(Catch2 REQUIRED)
add_executable (${name} ${sources}) add_executable (${name} ${sources})
target_link_libraries (${name} CORSIKA8 CONAN_PKG::catch2 CorsikaTestingCommon) target_link_libraries (${name} CORSIKA8 Catch2::Catch2WithMain CorsikaTestingCommon)
target_compile_options (${name} PRIVATE -g) # do not skip asserts target_compile_options (${name} PRIVATE -g) # do not skip asserts
target_include_directories (${name} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}) target_include_directories (${name} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR})
file (MAKE_DIRECTORY ${PROJECT_BINARY_DIR}/test_outputs/) file (MAKE_DIRECTORY ${PROJECT_BINARY_DIR}/test_outputs/)
......
#!/bin/bash
function show_usage(){
printf "\n\nUsage:"
printf "$0 options [parameters]\n"
printf "\n"
printf "Options:\n"
printf "\n -s or --source-directory:\n Corsika 8 download directory, which contains the 'conanfile.py' recipe. Default is the current directory."
printf "\n -d or --debug:\n Specify 'Debug' as build type for the installed dependences. This should be matched when building CORSIKA 8."
printf "\n -r or --release:\n Specify 'Release' as build type for the installed dependences. This should be matched when building CORSIKA 8."
printf "\n -rd or --release-with-debug:\n Specify 'RelWithDebInfo' as build type for the installed dependences. This should be matched when building CORSIKA 8."
printf "\n\nExample: ./conan-install.sh --source-directory /some_path/corsika --debug"
printf "\n -h or --help:\n Prints this message.\n"
exit 0
}
echo "|---------------------------------------------------|"
echo "|-----------------[ CORSIKA 8 ]---------------------|"
echo "|-----[ CONAN2 DEPENDENCIES INSTALL SCRIPT ]------- |"
echo "|---------------------------------------------------|"
echo "|-------------------- BEGIN ------------------------| "
echo " "
SOURCE=${BASH_SOURCE[0]}
while [ -L "$SOURCE" ]; do # resolve $SOURCE until the file is no longer a symlink
DIR=$( cd -P "$( dirname "$SOURCE" )" >/dev/null 2>&1 && pwd )
SOURCE=$(readlink "$SOURCE")
[[ $SOURCE != /* ]] && SOURCE=$DIR/$SOURCE # if $SOURCE was a relative symlink, we need to resolve it relative to the path where the symlink file was located
done
SCRIPT_DIR=$( cd -P "$( dirname "$SOURCE" )" >/dev/null 2>&1 && pwd )
BUILD_TYPE="RelWithDebInfo"
CORSIKA_DIR=${CURRENT_DIR}
CORSIKA_DIR_INFORMED=${CURRENT_DIR}
CONAN2_OUTPUT_FOLDER_NAME="conan_cmake"
printf "[ conan-install | info > This script is located at the directory: ${SCRIPT_DIR}\n"
if [[ "$1" == "--help" ]] || [[ "$1" == "-h" ]];then
show_usage
fi
while [ ! -z "$1" ]; do
case "$1" in
--source-directory|-s)
shift
CORSIKA_DIR=$(readlink -e ${1})
CORSIKA_DIR_INFORMED=${1}
;;
--debug|-d)
#shift
BUILD_TYPE="Debug"
;;
--release|-r)
#shift
BUILD_TYPE="Release"
;;
--release-with-debug|-rd)
#shift
BUILD_TYPE="RelWithDebInfo"
;;
--help|-h)
#shift
show_usage
;;
*)
show_usage
;;
esac
if [ $# -gt 0 ]; then
shift
fi
done
if [[ ${#CORSIKA_DIR} -eq 0 ]]; then
printf "[ conan-install | warning > Output folder for cmake scripts generated by conan2 not found.\n"
printf "[ conan-install | warning > Directory '${CORSIKA_DIR_INFORMED}' does not exist.\n"
exit 1
fi
if [ -d "${CORSIKA_DIR}/${CONAN2_OUTPUT_FOLDER_NAME}" ]; then
printf "[ conan-install | info > Output folder for cmake scripts generated by conan2 is: ${CORSIKA_DIR}/${CONAN2_OUTPUT_FOLDER_NAME}\n"
else
printf "[ conan-install | warning > Output folder for cmake scripts generated by conan2 not found.\n"
printf "[ conan-install | warning > Creating directory ${CORSIKA_DIR}/${CONAN2_OUTPUT_FOLDER_NAME}.\n"
mkdir -p "${CORSIKA_DIR}/${CONAN2_OUTPUT_FOLDER_NAME}"
fi
if [ -f "${CORSIKA_DIR}/conanfile.py" ]; then
printf "[ conan-install | info > conan2 recipe: ${CORSIKA_DIR}/conanfile.py\n"
else
printf "[ conan-install | warning > conan2 recipe not found.\n"
printf "[ conan-install | warning > File '${CORSIKA_DIR}/conanfile.py' does not exist.\n"
exit 1
fi
# Conan2 variabes
CONAN2_HOME=$(readlink -e `conan config home`)
CONAN2_PROFILE_NAME="corsika8"
printf "[ conan-install | info > conan2 home: ${CONAN2_HOME}\n"
# Conan2 commands
CONAN2_DEFAULT_PROFILE_COMMAND="conan profile detect --force"
CONAN2_PROFILE_COMMAND="conan profile detect --name ${CONAN2_PROFILE_NAME} --force"
CONAN2_INSTALL_COMMAND="conan install ${CORSIKA_DIR} --output-folder=${CORSIKA_DIR}/${CONAN2_OUTPUT_FOLDER_NAME} --build=missing --settings=build_type=${BUILD_TYPE} --profile=${CONAN2_PROFILE_NAME}"
CONAN2_SHOW_PROFLE_COMMAND="conan profile show -pr ${CONAN2_PROFILE_NAME}"
printf "[ conan-install | info > Creating default profile...\n\n"
eval $CONAN2_DEFAULT_PROFILE_COMMAND
if [ ! $? -eq 0 ]; then
printf "[ conan-install | error > Exit code 126 (Command invoked cannot execute):\n ${CONAN2_DEFAULT_PROFILE_COMMAND}.\n"
exit 126
fi
printf "[ conan-install | info > Creating '${CONAN2_PROFILE_NAME}' profile...\n\n"
eval $CONAN2_PROFILE_COMMAND
if [ ! $? -eq 0 ]; then
printf "[ conan-install | error > Exit code 126 (Command invoked cannot execute):\n ${CONAN2_PROFILE_COMMAND}.\n"
exit 126
fi
printf "\n\n[ conan-install | info > Editing '${CONAN2_PROFILE_NAME}' profile\n"
#========== cppstd setting ============
STD_NUMBER=`grep -n -m 1 "compiler.cppstd=" ${CONAN2_HOME}/profiles/${CONAN2_PROFILE_NAME} | cut -d: -f1`
SED_STD_COMMAND="sed -i '${STD_NUMBER}s/.*/compiler.cppstd=gnu17/' ${CONAN2_HOME}/profiles/${CONAN2_PROFILE_NAME}"
eval $SED_STD_COMMAND
if [ ! $? -eq 0 ]; then
printf "[ conan-install | error > Exit code 126 (Command invoked cannot execute):\n ${SED_STD_COMMAND}\n"
exit 126
fi
#========== libcxx setting ============
CXX_NUMBER=`grep -n -m 1 "compiler.libcxx=" ${CONAN2_HOME}/profiles/${CONAN2_PROFILE_NAME} | cut -d: -f1`
SED_CXX_COMMAND="sed -i '${CXX_NUMBER}s/.*/compiler.libcxx=libstdc++11/' ${CONAN2_HOME}/profiles/${CONAN2_PROFILE_NAME}"
eval $SED_CXX_COMMAND
if [ ! $? -eq 0 ]; then
printf "[ conan-install | error > Exit code 126 (Command invoked cannot execute):\n ${SED_CXX_COMMAND}\n"
exit 126
fi
eval $CONAN2_SHOW_PROFLE_COMMAND
if [ ! $? -eq 0 ]; then
printf "[ conan-install | error > Exit code 126 (Command invoked cannot execute):\n ${CONAN2_SHOW_PROFLE_COMMAND}\n"
exit 126
fi
printf "[ conan-install | info > ${CONAN2_INSTALL_COMMAND}\n"
eval ${CONAN2_INSTALL_COMMAND}
if [ ! $? -eq 0 ]; then
printf "[ conan-install | error > Exit code 126 (Command invoked cannot execute):\n ${CONAN2_INSTALL_COMMAND}\n"
exit 126
fi
CORSIKA_CMAKE_SCRIPT="#!/bin/bash
function show_usage(){
printf \"\n\nUsage:\"
printf \"\$0 options [parameters]\n\"
printf \"\n\"
printf \"Options:\n\"
printf \"\n -c or --cmake-flags:\n Additional flags and settings to cmake base command. Default is empty string.\"
printf \"\n\nExample: ./corsika-cmake.sh --cmake-flags '-DUSE_Pythia8_C8=C8' \"
printf \"\n\nNote: the source directory (the one containing CMakeLists.txt), CMAKE_BUILD_TYPE, CMAKE_POLICY_DEFAULT_CMP0091 and CMAKE_TOOLCHAIN_FILE are already set. Do not repeat them.\"
printf \"\n -h or --help:\n Prints this message.\n\"
exit 0
}
echo \"|---------------------------------------------------|\"
echo \"|-----------------[ CORSIKA 8 ]---------------------|\"
echo \"|----------[ CMAKE CONFIGURATION SCRIPT ]---------- |\"
echo \"|---------------------------------------------------|\"
echo \"|-------------------- BEGIN ------------------------|\"
CMAKE_BASE_SETTINGS=\"cmake -S ${CORSIKA_DIR} -D CONAN_CMAKE_DIR=${CORSIKA_DIR}/${CONAN2_OUTPUT_FOLDER_NAME} -D CMAKE_TOOLCHAIN_FILE=${CORSIKA_DIR}/${CONAN2_OUTPUT_FOLDER_NAME}/conan_toolchain.cmake -D CMAKE_POLICY_DEFAULT_CMP0091=NEW -DCMAKE_BUILD_TYPE=${BUILD_TYPE}\"
CMAKE_ADDITIONAL_SETTINGS=\"\"
if [[ \"\$1\" == \"--help\" ]] || [[ \"\$1\" == \"-h\" ]];then
show_usage
fi
while [ ! -z \"\$1\" ]; do
case \"\$1\" in
--cmake-flags|-c)
shift
CMAKE_ADDITIONAL_SETTINGS=\${1}
;;
--help|-h)
#shift
show_usage
;;
*)
show_usage
;;
esac
if [ \$# -gt 0 ]; then
shift
fi
done
printf \"[corsika-cmake | info > Issuing CMake command :\n\"
printf '[corsika-cmake | info > \e[1;36m%s\e[0m\n\n\n' \"\${CMAKE_BASE_SETTINGS} \${CMAKE_ADDITIONAL_SETTINGS}\"
eval \"\${CMAKE_BASE_SETTINGS} \${CMAKE_ADDITIONAL_SETTINGS}\"
if [ ! \$? -eq 0 ]; then
printf \"[ corsika | error > Project configuration (CMake) failed. \n\"
exit 1
fi
echo \"|---------------------------------------------------|\"
echo \"|-----------------[ CORSIKA 8 ]---------------------|\"
echo \"|----------[ CMAKE CONFIGURATION SCRIPT ]---------- |\"
echo \"|---------------------------------------------------|\"
echo \"|-------------------- END --------------------------|\"
"
printf "%s" "${CORSIKA_CMAKE_SCRIPT}" > ${SCRIPT_DIR}/corsika-cmake.sh
chmod +x ${SCRIPT_DIR}/corsika-cmake.sh
printf "\n\n[ conan-install | info > Copy and paste the commands below in the corsika-build directory:\n\n"
printf '\e[1;36m%s\e[0m\n' "> ${SCRIPT_DIR}/corsika-cmake.sh "
printf '\e[1;36m%s\e[0m\n' "> make -j8"
echo " "
echo "|---------------------------------------------------|"
echo "|-----------------[ CORSIKA 8 ]---------------------|"
echo "|-----[ CONAN2 DEPENDENCIES INSTALL SCRIPT ]------- |"
echo "|---------------------------------------------------|"
echo "|-------------------- END --------------------------| "
echo " "
exit 0
from conan import ConanFile
from conan.tools.cmake import cmake_layout,CMakeToolchain, CMakeDeps
class Pkg(ConanFile):
generators = "CMakeDeps", #"CMakeToolchain",
settings = "os", "arch", "compiler", "build_type"
default_options = {
'readline*:shared': 'True',
'arrow*:shared': 'False',
'arrow*:parquet': 'True',
'arrow*:fPIC': 'False',
'arrow*:with_re2': 'True',
'arrow*:with_protobuf': 'False',
'arrow*:with_openssl': 'False',
'arrow*:with_gflags': 'False',
'arrow*:with_glog': 'False',
'arrow*:with_grpc': 'False',
'arrow*:with_utf8proc': 'False',
'arrow*:with_zstd': 'False',
'arrow*:with_bz2': 'False',
'arrow*:with_lz4': 'True',
'arrow*:with_thrift': 'True',
'arrow*:with_boost': 'True',
'boost*:without_container': 'True',
'boost*:without_context': 'True',
'boost*:without_contract': 'True',
'boost*:without_coroutine': 'True',
'boost*:without_date_time': 'True',
'boost*:without_fiber': 'True',
'boost*:without_filesystem': 'False',
'boost*:without_graph': 'True',
'boost*:without_graph_parallel': 'True',
'boost*:without_iostreams': 'False',
'boost*:without_json': 'True',
'boost*:without_locale': 'True',
'boost*:without_log': 'True',
'boost*:without_math': 'False',
'boost*:without_mpi': 'True',
'boost*:without_nowide': 'True',
'boost*:without_program_options': 'True',
'boost*:without_python': 'True',
'boost*:without_serialization': 'False',
'boost*:without_stacktrace': 'True',
'boost*:without_system': 'False',
'boost*:without_test': 'True',
'boost*:without_thread': 'True',
'boost*:without_timer': 'True',
'boost*:without_type_erasure': 'True',
'boost*:without_wave': 'True'
}
def configure(self):
self.options['arrow'].with_boost = True
self.options['arrow'].parquet = True
self.options['arrow'].with_thrift = True
def requirements(self):
self.requires("spdlog/1.14.1", force=True)
self.requires("catch2/3.6.0")
self.requires("bzip2/1.0.8")
self.requires("boost/1.85.0", force=True)
self.requires("eigen/3.4.0")
self.requires("zlib/1.3.1")
self.requires("yaml-cpp/0.8.0")
self.requires("cli11/1.9.1")
self.requires("arrow/16.1.0")
self.requires("proposal/7.6.2")
def build_requirements(self):
self.tool_requires("readline/8.0")
self.tool_requires("bison/[>1.0]")
def generate(self):
tc = CMakeToolchain(self)
tc.absolute_paths = True
tc.generate()
[requires]
spdlog/1.9.2
catch2/2.13.8
eigen/3.3.8
boost/1.76.0
zlib/1.2.11
proposal/7.3.1
yaml-cpp/0.6.3
arrow/2.0.0
cli11/1.9.1
[build_requires]
readline/8.0 # this is only needed to fix a missing dependency in "bison" which is pulled-in by arrow
bison/[>1.0] # needed for arrow, and we HAVE to compile it
[generators]
cmake
[options]
readline:shared=True
arrow:shared=False
arrow:parquet=True
arrow:fPIC=False
arrow:with_re2=False
arrow:with_protobuf=False
arrow:with_openssl=False
arrow:with_gflags=False
arrow:with_glog=False
arrow:with_grpc=False
arrow:with_utf8proc=False
arrow:with_zstd=False
arrow:with_bz2=False
/* /*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
* *
* This software is distributed under the terms of the GNU General Public * This software is distributed under the terms of the 3-clause BSD license.
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of * See file LICENSE for a full version of the license.
* the license.
*/ */
// Another possibility: // Another possibility:
......
/* /*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
* *
* This software is distributed under the terms of the GNU General Public * This software is distributed under the terms of the 3-clause BSD license.
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of * See file LICENSE for a full version of the license.
* the license.
*/ */
#pragma once #pragma once
......
/* /*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
* *
* This software is distributed under the terms of the GNU General Public * This software is distributed under the terms of the 3-clause BSD license.
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of * See file LICENSE for a full version of the license.
* the license.
*/ */
#pragma once #pragma once
#include <corsika/framework/core/PhysicalUnits.hpp> #include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/core/Step.hpp>
#include <corsika/framework/process/ProcessReturn.hpp> #include <corsika/framework/process/ProcessReturn.hpp>
#include <corsika/framework/process/ContinuousProcessStepLength.hpp> #include <corsika/framework/process/ContinuousProcessStepLength.hpp>
...@@ -42,10 +42,15 @@ namespace corsika { ...@@ -42,10 +42,15 @@ namespace corsika {
, sequence_(pl) , sequence_(pl)
, output_(out) , output_(out)
, stack_(stack) , stack_(stack)
, forceInteraction_(false) { , forceInteraction_(false)
, forceDecay_(false) {
CORSIKA_LOG_INFO(c8_ascii_); CORSIKA_LOG_INFO(c8_ascii_);
CORSIKA_LOG_INFO("This is CORSIKA {}.{}.{}.{}", CORSIKA_RELEASE_NUMBER, CORSIKA_LOG_INFO("This is CORSIKA {}.{}.{}.{}", CORSIKA_RELEASE_NUMBER,
CORSIKA_MAJOR_NUMBER, CORSIKA_MINOR_NUMBER, CORSIKA_PATCH_NUMBER); CORSIKA_MAJOR_NUMBER, CORSIKA_MINOR_NUMBER, CORSIKA_PATCH_NUMBER);
CORSIKA_LOG_INFO(
"The C8 author list can be found at: "
"https://gitlab.iap.kit.edu/AirShowerPhysics/corsika/-/wikis/"
"Current-CORSIKA-8-author-list");
CORSIKA_LOG_INFO("Tracking algorithm: {} (version {})", TTracking::getName(), CORSIKA_LOG_INFO("Tracking algorithm: {} (version {})", TTracking::getName(),
TTracking::getVersion()); TTracking::getVersion());
if constexpr (stack_view_type::has_event) { if constexpr (stack_view_type::has_event) {
...@@ -92,6 +97,21 @@ namespace corsika { ...@@ -92,6 +97,21 @@ namespace corsika {
template <typename TTracking, typename TProcessList, typename TOutput, typename TStack> template <typename TTracking, typename TProcessList, typename TOutput, typename TStack>
inline void Cascade<TTracking, TProcessList, TOutput, TStack>::forceInteraction() { inline void Cascade<TTracking, TProcessList, TOutput, TStack>::forceInteraction() {
forceInteraction_ = true; forceInteraction_ = true;
if (forceDecay_) {
CORSIKA_LOG_ERROR("Cannot set forceInteraction when forceDecay is already set");
throw std::runtime_error(
"Cannot set forceInteraction when forceDecay is already set");
}
}
template <typename TTracking, typename TProcessList, typename TOutput, typename TStack>
inline void Cascade<TTracking, TProcessList, TOutput, TStack>::forceDecay() {
forceDecay_ = true;
if (forceInteraction_) {
CORSIKA_LOG_ERROR("Cannot set forceDecay when forceInteraction is already set");
throw std::runtime_error(
"Cannot set forceDecay when forceInteraction is already set");
}
} }
template <typename TTracking, typename TProcessList, typename TOutput, typename TStack> template <typename TTracking, typename TProcessList, typename TOutput, typename TStack>
...@@ -115,21 +135,37 @@ namespace corsika { ...@@ -115,21 +135,37 @@ namespace corsika {
FourMomentum const projectileP4{Elab, particle.getMomentum()}; FourMomentum const projectileP4{Elab, particle.getMomentum()};
// determine combined full inelastic cross section of the particles in the material // determine combined full inelastic cross section of the particles in the material
auto const targetMomentum = MomentumVector{
particle.getMomentum().getCoordinateSystem(), {0_GeV, 0_GeV, 0_GeV}};
auto const xs_function = [&](Code const targetId) -> CrossSectionType {
FourMomentum const targetP4{get_mass(targetId), targetMomentum};
return sequence_.getCrossSection(particle, targetId, targetP4);
};
CrossSectionType const total_cx = CrossSectionType const total_cx_pre = composition.getWeightedSum(xs_function);
composition.getWeightedSum([=](Code const targetId) -> CrossSectionType {
FourMomentum const targetP4(
get_mass(targetId),
MomentumVector(particle.getMomentum().getCoordinateSystem(),
{0_GeV, 0_GeV, 0_GeV}));
return sequence_.getCrossSection(particle, targetId, targetP4);
});
if (forceInteraction_) { if (forceInteraction_) {
CORSIKA_LOG_TRACE("forced interaction!"); CORSIKA_LOG_TRACE("forced interaction!");
forceInteraction_ = false; // just one (first) interaction forceInteraction_ = false; // just one (first) interaction
stack_view_type secondaries(particle); stack_view_type secondaries(particle);
interaction(secondaries, projectileP4, composition, total_cx); interaction(secondaries, projectileP4, composition, total_cx_pre);
sequence_.doSecondaries(secondaries);
particle.erase(); // primary particle is done
return;
}
if (forceDecay_) {
CORSIKA_LOG_TRACE("forced decay!");
forceDecay_ = false; // just one decay
stack_view_type secondaries(particle);
decay(secondaries, sequence_.getInverseLifetime(particle));
if (secondaries.getSize() == 1 && secondaries.getProjectile().getPID() ==
secondaries.getNextParticle().getPID()) {
throw std::runtime_error(
fmt::format("Particle {} decays into itself!",
get_name(secondaries.getProjectile().getPID())));
}
sequence_.doSecondaries(secondaries); sequence_.doSecondaries(secondaries);
particle.erase(); // primary particle is done particle.erase(); // primary particle is done
return; return;
...@@ -137,10 +173,10 @@ namespace corsika { ...@@ -137,10 +173,10 @@ namespace corsika {
// calculate interaction length in medium // calculate interaction length in medium
GrammageType const total_lambda = GrammageType const total_lambda =
(composition.getAverageMassNumber() * constants::u) / total_cx; (composition.getAverageMassNumber() * constants::u) / total_cx_pre;
// sample random exponential step length in grammage // sample random exponential step length in grammage
ExponentialDistribution expDist(total_lambda); ExponentialDistribution expDist{total_lambda};
GrammageType const next_interact = expDist(rng_); GrammageType const next_interact = expDist(rng_);
CORSIKA_LOG_DEBUG("total_lambda={} g/cm2, next_interact={} g/cm2", CORSIKA_LOG_DEBUG("total_lambda={} g/cm2, next_interact={} g/cm2",
...@@ -148,31 +184,31 @@ namespace corsika { ...@@ -148,31 +184,31 @@ namespace corsika {
double(next_interact / 1_g * 1_cm * 1_cm)); double(next_interact / 1_g * 1_cm * 1_cm));
// determine combined total inverse decay time // determine combined total inverse decay time
InverseTimeType const total_inv_lifetime = sequence_.getInverseLifetime(particle); InverseTimeType const total_inv_lifetime_pre = sequence_.getInverseLifetime(particle);
// sample random exponential decay time // sample random exponential decay time
ExponentialDistribution expDistDecay(1 / total_inv_lifetime); ExponentialDistribution expDistDecay(1 / total_inv_lifetime_pre);
TimeType const next_decay = expDistDecay(rng_); TimeType const next_decay = expDistDecay(rng_);
CORSIKA_LOG_DEBUG("total_lifetime={} ns, next_decay={} ns", CORSIKA_LOG_DEBUG("total_lifetime={} ns, next_decay={} ns",
(1 / total_inv_lifetime) / 1_ns, next_decay / 1_ns); (1 / total_inv_lifetime_pre) / 1_ns, next_decay / 1_ns);
// convert next_decay from time to length [m] // convert next_decay from time to length [m]
LengthType const distance_decay = next_decay * particle.getMomentum().getNorm() / LengthType const distance_decay = next_decay * particle.getMomentum().getNorm() /
particle.getEnergy() * constants::c; particle.getEnergy() * constants::c;
// determine geometric tracking // determine geometric tracking
auto [step, nextVol] = tracking_.getTrack(particle); auto [track, nextVol] = tracking_.getTrack(particle);
auto geomMaxLength = step.getLength(1); auto geomMaxLength = track.getLength(1);
// convert next_step from grammage to length // convert next_step from grammage to length
LengthType const distance_interact = LengthType const distance_interact =
currentLogicalNode->getModelProperties().getArclengthFromGrammage(step, currentLogicalNode->getModelProperties().getArclengthFromGrammage(track,
next_interact); next_interact);
// determine the maximum geometric step length // determine the maximum geometric step length
ContinuousProcessStepLength const continuousMaxStep = ContinuousProcessStepLength const continuousMaxStep =
sequence_.getMaxStepLength(particle, step); sequence_.getMaxStepLength(particle, track);
LengthType const continuous_max_dist = continuousMaxStep; LengthType const continuous_max_dist = continuousMaxStep;
// take minimum of geometry, interaction, decay for next step // take minimum of geometry, interaction, decay for next step
...@@ -180,16 +216,14 @@ namespace corsika { ...@@ -180,16 +216,14 @@ namespace corsika {
LengthType const min_non_continuous = std::min(min_discrete, geomMaxLength); LengthType const min_non_continuous = std::min(min_discrete, geomMaxLength);
LengthType const min_distance = std::min(min_non_continuous, continuous_max_dist); LengthType const min_distance = std::min(min_non_continuous, continuous_max_dist);
bool const isContinuous = continuous_max_dist < min_non_continuous;
// inform ContinuousProcesses (if applicable) that it is responsible for step-limit // inform ContinuousProcesses (if applicable) that it is responsible for step-limit
// this would become simpler if we follow the idea of Max to enumerate ALL types of // this would become simpler if we follow the idea of Max to enumerate ALL types of
// processes. Then non-continuous are included and no further logic is needed to // processes. Then non-continuous are included and no further logic is needed to
// distinguish between continuous and non-continuous limit. // distinguish between continuous and non-continuous limit.
ContinuousProcessIndex limitingId; auto const limitingId = isContinuous ? continuousMaxStep : ContinuousProcessIndex{};
bool const isContinuous = continuous_max_dist < min_non_continuous; // // the current step IS limited by a known continuous process
if (isContinuous) {
limitingId =
continuousMaxStep; // the current step IS limited by a known continuous process
}
CORSIKA_LOG_DEBUG( CORSIKA_LOG_DEBUG(
"transport particle by : {} m " "transport particle by : {} m "
...@@ -202,11 +236,12 @@ namespace corsika { ...@@ -202,11 +236,12 @@ namespace corsika {
// move particle along the trajectory to new position // move particle along the trajectory to new position
// also update momentum/direction/time // also update momentum/direction/time
step.setLength(min_distance); track.setLength(min_distance);
Step step{particle, track};
// apply all continuous processes on particle + track // apply all continuous processes on particle + track
if (sequence_.doContinuous(particle, step, limitingId) == if (sequence_.doContinuous(step, limitingId) == ProcessReturn::ParticleAbsorbed) {
ProcessReturn::ParticleAbsorbed) {
CORSIKA_LOG_DEBUG("Cascade: delete absorbed particle PID={} E={} GeV", CORSIKA_LOG_DEBUG("Cascade: delete absorbed particle PID={} E={} GeV",
particle.getPID(), particle.getEnergy() / 1_GeV); particle.getPID(), particle.getEnergy() / 1_GeV);
if (particle.isErased()) { if (particle.isErased()) {
...@@ -218,9 +253,10 @@ namespace corsika { ...@@ -218,9 +253,10 @@ namespace corsika {
} }
return; // particle is gone -> return return; // particle is gone -> return
} }
particle.setTime(particle.getTime() + step.getDuration()); particle.setTime(step.getTimePost());
particle.setPosition(step.getPosition(1)); particle.setPosition(step.getPositionPost());
particle.setDirection(step.getDirection(1)); particle.setDirection(step.getDirectionPost());
particle.setKineticEnergy(step.getEkinPost());
if (isContinuous) { if (isContinuous) {
return; // there is nothing further, step is finished return; // there is nothing further, step is finished
...@@ -267,7 +303,7 @@ namespace corsika { ...@@ -267,7 +303,7 @@ namespace corsika {
"Geometry check: numericalNodeAfterStep={} currentLogicalNode={}", "Geometry check: numericalNodeAfterStep={} currentLogicalNode={}",
fmt::ptr(numericalNodeAfterStep), fmt::ptr(currentLogicalNode)); fmt::ptr(numericalNodeAfterStep), fmt::ptr(currentLogicalNode));
if (numericalNodeAfterStep != currentLogicalNode) { if (numericalNodeAfterStep != currentLogicalNode) {
CORSIKA_LOG_ERROR( CORSIKA_LOG_DEBUG(
"expect to be in node currentLogicalNode={} but are in " "expect to be in node currentLogicalNode={} but are in "
"numericalNodeAfterStep={}. Continue, but without guarantee.", "numericalNodeAfterStep={}. Continue, but without guarantee.",
fmt::ptr(currentLogicalNode), fmt::ptr(numericalNodeAfterStep)); fmt::ptr(currentLogicalNode), fmt::ptr(numericalNodeAfterStep));
...@@ -284,7 +320,7 @@ namespace corsika { ...@@ -284,7 +320,7 @@ namespace corsika {
// secondaries, b) the projectile particle deleted (or // secondaries, b) the projectile particle deleted (or
// changed) // changed)
stack_view_type secondaries(particle); stack_view_type secondaries{particle};
/* /*
Create SecondaryView object on Stack. The data container Create SecondaryView object on Stack. The data container
...@@ -295,12 +331,19 @@ namespace corsika { ...@@ -295,12 +331,19 @@ namespace corsika {
important to use projectile/view (and not particle) for Interaction, important to use projectile/view (and not particle) for Interaction,
and Decay! and Decay!
*/ */
FourMomentum const projectileP4Post{particle.getEnergy(), particle.getMomentum()};
bool eraseParticle =
false; // only erase original particle if it decayed or interacted
if (distance_interact < distance_decay) { if (distance_interact < distance_decay) {
interaction(secondaries, projectileP4, composition, total_cx); eraseParticle = isInteracted(
interaction(secondaries, projectileP4Post, composition, total_cx_pre));
} else { } else {
[[maybe_unused]] auto projectile = secondaries.getProjectile(); [[maybe_unused]] auto projectile = secondaries.getProjectile();
if (decay(secondaries, total_inv_lifetime_pre) == ProcessReturn::Decayed) {
if (decay(secondaries, total_inv_lifetime) == ProcessReturn::Decayed) { eraseParticle = true;
if (secondaries.getSize() == 1 && if (secondaries.getSize() == 1 &&
projectile.getPID() == secondaries.getNextParticle().getPID()) { projectile.getPID() == secondaries.getNextParticle().getPID()) {
throw std::runtime_error(fmt::format("Particle {} decays into itself!", throw std::runtime_error(fmt::format("Particle {} decays into itself!",
...@@ -309,8 +352,11 @@ namespace corsika { ...@@ -309,8 +352,11 @@ namespace corsika {
} }
} }
sequence_.doSecondaries(secondaries); if (eraseParticle) {
particle.erase(); // doSecondaries() makes sense only if there was an actual event
sequence_.doSecondaries(secondaries);
particle.erase();
}
} }
template <typename TTracking, typename TProcessList, typename TOutput, typename TStack> template <typename TTracking, typename TProcessList, typename TOutput, typename TStack>
...@@ -327,7 +373,8 @@ namespace corsika { ...@@ -327,7 +373,8 @@ namespace corsika {
auto const returnCode = sequence_.selectDecay(view, sample_process); auto const returnCode = sequence_.selectDecay(view, sample_process);
if (returnCode != ProcessReturn::Decayed) { if (returnCode != ProcessReturn::Decayed) {
CORSIKA_LOG_DEBUG("Particle did not decay!"); CORSIKA_LOG_ERROR("Particle {} did not decay!",
get_name(view.getProjectile().getPID()));
} }
setEventType(view, history::EventType::Decay); setEventType(view, history::EventType::Decay);
return returnCode; return returnCode;
......
/* -*-c++-*- /*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
* *
* This software is distributed under the terms of the GNU General Public * This software is distributed under the terms of the 3-clause BSD license.
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of * See file LICENSE for a full version of the license.
* the license.
*/ */
#pragma once #pragma once
......
/* /*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
* *
* This software is distributed under the terms of the GNU General Public * This software is distributed under the terms of the 3-clause BSD license.
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of * See file LICENSE for a full version of the license.
* the license.
*/ */
#pragma once #pragma once
...@@ -14,14 +13,13 @@ ...@@ -14,14 +13,13 @@
namespace corsika { namespace corsika {
inline HEPEnergyType constexpr get_kinetic_energy_propagation_threshold( inline HEPEnergyType get_kinetic_energy_propagation_threshold(Code const code) {
Code const code) {
if (is_nucleus(code)) return particle::detail::threshold_nuclei; if (is_nucleus(code)) return particle::detail::threshold_nuclei;
return particle::detail::propagation_thresholds[static_cast<CodeIntType>(code)]; return particle::detail::propagation_thresholds[static_cast<CodeIntType>(code)];
} }
inline void constexpr set_kinetic_energy_propagation_threshold( inline void set_kinetic_energy_propagation_threshold(Code const code,
Code const code, HEPEnergyType const val) { HEPEnergyType const val) {
if (is_nucleus(code)) if (is_nucleus(code))
particle::detail::threshold_nuclei = val; particle::detail::threshold_nuclei = val;
else else
...@@ -33,12 +31,11 @@ namespace corsika { ...@@ -33,12 +31,11 @@ namespace corsika {
return particle::detail::masses[static_cast<CodeIntType>(code)]; return particle::detail::masses[static_cast<CodeIntType>(code)];
} }
inline HEPEnergyType constexpr get_energy_production_threshold(Code const p) { inline HEPEnergyType get_energy_production_threshold(Code const p) {
return particle::detail::production_thresholds[static_cast<CodeIntType>(p)]; return particle::detail::production_thresholds[static_cast<CodeIntType>(p)];
} }
inline void constexpr set_energy_production_threshold(Code const p, inline void set_energy_production_threshold(Code const p, HEPEnergyType const val) {
HEPEnergyType const val) {
particle::detail::production_thresholds[static_cast<CodeIntType>(p)] = val; particle::detail::production_thresholds[static_cast<CodeIntType>(p)] = val;
} }
...@@ -134,6 +131,11 @@ namespace corsika { ...@@ -134,6 +131,11 @@ namespace corsika {
if (std::abs(k) <= maxPDG) { if (std::abs(k) <= maxPDG) {
return particle::detail::conversionArray[k + maxPDG]; return particle::detail::conversionArray[k + maxPDG];
} else { } else {
if (1000000000 <= k && k <= 1009999990) { // nucleus (no L or I)
int const Z = (k - 1000000000) / 10000;
int const A = (k - 1000000000 - 10000 * Z) / 10;
return get_nucleus_code(A, Z);
}
return particle::detail::conversionMap.at(p); return particle::detail::conversionMap.at(p);
} }
} }
...@@ -149,7 +151,7 @@ namespace corsika { ...@@ -149,7 +151,7 @@ namespace corsika {
return get_mass(Code::Proton) * Z + (A - Z) * get_mass(Code::Neutron); return get_mass(Code::Proton) * Z + (A - Z) * get_mass(Code::Neutron);
} }
inline std::string_view get_nucleus_name(Code const code) { inline std::string get_nucleus_name(Code const code) {
size_t const A = get_nucleus_A(code); size_t const A = get_nucleus_A(code);
size_t const Z = get_nucleus_Z(code); size_t const Z = get_nucleus_Z(code);
return fmt::format("Nucleus_A{}_Z{}", A, Z); return fmt::format("Nucleus_A{}_Z{}", A, Z);
......
/*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the 3-clause BSD license.
* See file LICENSE for a full version of the license.
*/
#pragma once
#include <fmt/format.h>
#include <corsika/framework/core/ParticleProperties.hpp>
#include <boost/filesystem/path.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <Eigen/Dense>
//-----------------------------------
// STD
//-----------------------------------
namespace std {
auto inline format_as(std::_Put_time<char> const& arg) {
std::ostringstream os;
os << arg;
return os.str();
}
} // namespace std
//-----------------------------------
// CORSIKA
//-----------------------------------
namespace corsika {
// formatters for particle codes declared on ParticleProperties.hpp
auto inline format_as(Code code) { return get_name(code); }
auto inline format_as(PDGCode code) { return fmt::underlying(code); }
template <typename Type>
auto inline format_as(Type const& arg) {
std::ostringstream os;
os << arg;
return os.str();
}
} // namespace corsika
//-----------------------------------
// boost::filesystem
//-----------------------------------
namespace boost::filesystem {
auto inline format_as(path const& fname) { return fname.string().c_str(); }
} // namespace boost::filesystem
//-----------------------------------
// phys::units
//-----------------------------------
namespace phys::units {
template <typename Dimensions>
auto inline format_as(phys::units::quantity<Dimensions> const& arg) {
return io::to_string(arg);
}
} // namespace phys::units
//----------------------------------
// Eigen
//----------------------------------
namespace Eigen {
template <typename Scalar, int M, int N>
auto inline format_as(Matrix<Scalar, M, N> const& arg) {
std::ostringstream os;
os << arg;
return os.str();
}
template <typename Matrix1, typename Matrix2>
auto inline format_as(Product<Matrix1, Matrix2> const& arg) {
std::ostringstream os;
os << arg;
return os.str();
}
} // namespace Eigen
/* /*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
* *
* This software is distributed under the terms of the GNU General Public * This software is distributed under the terms of the 3-clause BSD license.
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of * See file LICENSE for a full version of the license.
* the license.
*/ */
#pragma once #pragma once
......
/* /*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
* *
* This software is distributed under the terms of the GNU General Public * This software is distributed under the terms of the 3-clause BSD license.
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of * See file LICENSE for a full version of the license.
* the license.
*/ */
#pragma once #pragma once
......
/* /*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
* *
* This software is distributed under the terms of the GNU General Public * This software is distributed under the terms of the 3-clause BSD license.
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of * See file LICENSE for a full version of the license.
* the license.
*/ */
#pragma once #pragma once
...@@ -81,7 +80,13 @@ namespace corsika { ...@@ -81,7 +80,13 @@ namespace corsika {
template <typename TDim> template <typename TDim>
inline CoordinateSystemPtr make_rotationToZ(CoordinateSystemPtr const& cs, inline CoordinateSystemPtr make_rotationToZ(CoordinateSystemPtr const& cs,
Vector<TDim> const& vVec) { Vector<TDim> const& vVec) {
auto const a = vVec.normalized().getComponents(cs).getEigenVector(); auto const vVecComp = vVec.getComponents(cs);
if (vVecComp.getX().magnitude() == 0 && vVecComp.getY().magnitude() == 0 &&
vVecComp.getZ().magnitude() == 0) {
return cs;
}
auto const a = vVecComp.normalized().getEigenVector();
auto const a1 = a(0), a2 = a(1), a3 = a(2); auto const a1 = a(0), a2 = a(1), a3 = a(2);
Eigen::Matrix3d A, B; Eigen::Matrix3d A, B;
......
/* /*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
* *
* This software is distributed under the terms of the GNU General Public * This software is distributed under the terms of the 3-clause BSD license.
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of * See file LICENSE for a full version of the license.
* the license.
*/ */
#pragma once #pragma once
......
/* /*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
* *
* This software is distributed under the terms of the GNU General Public * This software is distributed under the terms of the 3-clause BSD license.
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of * See file LICENSE for a full version of the license.
* the license.
*/ */
#pragma once #pragma once
......
/* /*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
* *
* This software is distributed under the terms of the GNU General Public * This software is distributed under the terms of the 3-clause BSD license.
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of * See file LICENSE for a full version of the license.
* the license.
*/ */
#pragma once #pragma once
...@@ -29,8 +28,9 @@ namespace corsika { ...@@ -29,8 +28,9 @@ namespace corsika {
} }
inline Point LeapFrogTrajectory::getPosition(double const u) const { inline Point LeapFrogTrajectory::getPosition(double const u) const {
Point position = initialPosition_ + initialVelocity_ * timeStep_ * u / 2; if (u == 0) return initialPosition_;
VelocityVector velocity = Point const position = initialPosition_ + initialVelocity_ * timeStep_ * u / 2;
VelocityVector const velocity =
initialVelocity_ + initialVelocity_.cross(magneticfield_) * timeStep_ * u * k_; initialVelocity_ + initialVelocity_.cross(magneticfield_) * timeStep_ * u * k_;
return position + velocity * timeStep_ * u / 2; return position + velocity * timeStep_ * u / 2;
} }
...@@ -40,6 +40,8 @@ namespace corsika { ...@@ -40,6 +40,8 @@ namespace corsika {
} }
inline DirectionVector LeapFrogTrajectory::getDirection(double const u) const { inline DirectionVector LeapFrogTrajectory::getDirection(double const u) const {
if (u == 0) return initialDirection_;
return (initialDirection_ + return (initialDirection_ +
initialDirection_.cross(magneticfield_) * timeStep_ * u * k_) initialDirection_.cross(magneticfield_) * timeStep_ * u * k_)
.normalized(); .normalized();
......
/* /*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
* *
* This software is distributed under the terms of the GNU General Public * This software is distributed under the terms of the 3-clause BSD license.
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of * See file LICENSE for a full version of the license.
* the license.
*/ */
#pragma once #pragma once
......