diff --git a/GeoLib/AnalyticalGeometry.cpp b/GeoLib/AnalyticalGeometry.cpp
index 4422f56d6f5bde8af99fe6790c86cd2ca27073c3..bf1a4f862646d737b05b1b499b63b1210cea028a 100644
--- a/GeoLib/AnalyticalGeometry.cpp
+++ b/GeoLib/AnalyticalGeometry.cpp
@@ -706,4 +706,51 @@ lineSegmentIntersect2d(MathLib::Point3d const& a, MathLib::Point3d const& b,
     }
 }
 
+void sortSegments(
+    MathLib::Point3d const& seg_beg_pnt,
+    std::vector<GeoLib::LineSegment>& sub_segments)
+{
+    double const eps(std::numeric_limits<double>::epsilon());
+
+    auto findNextSegment = [&eps](
+        MathLib::Point3d const& seg_beg_pnt,
+        std::vector<GeoLib::LineSegment>& sub_segments,
+        std::vector<GeoLib::LineSegment>::iterator& sub_seg_it)
+    {
+        if (sub_seg_it == sub_segments.end())
+            return;
+        // find appropriate segment for the given segment begin point
+        auto act_beg_seg_it = std::find_if(
+            sub_seg_it, sub_segments.end(),
+            [&seg_beg_pnt, &eps](GeoLib::LineSegment const& seg)
+            {
+                return MathLib::sqrDist(seg_beg_pnt, seg.getBeginPoint()) < eps ||
+                       MathLib::sqrDist(seg_beg_pnt, seg.getEndPoint()) < eps;
+            });
+        if (act_beg_seg_it == sub_segments.end())
+            return;
+        // if necessary correct orientation of segment, i.e. swap beg and end
+        if (MathLib::sqrDist(seg_beg_pnt, act_beg_seg_it->getEndPoint()) <
+            MathLib::sqrDist(seg_beg_pnt, act_beg_seg_it->getBeginPoint()))
+            std::swap(act_beg_seg_it->getBeginPoint(),
+                      act_beg_seg_it->getEndPoint());
+        assert(sub_seg_it != sub_segments.end());
+        // exchange segments within the container
+        if (sub_seg_it != act_beg_seg_it)
+            std::swap(*sub_seg_it, *act_beg_seg_it);
+    };
+
+    // find start segment
+    auto seg_it = sub_segments.begin();
+    findNextSegment(seg_beg_pnt, sub_segments, seg_it);
+
+    while (seg_it != sub_segments.end())
+    {
+        MathLib::Point3d & new_seg_beg_pnt(seg_it->getEndPoint());
+        seg_it++;
+        if (seg_it != sub_segments.end())
+            findNextSegment(new_seg_beg_pnt, sub_segments, seg_it);
+    }
+}
+
 } // end namespace GeoLib
diff --git a/GeoLib/AnalyticalGeometry.h b/GeoLib/AnalyticalGeometry.h
index ae7af19a7d98ad4adfd7bda8649d2853c6664f05..a952d7f502c70d17e88d2be51901d886b28e4e61 100644
--- a/GeoLib/AnalyticalGeometry.h
+++ b/GeoLib/AnalyticalGeometry.h
@@ -384,6 +384,15 @@ void computeAndInsertAllIntersectionPoints(GeoLib::PointVec &pnt_vec,
 GeoLib::Polygon rotatePolygonToXY(GeoLib::Polygon const& polygon_in,
     MathLib::Vector3 & plane_normal);
 
+/// Sorts the vector of segments such that the \f$i\f$-th segment is connected
+/// with the \f$i+1\f$st segment, i.e. the end point of the \f$i\f$-th segment
+/// is the start point of the \f$i+1\f$st segment.
+/// The current implementation requires that all segments have to be
+/// connectable. In order to obtain a unique result the segments are sorted such
+/// that the begin point of the first segment is \c seg_beg_pnt.
+void sortSegments(MathLib::Point3d const& seg_beg_pnt,
+    std::vector<GeoLib::LineSegment>& sub_segments);
+
 } // end namespace GeoLib
 
 #include "AnalyticalGeometry-impl.h"