From a2b61412b13893cc0d49cf48d609d139e945a64b Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@tu-dortmund.de>
Date: Mon, 11 Mar 2019 10:11:31 -0300
Subject: [PATCH] sliding planar exponential atmosphere
---
Environment/CMakeLists.txt | 1 +
Environment/SlidingPlanarExponential.h | 98 ++++++++++++++++++++++++++
2 files changed, 99 insertions(+)
create mode 100644 Environment/SlidingPlanarExponential.h
diff --git a/Environment/CMakeLists.txt b/Environment/CMakeLists.txt
index 47b973cb1..4b79788e6 100644
--- a/Environment/CMakeLists.txt
+++ b/Environment/CMakeLists.txt
@@ -10,6 +10,7 @@ set (
DensityFunction.h
Environment.h
FlatExponential.h
+ SlidingPlanarExponential.h
)
set (
diff --git a/Environment/SlidingPlanarExponential.h b/Environment/SlidingPlanarExponential.h
new file mode 100644
index 000000000..3dc059173
--- /dev/null
+++ b/Environment/SlidingPlanarExponential.h
@@ -0,0 +1,98 @@
+/*
+ * (c) Copyright 2019 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_Environment_SlidingPlanarExponential_h_
+#define _include_Environment_SlidingPlanarExponential_h_
+
+#include <corsika/environment/NuclearComposition.h>
+#include <corsika/geometry/Line.h>
+#include <corsika/environment/FlatExponential.h>
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/Trajectory.h>
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/random/RNGManager.h>
+#include <corsika/units/PhysicalUnits.h>
+
+#include <cassert>
+#include <limits>
+
+/**
+ *
+ */
+
+namespace corsika::environment {
+
+ template <class T>
+ class SlidingPlanarExponential : public T {
+ corsika::units::si::MassDensityType const fRho0;
+ units::si::LengthType const fLambda;
+ units::si::InverseLengthType const fInvLambda;
+ NuclearComposition const fNuclComp;
+ geometry::Vector<units::si::dimensionless_d> const fAxis;
+ geometry::Point const fP0;
+
+ public:
+ SlidingPlanarExponential(geometry::Point const& p0,
+ geometry::Vector<units::si::dimensionless_d> const& axis,
+ units::si::MassDensityType rho, units::si::LengthType lambda,
+ NuclearComposition pNuclComp)
+ : fRho0(rho)
+ , fLambda(lambda)
+ , fInvLambda(1 / lambda)
+ , fNuclComp(pNuclComp)
+ , fAxis(axis)
+ , fP0(p0) {}
+
+ corsika::units::si::MassDensityType GetMassDensity(
+ corsika::geometry::Point const& p) const override {
+ auto const height = (p - fP0).norm();
+ return fRho0 * exp(fInvLambda * height);
+ }
+
+ NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
+
+ corsika::units::si::GrammageType IntegratedGrammage(
+ corsika::geometry::Trajectory<corsika::geometry::Line> const& line,
+ corsika::units::si::LengthType pTo) const override {
+
+ auto const axis = (fP0 - line.GetR0()).normalized();
+ auto const vDotA = line.NormalizedDirection().dot(fAxis).magnitude();
+
+ if (vDotA == 0) {
+ return pTo * GetMassDensity(line.GetR0());
+ } else {
+ return GetMassDensity(line.GetR0()) * (fLambda / vDotA) *
+ (exp(vDotA * pTo / fLambda) - 1);
+ }
+ }
+
+ corsika::units::si::LengthType ArclengthFromGrammage(
+ corsika::geometry::Trajectory<corsika::geometry::Line> const& line,
+ corsika::units::si::GrammageType pGrammage) const override {
+ auto const axis = (fP0 - line.GetR0()).normalized();
+ auto const vDotA = line.NormalizedDirection().dot(fAxis).magnitude();
+
+ if (vDotA == 0) {
+ return pGrammage / GetMassDensity(line.GetR0());
+ } else {
+ auto const logArg = pGrammage * vDotA / (fRho0 * fLambda) + 1;
+ if (logArg > 0) {
+ return fLambda / vDotA * log(pGrammage * vDotA / (fRho0 * fLambda) + 1);
+ } else {
+ return std::numeric_limits<typename decltype(
+ pGrammage)::value_type>::infinity() *
+ corsika::units::si::meter;
+ }
+ }
+ }
+ };
+
+} // namespace corsika::environment
+#endif
--
GitLab