IAP GITLAB

Commit 0e397e07 authored by Andre Schmidt's avatar Andre Schmidt

Test

parent 7ddb4909
Pipeline #2007 failed with stages
This diff is collapsed.
/***************************************************************************
* 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;
}
/***************************************************************************
* 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
/*
* (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
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment