grammage (+inverse) conversion based on approximation of the integral

Environment/DensityFunction.h

+/**
+ * (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

Environment/LinearIntegrator.h

+/**
+ * (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