IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 22db345d authored by ralfulrich's avatar ralfulrich
Browse files

first small patches

parent 91194780
No related branches found
No related tags found
No related merge requests found
......@@ -24,13 +24,13 @@ namespace corsika {
}
template <typename TDerived>
GrammageType BaseExponential<TDerived>::integratedGrammage(
Trajectory<Line> const& line, units::si::LengthType vL,
Vector<units::si::dimensionless_d> const& axis) const {
if (vL == units::si::LengthType::zero()) { return units::si::GrammageType::zero(); }
GrammageType BaseExponential<TDerived>::getIntegratedGrammage(
Trajectory<Line> const& line, LengthType vL,
Vector<dimensionless_d> const& axis) const {
if (vL == LengthType::zero()) { return GrammageType::zero(); }
auto const uDotA = line.NormalizedDirection().dot(axis).magnitude();
auto const rhoStart = getImplementation().getMassDensity(line.GetR0());
auto const rhoStart = getImplementation().getMassDensity(line.getR0());
if (uDotA == 0) {
return vL * rhoStart;
......@@ -41,10 +41,10 @@ namespace corsika {
template <typename TDerived>
LengthType BaseExponential<TDerived>::getArclengthFromGrammage(
Trajectory<Line> const& line, units::si::GrammageType grammage,
Vector<units::si::dimensionless_d> const& axis) const {
Trajectory<Line> const& line, GrammageType grammage,
Vector<dimensionless_d> const& axis) const {
auto const uDotA = line.NormalizedDirection().dot(axis).magnitude();
auto const rhoStart = getImplementation().getMassDensity(line.GetR0());
auto const rhoStart = getImplementation().getMassDensity(line.getR0());
if (uDotA == 0) {
return grammage / rhoStart;
......@@ -54,15 +54,14 @@ namespace corsika {
return lambda_ / uDotA * log(logArg);
} else {
return std::numeric_limits<typename decltype(grammage)::value_type>::infinity() *
units::si::meter;
meter;
}
}
}
template <typename TDerived>
BaseExponential<TDerived>::BaseExponential(Point const& point,
units::si::MassDensityType rho0,
units::si::LengthType lambda)
BaseExponential<TDerived>::BaseExponential(Point const& point, MassDensityType rho0,
LengthType lambda)
: rho0_(rho0)
, lambda_(lambda)
, invLambda_(1 / lambda)
......
......@@ -21,17 +21,15 @@ namespace corsika {
template <typename T>
FlatExponential<T>::FlatExponential(Point const& point,
Vector<units::si::dimensionless_d> const& axis,
units::si::MassDensityType rho,
units::si::LengthType lambda,
Vector<dimensionless_d> const& axis,
MassDensityType rho, LengthType lambda,
NuclearComposition nuclComp)
: BaseExponential<FlatExponential<T>>(point, rho, lambda)
, axis_(axis)
, nuclComp_(nuclComp) {}
template <typename T>
units::si::MassDensityType FlatExponential<T>::getMassDensity(
Point const& point) const {
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_));
......@@ -43,16 +41,16 @@ namespace corsika {
}
template <typename T>
GrammageType FlatExponential<T>::integratedGrammage(
Trajectory<Line> const& line, units::si::LengthType to) const {
return BaseExponential<FlatExponential<T>>::integratedGrammage(line, to, axis_);
GrammageType FlatExponential<T>::getIntegratedGrammage(Trajectory<Line> const& line,
LengthType to) const {
return BaseExponential<FlatExponential<T>>::getIntegratedGrammage(line, to, axis_);
}
template <typename T>
LengthType FlatExponential<T>::getArclengthFromGrammage(
Trajectory<Line> const& line, units::si::GrammageType grammage) const {
return BaseExponential<FlatExponential<T>>::arclengthFromGrammage(line, grammage,
axis_);
LengthType FlatExponential<T>::getArclengthFromGrammage(Trajectory<Line> const& line,
GrammageType grammage) const {
return BaseExponential<FlatExponential<T>>::getArclengthFromGrammage(line, grammage,
axis_);
}
} // namespace corsika
......
......@@ -17,13 +17,13 @@
namespace corsika {
template <typename T>
HomogeneousMedium<T>::HomogeneousMedium(units::si::MassDensityType density,
HomogeneousMedium<T>::HomogeneousMedium(MassDensityType density,
NuclearComposition nuclComp)
: density_(density)
, nuclComp_(nuclComp) {}
template <typename T>
units::si::MassDensityType HomogeneousMedium<T>::getMassDensity(Point const&) const {
MassDensityType HomogeneousMedium<T>::getMassDensity(Point const&) const {
return density_;
}
template <typename T>
......@@ -33,13 +33,13 @@ namespace corsika {
template <typename T>
GrammageType HomogeneousMedium<T>::integratedGrammage(
Trajectory<Line> const&, units::si::LengthType to) const {
Trajectory<Line> const&, LengthType to) const {
return to * density_;
}
template <typename T>
LengthType HomogeneousMedium<T>::getArclengthFromGrammage(
Trajectory<Line> const&, units::si::GrammageType grammage) const {
Trajectory<Line> const&, GrammageType grammage) const {
return grammage / density_;
}
} // namespace corsika
......@@ -24,7 +24,7 @@ namespace corsika {
, densityFunction_(rhoTArgs...){}
template <typename T, typename TDensityFunction>
units::si::MassDensityType InhomogeneousMedium<T, TDensityFunction>::getMassDensity(
MassDensityType InhomogeneousMedium<T, TDensityFunction>::getMassDensity(
Point const& point) const {
return densityFunction_.evaluateAt(point);
}
......@@ -36,14 +36,14 @@ namespace corsika {
}
template <typename T, typename TDensityFunction>
GrammageType InhomogeneousMedium<T, TDensityFunction>::integratedGrammage(
Trajectory<Line> const& line, units::si::LengthType to) const {
return densityFunction_.integrateGrammage(line, to);
GrammageType InhomogeneousMedium<T, TDensityFunction>::getIntegratedGrammage(
Trajectory<Line> const& line, LengthType to) const {
return densityFunction_.getIntegrateGrammage(line, to);
}
template <typename T, typename TDensityFunction>
LengthType InhomogeneousMedium<T, TDensityFunction>::getArclengthFromGrammage(
Trajectory<Line> const& line, units::si::GrammageType grammage) const {
Trajectory<Line> const& line, GrammageType grammage) const {
return densityFunction_.getArclengthFromGrammage(line, grammage);
}
......
......@@ -21,16 +21,16 @@ namespace corsika {
template <typename TDerived>
auto LinearApproximationIntegrator<TDerived>::integrateGrammage(
Trajectory<Line> const& line, units::si::LengthType length) const {
auto const c0 = getImplementation().evaluateAt(line.GetPosition(0));
auto const c1 = getImplementation().rho_.FirstDerivative(line.GetPosition(0),
Trajectory<Line> const& line, LengthType length) const {
auto const c0 = getImplementation().evaluateAt(line.getPosition(0));
auto const c1 = getImplementation().rho_.FirstDerivative(line.getPosition(0),
line.NormalizedDirection());
return (c0 + 0.5 * c1 * length) * length;
}
template <typename TDerived>
auto LinearApproximationIntegrator<TDerived>::getArclengthFromGrammage(
Trajectory<Line> const& line, units::si::GrammageType grammage) const {
Trajectory<Line> const& line, GrammageType grammage) const {
auto const c0 = getImplementation().rho_(line.GetPosition(0));
auto const c1 = getImplementation().rho_.FirstDerivative(line.GetPosition(0),
line.NormalizedDirection());
......@@ -41,9 +41,8 @@ namespace corsika {
template <typename TDerived>
auto LinearApproximationIntegrator<TDerived>::getMaximumLength(
Trajectory<Line> const& line, [[maybe_unused]] double relError) const {
using namespace units::si;
[[maybe_unused]] auto const c1 = getImplementation().rho_.SecondDerivative(
line.GetPosition(0), line.NormalizedDirection());
line.getPosition(0), line.NormalizedDirection());
// todo: provide a real, working implementation
return 1_m * std::numeric_limits<double>::infinity();
......
......@@ -68,7 +68,7 @@ namespace corsika {
if (is_nucleus(compID)) {
return get_nucleus_A(compID) * fraction;
} else {
return get_mass(compID) / units::si::ConvertSIToHEP(constants::u) *
return get_mass(compID) / ConvertSIToHEP(constants::u) *
fraction;
}
})) {
......@@ -117,7 +117,7 @@ namespace corsika {
template <class TRNG>
corsika::Code NuclearComposition::sampleTarget(
std::vector<units::si::CrossSectionType> const& sigma, TRNG& randomStream) const {
std::vector<CrossSectionType> const& sigma, TRNG& randomStream) const {
using namespace units::si;
assert(sigma.size() == numberFractions_.size());
......
......@@ -14,10 +14,8 @@
namespace corsika {
template <typename TEnvModel>
ShowerAxis::ShowerAxis(Point const& pStart,
corsika::Vector<units::si::length_d> const& length,
Environment<TEnvModel> const& env,
int steps)
ShowerAxis::ShowerAxis(Point const& pStart, corsika::Vector<length_d> const& length,
Environment<TEnvModel> const& env, int steps)
: pointStart_(pStart)
, length_(length)
, max_length_(length_.norm())
......@@ -28,14 +26,14 @@ namespace corsika {
auto rho = [pStart, length, universe](double x) {
auto const p = pStart + length * x;
auto const* node = universe->GetContainingNode(p);
return node->GetModelProperties().getMassDensity(p).magnitude();
auto const* node = universe->getContainingNode(p);
return node->getModelProperties().getMassDensity(p).magnitude();
};
double error;
int k = 0;
X_[0] = units::si::GrammageType::zero();
auto sum = units::si::GrammageType::zero();
X_[0] = GrammageType::zero();
auto sum = GrammageType::zero();
for (int i = 1; i <= steps; ++i) {
auto const x_prev = (i - 1.) / steps;
......@@ -44,7 +42,7 @@ namespace corsika {
auto const r = boost::math::quadrature::gauss_kronrod<double, 15>::integrate(
rho, x_prev, x, 15, 1e-9, &error);
auto const result =
units::si::MassDensityType(phys::units::detail::magnitude_tag, r) * max_length_;
MassDensityType(phys::units::detail::magnitude_tag, r) * max_length_;
sum += result;
X_[i] = sum;
......@@ -66,7 +64,7 @@ namespace corsika {
if (lower < 0) {
CORSIKA_LOG_ERROR("cannot extrapolate to points behind point of injection l={} m",
l / 1_m);
l / 1_m);
throw std::runtime_error("cannot extrapolate to points behind point of injection");
}
......@@ -80,8 +78,8 @@ namespace corsika {
assert(0 <= lambda && lambda <= 1.);
CORSIKA_LOG_TRACE("ShowerAxis::X l={} m, lower={}, lambda={}, upper={}", l / 1_m, lower,
lambda, upper);
CORSIKA_LOG_TRACE("ShowerAxis::X l={} m, lower={}, lambda={}, upper={}", l / 1_m,
lower, lambda, upper);
// linear interpolation between X[lower] and X[upper]
return X_[upper] * lambda + X_[lower] * (1 - lambda);
......@@ -98,7 +96,7 @@ namespace corsika {
return X(projectedLength);
}
corsika::Vector<units::si::dimensionless_d> const& ShowerAxis::getDirection() const {
corsika::Vector<dimensionless_d> const& ShowerAxis::getDirection() const {
return axis_normalized_;
}
......
......@@ -16,14 +16,14 @@ namespace corsika {
template <typename T>
SlidingPlanarExponential<T>::SlidingPlanarExponential(
Point const& p0, units::si::MassDensityType rho0, units::si::LengthType lambda,
NuclearComposition nuclComp, units::si::LengthType referenceHeight)
Point const& p0, MassDensityType rho0, LengthType lambda,
NuclearComposition nuclComp, LengthType referenceHeight)
: BaseExponential<SlidingPlanarExponential<T>>(p0, rho0, lambda)
, nuclComp_(nuclComp)
, referenceHeight_(referenceHeight) {}
template <typename T>
units::si::MassDensityType SlidingPlanarExponential<T>::getMassDensity(
MassDensityType SlidingPlanarExponential<T>::getMassDensity(
Point const& point) const {
auto const height =
(point - BaseExponential<SlidingPlanarExponential<T>>::point_).norm() -
......@@ -38,19 +38,19 @@ namespace corsika {
}
template <typename T>
units::si::GrammageType SlidingPlanarExponential<T>::integratedGrammage(
Trajectory<Line> const& line, units::si::LengthType l) const {
GrammageType SlidingPlanarExponential<T>::integratedGrammage(
Trajectory<Line> const& line, LengthType l) const {
auto const axis =
(line.GetR0() - BaseExponential<SlidingPlanarExponential<T>>::point_).normalized();
(line.getR0() - BaseExponential<SlidingPlanarExponential<T>>::point_).normalized();
return BaseExponential<SlidingPlanarExponential<T>>::integratedGrammage(line, l,
axis);
}
template <typename T>
units::si::LengthType SlidingPlanarExponential<T>::arclengthFromGrammage(
Trajectory<Line> const& line, units::si::GrammageType const grammage) const {
LengthType SlidingPlanarExponential<T>::arclengthFromGrammage(
Trajectory<Line> const& line, GrammageType const grammage) const {
auto const axis =
(line.GetR0() - BaseExponential<SlidingPlanarExponential<T>>::point_).normalized();
(line.getR0() - BaseExponential<SlidingPlanarExponential<T>>::point_).normalized();
return BaseExponential<SlidingPlanarExponential<T>>::arclengthFromGrammage(
line, grammage, axis);
}
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment