diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt
index 621cab01d1facaa687c12f287c12ed2d54a6501f..925cdac6ef4a98b4dc63b2c01696c1d82aad9bc2 100644
--- a/Documentation/Examples/CMakeLists.txt
+++ b/Documentation/Examples/CMakeLists.txt
@@ -30,7 +30,6 @@ target_link_libraries (cascade_example SuperStupidStack CORSIKAunits CORSIKAlogg
ProcessStackInspector
ProcessTrackWriter
ProcessTrackingLine
- ProcessHadronicElasticModel
CORSIKAprocesses
CORSIKAparticles
CORSIKAgeometry
@@ -40,6 +39,23 @@ target_link_libraries (cascade_example SuperStupidStack CORSIKAunits CORSIKAlogg
install (TARGETS cascade_example DESTINATION share/examples)
CORSIKA_ADD_TEST (cascade_example)
+add_executable (boundary_example boundary_example.cc)
+target_compile_options(boundary_example PRIVATE -g) # do not skip asserts
+target_link_libraries (boundary_example SuperStupidStack CORSIKAunits CORSIKAlogging
+ CORSIKArandom
+ ProcessSibyll
+ CORSIKAcascade
+ ProcessTrackWriter
+ ProcessTrackingLine
+ CORSIKAprocesses
+ CORSIKAparticles
+ CORSIKAgeometry
+ CORSIKAenvironment
+ CORSIKAprocesssequence
+ )
+install (TARGETS boundary_example DESTINATION share/examples)
+CORSIKA_ADD_TEST (boundary_example)
+
add_executable (staticsequence_example staticsequence_example.cc)
target_compile_options(staticsequence_example PRIVATE -g) # do not skip asserts
target_link_libraries (staticsequence_example
diff --git a/Documentation/Examples/boundary_example.cc b/Documentation/Examples/boundary_example.cc
new file mode 100644
index 0000000000000000000000000000000000000000..0c767160e08123292379ad330941c2f37479b627
--- /dev/null
+++ b/Documentation/Examples/boundary_example.cc
@@ -0,0 +1,342 @@
+
+/*
+ * (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/tracking_line/TrackingLine.h>
+
+#include <corsika/setup/SetupEnvironment.h>
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+
+#include <corsika/environment/Environment.h>
+#include <corsika/environment/HomogeneousMedium.h>
+#include <corsika/environment/NuclearComposition.h>
+
+#include <corsika/geometry/Sphere.h>
+
+#include <corsika/process/sibyll/Decay.h>
+#include <corsika/process/sibyll/Interaction.h>
+#include <corsika/process/sibyll/NuclearInteraction.h>
+
+#include <corsika/process/track_writer/TrackWriter.h>
+
+#include <corsika/units/PhysicalUnits.h>
+
+#include <corsika/random/RNGManager.h>
+
+#include <corsika/utl/CorsikaFenv.h>
+
+#include <iostream>
+#include <limits>
+#include <typeinfo>
+
+using namespace corsika;
+using namespace corsika::process;
+using namespace corsika::units;
+using namespace corsika::particles;
+using namespace corsika::random;
+using namespace corsika::setup;
+using namespace corsika::geometry;
+using namespace corsika::environment;
+
+using namespace std;
+using namespace corsika::units::si;
+
+class ProcessCut : public process::ContinuousProcess<ProcessCut> {
+
+ HEPEnergyType fECut;
+
+ HEPEnergyType fEnergy = 0_GeV;
+ HEPEnergyType fEmEnergy = 0_GeV;
+ int fEmCount = 0;
+ HEPEnergyType fInvEnergy = 0_GeV;
+ int fInvCount = 0;
+
+public:
+ ProcessCut(const HEPEnergyType cut)
+ : fECut(cut) {}
+
+ template <typename Particle>
+ bool isBelowEnergyCut(Particle& p) const {
+ // nuclei
+ if (p.GetPID() == particles::Code::Nucleus) {
+ auto const ElabNuc = p.GetEnergy() / p.GetNuclearA();
+ auto const EcmNN = sqrt(2. * ElabNuc * 0.93827_GeV);
+ if (ElabNuc < fECut || EcmNN < 10_GeV)
+ return true;
+ else
+ return false;
+ } else {
+ // TODO: center-of-mass energy hard coded
+ const HEPEnergyType Ecm = sqrt(2. * p.GetEnergy() * 0.93827_GeV);
+ if (p.GetEnergy() < fECut || Ecm < 10_GeV)
+ return true;
+ else
+ return false;
+ }
+ }
+
+ bool isEmParticle(Code pCode) const {
+ bool is_em = false;
+ // FOR NOW: switch
+ switch (pCode) {
+ case Code::Electron:
+ is_em = true;
+ break;
+ case Code::Positron:
+ is_em = true;
+ break;
+ case Code::Gamma:
+ is_em = true;
+ break;
+ default:
+ break;
+ }
+ return is_em;
+ }
+
+ void defineEmParticles() const {
+ // create bool array identifying em particles
+ }
+
+ bool isInvisible(Code pCode) const {
+ bool is_inv = false;
+ // FOR NOW: switch
+ switch (pCode) {
+ case Code::NuE:
+ is_inv = true;
+ break;
+ case Code::NuEBar:
+ is_inv = true;
+ break;
+ case Code::NuMu:
+ is_inv = true;
+ break;
+ case Code::NuMuBar:
+ is_inv = true;
+ break;
+ case Code::MuPlus:
+ is_inv = true;
+ break;
+ case Code::MuMinus:
+ is_inv = true;
+ break;
+
+ case Code::Neutron:
+ is_inv = true;
+ break;
+
+ case Code::AntiNeutron:
+ is_inv = true;
+ break;
+
+ default:
+ break;
+ }
+ return is_inv;
+ }
+
+ template <typename Particle>
+ LengthType MaxStepLength(Particle& p, setup::Trajectory&) const {
+ cout << "ProcessCut: MinStep: pid: " << p.GetPID() << endl;
+ cout << "ProcessCut: MinStep: energy (GeV): " << p.GetEnergy() / 1_GeV << endl;
+ const Code pid = p.GetPID();
+ if (isEmParticle(pid) || isInvisible(pid) || isBelowEnergyCut(p)) {
+ cout << "ProcessCut: MinStep: next cut: " << 0. << endl;
+ return 0_m;
+ } else {
+ LengthType next_step = 1_m * std::numeric_limits<double>::infinity();
+ cout << "ProcessCut: MinStep: next cut: " << next_step << endl;
+ return next_step;
+ }
+ }
+
+ template <typename Particle, typename Stack>
+ EProcessReturn DoContinuous(Particle& p, setup::Trajectory&, Stack&) {
+ const Code pid = p.GetPID();
+ HEPEnergyType energy = p.GetEnergy();
+ cout << "ProcessCut: DoContinuous: " << pid << " E= " << energy
+ << ", EcutTot=" << (fEmEnergy + fInvEnergy + fEnergy) / 1_GeV << " GeV" << endl;
+ EProcessReturn ret = EProcessReturn::eOk;
+ if (isEmParticle(pid)) {
+ cout << "removing em. particle..." << endl;
+ fEmEnergy += energy;
+ fEmCount += 1;
+ // p.Delete();
+ ret = EProcessReturn::eParticleAbsorbed;
+ } else if (isInvisible(pid)) {
+ cout << "removing inv. particle..." << endl;
+ fInvEnergy += energy;
+ fInvCount += 1;
+ // p.Delete();
+ ret = EProcessReturn::eParticleAbsorbed;
+ } else if (isBelowEnergyCut(p)) {
+ cout << "removing low en. particle..." << endl;
+ fEnergy += energy;
+ // p.Delete();
+ ret = EProcessReturn::eParticleAbsorbed;
+ }
+ return ret;
+ }
+
+ void Init() {
+ fEmEnergy = 0. * 1_GeV;
+ fEmCount = 0;
+ fInvEnergy = 0. * 1_GeV;
+ fInvCount = 0;
+ fEnergy = 0. * 1_GeV;
+ // defineEmParticles();
+ }
+
+ void ShowResults() {
+ cout << " ******************************" << endl
+ << " ParticleCut: " << endl
+ << " energy in em. component (GeV): " << fEmEnergy / 1_GeV << endl
+ << " no. of em. particles injected: " << fEmCount << endl
+ << " energy in inv. component (GeV): " << fInvEnergy / 1_GeV << endl
+ << " no. of inv. particles injected: " << fInvCount << endl
+ << " energy below particle cut (GeV): " << fEnergy / 1_GeV << endl
+ << " ******************************" << endl;
+ }
+
+ HEPEnergyType GetInvEnergy() const { return fInvEnergy; }
+ HEPEnergyType GetCutEnergy() const { return fEnergy; }
+ HEPEnergyType GetEmEnergy() const { return fEmEnergy; }
+};
+
+template <bool deleteParticle>
+struct MyBoundaryCrossingProcess
+ : public BoundaryCrossingProcess<MyBoundaryCrossingProcess<deleteParticle>> {
+
+ MyBoundaryCrossingProcess(std::string const& filename) { fFile.open(filename); }
+
+ template <typename Particle>
+ EProcessReturn DoBoundaryCrossing(Particle& p,
+ typename Particle::BaseNodeType const& from,
+ typename Particle::BaseNodeType const& to) {
+ std::cout << "boundary crossing! from: " << &from << "; to: " << &to << std::endl;
+
+ auto const& name = particles::GetName(p.GetPID());
+ auto const start = p.GetPosition().GetCoordinates();
+
+ fFile << name << " " << start[0] / 1_m << ' ' << start[1] / 1_m << ' '
+ << start[2] / 1_m << '\n';
+
+ if constexpr (deleteParticle) { p.Delete(); }
+
+ return EProcessReturn::eOk;
+ }
+
+ void Init() {}
+
+private:
+ std::ofstream fFile;
+};
+
+//
+// The example main program for a particle cascade
+//
+int main() {
+ feenableexcept(FE_INVALID);
+ // initialize random number sequence(s)
+ random::RNGManager::GetInstance().RegisterRandomStream("cascade");
+
+ // setup environment, geometry
+ using EnvType = environment::Environment<setup::IEnvironmentModel>;
+ EnvType env;
+ auto& universe = *(env.GetUniverse());
+
+ const CoordinateSystem& rootCS = env.GetCoordinateSystem();
+
+ auto outerMedium = EnvType::CreateNode<Sphere>(
+ Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
+
+ // fraction of oxygen
+ auto const props =
+ outerMedium
+ ->SetModelProperties<environment::HomogeneousMedium<setup::IEnvironmentModel>>(
+ 1_kg / (1_m * 1_m * 1_m),
+ environment::NuclearComposition(
+ std::vector<particles::Code>{particles::Code::Proton},
+ std::vector<float>{1.f}));
+
+ auto innerMedium = EnvType::CreateNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 5_km);
+
+ innerMedium->SetModelProperties(props);
+
+ outerMedium->AddChild(std::move(innerMedium));
+
+ universe.AddChild(std::move(outerMedium));
+
+ // setup processes, decays and interactions
+ tracking_line::TrackingLine tracking;
+
+ random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
+ process::sibyll::Interaction sibyll;
+ process::sibyll::NuclearInteraction sibyllNuc(sibyll);
+ process::sibyll::Decay decay;
+ ProcessCut cut(20_GeV);
+
+ process::TrackWriter::TrackWriter trackWriter("tracks.dat");
+ MyBoundaryCrossingProcess<true> boundaryCrossing("crossings.dat");
+
+ // assemble all processes into an ordered process list
+ auto sequence = sibyll << sibyllNuc << decay << cut << boundaryCrossing << trackWriter;
+
+ // setup particle stack, and add primary particles
+ setup::Stack stack;
+ stack.Clear();
+ const Code beamCode = Code::Proton;
+ const HEPMassType mass = particles::GetMass(Code::Proton);
+ const HEPEnergyType E0 = 50_TeV;
+
+ std::uniform_real_distribution distTheta(0., 180.);
+ std::uniform_real_distribution distPhi(0., 360.);
+ std::mt19937 rng;
+
+ for (int i = 0; i < 100; ++i) {
+ auto const theta = distTheta(rng);
+ auto const phi = distPhi(rng);
+
+ auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
+ return sqrt((Elab - m) * (Elab + m));
+ };
+ HEPMomentumType P0 = elab2plab(E0, mass);
+ auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
+ return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
+ -ptot * cos(theta));
+ };
+ auto const [px, py, pz] =
+ momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0);
+ auto plab = corsika::stack::MomentumVector(rootCS, {px, py, pz});
+ cout << "input particle: " << beamCode << endl;
+ cout << "input angles: theta=" << theta << " phi=" << phi << endl;
+ cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
+ Point pos(rootCS, 0_m, 0_m, 0_m);
+ stack.AddParticle(
+ std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+ beamCode, E0, plab, pos, 0_ns});
+ }
+
+ // define air shower object, run simulation
+ cascade::Cascade EAS(env, tracking, sequence, stack);
+ EAS.Init();
+ EAS.Run();
+
+ cout << "Result: E0=" << E0 / 1_GeV << endl;
+ cut.ShowResults();
+ const HEPEnergyType Efinal =
+ cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy();
+ cout << "total energy (GeV): " << Efinal / 1_GeV << endl
+ << "relative difference (%): " << (Efinal / E0 - 1.) * 100 << endl;
+}
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 18c03158ab46f9909d06860f04cb2a173a88ae00..67ff9c8fa14173764044fc5b60749e6dfd9e6606 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -11,8 +11,6 @@
#include <corsika/cascade/Cascade.h>
#include <corsika/process/ProcessSequence.h>
-#include <corsika/process/hadronic_elastic_model/HadronicElasticModel.h>
-#include <corsika/process/stack_inspector/StackInspector.h>
#include <corsika/process/tracking_line/TrackingLine.h>
#include <corsika/setup/SetupEnvironment.h>
@@ -21,7 +19,6 @@
#include <corsika/environment/Environment.h>
#include <corsika/environment/HomogeneousMedium.h>
-#include <corsika/environment/NameModel.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/geometry/Sphere.h>
@@ -216,35 +213,6 @@ public:
HEPEnergyType GetEmEnergy() const { return fEmEnergy; }
};
-template <bool deleteParticle>
-struct MyBoundaryCrossingProcess
- : public BoundaryCrossingProcess<MyBoundaryCrossingProcess<deleteParticle>> {
-
- MyBoundaryCrossingProcess(std::string const& filename) { fFile.open(filename); }
-
- template <typename Particle>
- EProcessReturn DoBoundaryCrossing(Particle& p,
- typename Particle::BaseNodeType const& from,
- typename Particle::BaseNodeType const& to) {
- std::cout << "boundary crossing! from: " << &from << "; to: " << &to << std::endl;
-
- auto const& name = particles::GetName(p.GetPID());
- auto const start = p.GetPosition().GetCoordinates();
-
- fFile << name << " " << start[0] / 1_m << ' ' << start[1] / 1_m << ' '
- << start[2] / 1_m << '\n';
-
- if constexpr (deleteParticle) { p.Delete(); }
-
- return EProcessReturn::eOk;
- }
-
- void Init() {}
-
-private:
- std::ofstream fFile;
-};
-
//
// The example main program for a particle cascade
//
@@ -284,7 +252,6 @@ int main() {
// setup processes, decays and interactions
tracking_line::TrackingLine tracking;
- stack_inspector::StackInspector<setup::Stack> p0(true);
random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
process::sibyll::Interaction sibyll;
@@ -293,42 +260,42 @@ int main() {
ProcessCut cut(20_GeV);
process::TrackWriter::TrackWriter trackWriter("tracks.dat");
- MyBoundaryCrossingProcess<false> boundaryCrossing("crossings.dat");
// assemble all processes into an ordered process list
- auto sequence = sibyll << sibyllNuc << decay << cut << boundaryCrossing << trackWriter;
+ auto sequence = sibyll << sibyllNuc << decay << cut << trackWriter;
// setup particle stack, and add primary particle
setup::Stack stack;
stack.Clear();
- const Code beamCode = Code::Proton;
- const int nuclA = 56;
+ const Code beamCode = Code::Nucleus;
+ const int nuclA = 4;
const int nuclZ = int(nuclA / 2.15 + 0.7);
- const HEPMassType mass = particles::Proton::GetMass() * nuclZ +
- (nuclA - nuclZ) * particles::Neutron::GetMass();
- const HEPEnergyType E0 = nuclA * 100_TeV;
- double const phi = 0;
- double const theta = 0;
- auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
- return sqrt(Elab * Elab - m * m);
- };
- HEPMomentumType P0 = elab2plab(E0, mass);
- auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
- return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
- -ptot * cos(theta));
- };
- auto const [px, py, pz] =
- momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0);
- auto plab = corsika::stack::MomentumVector(rootCS, {px, py, pz});
- cout << "input particle: " << beamCode << endl;
- cout << "input angles: theta=" << theta << " phi=" << phi << endl;
- cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
- Point pos(rootCS, 0_m, 0_m, 0_m);
-
- stack.AddParticle(
- std::tuple<particles::Code, units::si::HEPEnergyType,
- corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
- beamCode, E0, plab, pos, 0_ns});
+ const HEPMassType mass = nuclA * particles::GetMass(Code::Proton);
+ const HEPEnergyType E0 = nuclA * 25_TeV;
+ double theta = 0.;
+ double phi = 0.;
+
+ {
+ auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
+ return sqrt((Elab - m) * (Elab + m));
+ };
+ HEPMomentumType P0 = elab2plab(E0, mass);
+ auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
+ return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
+ -ptot * cos(theta));
+ };
+ auto const [px, py, pz] =
+ momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0);
+ auto plab = corsika::stack::MomentumVector(rootCS, {px, py, pz});
+ cout << "input particle: " << beamCode << endl;
+ cout << "input angles: theta=" << theta << " phi=" << phi << endl;
+ cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
+ Point pos(rootCS, 0_m, 0_m, 0_m);
+ stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point,
+ units::si::TimeType, unsigned short, unsigned short>{
+ beamCode, E0, plab, pos, 0_ns, nuclA, nuclZ});
+ }
// define air shower object, run simulation
cascade::Cascade EAS(env, tracking, sequence, stack);
diff --git a/Environment/VolumeTreeNode.h b/Environment/VolumeTreeNode.h
index a9b790035817c7010cd85537cc71d85594e286ba..54bc0cc76d05a473fcdcd2a8cfa3d390f32f865b 100644
--- a/Environment/VolumeTreeNode.h
+++ b/Environment/VolumeTreeNode.h
@@ -90,6 +90,8 @@ namespace corsika::environment {
auto const& GetVolume() const { return *fGeoVolume; }
auto const& GetModelProperties() const { return *fModelProperties; }
+
+ auto const& HasModelProperties() const { return fModelProperties.get() != nullptr; }
template <typename ModelProperties, typename... Args>
auto SetModelProperties(Args&&... args) {
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index f47266e9e5012c35849dfa13f34c82071a0bebf6..4da4609aea684ec8dbefc4ca35230a527aec44a1 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -143,12 +143,10 @@ namespace corsika::cascade {
auto const* currentLogicalNode = particle.GetNode();
- auto const* currentNumericalNode =
- fEnvironment.GetUniverse()->GetContainingNode(particle.GetPosition());
-
- if (currentNumericalNode == &*fEnvironment.GetUniverse()) {
- throw std::runtime_error("particle entered void Universe");
- }
+ // assert that particle stays outside void Universe if it has no
+ // model properties set
+ assert(currentLogicalNode != &*fEnvironment.GetUniverse() ||
+ fEnvironment.GetUniverse()->HasModelProperties());
// convert next_step from grammage to length
LengthType const distance_interact =
@@ -230,7 +228,7 @@ namespace corsika::cascade {
auto const assertion = [&] {
auto const* numericalNodeAfterStep =
- fEnvironment.GetUniverse()->GetContainingNode(particle.GetPosition());
+ fEnvironment.GetUniverse()->GetContainingNode(particle.GetPosition());
return numericalNodeAfterStep == currentLogicalNode;
};
@@ -238,7 +236,7 @@ namespace corsika::cascade {
} else { // boundary crossing, step is limited by volume boundary
std::cout << "boundary crossing! next node = " << nextVol << std::endl;
particle.SetNode(nextVol);
- // DoBoundary may delete the particle (or not)
+ // DoBoundary may delete the particle (or not)
fProcessSequence.DoBoundaryCrossing(particle, *currentLogicalNode, *nextVol);
}
}
diff --git a/Tools/CMakeLists.txt b/Tools/CMakeLists.txt
index c0e4c12b9088e5384462b527158d1315fb47322d..c1347108d1d451651bb8a6a593c4acd7c8aab66d 100644
--- a/Tools/CMakeLists.txt
+++ b/Tools/CMakeLists.txt
@@ -1,4 +1,4 @@
-set (TOOLS_FILES plot_tracks.sh)
+set (TOOLS_FILES plot_tracks.sh plot_crossings.sh)
install (
FILES ${TOOLS_FILES}
diff --git a/Tools/plot_crossings.sh b/Tools/plot_crossings.sh
new file mode 100755
index 0000000000000000000000000000000000000000..d4a318e20ab7fe013ed89afae57b5a3c9220ea98
--- /dev/null
+++ b/Tools/plot_crossings.sh
@@ -0,0 +1,39 @@
+#!/bin/sh
+
+# (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.
+
+# with this script you can plot an animation of output of TrackWriter
+
+crossings_dat=$1
+if [ -z "$crossings_dat" ]; then
+ echo "usage: $0 <crossings.dat> [output.gif]" >&2
+ exit 1
+fi
+
+output=$2
+if [ -z "$output" ]; then
+ output="$crossings_dat.gif"
+fi
+
+cat <<EOF | gnuplot
+set term gif animate size 600,600
+set output "$output"
+
+set xlabel "x / m"
+set ylabel "y / m"
+set zlabel "z / m"
+set title "CORSIKA 8 preliminary"
+
+do for [t=0:360:1] {
+ set view 60, t
+ splot "$crossings_dat" u 2:3:4 w points t ""
+}
+EOF
+
+exit $?