IAP GITLAB

Commit cfb97fc2 authored by Ralf Ulrich's avatar Ralf Ulrich

magnetic fields, refactorred

parent 44d53229
......@@ -118,8 +118,7 @@ namespace corsika {
next_interact);
// determine the maximum geometric step length
LengthType const continuous_max_dist =
process_sequence_.getMaxStepLength(vParticle, step);
LengthType const continuous_max_dist = sequence_.getMaxStepLength(vParticle, step);
// take minimum of geometry, interaction, decay for next step
auto const min_distance =
......@@ -139,7 +138,7 @@ namespace corsika {
step.setLength(min_distance);
vParticle.setPosition(step.getPosition(1));
// assumption: tracking does not change absolute momentum:
vParticle.setMomentum(step.getDirection(1) * vParticle.getMomentum().norm());
vParticle.setMomentum(step.getDirection(1) * vParticle.getMomentum().getNorm());
vParticle.setTime(vParticle.getTime() + step.getDuration());
// apply all continuous processes on particle + track
......@@ -175,7 +174,7 @@ namespace corsika {
[[maybe_unused]] auto projectile = secondaries.getProjectile();
if (distance_interat < distance_decay) {
if (distance_interact < distance_decay) {
interaction(secondaries);
} else {
decay(secondaries);
......@@ -231,63 +230,63 @@ namespace corsika {
sequence_.doBoundaryCrossing(vParticle, *currentLogicalNode, *nextVol);
}
}
}
template <typename TTracking, typename TProcessList, typename TStack,
typename TStackView>
ProcessReturn Cascade<TTracking, TProcessList, TStack, TStackView>::decay(TStackView &
view) {
CORSIKA_LOG_DEBUG("decay");
InverseTimeType const actual_decay_time =
sequence_.getInverseLifetime(view.parent());
template <typename TTracking, typename TProcessList, typename TStack,
typename TStackView>
ProcessReturn Cascade<TTracking, TProcessList, TStack, TStackView>::decay(
TStackView& view) {
CORSIKA_LOG_DEBUG("decay");
InverseTimeType const actual_decay_time = sequence_.getInverseLifetime(view.parent());
UniformRealDistribution<InverseTimeType> uniDist(actual_decay_time);
const auto sample_process = uniDist(rng_);
UniformRealDistribution<InverseTimeType> uniDist(actual_decay_time);
const auto sample_process = uniDist(rng_);
auto const returnCode = sequence_.selectDecay(view, sample_process);
if (returnCode != ProcessReturn::Decayed) {
CORSIKA_LOG_WARN("Particle did not decay!");
}
setEventType(view, history::EventType::Decay);
return returnCode;
auto const returnCode = sequence_.selectDecay(view, sample_process);
if (returnCode != ProcessReturn::Decayed) {
CORSIKA_LOG_WARN("Particle did not decay!");
}
setEventType(view, history::EventType::Decay);
return returnCode;
}
template <typename TTracking, typename TProcessList, typename TStack,
typename TStackView>
ProcessReturn Cascade<TTracking, TProcessList, TStack, TStackView>::interaction(
TStackView & view) {
CORSIKA_LOG_DEBUG("collide");
template <typename TTracking, typename TProcessList, typename TStack,
typename TStackView>
ProcessReturn Cascade<TTracking, TProcessList, TStack, TStackView>::interaction(
TStackView& view) {
CORSIKA_LOG_DEBUG("collide");
InverseGrammageType const current_inv_length =
sequence_.getInverseInteractionLength(view.parent());
InverseGrammageType const current_inv_length =
sequence_.getInverseInteractionLength(view.parent());
UniformRealDistribution<InverseGrammageType> uniDist(current_inv_length);
UniformRealDistribution<InverseGrammageType> uniDist(current_inv_length);
const auto sample_process = uniDist(rng_);
auto const returnCode = sequence_.selectInteraction(view, sample_process);
if (returnCode != ProcessReturn::Interacted) {
CORSIKA_LOG_WARN("Particle did not interace!");
}
setEventType(view, history::EventType::Interaction);
return returnCode;
const auto sample_process = uniDist(rng_);
auto const returnCode = sequence_.selectInteraction(view, sample_process);
if (returnCode != ProcessReturn::Interacted) {
CORSIKA_LOG_WARN("Particle did not interace!");
}
setEventType(view, history::EventType::Interaction);
return returnCode;
}
template <typename TTracking, typename TProcessList, typename TStack,
typename TStackView>
void Cascade<TTracking, TProcessList, TStack, TStackView>::setNodes() {
std::for_each(stack_.begin(), stack_.end(), [&](auto& p) {
auto const* numericalNode =
environment_.getUniverse()->getContainingNode(p.getPosition());
p.setNode(numericalNode);
});
}
template <typename TTracking, typename TProcessList, typename TStack,
typename TStackView>
void Cascade<TTracking, TProcessList, TStack, TStackView>::setNodes() {
std::for_each(stack_.begin(), stack_.end(), [&](auto& p) {
auto const* numericalNode =
environment_.getUniverse()->getContainingNode(p.getPosition());
p.setNode(numericalNode);
});
}
template <typename TTracking, typename TProcessList, typename TStack,
typename TStackView>
void Cascade<TTracking, TProcessList, TStack, TStackView>::setEventType(
TStackView & view, [[maybe_unused]] history::EventType eventType) {
if constexpr (TStackView::has_event) {
for (auto&& sec : view) { sec.getEvent()->setEventType(eventType); }
}
template <typename TTracking, typename TProcessList, typename TStack,
typename TStackView>
void Cascade<TTracking, TProcessList, TStack, TStackView>::setEventType(
TStackView& view, [[maybe_unused]] history::EventType eventType) {
if constexpr (TStackView::has_event) {
for (auto&& sec : view) { sec.getEvent()->setEventType(eventType); }
}
}
} // namespace corsika
} // namespace corsika
/*
* (c) Copyright 2020 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
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
......@@ -22,19 +22,16 @@ 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));
}
inline Point Helix::getPositionFromArclength(LengthType const l) const {
Point Helix::getPositionFromArclength(LengthType const l) const {
return getPosition(getTimeFromArclength(l));
}
inline LengthType Helix::getArcLength(TimeType const t1, TimeType const t2) const {
LengthType Helix::getArcLength(TimeType const t1, TimeType const t2) const {
return (vPar_ + vPerp_).getNorm() * (t2 - t1);
}
inline TimeType Helix::getTimeFromArclength(LengthType const l) const {
TimeType Helix::getTimeFromArclength(LengthType const l) const {
return l / (vPar_ + vPerp_).getNorm();
}
......
......@@ -14,9 +14,13 @@
namespace corsika {
inline 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;
}
inline VelocityVector const& Line::getVelocity(TimeType const) const { return velocity_; }
inline VelocityVector const& Line::getVelocity(TimeType const) const {
return velocity_;
}
inline Point Line::getPositionFromArclength(LengthType const l) const {
return start_point_ + velocity_.normalized() * l;
......@@ -32,6 +36,8 @@ namespace corsika {
inline Point const& Line::getStartPoint() const { return start_point_; }
inline DirectionVector Line::getDirection() const { return velocity_.normalized(); }
inline VelocityVector const& Line::getVelocity() const { return velocity_; }
} // namespace corsika
......@@ -219,7 +219,7 @@ namespace corsika {
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 <
// decay_inv_sum / decay_inv_tot
// decay_inv_sum / decay_inv_tot
A_.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. *
***************************************************************************/
/*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
*
* 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.
*/
#pragma once
......@@ -30,7 +18,7 @@ namespace corsika::quartic_solver {
// 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) {
inline 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;
......@@ -70,7 +58,7 @@ namespace corsika::quartic_solver {
// 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) {
inline 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;
......
......@@ -10,7 +10,6 @@
#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/Line.hpp>
#include <corsika/framework/geometry/Point.hpp>
namespace corsika {
......@@ -22,12 +21,11 @@ namespace corsika {
template <typename TDerived>
GrammageType BaseExponential<TDerived>::getIntegratedGrammage(
setup::Trajectory const& line, LengthType vL,
DirectionVector const& axis) const {
setup::Trajectory const& traj, LengthType vL, DirectionVector const& axis) const {
if (vL == LengthType::zero()) { return GrammageType::zero(); }
auto const uDotA = line.getNormalizedDirection().dot(axis).magnitude();
auto const rhoStart = getImplementation().getMassDensity(line.getStartPoint());
auto const uDotA = traj.getDirection(0).dot(axis).magnitude();
auto const rhoStart = getImplementation().getMassDensity(traj.getPosition(0));
if (uDotA == 0) {
return vL * rhoStart;
......@@ -38,10 +36,10 @@ namespace corsika {
template <typename TDerived>
LengthType BaseExponential<TDerived>::getArclengthFromGrammage(
setup::Trajectory const& line, GrammageType grammage,
setup::Trajectory const& traj, GrammageType grammage,
DirectionVector const& axis) const {
auto const uDotA = line.getNormalizedDirection().dot(axis).magnitude();
auto const rhoStart = getImplementation().getMassDensity(line.getStartPoint());
auto const uDotA = traj.getDirection(0).dot(axis).magnitude();
auto const rhoStart = getImplementation().getMassDensity(traj.getPosition(0));
if (uDotA == 0) {
return grammage / rhoStart;
......
......@@ -27,9 +27,10 @@ namespace corsika {
template <typename T>
MassDensityType FlatExponential<T>::getMassDensity(Point const& point) const {
return BaseExponential<FlatExponential<T>>::rho0_ *
exp(BaseExponential<FlatExponential<T>>::invLambda_ *
(point - BaseExponential<FlatExponential<T>>::point_).dot(axis_));
return BaseExponential<FlatExponential<T>>::getRho0() *
exp(BaseExponential<FlatExponential<T>>::getInvLambda() *
(point - BaseExponential<FlatExponential<T>>::getAnchorPoint())
.dot(axis_));
}
template <typename T>
......
......@@ -15,9 +15,11 @@ namespace corsika {
template <typename TEnvModel>
ShowerAxis::ShowerAxis(Point const& pStart, Vector<length_d> const& length,
Environment<TEnvModel> const& env, int steps)
Environment<TEnvModel> const& env, bool const doThrow,
int const steps)
: pointStart_(pStart)
, length_(length)
, throw_(doThrow)
, max_length_(length_.getNorm())
, steplength_(max_length_ / steps)
, axis_normalized_(length / max_length_)
......@@ -69,7 +71,7 @@ namespace corsika {
throw std::runtime_error(
"cannot extrapolate to points behind point of injection");
}
return minimumX();
return getMinimumX();
}
if (upper >= X_.size()) {
......@@ -83,8 +85,8 @@ namespace corsika {
throw std::runtime_error(err.c_str());
}
}
CORSIKA_LOG_TRACE("showerAxis::X frac={}, fractionalBin={}, lower={}, upper={}", fraction,
fractionalBin, lower, upper);
CORSIKA_LOG_TRACE("showerAxis::X frac={}, fractionalBin={}, lower={}, upper={}",
fraction, fractionalBin, lower, upper);
assert(0 <= fraction && fraction <= 1.);
......
......@@ -23,10 +23,11 @@ namespace corsika {
template <typename T>
MassDensityType SlidingPlanarExponential<T>::getMassDensity(Point const& point) const {
auto const height =
(point - BaseExponential<SlidingPlanarExponential<T>>::point_).getNorm() -
(point - BaseExponential<SlidingPlanarExponential<T>>::getAnchorPoint())
.getNorm() -
referenceHeight_;
return BaseExponential<SlidingPlanarExponential<T>>::rho0_ *
exp(BaseExponential<SlidingPlanarExponential<T>>::invLambda_ * height);
return BaseExponential<SlidingPlanarExponential<T>>::getRho0() *
exp(BaseExponential<SlidingPlanarExponential<T>>::getInvLambda() * height);
}
template <typename T>
......@@ -36,22 +37,22 @@ namespace corsika {
template <typename T>
GrammageType SlidingPlanarExponential<T>::getIntegratedGrammage(
Trajectory<Line> const& line, LengthType l) const {
auto const axis =
(line.getStartPoint() - BaseExponential<SlidingPlanarExponential<T>>::point_)
.normalized();
return BaseExponential<SlidingPlanarExponential<T>>::getIntegratedGrammage(line, l,
setup::Trajectory const& traj, LengthType l) const {
auto const axis = (traj.getPosition(0) -
BaseExponential<SlidingPlanarExponential<T>>::getAnchorPoint())
.normalized();
return BaseExponential<SlidingPlanarExponential<T>>::getIntegratedGrammage(traj, l,
axis);
}
template <typename T>
LengthType SlidingPlanarExponential<T>::getArclengthFromGrammage(
Trajectory<Line> const& line, GrammageType const grammage) const {
auto const axis =
(line.getStartPoint() - BaseExponential<SlidingPlanarExponential<T>>::point_)
.normalized();
setup::Trajectory const& traj, GrammageType const grammage) const {
auto const axis = (traj.getPosition(0) -
BaseExponential<SlidingPlanarExponential<T>>::getAnchorPoint())
.normalized();
return BaseExponential<SlidingPlanarExponential<T>>::getArclengthFromGrammage(
line, grammage, axis);
traj, grammage, axis);
}
} // namespace corsika
......@@ -80,13 +80,14 @@ namespace corsika {
excludedNodes_.push_back(pNode.get());
}
template <typename IModelProperties>
template <class MediumType, typename... Args>
auto VolumeTreeNode<IModelProperties>::createMedium(Args&&... args) {
static_assert(std::is_base_of_v<IMediumModel, MediumType>,
"unusable type provided, needs to be derived from \"IMediumModel\"");
return std::make_shared<MediumType>(std::forward<Args>(args)...);
}
/*
template <typename IModelProperties>
template <class MediumType, typename... Args>
auto VolumeTreeNode<IModelProperties>::createMedium(Args&&... args) {
static_assert(std::is_base_of_v<IMediumModel, MediumType>,
"unusable type provided, needs to be derived from \"IMediumModel\"");
return std::make_shared<MediumType>(std::forward<Args>(args)...);
}
*/
} // namespace corsika
......@@ -29,12 +29,12 @@ namespace corsika {
corsika::setup::Stack::particle_type& particle,
corsika::setup::Trajectory& trajectory) {
TimeType const timeOfIntersection =
(plane_.getCenter() - trajectory.getStartPoint()).dot(plane_.getNormal()) /
trajectory.getVelocity().dot(plane_.getNormal());
(plane_.getCenter() - trajectory.getPosition(0)).dot(plane_.getNormal()) /
trajectory.getVelocity(0).dot(plane_.getNormal());
if (timeOfIntersection < TimeType::zero()) { return ProcessReturn::Ok; }
if (plane_.isAbove(trajectory.getStartPoint()) ==
if (plane_.isAbove(trajectory.getPosition(0)) ==
plane_.isAbove(trajectory.getPosition(1))) {
return ProcessReturn::Ok;
}
......@@ -71,10 +71,11 @@ namespace corsika {
auto const* currentLogicalVolumeNode = vParticle.getNode();
auto magneticfield = currentLogicalVolumeNode->getModelProperties().getMagneticField(
vParticle.getPosition());
auto direction = trajectory.getLine().getV0().normalized();
auto direction = trajectory.getVelocity(0).normalized();
if (chargeNumber != 0 &&
abs(plane_.getNormal().dot(trajectory.getLine().getV0().cross(magneticfield))) *
abs(plane_.getNormal().dot(
trajectory.getLine().getVelocity().cross(magneticfield))) *
1_s / 1_m / 1_T >
1e-6) {
auto const* currentLogicalVolumeNode = vParticle.getNode();
......@@ -82,10 +83,10 @@ namespace corsika {
currentLogicalVolumeNode->getModelProperties().getMagneticField(
vParticle.getPosition());
auto k =
chargeNumber * constants::c * 1_eV / (vParticle.getMomentum().norm() * 1_V);
chargeNumber * constants::c * 1_eV / (vParticle.getMomentum().getNorm() * 1_V);
if (direction.dot(plane_.getNormal()) * direction.dot(plane_.getNormal()) -
(plane_.getNormal().dot(trajectory.getLine().getR0() - plane_.getCenter()) *
(plane_.getNormal().dot(trajectory.getPosition(0) - plane_.getCenter()) *
plane_.getNormal().dot(direction.cross(magneticfield)) * 2 * k) <
0) {
return std::numeric_limits<double>::infinity() * 1_m;
......@@ -93,16 +94,14 @@ namespace corsika {
LengthType MaxStepLength1 =
(sqrt(direction.dot(plane_.getNormal()) * direction.dot(plane_.getNormal()) -
(plane_.getNormal().dot(trajectory.getLine().getR0() -
plane_.getCenter()) *
(plane_.getNormal().dot(trajectory.getPosition(0) - plane_.getCenter()) *
plane_.getNormal().dot(direction.cross(magneticfield)) * 2 * k)) -
direction.dot(plane_.getNormal()) / direction.getNorm()) /
(plane_.getNormal().dot(direction.cross(magneticfield)) * k);
LengthType MaxStepLength2 =
(-sqrt(direction.dot(plane_.getNormal()) * direction.dot(plane_.getNormal()) -
(plane_.getNormal().dot(trajectory.getLine().getR0() -
plane_.getCenter()) *
(plane_.getNormal().dot(trajectory.getPosition(0) - plane_.getCenter()) *
plane_.getNormal().dot(direction.cross(magneticfield)) * 2 * k)) -
direction.dot(plane_.getNormal()) / direction.getNorm()) /
(plane_.getNormal().dot(direction.cross(magneticfield)) * k);
......@@ -113,27 +112,29 @@ namespace corsika {
std::cout << " steplength to obs plane 2: " << MaxStepLength2 << std::endl;
return MaxStepLength2 *
(direction + direction.cross(magneticfield) * MaxStepLength2 * k / 2)
.norm() *
.getNorm() *
1.001;
} else if (MaxStepLength2 <= 0_m || MaxStepLength1 < MaxStepLength2) {
std::cout << " steplength to obs plane 1: " << MaxStepLength1 << std::endl;
return MaxStepLength1 *
(direction + direction.cross(magneticfield) * MaxStepLength2 * k / 2)
.norm() *
.getNorm() *
1.001;
}
}
TimeType const timeOfIntersection =
(plane_.getCenter() - trajectory.getStartPoint()).dot(plane_.getNormal()) /
trajectory.getVelocity().dot(plane_.getNormal());
(plane_.getCenter() - trajectory.getPosition(0)).dot(plane_.getNormal()) /
trajectory.getVelocity(0).dot(plane_.getNormal());
if (timeOfIntersection < TimeType::zero()) {
return std::numeric_limits<double>::infinity() * 1_m;
}
auto const pointOfIntersection = trajectory.getPosition(timeOfIntersection);
auto dist = (trajectory.getStartPoint() - pointOfIntersection).getNorm() * 1.0001;
double const fractionOfIntersection = timeOfIntersection / trajectory.getDuration();
auto const pointOfIntersection = trajectory.getPosition(fractionOfIntersection);
auto dist = (trajectory.getPosition(0) - pointOfIntersection).getNorm() * 1.0001;
CORSIKA_LOG_TRACE("ObservationPlane w/o magnetic field: getMaxStepLength l={} m",
dist / 1_m);
return dist;
......
......@@ -44,8 +44,8 @@ namespace corsika {
<< std::setw(width_) << std::scientific << std::setprecision(precision_) << start[2] / 1_m
<< std::setw(width_) << std::scientific << std::setprecision(precision_) << delta[0] / 1_m
<< std::setw(width_) << std::scientific << std::setprecision(precision_) << delta[1] / 1_m
<< std::setw(width) << std::scientific << std::setprecision(precision) << delta[2] / 1_m
<< std::setw(width) << std::scientific << std::setprecision(precision) << delta.norm() / 1_m
<< std::setw(width_) << std::scientific << std::setprecision(precision_) << delta[2] / 1_m
<< std::setw(width_) << std::scientific << std::setprecision(precision_) << delta.getNorm() / 1_m
<< '\n';
// clang-format on
......
......@@ -8,8 +8,6 @@
#pragma once
#include <corsika/modules/TrackingLine.hpp>
#include <corsika/framework/geometry/Point.hpp>
#include <corsika/framework/geometry/QuantityVector.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
......
/*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
*
* 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.
*/
#pragma once
#include <boost/histogram.hpp>
......
......@@ -82,7 +82,7 @@ namespace corsika {
Cascade(Cascade const&) = default;
Cascade(Cascade&&) = default;
~Cascade() = default;
Cascade & operator=(Cascade const&)) = default;
Cascade& operator=(Cascade const&) = default;
Cascade(Environment<MediumInterface> const& env, TTracking& tr, TProcessList& pl,
TStack& stack)
: environment_(env)
......
......@@ -48,7 +48,7 @@ namespace corsika {
inline Point getPosition(TimeType const t) const;
VelocityVec getVelocity(TimeType const t) const;
VelocityVec getVelocity(TimeType const t) const;
inline Point getPositionFromArclength(LengthType const l) const;
......
......@@ -37,8 +37,8 @@ namespace corsika {
Intersections(TimeType&& t)
: has_intersections_(true)
, intersections_(
std::make_pair(t, std::numeric_limits<TimeType::value_type>::infinity() * second)) {}
, intersections_(std::make_pair(
t, std::numeric_limits<TimeType::value_type>::infinity() * second)) {}
bool hasIntersections() const { return has_intersections_; }
......
......@@ -43,6 +43,8 @@ namespace corsika {
inline Point const& getStartPoint() const;
inline DirectionVector getDirection() const;
inline VelocityVector const& getVelocity() const;
private:
......
......@@ -37,7 +37,7 @@ namespace corsika {
LengthType getRadius() const;
void setRadius(LengthType const);
private:
Point center_;
LengthType radius_;
......
/***************************************************************************
* 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. *
***************************************************************************/
/*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
*
* 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.
*/
#pragma once
......
......@@ -13,7 +13,6 @@
#include <corsika/framework/geometry/Line.hpp>
#include <corsika/framework/geometry/Point.hpp>
#include <corsika/setup/SetupTrajectory.hpp>
#include <limits>
namespace corsika {
......@@ -25,7 +24,6 @@ namespace corsika {
template <typename TDerived>
class BaseExponential {
protected:
auto const& getImplementation() const;
// clang-format off
......@@ -70,12 +68,16 @@ namespace corsika {
public:
BaseExponential(Point const& point, MassDensityType rho0, LengthType lambda);
Point const& getAnchorPoint() const { return point_; }
MassDensityType getRho0() const { return rho0_; }
InverseLengthType getInvLambda() const { return invLambda_; }