IAP GITLAB

Commit 44d53229 authored by Ralf Ulrich's avatar Ralf Ulrich

first part of magnetic field changes

parent ff9688fc
......@@ -115,7 +115,7 @@ endif (NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
#+++++++++++++++++++++++++++++
# Setup external dependencies.
#
include (conan.cmake)
include (conan) # self-provided in 'cmake' directory
#
# download and build all dependencies
conan_cmake_run (CONANFILE conanfile.txt
......
This diff is collapsed.
......@@ -22,16 +22,19 @@ namespace corsika {
(vPerp_ * (std::cos(omegaC_ * t) - 1) + uPerp_ * std::sin(omegaC_ * t)) /
omegaC_;
}
inline VelocityVector Helix::getVelocity(TimeType const t) const {
return vPar_ + (vPerp_ * (cos(omegaC_ * t) - 1) + uPerp_ * sin(omegaC_ * t));
}
Point Helix::getPositionFromArclength(LengthType const l) const {
inline Point Helix::getPositionFromArclength(LengthType const l) const {
return getPosition(getTimeFromArclength(l));
}
LengthType Helix::getArcLength(TimeType const t1, TimeType const t2) const {
inline LengthType Helix::getArcLength(TimeType const t1, TimeType const t2) const {
return (vPar_ + vPerp_).getNorm() * (t2 - t1);
}
TimeType Helix::getTimeFromArclength(LengthType const l) const {
inline TimeType Helix::getTimeFromArclength(LengthType const l) const {
return l / (vPar_ + vPerp_).getNorm();
}
......
......@@ -14,22 +14,24 @@
namespace corsika {
Point Line::getPosition(TimeType const t) const { return start_point_ + velocity_ * t; }
inline Point Line::getPosition(TimeType const t) const { return start_point_ + velocity_ * t; }
Point Line::getPositionFromArclength(LengthType const l) const {
inline VelocityVector const& Line::getVelocity(TimeType const) const { return velocity_; }
inline Point Line::getPositionFromArclength(LengthType const l) const {
return start_point_ + velocity_.normalized() * l;
}
LengthType Line::getArcLength(TimeType const t1, TimeType const t2) const {
inline LengthType Line::getArcLength(TimeType const t1, TimeType const t2) const {
return velocity_.getNorm() * (t2 - t1);
}
TimeType Line::getTimeFromArclength(LengthType const t) const {
inline TimeType Line::getTimeFromArclength(LengthType const t) const {
return t / velocity_.getNorm();
}
Point const& Line::getStartPoint() const { return start_point_; }
inline Point const& Line::getStartPoint() const { return start_point_; }
Line::VelocityVec const& Line::getVelocity() const { return velocity_; }
inline VelocityVector const& Line::getVelocity() const { return velocity_; }
} // namespace corsika
......@@ -13,12 +13,16 @@
namespace corsika {
bool Sphere::isInside(Point const& p) const {
inline bool Sphere::isInside(Point const& p) const {
return radius_ * radius_ > (center_ - p).getSquaredNorm();
}
Point const& Sphere::getCenter() const { return center_; }
inline Point const& Sphere::getCenter() const { return center_; }
LengthType Sphere::getRadius() const { return radius_; }
inline void Sphere::setCenter(Point const& p) { center_ = p; }
inline LengthType Sphere::getRadius() const { return radius_; }
inline void Sphere::setRadius(LengthType const r) { radius_ = r; }
} // namespace corsika
......@@ -8,43 +8,149 @@
#pragma once
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/Line.hpp>
#include <corsika/framework/geometry/Point.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/PhysicalGeometry.hpp>
namespace corsika {
template <typename TType>
Point Trajectory<TType>::getPosition(double const u) const {
return TType::getPosition(timeLength_ * u);
inline VelocityVector LineTrajectory::getVelocity(double const u) const {
return initialVelocity_ * (1 - u) + finalVelocity_ * u;
}
inline TimeType LineTrajectory::getDuration(double const u) const {
return u * timeStep_;
}
inline LengthType LineTrajectory::getLength(double const u) const {
if (timeLength_ == 0_s) return 0_m;
if (timeStep_ == std::numeric_limits<TimeType::value_type>::infinity() * 1_s)
return std::numeric_limits<LengthType::value_type>::infinity() * 1_m;
return getDistance(u) * timeStep_ / timeLength_;
}
///! set new duration along potentially bend trajectory.
inline void LineTrajectory::setLength(LengthType const limit) {
setDuration(line_.getTimeFromArclength(limit));
}
///! set new duration along potentially bend trajectory.
// Scale other properties by "limit/timeLength_"
inline void LineTrajectory::setDuration(TimeType const limit) {
if (timeStep_ == 0_s) {
timeLength_ = 0_s;
setFinalVelocity(getVelocity(0));
timeStep_ = limit;
} else {
// for infinite steps there can't be a difference between
// curved and straight trajectory, this is fundamentally
// undefined: assume they are the same (which, i.e. is always correct for a
// straight line trajectory).
//
// Final note: only straight-line trajectories should have
// infinite steps! Everything else is ill-defined.
if (timeStep_ == std::numeric_limits<TimeType::value_type>::infinity() * 1_s ||
timeLength_ == std::numeric_limits<TimeType::value_type>::infinity() * 1_s) {
timeLength_ = limit;
timeStep_ = limit;
// ...and don't touch velocity
} else {
const double scale = limit / timeStep_;
timeLength_ *= scale;
setFinalVelocity(getVelocity(scale));
timeStep_ = limit;
}
}
}
inline LengthType LineTrajectory::getDistance(double const u) const {
assert(u <= 1);
assert(u >= 0);
return line_.getArcLength(0 * second, u * timeLength_);
}
template <typename TType>
TimeType Trajectory<TType>::getDuration() const {
return timeLength_;
inline Line const LeapFrogTrajectory::getLine() const {
auto D = getPosition(1) - getPosition(0);
auto d = D.getNorm();
auto v = initialVelocity_;
if (d > 1_um) { // if trajectory is ultra-short, we do not
// re-calculate velocity, just use initial
// value. Otherwise, this is numerically unstable
v = D / d * getVelocity(0).getNorm();
}
return Line(getPosition(0), v);
}
template <typename TType>
LengthType Trajectory<TType>::getLength() const {
return getDistance(timeLength_);
inline Point LeapFrogTrajectory::getPosition(double const u) const {
Point position = initialPosition_ + initialVelocity_ * timeStep_ * u / 2;
VelocityVector velocity =
initialVelocity_ + initialVelocity_.cross(magneticfield_) * timeStep_ * u * k_;
return position + velocity * timeStep_ * u / 2;
}
template <typename TType>
LengthType Trajectory<TType>::getDistance(TimeType const t) const {
assert(t <= timeLength_);
assert(t >= 0 * second);
return TType::getArcLength(0 * second, t);
inline VelocityVector LeapFrogTrajectory::getVelocity(double const u) const {
return initialVelocity_ + initialVelocity_.cross(magneticfield_) * timeStep_ * u * k_;
}
template <typename TType>
void Trajectory<TType>::getLimitEndTo(LengthType const limit) {
timeLength_ = TType::getTimeFromArclength(limit);
inline DirectionVector LeapFrogTrajectory::getDirection(double const u) const {
return getVelocity(u).normalized();
}
template <typename TType>
auto Trajectory<TType>::getNormalizedDirection() const {
static_assert(std::is_same_v<TType, corsika::Line>);
return TType::getVelocity().normalized();
///! duration along potentially bend trajectory
inline TimeType LeapFrogTrajectory::getDuration(double const u) const {
return u * timeStep_ *
(double(getVelocity(u).getNorm() / initialVelocity_.getNorm()) + 1.0) / 2;
}
///! total length along potentially bend trajectory
inline LengthType LeapFrogTrajectory::getLength(double const u) const {
return timeStep_ * initialVelocity_.getNorm() * u;
}
///! set new duration along potentially bend trajectory.
inline void LeapFrogTrajectory::setLength(LengthType const limit) {
if (initialVelocity_.getNorm() == 0_m / 1_s) setDuration(0_s);
setDuration(limit / initialVelocity_.getNorm());
}
///! set new duration along potentially bend trajectory.
// Scale other properties by "limit/timeLength_"
inline void LeapFrogTrajectory::setDuration(TimeType const limit) { timeStep_ = limit; }
/*
template <typename TType>
Point Trajectory<TType>::getPosition(double const u) const {
return TType::getPosition(timeLength_ * u);
}
template <typename TType>
TimeType Trajectory<TType>::getDuration() const {
return timeLength_;
}
template <typename TType>
LengthType Trajectory<TType>::getLength() const {
return getDistance(timeLength_);
}
template <typename TType>
LengthType Trajectory<TType>::getDistance(TimeType const t) const {
assert(t <= timeLength_);
assert(t >= 0 * second);
return TType::getArcLength(0 * second, t);
}
template <typename TType>
void Trajectory<TType>::getLimitEndTo(LengthType const limit) {
timeLength_ = TType::getTimeFromArclength(limit);
}
template <typename TType>
auto Trajectory<TType>::getNormalizedDirection() const {
static_assert(std::is_same_v<TType, corsika::Line>);
return TType::getVelocity().normalized();
}
*/
} // namespace corsika
......@@ -159,7 +159,7 @@ namespace corsika {
// if this is not a ContinuousProcess --> evaluate probability
lambda_inv_sum += A_.getInverseInteractionLength(view.parent());
// check if we should execute THIS process and then EXIT
if (lambda_inv_select < lambda_inv_sum) {
if (lambda_inv_select <= lambda_inv_sum) {
A_.doInteraction(view);
return ProcessReturn::Interacted;
}
......@@ -174,7 +174,7 @@ namespace corsika {
lambda_inv_sum += B_.getInverseInteractionLength(view.parent());
// soon as SecondaryView::parent() is migrated!
// check if we should execute THIS process and then EXIT
if (lambda_inv_select < lambda_inv_sum) {
if (lambda_inv_select <= lambda_inv_sum) {
B_.doInteraction(view);
return ProcessReturn::Interacted;
}
......@@ -218,7 +218,7 @@ namespace corsika {
// if this is not a ContinuousProcess --> evaluate probability
decay_inv_sum += A_.getInverseLifetime(view.parent());
// check if we should execute THIS process and then EXIT
if (decay_inv_select < decay_inv_sum) { // more pedagogical: rndm_select <
if (decay_inv_select <= decay_inv_sum) { // more pedagogical: rndm_select <
// decay_inv_sum / decay_inv_tot
A_.doDecay(view);
return ProcessReturn::Decayed;
......@@ -232,7 +232,7 @@ namespace corsika {
// if this is not a ContinuousProcess --> evaluate probability
decay_inv_sum += B_.getInverseLifetime(view.parent());
// check if we should execute THIS process and then EXIT
if (decay_inv_select < decay_inv_sum) {
if (decay_inv_select <= decay_inv_sum) {
B_.doDecay(view);
return ProcessReturn::Decayed;
}
......
/***************************************************************************
* 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. *
***************************************************************************/
#pragma once
#include <cmath>
namespace corsika::quartic_solver {
//---------------------------------------------------------------------------
// solve cubic equation x^3 + a*x^2 + b*x + c = 0
// 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;
}
} // namespace corsika::quartic_solver
\ No newline at end of file
......@@ -12,7 +12,6 @@
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/Line.hpp>
#include <corsika/framework/geometry/Point.hpp>
#include <corsika/framework/geometry/Trajectory.hpp>
namespace corsika {
......@@ -23,8 +22,8 @@ namespace corsika {
template <typename TDerived>
GrammageType BaseExponential<TDerived>::getIntegratedGrammage(
Trajectory<Line> const& line, LengthType vL,
Vector<dimensionless_d> const& axis) const {
setup::Trajectory const& line, LengthType vL,
DirectionVector const& axis) const {
if (vL == LengthType::zero()) { return GrammageType::zero(); }
auto const uDotA = line.getNormalizedDirection().dot(axis).magnitude();
......@@ -39,8 +38,8 @@ namespace corsika {
template <typename TDerived>
LengthType BaseExponential<TDerived>::getArclengthFromGrammage(
Trajectory<Line> const& line, GrammageType grammage,
Vector<dimensionless_d> const& axis) const {
setup::Trajectory const& line, GrammageType grammage,
DirectionVector const& axis) const {
auto const uDotA = line.getNormalizedDirection().dot(axis).magnitude();
auto const rhoStart = getImplementation().getMassDensity(line.getStartPoint());
......
......@@ -11,7 +11,6 @@
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/Line.hpp>
#include <corsika/framework/geometry/Point.hpp>
#include <corsika/framework/geometry/Trajectory.hpp>
#include <corsika/media/BaseExponential.hpp>
#include <corsika/media/NuclearComposition.hpp>
......@@ -39,13 +38,13 @@ namespace corsika {
}
template <typename T>
GrammageType FlatExponential<T>::getIntegratedGrammage(Trajectory<Line> const& line,
GrammageType FlatExponential<T>::getIntegratedGrammage(setup::Trajectory const& line,
LengthType to) const {
return BaseExponential<FlatExponential<T>>::getIntegratedGrammage(line, to, axis_);
}
template <typename T>
LengthType FlatExponential<T>::getArclengthFromGrammage(Trajectory<Line> const& line,
LengthType FlatExponential<T>::getArclengthFromGrammage(setup::Trajectory const& line,
GrammageType grammage) const {
return BaseExponential<FlatExponential<T>>::getArclengthFromGrammage(line, grammage,
axis_);
......
......@@ -11,7 +11,6 @@
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/Line.hpp>
#include <corsika/framework/geometry/Point.hpp>
#include <corsika/framework/geometry/Trajectory.hpp>
#include <corsika/media/NuclearComposition.hpp>
namespace corsika {
......@@ -32,13 +31,13 @@ namespace corsika {
}
template <typename T>
GrammageType HomogeneousMedium<T>::getIntegratedGrammage(Trajectory<Line> const&,
GrammageType HomogeneousMedium<T>::getIntegratedGrammage(setup::Trajectory const&,
LengthType to) const {
return to * density_;
}
template <typename T>
LengthType HomogeneousMedium<T>::getArclengthFromGrammage(Trajectory<Line> const&,
LengthType HomogeneousMedium<T>::getArclengthFromGrammage(setup::Trajectory const&,
GrammageType grammage) const {
return grammage / density_;
}
......
......@@ -11,7 +11,6 @@
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/Line.hpp>
#include <corsika/framework/geometry/Point.hpp>
#include <corsika/framework/geometry/Trajectory.hpp>
#include <corsika/media/NuclearComposition.hpp>
namespace corsika {
......@@ -37,13 +36,13 @@ namespace corsika {
template <typename T, typename TDensityFunction>
GrammageType InhomogeneousMedium<T, TDensityFunction>::getIntegratedGrammage(
Trajectory<Line> const& line, LengthType to) const {
setup::Trajectory const& line, LengthType to) const {
return densityFunction_.getIntegrateGrammage(line, to);
}
template <typename T, typename TDensityFunction>
LengthType InhomogeneousMedium<T, TDensityFunction>::getArclengthFromGrammage(
Trajectory<Line> const& line, GrammageType grammage) const {
setup::Trajectory const& line, GrammageType grammage) const {
return densityFunction_.getArclengthFromGrammage(line, grammage);
}
......
......@@ -19,28 +19,28 @@ namespace corsika {
template <typename TDerived>
auto LinearApproximationIntegrator<TDerived>::getIntegrateGrammage(
Trajectory<Line> const& line, LengthType length) const {
setup::Trajectory const& line, LengthType length) const {
auto const c0 = getImplementation().evaluateAt(line.getPosition(0));
auto const c1 = getImplementation().rho_.getFirstDerivative(
line.getPosition(0), line.getNormalizedDirection());
auto const c1 = getImplementation().rho_.getFirstDerivative(line.getPosition(0),
line.getDirection(0));
return (c0 + 0.5 * c1 * length) * length;
}
template <typename TDerived>
auto LinearApproximationIntegrator<TDerived>::getArclengthFromGrammage(
Trajectory<Line> const& line, GrammageType grammage) const {
setup::Trajectory const& line, GrammageType grammage) const {
auto const c0 = getImplementation().rho_(line.getPosition(0));
auto const c1 = getImplementation().rho_.getFirstDerivative(
line.getPosition(0), line.getNormalizedDirection());
auto const c1 = getImplementation().rho_.getFirstDerivative(line.getPosition(0),
line.getDirection(0));
return (1 - 0.5 * grammage * c1 / (c0 * c0)) * grammage / c0;
}
template <typename TDerived>
auto LinearApproximationIntegrator<TDerived>::getMaximumLength(
Trajectory<Line> const& line, [[maybe_unused]] double relError) const {
setup::Trajectory const& line, [[maybe_unused]] double relError) const {
[[maybe_unused]] auto const c1 = getImplementation().rho_.getSecondDerivative(
line.getPosition(0), line.getNormalizedDirection());
line.getPosition(0), line.getDirection(0));
// todo: provide a real, working implementation
return 1_m * std::numeric_limits<double>::infinity();
......
......@@ -57,32 +57,42 @@ namespace corsika {
}
GrammageType ShowerAxis::getX(LengthType l) const {
auto const fractionalBin = l / steplength_;
double const fractionalBin = l / steplength_;
int const lower = fractionalBin; // indices of nearest X support points
auto const lambda = fractionalBin - lower;
decltype(X_.size()) const upper = lower + 1;
double const fraction = fractionalBin - lower;
unsigned int const upper = lower + 1;
if (lower < 0) {
if (fractionalBin < 0) {
CORSIKA_LOG_ERROR("cannot extrapolate to points behind point of injection l={} m",
l / 1_m);
throw std::runtime_error("cannot extrapolate to points behind point of injection");
if (throw_) {
throw std::runtime_error(
"cannot extrapolate to points behind point of injection");
}
return minimumX();
}
if (upper >= X_.size()) {
const std::string err =
fmt::format("shower axis too short, cannot extrapolate (l / max_length_ = {} )",
l / max_length_);
CORSIKA_LOG_ERROR(err);
throw std::runtime_error(err.c_str());
CORSIKA_LOG_ERROR(
"shower axis too short, cannot extrapolate (l / max_length_ = {} )",
l / max_length_);
if (throw_) {
const std::string err = fmt::format(
"shower axis too short, cannot extrapolate (l / max_length_ = {} )",
l / max_length_);
throw std::runtime_error(err.c_str());
}
}
CORSIKA_LOG_TRACE("showerAxis::X frac={}, fractionalBin={}, lower={}, upper={}", fraction,
fractionalBin, lower, upper);
assert(0 <= lambda && lambda <= 1.);
assert(0 <= fraction && fraction <= 1.);
CORSIKA_LOG_TRACE("ShowerAxis::getX l={} m, lower={}, lambda={}, upper={}", l / 1_m,
lower, lambda, upper);
CORSIKA_LOG_TRACE("ShowerAxis::getX l={} m, lower={}, fraction={}, upper={}", l / 1_m,
lower, fraction, upper);
// linear interpolation between getX[lower] and X[upper]
return X_[upper] * lambda + X_[lower] * (1 - lambda);
return X_[upper] * fraction + X_[lower] * (1 - fraction);
}
LengthType ShowerAxis::getSteplength() const { return steplength_; }
......
......@@ -12,14 +12,17 @@
namespace corsika {