Merge branch 'andre-master' into 'master'

Andre master

See merge request !6

Documentation/Examples/vertical_EAS.cc

@@ -124,7 +124,7 @@ int main(int argc, char** argv) {
   node4->AddChild(std::move(node5));
   node3->AddChild(std::move(node4));
   node2->AddChild(std::move(node3));
-  node1->AddChild(std::move(node2)); 
+  node1->AddChild(std::move(node2));
   env.GetUniverse()->AddChild(std::move(node1));
 
   // setup particle stack, and add primary particle

Framework/Cascade/Cascade.h

@@ -208,7 +208,7 @@ namespace corsika::cascade {
       LengthType const distance_max = fProcessSequence.MaxStepLength(vParticle, step);
       std::cout << "distance_max=" << distance_max << std::endl;
 
-      // take minimum of geometry, interaction, decay, magnetic field for next step
+      // take minimum of geometry, interaction, decay for next step
       auto const min_distance = std::min(
           {distance_interact, distance_decay, distance_max, geomMaxLength});
 

Framework/Utilities/quartic.cpp

+/***************************************************************************
+ *   Copyright (C) 2016 by Саша Миленковић                                 *
+ *   sasa.milenkovic.xyz@gmail.com                                         *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *   This program is distributed in the hope that it will be useful,       *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
+ *   GNU General Public License for more details.                          *
+ *   ( http://www.gnu.org/licenses/gpl-3.0.en.html )                       *
+ *									   *
+ *   You should have received a copy of the GNU General Public License     *
+ *   along with this program; if not, write to the                         *
+ *   Free Software Foundation, Inc.,                                       *
+ *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
+ ***************************************************************************/
+
+#include <cmath>
+#include "quartic.h"
+
+//---------------------------------------------------------------------------
+// solve cubic equation x^3 + a*x^2 + b*x + c
+// x - array of size 3
+// In case 3 real roots: => x[0], x[1], x[2], return 3
+//         2 real roots: x[0], x[1],          return 2
+//         1 real root : x[0], x[1] ± i*x[2], return 1
+unsigned int solveP3(double *x,double a,double b,double c) {
+	double a2 = a*a;
+    	double q  = (a2 - 3*b)/9;
+	double r  = (a*(2*a2-9*b) + 27*c)/54;
+    	double r2 = r*r;
+	double q3 = q*q*q;
+	double A,B;
+    	if(r2<q3) 
+    	{
+    		double t=r/sqrt(q3);
+    		if( t<-1) t=-1;
+    		if( t> 1) t= 1;
+    		t=acos(t);
+    		a/=3; q=-2*sqrt(q);
+    		x[0]=q*cos(t/3)-a;
+    		x[1]=q*cos((t+M_2PI)/3)-a;
+    		x[2]=q*cos((t-M_2PI)/3)-a;
+    		return 3;
+    	} 
+    	else 
+    	{
+    		A =-pow(fabs(r)+sqrt(r2-q3),1./3);
+    		if( r<0 ) A=-A;
+    		B = (0==A ? 0 : q/A);
+          
+		a/=3;
+		x[0] =(A+B)-a;
+		x[1] =-0.5*(A+B)-a;
+		x[2] = 0.5*sqrt(3.)*(A-B);
+		if(fabs(x[2])<eps) { x[2]=x[1]; return 2; }
+		
+		return 1;
+        }
+}
+
+//---------------------------------------------------------------------------
+// solve quartic equation x^4 + a*x^3 + b*x^2 + c*x + d
+// Attention - this function returns dynamically allocated array. It has to be released afterwards.
+DComplex* solve_quartic(double a, double b, double c, double d)
+{
+	double a3 = -b;
+	double b3 =  a*c -4.*d;
+	double c3 = -a*a*d - c*c + 4.*b*d;
+
+	// cubic resolvent
+	// y^3 − b*y^2 + (ac−4d)*y − a^2*d−c^2+4*b*d = 0
+	
+	double x3[3];
+	unsigned int iZeroes = solveP3(x3, a3, b3, c3);
+	
+	double q1, q2, p1, p2, D, sqD, y;
+
+	y = x3[0];
+	// The essence - choosing Y with maximal absolute value.
+	if(iZeroes != 1)
+	{
+		if(fabs(x3[1]) > fabs(y)) y = x3[1];
+		if(fabs(x3[2]) > fabs(y)) y = x3[2];
+	}
+
+	// h1+h2 = y && h1*h2 = d  <=>  h^2 -y*h + d = 0    (h === q)
+
+	D = y*y - 4*d;
+	if(fabs(D) < eps) //in other words - D==0
+	{
+		q1 = q2 = y * 0.5;
+		// g1+g2 = a && g1+g2 = b-y   <=>   g^2 - a*g + b-y = 0    (p === g)
+		D = a*a - 4*(b-y);
+		if(fabs(D) < eps) //in other words - D==0
+			p1 = p2 = a * 0.5;
+
+		else
+		{
+			sqD = sqrt(D);
+			p1 = (a + sqD) * 0.5;
+			p2 = (a - sqD) * 0.5;
+		}
+	}
+	else
+	{
+		sqD = sqrt(D);
+		q1 = (y + sqD) * 0.5;
+		q2 = (y - sqD) * 0.5;
+		// g1+g2 = a && g1*h2 + g2*h1 = c       ( && g === p )  Krammer
+		p1 = (a*q1-c)/(q1-q2);
+		p2 = (c-a*q2)/(q1-q2);
+	}
+
+	DComplex* retval = new DComplex[4];
+
+	// solving quadratic eq. - x^2 + p1*x + q1 = 0
+	D = p1*p1 - 4*q1;
+	if(D < 0.0)
+	{
+		retval[0].real( -p1 * 0.5 );
+		retval[0].imag( sqrt(-D) * 0.5 );
+		retval[1] = std::conj(retval[0]);
+	}
+	else
+	{
+		sqD = sqrt(D);
+		retval[0].real( (-p1 + sqD) * 0.5 );
+		retval[1].real( (-p1 - sqD) * 0.5 );
+	}
+
+	// solving quadratic eq. - x^2 + p2*x + q2 = 0
+	D = p2*p2 - 4*q2;
+	if(D < 0.0)
+	{
+		retval[2].real( -p2 * 0.5 );
+		retval[2].imag( sqrt(-D) * 0.5 );
+		retval[3] = std::conj(retval[2]);
+	}
+	else
+	{
+		sqD = sqrt(D);
+		retval[2].real( (-p2 + sqD) * 0.5 );
+		retval[3].real( (-p2 - sqD) * 0.5 );
+	}
+
+    return retval;
+}
+

Framework/Utilities/quartic.h

+/***************************************************************************
+ *   Copyright (C) 2016 by Саша Миленковић                                 *
+ *   sasa.milenkovic.xyz@gmail.com                                         *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *   This program is distributed in the hope that it will be useful,       *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
+ *   GNU General Public License for more details.                          *
+ *   ( http://www.gnu.org/licenses/gpl-3.0.en.html )                       *
+ *									   *
+ *   You should have received a copy of the GNU General Public License     *
+ *   along with this program; if not, write to the                         *
+ *   Free Software Foundation, Inc.,                                       *
+ *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
+ ***************************************************************************/
+ 
+#ifndef QUARTIC_H_INCLUDED
+#define QUARTIC_H_INCLUDED
+
+#include <complex>
+
+const double PI = 3.141592653589793238463L;
+const double M_2PI = 2*PI;
+const double eps=1e-12;
+
+typedef std::complex<double> DComplex;
+
+//---------------------------------------------------------------------------
+// useful for testing
+ inline DComplex polinom_2(DComplex x, double a, double b)
+ {
+	 //Horner's scheme for x*x + a*x + b
+	 return x * (x + a) + b;
+ }
+
+//---------------------------------------------------------------------------
+// useful for testing
+ inline DComplex polinom_3(DComplex x, double a, double b, double c)
+ {
+	 //Horner's scheme for x*x*x + a*x*x + b*x + c;
+	 return x * (x * (x + a) + b) + c;
+ }
+
+//---------------------------------------------------------------------------
+// useful for testing
+ inline DComplex polinom_4(DComplex x, double a, double b, double c, double d)
+ {
+	 //Horner's scheme for x*x*x*x + a*x*x*x + b*x*x + c*x + d;
+	 return x * (x * (x * (x + a) + b) + c) + d;
+ }
+
+//---------------------------------------------------------------------------
+// x - array of size 3
+// In case 3 real roots: => x[0], x[1], x[2], return 3
+//         2 real roots: x[0], x[1],          return 2
+//         1 real root : x[0], x[1] ± i*x[2], return 1
+unsigned int solveP3(double* x, double a, double b, double c);
+
+//---------------------------------------------------------------------------
+// solve quartic equation x^4 + a*x^3 + b*x^2 + c*x + d
+// Attention - this function returns dynamically allocated array. It has to be released afterwards.
+DComplex* solve_quartic(double a, double b, double c, double d);
+
+
+#endif // QUARTIC_H_INCLUDED

Processes/TrackingLine/TrackingLine_interpolation.h

+/*
+ * (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_processes_TrackingLine_h_
+#define _include_corsika_processes_TrackingLine_h_
+
+#include <corsika/geometry/Line.h>
+#include <corsika/geometry/Plane.h>
+#include <corsika/geometry/Sphere.h>
+#include <corsika/geometry/Trajectory.h>
+#include <corsika/geometry/Vector.h>
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/units/PhysicalUnits.h>
+#include <corsika/utl/quartic.h>
+#include <cmath>
+#include <optional>
+#include <type_traits>
+#include <utility>
+
+namespace corsika::environment {
+  template <typename IEnvironmentModel>
+  class Environment;
+  template <typename IEnvironmentModel>
+  class VolumeTreeNode;
+} // namespace corsika::environment
+
+namespace corsika::process {
+
+  namespace tracking_line {
+
+    std::optional<std::pair<corsika::units::si::TimeType, corsika::units::si::TimeType>>
+    TimeOfIntersection(geometry::Line const&, geometry::Sphere const&);
+
+    corsika::units::si::TimeType TimeOfIntersection(geometry::Line const&,
+                                                    geometry::Plane const&);
+
+    class TrackingLine {
+
+    public:
+      TrackingLine(){};
+
+      template <typename Particle> // was Stack previously, and argument was
+                                   // Stack::StackIterator
+      auto GetTrack(Particle const& p) {
+        using namespace corsika::units::si;
+        using namespace corsika::geometry;
+        geometry::Vector<SpeedType::dimension_type> velocity =
+            p.GetMomentum() / p.GetEnergy() * corsika::units::constants::c;
+            
+            std::complex<double>*  solutions = solve_quartic(1, 0, 1, -20);
+            std::vector<double> tmp;
+            for(int i = 0; i < 4; i++) {
+              if(solutions[i].imag() == 0 && solutions[i].real() >= 0) {
+                tmp.push_back(solutions[i].real());
+                }
+            }
+            double s = *std::min_element(tmp.begin(),tmp.end());
+            std::cout << "s = " << s << std::endl;
+            std::cout << "x1 = " << (solutions[0].real()>=0. ? " " : "") << solutions[0].real(); if(solutions[0].imag()!=0.0) std::cout << "   +   i * " <<  solutions[0].imag(); std::cout << std::endl;
+		std::cout << "x2 = " << (solutions[1].real()>=0. ? " " : "") << solutions[1].real(); if(solutions[1].imag()!=0.0) std::cout << "   -   i * " << -solutions[1].imag(); std::cout << std::endl;
+		std::cout << "x3 = " << (solutions[2].real()>=0. ? " " : "") << solutions[2].real(); if(solutions[2].imag()!=0.0) std::cout << "   +   i * " <<  solutions[2].imag(); std::cout << std::endl;
+		std::cout << "x4 = " << (solutions[3].real()>=0. ? " " : "") << solutions[3].real(); if(solutions[3].imag()!=0.0) std::cout << "   -   i * " << -solutions[3].imag(); std::cout << std::endl;
+
+        auto const currentPosition = p.GetPosition();
+        std::cout << "TrackingLine pid: " << p.GetPID()
+                  << " , E = " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
+        std::cout << "TrackingLine pos: "
+                  << currentPosition.GetCoordinates()
+                  // << " [" << p.GetNode()->GetModelProperties().GetName() << "]"
+                  << std::endl;
+        std::cout << "TrackingLine   E: " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
+        
+        // determine velocity after adding magnetic field
+        auto const* currentLogicalVolumeNode = p.GetNode();
+        int chargeNumber;
+        if(corsika::particles::IsNucleus(p.GetPID())) {
+        	chargeNumber = p.GetNuclearZ();
+        } else {
+        	chargeNumber = corsika::particles::GetChargeNumber(p.GetPID());
+        }
+        geometry::Vector<dimensionless_d> const directionBefore = velocity.normalized();
+        auto magMaxLength = 1_m/0;
+        auto directionAfter = directionBefore;
+        if(chargeNumber != 0) {
+        	auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(currentPosition);
+        	std::cout << "TrackingLine   B: " << magneticfield.GetComponents() / 1_uT << " uT " << std::endl;
+        	geometry::Vector<SpeedType::dimension_type> const velocityVerticalMag = velocity -
+        		velocity.parallelProjectionOnto(magneticfield);
+		LengthType const gyroradius = p.GetEnergy() * velocityVerticalMag.GetNorm() * 1_V / 
+					      (corsika::units::constants::cSquared * abs(chargeNumber) * 
+					      magneticfield.GetNorm() * 1_eV);
+ 		//steplength should consider more things than just gyroradius
+    double maxAngle = 0.1; 
+		LengthType const Steplength = 2 * sin(maxAngle * M_PI / 180) * gyroradius;
+		// First Movement
+		auto position = currentPosition + directionBefore * Steplength / 2;
+		// Change of direction by magnetic field at position
+		magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(position);
+		directionAfter = directionBefore + directionBefore.cross(magneticfield) * chargeNumber * 
+				 Steplength * corsika::units::constants::cSquared * 1_eV / 
+				 (p.GetEnergy() * velocity.GetNorm() * 1_V); 
+		// Second Movement
+		position = position + directionAfter * Steplength / 2;
+		magMaxLength = (position - currentPosition).GetNorm();
+		geometry::Vector<dimensionless_d> const direction = (position - currentPosition) / 
+									magMaxLength;
+		velocity = direction * velocity.GetNorm();
+		std::cout << "TrackingLine   p: " << (direction * p.GetMomentum().GetNorm()).GetComponents() / 1_GeV
+                  << " GeV " << std::endl;
+
+        } else {
+        	std::cout << "TrackingLine   p: " << p.GetMomentum().GetComponents() / 1_GeV
+                  << " GeV " << std::endl;
+        }
+        
+        std::cout << "TrackingLine   v: " << velocity.GetComponents() << std::endl;
+        
+        geometry::Line line(currentPosition, velocity);
+
+        //auto const* currentLogicalVolumeNode = p.GetNode();
+        //~ auto const* currentNumericalVolumeNode =
+        //~ fEnvironment.GetUniverse()->GetContainingNode(currentPosition);
+        auto const numericallyInside =
+            currentLogicalVolumeNode->GetVolume().Contains(currentPosition);
+
+        std::cout << "numericallyInside = " << (numericallyInside ? "true" : "false");
+
+        auto const& children = currentLogicalVolumeNode->GetChildNodes();
+        auto const& excluded = currentLogicalVolumeNode->GetExcludedNodes();
+
+        std::vector<std::pair<TimeType, decltype(p.GetNode())>> intersections;
+
+        // for entering from outside
+        auto addIfIntersects = [&](auto const& vtn) {
+          auto const& volume = vtn.GetVolume();
+          auto const& sphere = dynamic_cast<geometry::Sphere const&>(
+              volume); // for the moment we are a bit bold here and assume
+          // everything is a sphere, crashes with exception if not
+
+          if (auto opt = TimeOfIntersection(line, sphere); opt.has_value()) {
+            auto const [t1, t2] = *opt;
+            std::cout << "intersection times: " << t1 / 1_s << "; "
+                      << t2 / 1_s
+                      // << " " << vtn.GetModelProperties().GetName()
+                      << std::endl;
+            if (t1.magnitude() > 0)
+              intersections.emplace_back(t1, &vtn);
+            else if (t2.magnitude() > 0)
+              std::cout << "inside other volume" << std::endl;
+          }
+        };
+
+        for (auto const& child : children) { addIfIntersects(*child); }
+        for (auto const* ex : excluded) { addIfIntersects(*ex); }
+
+        {
+          auto const& sphere = dynamic_cast<geometry::Sphere const&>(
+              currentLogicalVolumeNode->GetVolume());
+          // for the moment we are a bit bold here and assume
+          // everything is a sphere, crashes with exception if not
+          [[maybe_unused]] auto const [t1, t2] = *TimeOfIntersection(line, sphere);
+          [[maybe_unused]] auto dummy_t1 = t1;
+          intersections.emplace_back(t2, currentLogicalVolumeNode->GetParent());
+        }
+
+        auto const minIter = std::min_element(
+            intersections.cbegin(), intersections.cend(),
+            [](auto const& a, auto const& b) { return a.first < b.first; });
+
+        TimeType min;
+
+        if (minIter == intersections.cend()) {
+          min = 1_s; // todo: do sth. more reasonable as soon as tracking is able
+          // to handle the numerics properly
+          throw std::runtime_error("no intersection with anything!");
+        } else {
+          min = minIter->first;
+        }
+
+        std::cout << " t-intersect: "
+                  << min
+                  // << " " << minIter->second->GetModelProperties().GetName()
+                  << std::endl;
+
+        return std::make_tuple(geometry::Trajectory<geometry::Line>(line, min),
+                               velocity.norm() * min, minIter->second, magMaxLength, 
+                               directionBefore, directionAfter);
+      }
+    };
+
+  } // namespace tracking_line
+
+} // namespace corsika::process
+
+#endif