Index: DamEngine/trunk/src/Deltares.DamEngine.Data/Geometry/GeometryPointString.cs
===================================================================
diff -u -r4000 -r4052
--- DamEngine/trunk/src/Deltares.DamEngine.Data/Geometry/GeometryPointString.cs	(.../GeometryPointString.cs)	(revision 4000)
+++ DamEngine/trunk/src/Deltares.DamEngine.Data/Geometry/GeometryPointString.cs	(.../GeometryPointString.cs)	(revision 4052)
@@ -24,970 +24,969 @@
 using System.Linq;
 using Deltares.DamEngine.Data.Standard;
 
-namespace Deltares.DamEngine.Data.Geometry
+namespace Deltares.DamEngine.Data.Geometry;
+
+/// <summary>
+/// Classifications of an objected with regards to a geometric line in the XZ-plane.
+/// </summary>
+public enum RelativeXzPosition
 {
     /// <summary>
-    /// Classifications of an objected with regards to a geometric line in the XZ-plane.
+    /// Indicates that the object is considered 'on' the geometric line.
     /// </summary>
-    public enum RelativeXzPosition
-    {
-        /// <summary>
-        /// Indicates that the object is considered 'on' the geometric line.
-        /// </summary>
-        OnGeometricLine,
+    OnGeometricLine,
 
-        /// <summary>
-        /// Indicates that the object is considered 'above' the geometric line (object Z
-        /// coordinate higher than that of the line).
-        /// </summary>
-        AboveGeometricLine,
+    /// <summary>
+    /// Indicates that the object is considered 'above' the geometric line (object Z
+    /// coordinate higher than that of the line).
+    /// </summary>
+    AboveGeometricLine,
 
-        /// <summary>
-        /// Indicates that the object is considered 'below' the geometric line (object Z
-        /// coordinates lower than that of the line).
-        /// </summary>
-        BelowGeometricLine,
+    /// <summary>
+    /// Indicates that the object is considered 'below' the geometric line (object Z
+    /// coordinates lower than that of the line).
+    /// </summary>
+    BelowGeometricLine,
 
-        /// <summary>
-        /// Indicates that the object is considered 'outside the scope' of the geometric line.
-        /// </summary>
-        BeyondGeometricLine
-    }
+    /// <summary>
+    /// Indicates that the object is considered 'outside the scope' of the geometric line.
+    /// </summary>
+    BeyondGeometricLine
+}
 
+/// <summary>
+/// Type of extrapolation
+/// </summary>
+public enum ExtraPolationMode
+{
     /// <summary>
-    /// Type of extrapolation
+    /// No extrapolation should be used.
     /// </summary>
-    public enum ExtraPolationMode
-    {
-        /// <summary>
-        /// No extrapolation should be used.
-        /// </summary>
-        Beyond,
+    Beyond,
 
-        /// <summary>
-        /// Used constant (0th order) extrapolation at the extremes.
-        /// </summary>
-        Horizontal
-    }
+    /// <summary>
+    /// Used constant (0th order) extrapolation at the extremes.
+    /// </summary>
+    Horizontal
+}
 
+/// <summary>
+/// Collection of points (X,Z) in world coordinates used for describing Surface lines
+/// </summary>
+public class GeometryPointString : GeometryObject
+{
     /// <summary>
-    /// Collection of points (X,Z) in world coordinates used for describing Surface lines
+    /// Matching distance where a point within this range is considered on the same point.
     /// </summary>
-    public class GeometryPointString : GeometryObject
-    {
-        /// <summary>
-        /// Matching distance where a point within this range is considered on the same point.
-        /// </summary>
-        private const double epsilon = GeometryConstants.Accuracy;
+    private const double epsilon = GeometryConstants.Accuracy;
 
-        /// <summary>
-        /// The calculate points as protected field (to be able to prevent recursive calls to CalcPoints)
-        /// </summary>
-        protected readonly List<Point2D> calcPoints = new List<Point2D>();
+    /// <summary>
+    /// The calculate points as protected field (to be able to prevent recursive calls to CalcPoints)
+    /// </summary>
+    protected readonly List<Point2D> calcPoints = new List<Point2D>();
 
-        private readonly List<GeometryPoint> points = new List<GeometryPoint>();
-        private bool isFrozen;
-        private bool hasNaNx;
-        private double frozenMaxZ = double.NaN;
+    private readonly List<GeometryPoint> points = new List<GeometryPoint>();
+    private bool isFrozen;
+    private bool hasNaNx;
+    private double frozenMaxZ = double.NaN;
 
-        // sortedPoints must never be used outside this class. Either the GPS concerned must have sorted points but then they already are 
-        // (eg. surfaceline, headline) or they may be unsorted in which case using the sorted list in other classes leads to errors (eg.
-        // geometrysurfaces, waternetlines) 
-        private List<Point2D> sortedPoints;
+    // sortedPoints must never be used outside this class. Either the GPS concerned must have sorted points but then they already are 
+    // (eg. surfaceline, headline) or they may be unsorted in which case using the sorted list in other classes leads to errors (eg.
+    // geometrysurfaces, waternetlines) 
+    private List<Point2D> sortedPoints;
 
-        /// <summary>
-        /// Gets the <see cref="GeometryPoint"/> at the specified index.
-        /// </summary>
-        /// <value>
-        /// The <see cref="GeometryPoint"/>.
-        /// </value>
-        /// <param name="index">The index.</param>
-        /// <exception cref="ArgumentOutOfRangeException">When <paramref name="index"/> less
-        /// than zero or is greater or equals to <see cref="Count"/>.</exception>
-        public Point2D this[int index]
+    /// <summary>
+    /// Gets the <see cref="GeometryPoint"/> at the specified index.
+    /// </summary>
+    /// <value>
+    /// The <see cref="GeometryPoint"/>.
+    /// </value>
+    /// <param name="index">The index.</param>
+    /// <exception cref="ArgumentOutOfRangeException">When <paramref name="index"/> less
+    /// than zero or is greater or equals to <see cref="Count"/>.</exception>
+    public Point2D this[int index]
+    {
+        get
         {
-            get
-            {
-                return calcPoints[index];
-            }
+            return calcPoints[index];
         }
+    }
 
-        /// <summary>
-        /// The calculate points (to be used in calcualtion instead of Points for better performance)
-        /// </summary>
-        public virtual List<Point2D> CalcPoints
+    /// <summary>
+    /// The calculate points (to be used in calcualtion instead of Points for better performance)
+    /// </summary>
+    public virtual List<Point2D> CalcPoints
+    {
+        get
         {
-            get
-            {
-                return calcPoints;
-            }
+            return calcPoints;
         }
+    }
 
-        /// <summary>
-        /// List of points that describe the physical surface line or surface.
-        /// </summary>
-        /// <value>
-        /// The points.
-        /// </value>
-        public virtual IList<GeometryPoint> Points
+    /// <summary>
+    /// List of points that describe the physical surface line or surface.
+    /// </summary>
+    /// <value>
+    /// The points.
+    /// </value>
+    public virtual IList<GeometryPoint> Points
+    {
+        get
         {
-            get
-            {
-                return points;
-            }
+            return points;
         }
+    }
 
-        /// <summary>
-        /// Gets the count of the points.
-        /// </summary>
-        /// <value>
-        /// The count.
-        /// </value>
-        public int Count
+    /// <summary>
+    /// Gets the count of the points.
+    /// </summary>
+    /// <value>
+    /// The count.
+    /// </value>
+    public int Count
+    {
+        get
         {
-            get
-            {
-                return calcPoints.Count;
-            }
+            return calcPoints.Count;
         }
+    }
 
-        /// <summary>
-        /// Freezes this instance.
-        /// </summary>
-        public void Freeze()
+    /// <summary>
+    /// Freezes this instance.
+    /// </summary>
+    public void Freeze()
+    {
+        if (!isFrozen)
         {
-            if (!isFrozen)
+            sortedPoints = new List<Point2D>(calcPoints.Count);
+            foreach (Point2D point in calcPoints)
             {
-                sortedPoints = new List<Point2D>(calcPoints.Count);
-                foreach (Point2D point in calcPoints)
-                {
-                    sortedPoints.Add(point);
-                    hasNaNx = hasNaNx || double.IsNaN(point.X);
-                    frozenMaxZ = Math.Max(frozenMaxZ, point.Z);
-                }
-
-                sortedPoints.Sort();
+                sortedPoints.Add(point);
+                hasNaNx = hasNaNx || double.IsNaN(point.X);
+                frozenMaxZ = Math.Max(frozenMaxZ, point.Z);
             }
 
-            isFrozen = true;
+            sortedPoints.Sort();
         }
 
-        /// <summary>
-        /// Clones this object except the points
-        /// </summary>
-        /// <returns></returns>
-        public virtual GeometryPointString Clone()
-        {
-            var clone = new GeometryPointString();
-            this.CloneProperties(clone); // exludes the points !
-            clone.Points.Clear();
-            foreach (GeometryPoint point in Points)
-            {
-                clone.Points.Add(new GeometryPoint(point));
-            }
+        isFrozen = true;
+    }
 
-            clone.SyncCalcPoints();
-            return clone;
+    /// <summary>
+    /// Clones this object except the points
+    /// </summary>
+    /// <returns></returns>
+    public virtual GeometryPointString Clone()
+    {
+        var clone = new GeometryPointString();
+        this.CloneProperties(clone); // exludes the points !
+        clone.Points.Clear();
+        foreach (GeometryPoint point in Points)
+        {
+            clone.Points.Add(new GeometryPoint(point));
         }
 
-        /// <summary>
-        /// Gets the minimum z.
-        /// Note: CalcPoints must be used instead of calcPoints as otherwise unclear behaviour of Linq spoils the result
-        /// Change back to calcPoints and Benchmark4_01l in the WaternetCreatorTests will fail!
-        /// </summary>
-        /// <returns></returns>
-        public double GetMinZ()
+        clone.SyncCalcPoints();
+        return clone;
+    }
+
+    /// <summary>
+    /// Gets the minimum z.
+    /// Note: CalcPoints must be used instead of calcPoints as otherwise unclear behaviour of Linq spoils the result
+    /// Change back to calcPoints and Benchmark4_01l in the WaternetCreatorTests will fail!
+    /// </summary>
+    /// <returns></returns>
+    public double GetMinZ()
+    {
+        return CalcPoints.Any()
+                   ? CalcPoints.Min(p => p.Z)
+                   : double.NaN;
+    }
+
+    /// <summary>
+    /// Gets the maximum z.
+    /// Note: CalcPoints must be used instead of calcPoints as otherwise unclear behaviour of Linq spoils the result
+    /// Change back to calcPoints and Benchmark4_01l in the WaternetCreatorTests will fail!
+    /// </summary>
+    /// <returns></returns>
+    public double GetMaxZ()
+    {
+        if (isFrozen)
         {
-            return CalcPoints.Any()
-                       ? CalcPoints.Min(p => p.Z)
-                       : double.NaN;
+            return frozenMaxZ;
         }
 
-        /// <summary>
-        /// Gets the maximum z.
-        /// Note: CalcPoints must be used instead of calcPoints as otherwise unclear behaviour of Linq spoils the result
-        /// Change back to calcPoints and Benchmark4_01l in the WaternetCreatorTests will fail!
-        /// </summary>
-        /// <returns></returns>
-        public double GetMaxZ()
+        return CalcPoints.Any()
+                   ? CalcPoints.Max(p => p.Z)
+                   : double.NaN;
+    }
+
+    /// <summary>
+    /// The minimal X value among <see cref="CalcPoints"/>.
+    /// Note: CalcPoints must be used instead of calcPoints as otherwise unclear behaviour of Linq spoils the result
+    /// Change back to calcPoints and Benchmark4_01l in the WaternetCreatorTests will fail!
+    /// </summary>
+    /// <returns>The minimal X value or <see cref="double.NaN"/> in case there are no points.</returns>
+    public double GetMinX()
+    {
+        if (isFrozen && !hasNaNx)
         {
-            if (isFrozen)
-            {
-                return frozenMaxZ;
-            }
+            return sortedPoints[0].X;
+        }
 
-            return CalcPoints.Any()
-                       ? CalcPoints.Max(p => p.Z)
-                       : double.NaN;
+        return CalcPoints.Any(p => !double.IsNaN(p.X))
+                   ? CalcPoints.Where(p => !double.IsNaN(p.X)).Min(p => p.X)
+                   : double.NaN;
+    }
+
+    /// <summary>
+    /// Gets the maximum x.
+    /// Note: CalcPoints must be used instead of calcPoints as otherwise unclear behaviour of Linq spoils the result
+    /// Change back to calcPoints and Benchmark4_01l in the WaternetCreatorTests will fail!
+    /// </summary>
+    /// <returns></returns>
+    public double GetMaxX()
+    {
+        if (isFrozen && !hasNaNx)
+        {
+            return sortedPoints[sortedPoints.Count - 1].X;
         }
 
-        /// <summary>
-        /// The minimal X value among <see cref="CalcPoints"/>.
-        /// Note: CalcPoints must be used instead of calcPoints as otherwise unclear behaviour of Linq spoils the result
-        /// Change back to calcPoints and Benchmark4_01l in the WaternetCreatorTests will fail!
-        /// </summary>
-        /// <returns>The minimal X value or <see cref="double.NaN"/> in case there are no points.</returns>
-        public double GetMinX()
+        return CalcPoints.Any(p => !double.IsNaN(p.X))
+                   ? CalcPoints.Where(p => !double.IsNaN(p.X)).Max(p => p.X)
+                   : double.NaN;
+    }
+
+    /// <summary>
+    /// Finds all vertical XZ intersections with this geometric point string, returning
+    /// the highest value among the intersections.
+    /// </summary>
+    /// <param name="x">X coordinate</param>
+    /// <returns>
+    /// The z value determined by <see cref="CalcPoints" /> or <see cref="double.NaN" />
+    /// when <see cref="CalcPoints" /> is empty.
+    /// </returns>
+    /// <remarks>
+    /// Uses constant extrapolation and linear interpolation.
+    /// </remarks>        
+    public double GetZAtUnsortedX(double x)
+    {
+        if (calcPoints.Count > 0)
         {
-            if (isFrozen && !hasNaNx)
+            var verticalLineAtX = new Line(new Point2D
+                                           {
+                                               X = x,
+                                               Z = GetMaxZ()
+                                           },
+                                           new Point2D
+                                           {
+                                               X = x,
+                                               Z = GetMinZ()
+                                           });
+            if (Math.Abs(verticalLineAtX.BeginPoint.Z - verticalLineAtX.EndPoint.Z) < GeometryConstants.Accuracy)
             {
-                return sortedPoints[0].X;
+                verticalLineAtX.BeginPoint.Z += 1.0;
+                verticalLineAtX.EndPoint.Z -= 1.0;
             }
 
-            return CalcPoints.Any(p => !double.IsNaN(p.X))
-                       ? CalcPoints.Where(p => !double.IsNaN(p.X)).Min(p => p.X)
-                       : double.NaN;
+            IList<Point2D> intersectionPoints = IntersectionPointsXzWithLineXz(verticalLineAtX);
+            if (intersectionPoints.Count != 0)
+            {
+                return intersectionPoints.Max(gp => gp.Z);
+            }
+
+            // Remain consistent with GetZAtX, do constant extrapolation:
+            double zAtX;
+            if (DoConstantExtrapolationXz(x, out zAtX))
+            {
+                return zAtX;
+            }
         }
 
-        /// <summary>
-        /// Gets the maximum x.
-        /// Note: CalcPoints must be used instead of calcPoints as otherwise unclear behaviour of Linq spoils the result
-        /// Change back to calcPoints and Benchmark4_01l in the WaternetCreatorTests will fail!
-        /// </summary>
-        /// <returns></returns>
-        public double GetMaxX()
+        return double.NaN;
+    }
+
+    /// <summary>
+    /// Retrieves the Z value for the given x.
+    /// </summary>
+    /// <param name="x">X coordinate</param>
+    /// <returns>
+    /// The z value determined by <see cref="CalcPoints"/> or <see cref="double.NaN"/>
+    /// when <see cref="CalcPoints"/> is empty.
+    /// </returns>
+    /// <remarks>
+    /// <para>Uses constant extrapolation and linear interpolation.</para>
+    /// </remarks>
+    public virtual double GetZatX(double x)
+    {
+        if (calcPoints.Any())
         {
-            if (isFrozen && !hasNaNx)
+            double zAtX;
+            if (DoConstantExtrapolationXz(x, out zAtX))
             {
-                return sortedPoints[sortedPoints.Count - 1].X;
+                return zAtX;
             }
-
-            return CalcPoints.Any(p => !double.IsNaN(p.X))
-                       ? CalcPoints.Where(p => !double.IsNaN(p.X)).Max(p => p.X)
-                       : double.NaN;
         }
 
-        /// <summary>
-        /// Finds all vertical XZ intersections with this geometric point string, returning
-        /// the highest value among the intersections.
-        /// </summary>
-        /// <param name="x">X coordinate</param>
-        /// <returns>
-        /// The z value determined by <see cref="CalcPoints" /> or <see cref="double.NaN" />
-        /// when <see cref="CalcPoints" /> is empty.
-        /// </returns>
-        /// <remarks>
-        /// Uses constant extrapolation and linear interpolation.
-        /// </remarks>        
-        public double GetZAtUnsortedX(double x)
+        for (var i = 0; i < calcPoints.Count - 1; i++)
         {
-            if (calcPoints.Count > 0)
-            {
-                var verticalLineAtX = new Line(new Point2D
-                                               {
-                                                   X = x,
-                                                   Z = GetMaxZ()
-                                               },
-                                               new Point2D
-                                               {
-                                                   X = x,
-                                                   Z = GetMinZ()
-                                               });
-                if (Math.Abs(verticalLineAtX.BeginPoint.Z - verticalLineAtX.EndPoint.Z) < GeometryConstants.Accuracy)
-                {
-                    verticalLineAtX.BeginPoint.Z += 1.0;
-                    verticalLineAtX.EndPoint.Z -= 1.0;
-                }
+            Point2D current = calcPoints[i];
+            Point2D next = calcPoints[i + 1];
 
-                IList<Point2D> intersectionPoints = IntersectionPointsXzWithLineXz(verticalLineAtX);
-                if (intersectionPoints.Count != 0)
-                {
-                    return intersectionPoints.Max(gp => gp.Z);
-                }
+            double leftOffset = x - current.X;
+            double rightOffset = next.X - x;
 
-                // Remain consistent with GetZAtX, do constant extrapolation:
-                double zAtX;
-                if (DoConstantExtrapolationXz(x, out zAtX))
-                {
-                    return zAtX;
-                }
+            if (Math.Abs(leftOffset) < epsilon)
+            {
+                return current.Z;
             }
 
-            return double.NaN;
-        }
-
-        /// <summary>
-        /// Retrieves the Z value for the given x.
-        /// </summary>
-        /// <param name="x">X coordinate</param>
-        /// <returns>
-        /// The z value determined by <see cref="CalcPoints"/> or <see cref="double.NaN"/>
-        /// when <see cref="CalcPoints"/> is empty.
-        /// </returns>
-        /// <remarks>
-        /// <para>Uses constant extrapolation and linear interpolation.</para>
-        /// </remarks>
-        public virtual double GetZatX(double x)
-        {
-            if (calcPoints.Any())
+            if (Math.Abs(rightOffset) < epsilon)
             {
-                double zAtX;
-                if (DoConstantExtrapolationXz(x, out zAtX))
-                {
-                    return zAtX;
-                }
+                return next.Z;
             }
 
-            for (var i = 0; i < calcPoints.Count - 1; i++)
+            if (leftOffset >= 0 && rightOffset >= 0)
             {
-                Point2D current = calcPoints[i];
-                Point2D next = calcPoints[i + 1];
+                double fraction = leftOffset / (leftOffset + rightOffset);
 
-                double leftOffset = x - current.X;
-                double rightOffset = next.X - x;
+                return (1.0 - fraction) * current.Z + fraction * next.Z;
+            }
+        }
 
-                if (Math.Abs(leftOffset) < epsilon)
-                {
-                    return current.Z;
-                }
+        return Double.NaN;
+    }
 
-                if (Math.Abs(rightOffset) < epsilon)
-                {
-                    return next.Z;
-                }
+    /// <summary>
+    /// Gets the z at x starting from index i in the point list.
+    /// This is only meant as a fast(er) version for GetZAtX for lists that are sorted by X, so use with care.
+    /// </summary>
+    /// <param name="x">The x.</param>
+    /// <param name="i">The i.</param>
+    /// <returns></returns>
+    public virtual double GetZatX(double x, ref int i)
+    {
+        for (; i < calcPoints.Count - 1; i++)
+        {
+            Point2D current = calcPoints[i];
+            double leftOffset = x - current.X;
+            Point2D next = calcPoints[i + 1];
+            double rightOffset = next.X - x;
 
-                if (leftOffset >= 0 && rightOffset >= 0)
-                {
-                    double fraction = leftOffset / (leftOffset + rightOffset);
+            if (leftOffset < epsilon)
+            {
+                return current.Z;
+            }
 
-                    return (1.0 - fraction) * current.Z + fraction * next.Z;
-                }
+            if (rightOffset >= epsilon)
+            {
+                double fraction = leftOffset / (leftOffset + rightOffset);
+                return (1.0 - fraction) * current.Z + fraction * next.Z;
             }
 
-            return Double.NaN;
+            if (i + 1 == calcPoints.Count - 1)
+            {
+                return next.Z;
+            }
         }
 
-        /// <summary>
-        /// Gets the z at x starting from index i in the point list.
-        /// This is only meant as a fast(er) version for GetZAtX for lists that are sorted by X, so use with care.
-        /// </summary>
-        /// <param name="x">The x.</param>
-        /// <param name="i">The i.</param>
-        /// <returns></returns>
-        public virtual double GetZatX(double x, ref int i)
-        {
-            for (; i < calcPoints.Count - 1; i++)
-            {
-                Point2D current = calcPoints[i];
-                double leftOffset = x - current.X;
-                Point2D next = calcPoints[i + 1];
-                double rightOffset = next.X - x;
+        return Double.NaN;
+    }
 
-                if (leftOffset < epsilon)
-                {
-                    return current.Z;
-                }
+    /// <summary>
+    /// Gets all z-values at x for a line.
+    /// </summary>
+    /// <param name="x">The x.</param>
+    /// <returns></returns>
+    public List<double> GetAllZatXForLine(double x)
+    {
+        return GetAllZatX(x, false);
+    }
 
-                if (rightOffset >= epsilon)
-                {
-                    double fraction = leftOffset / (leftOffset + rightOffset);
-                    return (1.0 - fraction) * current.Z + fraction * next.Z;
-                }
+    /// <summary>
+    /// Gets all z values at x for a surface.
+    /// </summary>
+    /// <param name="x">The x.</param>
+    /// <returns></returns>
+    public List<double> GetAllZatXForSurface(double x)
+    {
+        return GetAllZatX(x, true);
+    }
 
-                if (i + 1 == calcPoints.Count - 1)
-                {
-                    return next.Z;
-                }
+    /// <summary>
+    /// Finds the first intersection of the line with a given horizontal line.
+    /// </summary>
+    /// <param name="z">The height level of the horizontal line.</param>
+    /// <param name="xmin">Optional: Discard intersections whose X value is less then 
+    /// or equal to this value.</param>
+    /// <returns>The first intersection point matching the criteria.</returns>
+    /// <seealso cref="GetXatZStartingAt"/>
+    public double GetXatZ(double z, double xmin = Double.MinValue)
+    {
+        return GetXatZStartingAt(0, z, xmin);
+    }
+
+    /// <summary>
+    /// Finds the first intersection of the line with a given horizontal line.
+    /// </summary>
+    /// <param name="start">Start considering <see cref="GeometryPoint"/> instances from
+    /// <see cref="CalcPoints"/> starting from this index.</param>
+    /// <param name="z">The height level of the horizontal line.</param>
+    /// <param name="xmin">Optional: Discard intersections whose X value is less then 
+    /// or equal to this value.</param>
+    /// <returns>The first intersection point matching the criteria.</returns>
+    /// <seealso cref="GetXatZ"/>
+    public double GetXatZStartingAt(int start, double z, double xmin = Double.MinValue)
+    {
+        for (int i = start; i < calcPoints.Count - 1; i++)
+        {
+            Point2D current = calcPoints[i];
+            Point2D next = calcPoints[i + 1];
+
+            if (IsSegmentNotCrossingZ(z, current, next)) // Performance micro-optimization
+            {
+                continue;
             }
 
-            return Double.NaN;
+            double x = GetXIntersectingZ(z, current, next);
+            if (x > xmin)
+            {
+                return x;
+            }
         }
 
-        /// <summary>
-        /// Gets all z-values at x for a line.
-        /// </summary>
-        /// <param name="x">The x.</param>
-        /// <returns></returns>
-        public List<double> GetAllZatXForLine(double x)
+        return Double.NaN;
+    }
+
+    /// <summary>
+    /// Gets the points at x.
+    /// </summary>
+    /// <param name="x">The x.</param>
+    /// <returns></returns>
+    public IEnumerable<Point2D> GetPointsAtX(double x)
+    {
+        return (from Point2D point in calcPoints
+                where point.X.AlmostEquals(x, GeometryPoint.Precision)
+                select point);
+    }
+
+    /// <summary>
+    /// Gets the Point2D at the specified coordinates.
+    /// </summary>
+    /// <param name="x">The x.</param>
+    /// <param name="z">The z.</param>
+    /// <returns></returns>
+    public Point2D GetPointAt(double x, double z)
+    {
+        return calcPoints.FirstOrDefault(
+            point => point.X.AlmostEquals(x, GeometryPoint.Precision) &&
+                     point.Z.AlmostEquals(z, GeometryPoint.Precision));
+    }
+
+    /// <summary>
+    /// Gets the GeometryPoint at the specified coordinates.
+    /// </summary>
+    /// <param name="x">The x.</param>
+    /// <param name="z">The z.</param>
+    /// <returns></returns>
+    public GeometryPoint GetGeometryPointAt(double x, double z)
+    {
+        return Points.FirstOrDefault(
+            point => point.X.AlmostEquals(x, GeometryPoint.Precision) &&
+                     point.Z.AlmostEquals(z, GeometryPoint.Precision));
+    }
+
+    //#Bka: this should never be 1 method! Split it into 2 methods
+    // It is used for sliplanes only and therefor uses Points instead of calcPoints
+    /// <summary>
+    /// Adds the or remove list object.
+    /// </summary>
+    /// <param name="aPoint">a point.</param>
+    /// <param name="aToAdd">if set to <c>true</c> [a to add].</param>
+    /// <param name="aIndex">a index.</param>
+    public virtual void AddOrRemoveListObject(object aPoint, bool aToAdd, int aIndex)
+    {
+        var point = aPoint as GeometryPoint;
+        if (point == null)
         {
-            return GetAllZatX(x, false);
+            return;
         }
 
-        /// <summary>
-        /// Gets all z values at x for a surface.
-        /// </summary>
-        /// <param name="x">The x.</param>
-        /// <returns></returns>
-        public List<double> GetAllZatXForSurface(double x)
+        if (aToAdd)
         {
-            return GetAllZatX(x, true);
+            Points.Insert(aIndex, point); //#bka what happens if index already exists? what if point exists? What if both???                
         }
-
-        /// <summary>
-        /// Finds the first intersection of the line with a given horizontal line.
-        /// </summary>
-        /// <param name="z">The height level of the horizontal line.</param>
-        /// <param name="xmin">Optional: Discard intersections whose X value is less then 
-        /// or equal to this value.</param>
-        /// <returns>The first intersection point matching the criteria.</returns>
-        /// <seealso cref="GetXatZStartingAt"/>
-        public double GetXatZ(double z, double xmin = Double.MinValue)
+        else
         {
-            return GetXatZStartingAt(0, z, xmin);
+            if (Points.Contains(point))
+            {
+                Points.Remove(point);
+            }
         }
+    }
 
-        /// <summary>
-        /// Finds the first intersection of the line with a given horizontal line.
-        /// </summary>
-        /// <param name="start">Start considering <see cref="GeometryPoint"/> instances from
-        /// <see cref="CalcPoints"/> starting from this index.</param>
-        /// <param name="z">The height level of the horizontal line.</param>
-        /// <param name="xmin">Optional: Discard intersections whose X value is less then 
-        /// or equal to this value.</param>
-        /// <returns>The first intersection point matching the criteria.</returns>
-        /// <seealso cref="GetXatZ"/>
-        public double GetXatZStartingAt(int start, double z, double xmin = Double.MinValue)
+    /// <summary>
+    /// Returns ALL intersection points that are found for a line costructed at level z, 
+    /// including the double ones as the number of points can be relevant!
+    /// </summary>
+    /// <param name="z">the level at which the line for intersections is constructed.</param>
+    /// <returns></returns>
+    public IList<double> IntersectionsXAtZ(double z)
+    {
+        var intersectionsX = new List<double>();
+
+        var lineAtZ = new Line
         {
-            for (int i = start; i < calcPoints.Count - 1; i++)
-            {
-                Point2D current = calcPoints[i];
-                Point2D next = calcPoints[i + 1];
+            BeginPoint = new Point2D(GetMinX(), z),
+            EndPoint = new Point2D(GetMaxX(), z)
+        };
 
-                if (IsSegmentNotCrossingZ(z, current, next)) // Performance micro-optimization
-                {
-                    continue;
-                }
+        intersectionsX.AddRange(IntersectionPointsXzWithLineXzWithAllPoints(lineAtZ).Select(gp => gp.X));
+        return intersectionsX;
+    }
 
-                double x = GetXIntersectingZ(z, current, next);
-                if (x > xmin)
-                {
-                    return x;
-                }
-            }
+    /// <summary>
+    /// Returns ALL intersection points that are found, including the double ones as the number of points can be relevant!
+    /// </summary>
+    /// <param name="line"></param>
+    /// <returns></returns>
+    public IList<Point2D> IntersectionPointsXzWithLineXzWithAllPoints(Line line)
+    {
+        return IntersectionPointsWithLineCore(line, true);
+    }
 
-            return Double.NaN;
-        }
+    /// <summary>
+    /// Returns the UNIQUE intersectionpoints (so can not be used where number of interscections is relevant)
+    /// </summary>
+    /// <param name="line"></param>
+    /// <returns></returns>
+    public IList<Point2D> IntersectionPointsXzWithLineXz(Line line)
+    {
+        return IntersectionPointsWithLineCore(line, false);
+    }
 
-        /// <summary>
-        /// Gets the points at x.
-        /// </summary>
-        /// <param name="x">The x.</param>
-        /// <returns></returns>
-        public IEnumerable<Point2D> GetPointsAtX(double x)
-        {
-            return (from Point2D point in calcPoints
-                    where point.X.AlmostEquals(x, GeometryPoint.Precision)
-                    select point);
-        }
+    /// <summary>
+    /// Find intersection (xz-plane) from this surface with another surface
+    /// or phratic line 
+    /// </summary>
+    /// <param name="list"></param>
+    /// <returns></returns>
+    /// <remarks>Considers all <see cref="GeometryPoint"/> in <see cref="CalcPoints"/> to
+    /// such that they describe a closed loop.</remarks>
+    /// <seealso cref="IntersectionXzPointsWithGeometryPointList"/>
+    public IList<Point2D> ClosedGeometryIntersectionXzPointsWithGeometryPointList(IList<Point2D> list)
+    {
+        return IntersectWithPointsListCore(list, true);
+    }
 
-        /// <summary>
-        /// Gets the Point2D at the specified coordinates.
-        /// </summary>
-        /// <param name="x">The x.</param>
-        /// <param name="z">The z.</param>
-        /// <returns></returns>
-        public Point2D GetPointAt(double x, double z)
-        {
-            return calcPoints.FirstOrDefault(
-                point => point.X.AlmostEquals(x, GeometryPoint.Precision) &&
-                         point.Z.AlmostEquals(z, GeometryPoint.Precision));
-        }
+    /// <summary>
+    /// Finds all intersections in the XZ-plane the given <see cref="GeometryPointString"/>.
+    /// </summary>
+    /// <param name="externalSurface">The geometry point string.</param>
+    /// <returns></returns>
+    /// <seealso cref="IntersectionXzPointsWithGeometryPointList"/>
+    public List<Point2D> IntersectionXzPointsWithGeometryString(GeometryPointString externalSurface)
+    {
+        return IntersectionXzPointsWithGeometryPointList(externalSurface.CalcPoints);
+    }
 
-        /// <summary>
-        /// Gets the GeometryPoint at the specified coordinates.
-        /// </summary>
-        /// <param name="x">The x.</param>
-        /// <param name="z">The z.</param>
-        /// <returns></returns>
-        public GeometryPoint GetGeometryPointAt(double x, double z)
-        {
-            return Points.FirstOrDefault(
-                point => point.X.AlmostEquals(x, GeometryPoint.Precision) &&
-                         point.Z.AlmostEquals(z, GeometryPoint.Precision));
-        }
+    /// <summary>
+    /// Sorts the points by x ascending (only to be used for Surface lines).
+    /// </summary>
+    public virtual void SortPointsByXAscending()
+    {
+        calcPoints.Sort();
+        points.Sort();
+    }
 
-        //#Bka: this should never be 1 method! Split it into 2 methods
-        // It is used for sliplanes only and therefor uses Points instead of calcPoints
-        /// <summary>
-        /// Adds the or remove list object.
-        /// </summary>
-        /// <param name="aPoint">a point.</param>
-        /// <param name="aToAdd">if set to <c>true</c> [a to add].</param>
-        /// <param name="aIndex">a index.</param>
-        public virtual void AddOrRemoveListObject(object aPoint, bool aToAdd, int aIndex)
+    /// <summary>
+    /// Removes the non ascending points on surface line.
+    /// </summary>
+    public void RemoveNonAscendingPointsOnSurfaceLine()
+    {
+        for (int i = points.Count - 1; i > 0; i--)
         {
-            var point = aPoint as GeometryPoint;
-            if (point == null)
+            if (Math.Abs(points[i].X - points[i - 1].X) < GeometryConstants.Accuracy)
             {
-                return;
+                if (Math.Abs(points[i].Z - points[i - 1].Z) > GeometryConstants.Accuracy)
+                {
+                    points.Remove(points[i]);
+                }
             }
+        }
+    }
 
-            if (aToAdd)
+    /// <summary>
+    /// Gets a part of the geometry point string defined by a begin and end point
+    /// </summary>
+    /// <param name="begin"></param>
+    /// <param name="end"></param>
+    /// <returns></returns>
+    public GeometryPointString GetPart(double begin, double end)
+    {
+        var part = new GeometryPointString();
+        var filling = false;
+        var filled = false;
+
+        for (var i = 0; i < calcPoints.Count; i++)
+        {
+            if (!filling && !filled)
             {
-                Points.Insert(aIndex, point); //#bka what happens if index already exists? what if point exists? What if both???                
+                filling = calcPoints[i].X >= begin - epsilon;
             }
             else
             {
-                if (Points.Contains(point))
+                filled = calcPoints[i].X >= end - epsilon;
+                if (filled)
                 {
-                    Points.Remove(point);
+                    filling = false;
                 }
             }
-        }
 
-        /// <summary>
-        /// Returns ALL intersection points that are found for a line costructed at level z, 
-        /// including the double ones as the number of points can be relevant!
-        /// </summary>
-        /// <param name="z">the level at which the line for intersections is constructed.</param>
-        /// <returns></returns>
-        public IList<double> IntersectionsXAtZ(double z)
-        {
-            var intersectionsX = new List<double>();
-
-            var lineAtZ = new Line
+            if (filling)
             {
-                BeginPoint = new Point2D(GetMinX(), z),
-                EndPoint = new Point2D(GetMaxX(), z)
-            };
-
-            intersectionsX.AddRange(IntersectionPointsXzWithLineXzWithAllPoints(lineAtZ).Select(gp => gp.X));
-            return intersectionsX;
+                part.calcPoints.Add(calcPoints[i]);
+            }
         }
 
-        /// <summary>
-        /// Returns ALL intersection points that are found, including the double ones as the number of points can be relevant!
-        /// </summary>
-        /// <param name="line"></param>
-        /// <returns></returns>
-        public IList<Point2D> IntersectionPointsXzWithLineXzWithAllPoints(Line line)
-        {
-            return IntersectionPointsWithLineCore(line, true);
-        }
+        var beginPoint = new Point2D(begin, GetZatX(begin));
+        var endPoint = new Point2D(end, GetZatX(end));
 
-        /// <summary>
-        /// Returns the UNIQUE intersectionpoints (so can not be used where number of interscections is relevant)
-        /// </summary>
-        /// <param name="line"></param>
-        /// <returns></returns>
-        public IList<Point2D> IntersectionPointsXzWithLineXz(Line line)
+        if (part.calcPoints.Count == 0)
         {
-            return IntersectionPointsWithLineCore(line, false);
+            part.calcPoints.Add(beginPoint);
+            part.calcPoints.Add(endPoint);
         }
 
-        /// <summary>
-        /// Find intersection (xz-plane) from this surface with another surface
-        /// or phratic line 
-        /// </summary>
-        /// <param name="list"></param>
-        /// <returns></returns>
-        /// <remarks>Considers all <see cref="GeometryPoint"/> in <see cref="CalcPoints"/> to
-        /// such that they describe a closed loop.</remarks>
-        /// <seealso cref="IntersectionXzPointsWithGeometryPointList"/>
-        public IList<Point2D> ClosedGeometryIntersectionXzPointsWithGeometryPointList(IList<Point2D> list)
+        if (!part.calcPoints[0].LocationEquals(beginPoint))
         {
-            return IntersectWithPointsListCore(list, true);
+            part.calcPoints.Insert(0, beginPoint);
         }
 
-        /// <summary>
-        /// Finds all intersections in the XZ-plane the given <see cref="GeometryPointString"/>.
-        /// </summary>
-        /// <param name="externalSurface">The geometry point string.</param>
-        /// <returns></returns>
-        /// <seealso cref="IntersectionXzPointsWithGeometryPointList"/>
-        public List<Point2D> IntersectionXzPointsWithGeometryString(GeometryPointString externalSurface)
+        if (!part.calcPoints[part.calcPoints.Count - 1].LocationEquals(endPoint))
         {
-            return IntersectionXzPointsWithGeometryPointList(externalSurface.CalcPoints);
+            part.calcPoints.Add(endPoint);
         }
 
-        /// <summary>
-        /// Sorts the points by x ascending (only to be used for Surface lines).
-        /// </summary>
-        public virtual void SortPointsByXAscending()
-        {
-            calcPoints.Sort();
-            points.Sort();
-        }
+        SyncPoints();
+        return part;
+    }
 
-        /// <summary>
-        /// Removes the non ascending points on surface line.
-        /// </summary>
-        public void RemoveNonAscendingPointsOnSurfaceLine()
+    /// <summary>
+    /// Removes all double points at a location, if consecutive
+    /// </summary>
+    public void CondensePoints()
+    {
+        for (int i = calcPoints.Count - 1; i > 0; i--)
         {
-            for (int i = points.Count - 1; i > 0; i--)
+            if (calcPoints[i].LocationEquals(calcPoints[i - 1]))
             {
-                if (Math.Abs(points[i].X - points[i - 1].X) < GeometryConstants.Accuracy)
-                {
-                    if (Math.Abs(points[i].Z - points[i - 1].Z) > GeometryConstants.Accuracy)
-                    {
-                        points.Remove(points[i]);
-                    }
-                }
+                calcPoints.RemoveAt(i);
             }
         }
+    }
 
-        /// <summary>
-        /// Gets a part of the geometry point string defined by a begin and end point
-        /// </summary>
-        /// <param name="begin"></param>
-        /// <param name="end"></param>
-        /// <returns></returns>
-        public GeometryPointString GetPart(double begin, double end)
-        {
-            var part = new GeometryPointString();
-            var filling = false;
-            var filled = false;
+    /// <summary>
+    /// Removes all points which don't influence the shape of the line
+    /// </summary>
+    public void RemoveUnnecessaryPoints()
+    {
+        const double slopeTolerance = 1E-10;
 
-            for (var i = 0; i < calcPoints.Count; i++)
-            {
-                if (!filling && !filled)
-                {
-                    filling = calcPoints[i].X >= begin - epsilon;
-                }
-                else
-                {
-                    filled = calcPoints[i].X >= end - epsilon;
-                    if (filled)
-                    {
-                        filling = false;
-                    }
-                }
+        CondensePoints();
 
-                if (filling)
-                {
-                    part.calcPoints.Add(calcPoints[i]);
-                }
-            }
+        for (int i = calcPoints.Count - 2; i > 1; i--)
+        {
+            // if the slope of the line before the point is equal to the slope after the point, the point can be removed
+            double slopeBefore = (calcPoints[i].Z - calcPoints[i - 1].Z) / (calcPoints[i].X - calcPoints[i - 1].X);
+            double slopeAfter = (calcPoints[i + 1].Z - calcPoints[i].Z) / (calcPoints[i + 1].X - calcPoints[i].X);
 
-            var beginPoint = new Point2D(begin, GetZatX(begin));
-            var endPoint = new Point2D(end, GetZatX(end));
-
-            if (part.calcPoints.Count == 0)
+            if (Routines2D.AreEqual(slopeBefore, slopeAfter, slopeTolerance))
             {
-                part.calcPoints.Add(beginPoint);
-                part.calcPoints.Add(endPoint);
+                calcPoints.RemoveAt(i);
             }
+        }
+    }
 
-            if (!part.calcPoints[0].LocationEquals(beginPoint))
+    /// <summary>
+    /// Synchronizes the calculation points.
+    /// </summary>
+    public void SyncCalcPoints()
+    {
+        calcPoints.Clear();
+        foreach (GeometryPoint geometryPoint in Points)
+        {
+            var p2D = new Point2D
             {
-                part.calcPoints.Insert(0, beginPoint);
-            }
+                X = geometryPoint.X,
+                Z = geometryPoint.Z
+            };
+            calcPoints.Add(p2D);
+        }
+    }
 
-            if (!part.calcPoints[part.calcPoints.Count - 1].LocationEquals(endPoint))
+    /// <summary>
+    /// Synchronizes the points.
+    /// </summary>
+    public void SyncPoints()
+    {
+        points.Clear();
+        foreach (Point2D p2D in calcPoints)
+        {
+            var geometryPoint = new GeometryPoint
             {
-                part.calcPoints.Add(endPoint);
-            }
-
-            SyncPoints();
-            return part;
+                X = p2D.X,
+                Z = p2D.Z
+            };
+            points.Add(geometryPoint);
         }
+    }
 
-        /// <summary>
-        /// Removes all double points at a location, if consecutive
-        /// </summary>
-        public void CondensePoints()
+    /// <summary>
+    /// Gets the surrounding rectangle around the geometry point string
+    /// </summary>
+    /// <returns></returns>
+    public override GeometryBounds GetGeometryBounds()
+    {
+        if (!Points.Any())
         {
-            for (int i = calcPoints.Count - 1; i > 0; i--)
+            // Sync with calcPoints
+            SyncPoints();
+            // if still no points, then return null
+            if (!Points.Any())
             {
-                if (calcPoints[i].LocationEquals(calcPoints[i - 1]))
-                {
-                    calcPoints.RemoveAt(i);
-                }
+                return null;
             }
         }
 
-        /// <summary>
-        /// Removes all points which don't influence the shape of the line
-        /// </summary>
-        public void RemoveUnnecessaryPoints()
+        GeometryBounds bounds = Points[0].GetGeometryBounds();
+        for (var i = 1; i < Points.Count; i++)
         {
-            const double slopeTolerance = 1E-10;
+            GeometryPoint point = Points[i];
 
-            CondensePoints();
+            bounds.Left = Math.Min(bounds.Left, point.X);
+            bounds.Right = Math.Max(bounds.Right, point.X);
+            bounds.Top = Math.Max(bounds.Top, point.Z);
+            bounds.Bottom = Math.Min(bounds.Bottom, point.Z);
+        }
 
-            for (int i = calcPoints.Count - 2; i > 1; i--)
-            {
-                // if the slope of the line before the point is equal to the slope after the point, the point can be removed
-                double slopeBefore = (calcPoints[i].Z - calcPoints[i - 1].Z) / (calcPoints[i].X - calcPoints[i - 1].X);
-                double slopeAfter = (calcPoints[i + 1].Z - calcPoints[i].Z) / (calcPoints[i + 1].X - calcPoints[i].X);
+        return bounds;
+    }
 
-                if (Routines2D.AreEqual(slopeBefore, slopeAfter, slopeTolerance))
-                {
-                    calcPoints.RemoveAt(i);
-                }
-            }
-        }
+    private IList<Point2D> IntersectionPointsXzClosedStringWithLineXz(Line line)
+    {
+        IList<Point2D> intersectionPointsWithLine = IntersectionPointsXzWithLineXz(line);
 
-        /// <summary>
-        /// Synchronizes the calculation points.
-        /// </summary>
-        public void SyncCalcPoints()
+        // close the poly line 
+        if (calcPoints.Count > 0)
         {
-            calcPoints.Clear();
-            foreach (GeometryPoint geometryPoint in Points)
-            {
-                var p2D = new Point2D
-                {
-                    X = geometryPoint.X,
-                    Z = geometryPoint.Z
-                };
-                calcPoints.Add(p2D);
-            }
+            DoIntersectAndAddToCollection(line, calcPoints[calcPoints.Count - 1], calcPoints[0],
+                                          intersectionPointsWithLine, false);
         }
 
-        /// <summary>
-        /// Synchronizes the points.
-        /// </summary>
-        public void SyncPoints()
-        {
-            points.Clear();
-            foreach (Point2D p2D in calcPoints)
-            {
-                var geometryPoint = new GeometryPoint
-                {
-                    X = p2D.X,
-                    Z = p2D.Z
-                };
-                points.Add(geometryPoint);
-            }
-        }
+        return intersectionPointsWithLine;
+    }
 
-        /// <summary>
-        /// Gets the surrounding rectangle around the geometry point string
-        /// </summary>
-        /// <returns></returns>
-        public override GeometryBounds GetGeometryBounds()
+    /// <summary>
+    /// Finds all intersections in the XZ-plane the given <see cref="GeometryPoint"/> list.
+    /// </summary>
+    /// <param name="list">The list.</param>
+    /// <returns></returns>
+    /// <seealso cref="ClosedGeometryIntersectionXzPointsWithGeometryPointList"/>
+    /// <seealso cref="IntersectionXzPointsWithGeometryString"/>
+    private List<Point2D> IntersectionXzPointsWithGeometryPointList(IList<Point2D> list)
+    {
+        return IntersectWithPointsListCore(list, false);
+    }
+
+    /// <summary>
+    /// Checks if constant extrapolation can be applied, and if so set the Z value.
+    /// </summary>
+    /// <param name="x">The evaluated X coordinate.</param>
+    /// <param name="z">Output param: Extrapolated Z value, or <see cref="double.NaN"/>
+    /// when no extrapolation is possible.</param>
+    /// <returns>True if extrapolation possible; false otherwise.</returns>
+    private bool DoConstantExtrapolationXz(double x, out double z)
+    {
+        if (calcPoints.Count > 0)
         {
-            if (!Points.Any())
+            Point2D first = calcPoints[0];
+            if (x < first.X || Math.Abs(x - first.X) < epsilon)
             {
-                // Sync with calcPoints
-                SyncPoints();
-                // if still no points, then return null
-                if (!Points.Any())
-                {
-                    return null;
-                }
+                z = first.Z;
+                return true;
             }
 
-            GeometryBounds bounds = Points[0].GetGeometryBounds();
-            for (var i = 1; i < Points.Count; i++)
+            Point2D last = calcPoints[calcPoints.Count - 1];
+            if (x > last.X)
             {
-                GeometryPoint point = Points[i];
-
-                bounds.Left = Math.Min(bounds.Left, point.X);
-                bounds.Right = Math.Max(bounds.Right, point.X);
-                bounds.Top = Math.Max(bounds.Top, point.Z);
-                bounds.Bottom = Math.Min(bounds.Bottom, point.Z);
+                z = last.Z;
+                return true;
             }
-
-            return bounds;
         }
 
-        private IList<Point2D> IntersectionPointsXzClosedStringWithLineXz(Line line)
-        {
-            IList<Point2D> intersectionPointsWithLine = IntersectionPointsXzWithLineXz(line);
+        z = double.NaN;
+        return false;
+    }
 
-            // close the poly line 
-            if (calcPoints.Count > 0)
-            {
-                DoIntersectAndAddToCollection(line, calcPoints[calcPoints.Count - 1], calcPoints[0],
-                                              intersectionPointsWithLine, false);
-            }
+    /// <summary>
+    /// Can be used for a Line or for a Surface where a surface is supposed to closed.
+    /// In case of a waternet line it is possible
+    /// to have more z values at a give x coor
+    /// Furthermore a z is not needed at al x values
+    /// </summary>
+    /// <param name="x"></param>
+    /// <param name="asSurface"></param>
+    /// <returns></returns>
+    private List<double> GetAllZatX(double x, bool asSurface)
+    {
+        var result = new List<double>();
 
-            return intersectionPointsWithLine;
+        if (calcPoints.Count == 1 && Math.Abs(calcPoints[0].X - x) < epsilon)
+        {
+            result.Add(calcPoints[0].Z);
         }
 
-        /// <summary>
-        /// Finds all intersections in the XZ-plane the given <see cref="GeometryPoint"/> list.
-        /// </summary>
-        /// <param name="list">The list.</param>
-        /// <returns></returns>
-        /// <seealso cref="ClosedGeometryIntersectionXzPointsWithGeometryPointList"/>
-        /// <seealso cref="IntersectionXzPointsWithGeometryString"/>
-        private List<Point2D> IntersectionXzPointsWithGeometryPointList(IList<Point2D> list)
+        int pointsCount = calcPoints.Count - 1;
+        if (asSurface)
         {
-            return IntersectWithPointsListCore(list, false);
+            pointsCount = calcPoints.Count;
         }
 
-        /// <summary>
-        /// Checks if constant extrapolation can be applied, and if so set the Z value.
-        /// </summary>
-        /// <param name="x">The evaluated X coordinate.</param>
-        /// <param name="z">Output param: Extrapolated Z value, or <see cref="double.NaN"/>
-        /// when no extrapolation is possible.</param>
-        /// <returns>True if extrapolation possible; false otherwise.</returns>
-        private bool DoConstantExtrapolationXz(double x, out double z)
+        for (var i = 0; i < pointsCount; i++)
         {
-            if (calcPoints.Count > 0)
+            Point2D current;
+            Point2D next;
+            if (i == calcPoints.Count - 1)
             {
-                Point2D first = calcPoints[0];
-                if (x < first.X || Math.Abs(x - first.X) < epsilon)
-                {
-                    z = first.Z;
-                    return true;
-                }
-
-                Point2D last = calcPoints[calcPoints.Count - 1];
-                if (x > last.X)
-                {
-                    z = last.Z;
-                    return true;
-                }
+                current = calcPoints[i];
+                next = calcPoints[0];
             }
+            else
+            {
+                current = calcPoints[i];
+                next = calcPoints[i + 1];
+            }
 
-            z = double.NaN;
-            return false;
-        }
+            double leftOffset = x - current.X;
+            double rightOffset = next.X - x;
 
-        /// <summary>
-        /// Can be used for a Line or for a Surface where a surface is supposed to closed.
-        /// In case of a waternet line it is possible
-        /// to have more z values at a give x coor
-        /// Furthermore a z is not needed at al x values
-        /// </summary>
-        /// <param name="x"></param>
-        /// <param name="asSurface"></param>
-        /// <returns></returns>
-        private List<double> GetAllZatX(double x, bool asSurface)
-        {
-            var result = new List<double>();
-
-            if (calcPoints.Count == 1 && Math.Abs(calcPoints[0].X - x) < epsilon)
+            var matchedWithAPoint = false;
+            if (Math.Abs(leftOffset) < epsilon)
             {
-                result.Add(calcPoints[0].Z);
+                result.Add(current.Z);
+                matchedWithAPoint = true;
             }
 
-            int pointsCount = calcPoints.Count - 1;
-            if (asSurface)
+            if (Math.Abs(rightOffset) < epsilon)
             {
-                pointsCount = calcPoints.Count;
+                result.Add(next.Z);
+                matchedWithAPoint = true;
             }
 
-            for (var i = 0; i < pointsCount; i++)
+            if (!matchedWithAPoint)
             {
-                Point2D current;
-                Point2D next;
-                if (i == calcPoints.Count - 1)
+                if (leftOffset > 0 && rightOffset > 0)
                 {
-                    current = calcPoints[i];
-                    next = calcPoints[0];
-                }
-                else
-                {
-                    current = calcPoints[i];
-                    next = calcPoints[i + 1];
-                }
+                    double fraction = leftOffset / (leftOffset + rightOffset);
 
-                double leftOffset = x - current.X;
-                double rightOffset = next.X - x;
-
-                var matchedWithAPoint = false;
-                if (Math.Abs(leftOffset) < epsilon)
-                {
-                    result.Add(current.Z);
-                    matchedWithAPoint = true;
+                    result.Add((1.0 - fraction) * current.Z + fraction * next.Z);
                 }
 
-                if (Math.Abs(rightOffset) < epsilon)
+                // if both ofsets are negative the waterline goes back
+                if ((leftOffset < 0) && (rightOffset < 0))
                 {
-                    result.Add(next.Z);
-                    matchedWithAPoint = true;
-                }
+                    double fraction = rightOffset / (rightOffset + leftOffset);
 
-                if (!matchedWithAPoint)
-                {
-                    if (leftOffset > 0 && rightOffset > 0)
-                    {
-                        double fraction = leftOffset / (leftOffset + rightOffset);
-
-                        result.Add((1.0 - fraction) * current.Z + fraction * next.Z);
-                    }
-
-                    // if both ofsets are negative the waterline goes back
-                    if ((leftOffset < 0) && (rightOffset < 0))
-                    {
-                        double fraction = rightOffset / (rightOffset + leftOffset);
-
-                        result.Add((1.0 - fraction) * next.Z + fraction * current.Z);
-                    }
+                    result.Add((1.0 - fraction) * next.Z + fraction * current.Z);
                 }
             }
-
-            return result.Distinct().ToList();
         }
 
-        private static bool IsSegmentNotCrossingZ(double z, Point2D current, Point2D next)
-        {
-            if (double.IsNaN(z))
-            {
-                return true;
-            }
+        return result.Distinct().ToList();
+    }
 
-            double leftOffset = Math.Abs(current.Z - z);
-            double rightOffset = Math.Abs(next.Z - z);
-
-            int currentSign = leftOffset < epsilon ? 0 : Math.Sign(current.Z - z);
-            int nextSign = rightOffset < epsilon ? 0 : Math.Sign(next.Z - z);
-            return currentSign == nextSign && currentSign != 0;
+    private static bool IsSegmentNotCrossingZ(double z, Point2D current, Point2D next)
+    {
+        if (double.IsNaN(z))
+        {
+            return true;
         }
 
-        private static double GetXIntersectingZ(double z, Point2D current, Point2D next)
-        {
-            double leftOffset = Math.Abs(current.Z - z);
-            if (leftOffset < epsilon)
-            {
-                return current.X;
-            }
+        double leftOffset = Math.Abs(current.Z - z);
+        double rightOffset = Math.Abs(next.Z - z);
 
-            double rightOffset = Math.Abs(next.Z - z);
-            if (rightOffset < epsilon)
-            {
-                return next.X;
-            }
+        int currentSign = leftOffset < epsilon ? 0 : Math.Sign(current.Z - z);
+        int nextSign = rightOffset < epsilon ? 0 : Math.Sign(next.Z - z);
+        return currentSign == nextSign && currentSign != 0;
+    }
 
-            double fraction = leftOffset / (leftOffset + rightOffset);
-            return GeneralMathRoutines.LinearInterpolate(current.X, next.X, fraction);
+    private static double GetXIntersectingZ(double z, Point2D current, Point2D next)
+    {
+        double leftOffset = Math.Abs(current.Z - z);
+        if (leftOffset < epsilon)
+        {
+            return current.X;
         }
 
-        private IList<Point2D> IntersectionPointsWithLineCore(Line line, bool allowDuplicates)
+        double rightOffset = Math.Abs(next.Z - z);
+        if (rightOffset < epsilon)
         {
-            var intersectionPointsWithLine = new List<Point2D>();
+            return next.X;
+        }
 
-            for (var pointIndex = 0; pointIndex < calcPoints.Count - 1; pointIndex++)
-            {
-                DoIntersectAndAddToCollection(line, calcPoints[pointIndex], calcPoints[pointIndex + 1],
-                                              intersectionPointsWithLine, allowDuplicates);
-            }
+        double fraction = leftOffset / (leftOffset + rightOffset);
+        return GeneralMathRoutines.LinearInterpolate(current.X, next.X, fraction);
+    }
 
-            return intersectionPointsWithLine;
-        }
+    private IList<Point2D> IntersectionPointsWithLineCore(Line line, bool allowDuplicates)
+    {
+        var intersectionPointsWithLine = new List<Point2D>();
 
-        private static void DoIntersectAndAddToCollection(Line line, Point2D begin, Point2D end, ICollection<Point2D> intersectionPointsWithLine, bool allowDuplicates)
+        for (var pointIndex = 0; pointIndex < calcPoints.Count - 1; pointIndex++)
         {
-            var lineInPoly = new Line
-            {
-                BeginPoint = begin,
-                EndPoint = end
-            };
-            Point2D intersectionPoint = lineInPoly.GetIntersectPointXz(line);
-            if (intersectionPoint != null && (allowDuplicates || NoPointSameXzLocation(intersectionPointsWithLine, intersectionPoint)))
-            {
-                intersectionPointsWithLine.Add(intersectionPoint);
-            }
+            DoIntersectAndAddToCollection(line, calcPoints[pointIndex], calcPoints[pointIndex + 1],
+                                          intersectionPointsWithLine, allowDuplicates);
         }
 
-        private static bool NoPointSameXzLocation(IEnumerable<Point2D> collection, Point2D point)
+        return intersectionPointsWithLine;
+    }
+
+    private static void DoIntersectAndAddToCollection(Line line, Point2D begin, Point2D end, ICollection<Point2D> intersectionPointsWithLine, bool allowDuplicates)
+    {
+        var lineInPoly = new Line
         {
-            return !collection.Any(
-                       p => Math.Abs(p.X - point.X) < GeometryConstants.Accuracy &&
-                            Math.Abs(p.Z - point.Z) < GeometryConstants.Accuracy);
+            BeginPoint = begin,
+            EndPoint = end
+        };
+        Point2D intersectionPoint = lineInPoly.GetIntersectPointXz(line);
+        if (intersectionPoint != null && (allowDuplicates || NoPointSameXzLocation(intersectionPointsWithLine, intersectionPoint)))
+        {
+            intersectionPointsWithLine.Add(intersectionPoint);
         }
+    }
 
-        private List<Point2D> IntersectWithPointsListCore(IList<Point2D> list, bool closePointString)
+    private static bool NoPointSameXzLocation(IEnumerable<Point2D> collection, Point2D point)
+    {
+        return !collection.Any(
+                   p => Math.Abs(p.X - point.X) < GeometryConstants.Accuracy &&
+                        Math.Abs(p.Z - point.Z) < GeometryConstants.Accuracy);
+    }
+
+    private List<Point2D> IntersectWithPointsListCore(IList<Point2D> list, bool closePointString)
+    {
+        var result = new List<Point2D>();
+        var line = new Line();
+        for (var externalPointIndex = 0; externalPointIndex < list.Count - 1; externalPointIndex++)
         {
-            var result = new List<Point2D>();
-            var line = new Line();
-            for (var externalPointIndex = 0; externalPointIndex < list.Count - 1; externalPointIndex++)
-            {
-                line.BeginPoint = list[externalPointIndex];
-                line.EndPoint = list[externalPointIndex + 1];
+            line.BeginPoint = list[externalPointIndex];
+            line.EndPoint = list[externalPointIndex + 1];
 
-                result.AddRange(!closePointString
-                                    ? IntersectionPointsXzWithLineXzWithAllPoints(line)
-                                    : IntersectionPointsXzClosedStringWithLineXz(line));
-            }
-
-            return result;
+            result.AddRange(!closePointString
+                                ? IntersectionPointsXzWithLineXzWithAllPoints(line)
+                                : IntersectionPointsXzClosedStringWithLineXz(line));
         }
+
+        return result;
     }
 }
\ No newline at end of file