diff --git a/Applications/Utils/MeshGeoTools/VerticalSliceFromLayers.cpp b/Applications/Utils/MeshGeoTools/VerticalSliceFromLayers.cpp
index 8243043a04bb3e5dc3ff020e973d9607be31c40d..e6a2d691783c5a6135939a741b8e3cbc96166a74 100644
--- a/Applications/Utils/MeshGeoTools/VerticalSliceFromLayers.cpp
+++ b/Applications/Utils/MeshGeoTools/VerticalSliceFromLayers.cpp
@@ -193,10 +193,9 @@ void rotateGeometryToXY(std::vector<GeoLib::Point*>& points,
                         MathLib::DenseMatrix<double>& rotation_matrix,
                         double& z_shift)
 {
-    MathLib::Vector3 plane_normal;
-    double d;
     // compute the plane normal
-    GeoLib::getNewellPlane(points.begin(), points.end(), plane_normal, d);
+    auto const [plane_normal, d] =
+        GeoLib::getNewellPlane(points.begin(), points.end());
     // rotate points into x-y-plane
     GeoLib::computeRotationMatrixToXY(plane_normal, rotation_matrix);
     GeoLib::rotatePoints(rotation_matrix, points.begin(), points.end());
diff --git a/GeoLib/AnalyticalGeometry-impl.h b/GeoLib/AnalyticalGeometry-impl.h
index 6ad53498b1970e563c9f03d3e8bf78b1f5b33ad2..12086e2cf98f499053ae3a5ba82e47dfc94b51b7 100644
--- a/GeoLib/AnalyticalGeometry-impl.h
+++ b/GeoLib/AnalyticalGeometry-impl.h
@@ -10,12 +10,11 @@
 
 namespace GeoLib
 {
-
 template <typename InputIterator>
-void getNewellPlane (InputIterator pnts_begin, InputIterator pnts_end,
-                     MathLib::Vector3 &plane_normal,
-                     double& d)
+std::pair<MathLib::Vector3, double> getNewellPlane(InputIterator pnts_begin,
+                                                   InputIterator pnts_end)
 {
+    MathLib::Vector3 plane_normal;
     MathLib::Vector3 centroid;
     for (auto i=std::prev(pnts_end), j=pnts_begin; j!=pnts_end; i = j, ++j) {
         auto &pt_i = *(*i);
@@ -33,18 +32,19 @@ void getNewellPlane (InputIterator pnts_begin, InputIterator pnts_end,
     plane_normal.normalize();
     auto n_pnts(std::distance(pnts_begin, pnts_end));
     assert(n_pnts > 2);
+    double d = 0.0;
     if (n_pnts > 0)
     {
         d = MathLib::scalarProduct(centroid, plane_normal) / n_pnts;
     }
+    return std::make_pair(plane_normal, d);
 }
 
 template <class T_POINT>
-void getNewellPlane (const std::vector<T_POINT*>& pnts,
-                     MathLib::Vector3 &plane_normal,
-                     double& d)
+std::pair<MathLib::Vector3, double> getNewellPlane(
+    const std::vector<T_POINT*>& pnts)
 {
-    getNewellPlane(pnts.begin(), pnts.end(), plane_normal, d);
+    return getNewellPlane(pnts.begin(), pnts.end());
 }
 
 template <class T_POINT>
@@ -216,11 +216,9 @@ MathLib::DenseMatrix<double> rotatePointsToXY(InputIterator1 p_pnts_begin,
 {
     assert(std::distance(p_pnts_begin, p_pnts_end) > 2);
 
-    // calculate supporting plane
-    MathLib::Vector3 plane_normal;
-    double d;
     // compute the plane normal
-    GeoLib::getNewellPlane(p_pnts_begin, p_pnts_end, plane_normal, d);
+    auto const [plane_normal, d] =
+        GeoLib::getNewellPlane(p_pnts_begin, p_pnts_end);
 
     // rotate points into x-y-plane
     MathLib::DenseMatrix<double> rot_mat(3, 3);
diff --git a/GeoLib/AnalyticalGeometry.cpp b/GeoLib/AnalyticalGeometry.cpp
index 058f67044b7f686cd62e87e7173237440db1b970..9ee2d4b28bd48cdd2d2e5eb4075b4455d868b60f 100644
--- a/GeoLib/AnalyticalGeometry.cpp
+++ b/GeoLib/AnalyticalGeometry.cpp
@@ -374,7 +374,7 @@ GeoLib::Polygon rotatePolygonToXY(GeoLib::Polygon const& polygon_in,
 
     // 2 rotate points
     double d_polygon (0.0);
-    GeoLib::getNewellPlane (*polygon_pnts, plane_normal, d_polygon);
+    std::tie(plane_normal, d_polygon) = GeoLib::getNewellPlane(*polygon_pnts);
     MathLib::DenseMatrix<double> rot_mat(3,3);
     GeoLib::computeRotationMatrixToXY(plane_normal, rot_mat);
     GeoLib::rotatePoints(rot_mat, *polygon_pnts);
diff --git a/GeoLib/AnalyticalGeometry.h b/GeoLib/AnalyticalGeometry.h
index 1a8156f4b1d15e9b8c7809ef53de5bc8c64e4be1..d091d1e53073d9016bbce41bd9ec1b4345d7d7c3 100644
--- a/GeoLib/AnalyticalGeometry.h
+++ b/GeoLib/AnalyticalGeometry.h
@@ -54,13 +54,12 @@ Orientation getOrientationFast(MathLib::Point3d const& p0,
  * \cite Ericson:2004:RCD:1121584 .
  * @param pnts_begin Iterator pointing to the initial point of a closed polyline describing a polygon
  * @param pnts_end Iterator pointing to the element following the last point of a closed polyline describing a polygon
- * @param plane_normal the normal of the plane the polygon is located in
- * @param d parameter from the plane equation
+ * @return pair of plane_normal and the parameter d: the normal of the plane the
+ * polygon is located in, d parameter from the plane equation
  */
 template <typename InputIterator>
-void getNewellPlane (InputIterator pnts_begin, InputIterator pnts_end,
-                     MathLib::Vector3 &plane_normal,
-                     double& d);
+std::pair<MathLib::Vector3, double> getNewellPlane(InputIterator pnts_begin,
+                                                   InputIterator pnts_end);
 
 /**
  * compute a supporting plane (represented by plane_normal and the value d) for the polygon
@@ -68,18 +67,18 @@ void getNewellPlane (InputIterator pnts_begin, InputIterator pnts_end,
  * it holds \f$ n \cdot p + d = 0\f$. The Newell algorithm is described in
  * \cite Ericson:2004:RCD:1121584 .
  * @param pnts points of a closed polyline describing a polygon
- * @param plane_normal the normal of the plane the polygon is located in
- * @param d parameter from the plane equation
+ * @return pair of plane_normal and the parameter d: the normal of the plane the
+ * polygon is located in, parameter d from the plane equation
  */
 template <class T_POINT>
-void getNewellPlane (const std::vector<T_POINT*>& pnts,
-                     MathLib::Vector3 &plane_normal,
-                     double& d);
+std::pair<MathLib::Vector3, double> getNewellPlane(
+    const std::vector<T_POINT*>& pnts);
 
 /** Same as getNewellPlane(pnts, plane_normal, d).
  */
 template <class T_POINT>
-std::pair<MathLib::Vector3, double> getNewellPlane(const std::vector<T_POINT>& pnts);
+std::pair<MathLib::Vector3, double> getNewellPlane(
+    const std::vector<T_POINT>& pnts);
 
 /**
  * Computes a rotation matrix that rotates the given 2D normal vector parallel to X-axis