From 4a956d6b31b25d645196d9d957d4f1c960d35f69 Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@kit.edu>
Date: Wed, 30 Jan 2019 20:14:33 +0100
Subject: [PATCH] grammage (+inverse) conversion based on approximation of the
 integral

---
 Environment/DensityFunction.h  | 39 ++++++++++++++++++++++++++++
 Environment/LinearIntegrator.h | 46 ++++++++++++++++++++++++++++++++++
 2 files changed, 85 insertions(+)
 create mode 100644 Environment/DensityFunction.h
 create mode 100644 Environment/LinearIntegrator.h

diff --git a/Environment/DensityFunction.h b/Environment/DensityFunction.h
new file mode 100644
index 000000000..a374c3479
--- /dev/null
+++ b/Environment/DensityFunction.h
@@ -0,0 +1,39 @@
+/**
+ * (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_environment_DensityFunction_h_
+#define _include_environment_DensityFunction_h_
+
+#include <corsika/environment/LinearIntegrator.h>
+#include <corsika/geometry/Line.h>
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/Trajectory.h>
+
+namespace corsika::environment {
+
+  template <class TDerivableRho>
+  // TODO make LinearApproximator a template parameter to facilitate exchangability
+  class DensityFunction : public LinearApproximator<DensityFunction<TDerivableRho>> {
+    friend class LinearApproximator<DensityFunction<TDerivableRho>>;
+
+    TDerivableRho fRho; //!< functor for density
+
+  public:
+    DensityFunction(TDerivableRho rho)
+        : fRho(rho) {}
+
+    corsika::units::si::MassDensityType EvaluateAt(
+        corsika::geometry::Point const& p) const {
+      return fRho(p);
+    }
+  };
+} // namespace corsika::environment
+
+#endif
diff --git a/Environment/LinearIntegrator.h b/Environment/LinearIntegrator.h
new file mode 100644
index 000000000..1bea44a29
--- /dev/null
+++ b/Environment/LinearIntegrator.h
@@ -0,0 +1,46 @@
+/**
+ * (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_environment_LinearIntegrator_h_
+#define _include_environment_LinearIntegrator_h_
+
+#include <corsika/geometry/Line.h>
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/Trajectory.h>
+
+namespace corsika::environment {
+  template <class TDerived>
+  class LinearApproximator {
+    auto const& GetImplementation() const { return *static_cast<TDerived const*>(this); }
+
+  public:
+    auto IntegrateGrammage(
+        corsika::geometry::Trajectory<corsika::geometry::Line> const& line,
+        corsika::units::si::LengthType length) const {
+      auto const c0 = GetImplementation().fRho(line.GetPosition(0));
+      auto const c1 = GetImplementation().fRho.Derivative<1>(line.GetPosition(0),
+                                                             line.NormalizedDirection());
+
+      return (c0 + 0.5 * c1 * length) * length;
+    }
+
+    auto ArclengthFromGrammage(
+        corsika::geometry::Trajectory<corsika::geometry::Line> const& line,
+        corsika::units::si::GrammageType grammage) const {
+      auto const c0 = GetImplementation().fRho(line.GetPosition(0));
+      auto const c1 = GetImplementation().fRho.Derivative<1>(line.GetPosition(0),
+                                                             line.NormalizedDirection());
+
+      return (1 - 0.5 * grammage * c1 / (c0 * c0)) * grammage / c0;
+    }
+  };
+} // namespace corsika::environment
+
+#endif
-- 
GitLab