Index: DamEngine/trunk/src/Deltares.DamEngine.Data/Geotechnics/SoilProfile2D.cs
===================================================================
diff -u -r6404 -r6557
--- DamEngine/trunk/src/Deltares.DamEngine.Data/Geotechnics/SoilProfile2D.cs	(.../SoilProfile2D.cs)	(revision 6404)
+++ DamEngine/trunk/src/Deltares.DamEngine.Data/Geotechnics/SoilProfile2D.cs	(.../SoilProfile2D.cs)	(revision 6557)
@@ -34,7 +34,6 @@
 /// </summary>   
 public class SoilProfile2D : SoilProfile
 {
-    private const double deviation = 1E-05;
     private const double toleranceAlmostEqual = 1e-07;
     private readonly List<SoilLayer2D> surfaces = new List<SoilLayer2D>();
 
@@ -163,33 +162,35 @@
     public static SoilLayer2D DetermineOriginalLayerFromOldSurfaces(GeometrySurface geometrySurface,
                                                                     IEnumerable<SoilLayer2D> oldSurfaces)
     {
-        var point = new Point2D(0.0, 0.0);
-        var isPointInOuterLoopAndOldSurface = false;
-        point = IsPointInOuterLoopAndOldSurface(geometrySurface, oldSurfaces, point, ref isPointInOuterLoopAndOldSurface);
-
-        if (!isPointInOuterLoopAndOldSurface)
+        Point2D point = geometrySurface.DetermineValidTestPointBasedOnNewSurface();
+        // If the test point is beyond the x-limits of the old surfaces, there cannot be an old layer.
+        if (IsSurfaceLeftOrRightOfOldSurfaces(point, oldSurfaces))
         {
-            isPointInOuterLoopAndOldSurface = IsPointInOldSurfaceJustBelowTopOfNewGeometryWithinItsLimits(geometrySurface,
-                                                                                                          oldSurfaces, point,
-                                                                                                          isPointInOuterLoopAndOldSurface);
+            return null;
         }
+        
+        return DetermineOldSurfaceFromTestPoint(oldSurfaces, point);
+    }
 
-        if (isPointInOuterLoopAndOldSurface)
+    /// <summary>
+    /// Determines whether the given point is left or right of the old surfaces.
+    /// </summary>
+    /// <param name="point"></param>
+    /// <param name="oldSurfaces"></param>
+    /// <returns>true when the given point is left or right of the old surfaces else false.</returns>
+    public static bool IsSurfaceLeftOrRightOfOldSurfaces(Point2D point, IEnumerable<SoilLayer2D> oldSurfaces)
+    {
+        // If the test point is beyound the x-limits of the old surfaces, there can not be an old layer
+        double xminFromSurfaces = DetermineXminFromSurfaces(oldSurfaces);
+        double xmaxFromSurfaces = DetermineXmaxFromSurfaces(oldSurfaces);
+        if (point.X > xmaxFromSurfaces || point.X < xminFromSurfaces)
         {
-            if (IsPointInPreviousOuterLoopOfOldSurface(oldSurfaces, point, out SoilLayer2D soilLayer2D))
-            {
-                return soilLayer2D;
-            }
-
-            if (IsPointInOuterLoopOfOldSurface(oldSurfaces, point, out SoilLayer2D originalLayerFromOldSurfaces1))
-            {
-                return originalLayerFromOldSurfaces1;
-            }
+            return true;
         }
 
-        return null;
+        return false;
     }
-
+    
     /// <summary>
     /// Returns a <see cref="System.String" /> that represents this instance.
     /// </summary>
@@ -201,17 +202,22 @@
         return Name;
     }
 
-    private static bool IsPointInOuterLoopOfOldSurface(IEnumerable<SoilLayer2D> oldSurfaces, Point2D point,
-                                                       out SoilLayer2D originalLayerFromOldSurfaces1)
+    /// <summary>
+    /// Determine the old surface from the test point.
+    /// </summary>
+    /// <param name="oldSurfaces"></param>
+    /// <param name="point"></param>
+    /// <returns>The found old surface/layer else null.</returns>
+    public static SoilLayer2D DetermineOldSurfaceFromTestPoint(IEnumerable<SoilLayer2D> oldSurfaces, Point2D point)
     {
-        originalLayerFromOldSurfaces1 = null;
         foreach (SoilLayer2D oldSurface in oldSurfaces)
         {
             GeometryLoop outerLoop = oldSurface.GeometrySurface.OuterLoop;
-            if (outerLoop != null && outerLoop.CurveList.Count > 2 && Routines2D.CheckIfPointIsInPolygon(outerLoop, point.X,
-                                                                                                         point.Z) == PointInPolygon.InsidePolygon)
+            if (outerLoop != null && outerLoop.CurveList.Count > 2 && 
+                Routines2D.CheckIfPointIsInPolygon(outerLoop, point.X, point.Z) == PointInPolygon.InsidePolygon)
             {
                 var isPointInOuterLoopOfOldSurface = true;
+                // Make sure point is NOT in (one of) the inner loop(s). 
                 foreach (GeometryLoop innerLoop in oldSurface.GeometrySurface.InnerLoops)
                 {
                     if (Routines2D.CheckIfPointIsInPolygon(innerLoop, point.X, point.Z) == PointInPolygon.InsidePolygon)
@@ -222,117 +228,13 @@
 
                 if (isPointInOuterLoopOfOldSurface)
                 {
-                    originalLayerFromOldSurfaces1 = oldSurface;
-                    return true;
+                    return oldSurface;
                 }
             }
         }
-
-        return false;
+        return null;
     }
-
-    private static bool IsPointInPreviousOuterLoopOfOldSurface(IEnumerable<SoilLayer2D> oldSurfaces, Point2D point,
-                                                               out SoilLayer2D soilLayer2D)
-    {
-        soilLayer2D = null;
-        foreach (SoilLayer2D oldSurface in oldSurfaces)
-        {
-            GeometryLoop previousOuterLoop = oldSurface.GeometrySurface.PreviousOuterLoop;
-            if (previousOuterLoop != null && previousOuterLoop.CurveList.Count > 2 &&
-                Routines2D.CheckIfPointIsInPolygon(previousOuterLoop, point.X, point.Z) == PointInPolygon.InsidePolygon)
-            {
-                var isPointInPreviousOuterLoopOfOldSurface = true;
-                foreach (GeometryLoop previousInnerLoop in oldSurface.GeometrySurface.PreviousInnerLoops)
-                {
-                    if (Routines2D.CheckIfPointIsInPolygon(previousInnerLoop, point.X, point.Z) == PointInPolygon.InsidePolygon)
-                    {
-                        isPointInPreviousOuterLoopOfOldSurface = false;
-                    }
-                }
-
-                if (isPointInPreviousOuterLoopOfOldSurface)
-                {
-                    soilLayer2D = oldSurface;
-                    return true;
-                }
-            }
-        }
-
-        return false;
-    }
-
-    private static bool IsPointInOldSurfaceJustBelowTopOfNewGeometryWithinItsLimits(GeometrySurface geometrySurface,
-                                                                                    IEnumerable<SoilLayer2D> oldSurfaces, Point2D point, bool isPointInOuterLoopAndOldSurface)
-    {
-        GeometryPointString topGeometrySurface = geometrySurface.DetermineTopGeometrySurface();
-        topGeometrySurface.SortPointsByXAscending();
-        Point2D geometryPoint1 = topGeometrySurface[0];
-        geometryPoint1.X -= deviation;
-        geometryPoint1.Z -= deviation;
-        Point2D geometryPoint2 = topGeometrySurface[checked(topGeometrySurface.Count - 1)];
-        geometryPoint2.X += deviation;
-        geometryPoint2.Z -= deviation;
-        bool isPoint1WithinOldSurfaces = IsPointWithinOldSurfaces(geometryPoint1, oldSurfaces, -deviation);
-        bool isPoint2WithinOldSurfaces = IsPointWithinOldSurfaces(geometryPoint2, oldSurfaces, -deviation);
-        double d = double.NaN;
-        if (isPoint1WithinOldSurfaces && !isPoint2WithinOldSurfaces)
-        {
-            point.X = geometryPoint1.X;
-            point.Z = geometryPoint1.Z;
-            isPointInOuterLoopAndOldSurface = true;
-            d = geometryPoint1.X + deviation;
-        }
-
-        if (!isPoint1WithinOldSurfaces && isPoint2WithinOldSurfaces)
-        {
-            point.X = geometryPoint2.X;
-            point.Z = geometryPoint2.Z;
-            isPointInOuterLoopAndOldSurface = true;
-            d = geometryPoint2.X - deviation;
-        }
-
-        if (!double.IsNaN(d))
-        {
-            double xminFromSurfaces = DetermineXminFromSurfaces(oldSurfaces);
-            double xmaxFromSurfaces = DetermineXmaxFromSurfaces(oldSurfaces);
-            if (d <= xmaxFromSurfaces && d >= xminFromSurfaces)
-            {
-                isPointInOuterLoopAndOldSurface = false;
-            }
-        }
-
-        return isPointInOuterLoopAndOldSurface;
-    }
-
-    private static Point2D IsPointInOuterLoopAndOldSurface(GeometrySurface geometrySurface, IEnumerable<SoilLayer2D> oldSurfaces, Point2D point, ref bool isPointInOuterLoopAndOldSurface)
-    {
-        foreach (GeometryCurve curve in geometrySurface.OuterLoop.CurveList)
-        {
-            point = new Point2D((curve.HeadPoint.X + curve.EndPoint.X) / 2.0, (curve.HeadPoint.Z + curve.EndPoint.Z) / 2.0);
-            if (IsPointWithinOldSurfaces(point, oldSurfaces, deviation) || IsPointWithinOldSurfaces(point, oldSurfaces,
-                                                                                                    -deviation))
-            {
-                point.Z += deviation;
-                if (Routines2D.CheckIfPointIsInPolygon(geometrySurface.OuterLoop, point.X, point.Z) ==
-                    PointInPolygon.InsidePolygon)
-                {
-                    isPointInOuterLoopAndOldSurface = true;
-                    break;
-                }
-
-                point.Z -= 2 * deviation;
-                if (Routines2D.CheckIfPointIsInPolygon(geometrySurface.OuterLoop, point.X, point.Z) ==
-                    PointInPolygon.InsidePolygon)
-                {
-                    isPointInOuterLoopAndOldSurface = true;
-                    break;
-                }
-            }
-        }
-
-        return point;
-    }
-
+    
     private static double DetermineXminFromSurfaces(IEnumerable<SoilLayer2D> oldSurfaces)
     {
         var xminFromSurfaces = double.MaxValue;
@@ -355,51 +257,6 @@
         return xmaxFromSurfaces;
     }
 
-    private static bool IsPointWithinOldSurfaces(Point2D point, IEnumerable<SoilLayer2D> oldSurfaces, double verticalShift)
-    {
-        point.Z += verticalShift;
-        var shiftedPoint = new Point2D(point.X, point.Z);
-        foreach (SoilLayer2D oldSurface in oldSurfaces)
-        {
-            GeometryLoop outerLoop = oldSurface.GeometrySurface.OuterLoop;
-            List<GeometryLoop> innerLoops = oldSurface.GeometrySurface.InnerLoops;
-            bool isPointInSurface = IsPointInGivenOuterLoopOfOldSurface(shiftedPoint, outerLoop, innerLoops);
-            if (!isPointInSurface)
-            {
-                GeometryLoop previousOuterLoop = oldSurface.GeometrySurface.PreviousOuterLoop;
-                List<GeometryLoop> previousInnerLoops = oldSurface.GeometrySurface.PreviousInnerLoops;
-                isPointInSurface = IsPointInGivenOuterLoopOfOldSurface(shiftedPoint, previousOuterLoop, previousInnerLoops);
-            }
-
-            if (isPointInSurface)
-            {
-                return true;
-            }
-        }
-
-        return false;
-    }
-
-    private static bool IsPointInGivenOuterLoopOfOldSurface(Point2D point, GeometryLoop outerLoop,
-                                                            List<GeometryLoop> innerLoops)
-    {
-        var isPointInSurface = false;
-        if (outerLoop != null && Routines2D.CheckIfPointIsInPolygon(outerLoop, point.X, point.Z) == PointInPolygon.InsidePolygon)
-        {
-            isPointInSurface = true;
-            // Make sure the point is in the outer loop, not in the inner loop(s).
-            foreach (GeometryLoop innerLoop in innerLoops)
-            {
-                if (Routines2D.CheckIfPointIsInPolygon(innerLoop, point.X, point.Z) == PointInPolygon.InsidePolygon)
-                {
-                    return false;
-                }
-            }
-        }
-
-        return isPointInSurface;
-    }
-
     private SoilProfile1D DetermineSoilProfile1DAtX(double x)
     {
         var soilProfile = new SoilProfile1D