using System; using System.Collections.Generic; using System.Diagnostics; using Deltares.Mathematics; using Deltares.Standard; // ======================================================================================================================= // Class name: // // Description: // // Copyright (c) 2005-2009 Deltares // // Date ID Modification // 2008-03-31 Best Created // ======================================================================================================================= namespace Deltares.Stability.Calculation.Inner { public struct PolyPoint { public double X; public double Y; } // end PolyPoint } namespace Deltares.Stability.Calculation.Inner { public class MStabDatafunctions { public static double NormalDistribution(double Argm) { double result; const double Cnst = 2.506628274631; // sqrt(2*Pi) double Arg2; double ArgA; int Term; int Sign; result = 0; Arg2 = Argm*Argm; ArgA = Math.Abs(Argm); Term = 15; // supplies accuracy of 4.7(-10) in X=2.6 Sign = 2*(Term%2) - 1; if (ArgA > 2.6) { do { result = Term/(ArgA + result); Term -= 1; } while (!(Term == 0)); result = Math.Exp(-Arg2/2)/Cnst/(ArgA + result); if (Argm > 0) { result = 1 - result; } } else { do { result = Arg2/(2 + (1 + Sign*result)/Term); Term -= 1; Sign = -Sign; } while (!(Term == 0)); result = 0.5 + Math.Exp(-Arg2/2)/Cnst*Argm/(1 - result); } return result; } // ======================================================================================================================= // Method : NormalDistrInverse(Argm: Double) : Double; // // Member of class : // // Access : // // Description: Returns Inverse Gaussian Probability Function; // accuracy 7.2(-10) // // Pre-condition : Argument Range between 0 and 1 // // Post-condition : For Outcome �MaxDouble, check Argument // // Comment : Abramowitz and Stegun, 26.2.23; // Taylor Series correction with two terms; // Outcome is between -MaxDouble and +MaxDouble, // which practically means between -10 and +10; // MaxDouble = 1.7E+308; If this value is obtained, // it probably indicates an invalid argument range. // // Method History : // Date ID Modification // 1996-12-09 Sel Created // 1997-03-06 Sel MaxDouble added to indicate invalid argument range // 1999-04-14 Tuk Solved compiler warning by initializing Result. // ======================================================================================================================= public static double NormalDistrInverse(double Argm) { double result; double Temp; result = 0; if ((Argm <= 0) || (Argm >= 1)) { if (Argm <= 0) { result = Constants.CMaxDouble; } if (Argm >= 1) { result = Constants.CMaxDouble; } } else { if (Argm < 0.5) { Temp = Math.Sqrt(-2*Math.Log(Argm)); } else { Temp = Math.Sqrt(-2*Math.Log(1 - Argm)); } result = Temp - (2.515517 + Temp*(0.802853 + Temp*0.010328))/(1 + Temp*(1.432788 + Temp*(0.189269 + Temp*0.001308))); if (Argm < 0.5) { result = -result; } Temp = 2*result*(Argm - NormalDistribution(result)); Temp = 1 - Math.Sqrt(1 - Temp*Math.Sqrt(2*Math.PI)*Math.Exp(result*result/2)); result = result + Temp/result; } return result; } // ======================================================================================================================= // Description Get the Curve number of surface // // Date ID Modification // 2008-08-13 Best Created // ======================================================================================================================= public static int GetCurveNumberFromPoint(TPoints APoint, TPoints[] ASurface) { int result; int I; result = -1; for (I = ASurface.GetLowerBound(0) + 1; I <= ASurface.GetUpperBound(0); I ++) { if ((APoint.XCoor >= ASurface[I - 1].XCoor - 0.1*Constants.CGeoAccu) && (APoint.XCoor <= ASurface[I].XCoor + 0.1*Constants.CGeoAccu)) { if (((APoint.ZCoor >= Math.Min(ASurface[I - 1].ZCoor, ASurface[I].ZCoor) - 0.1*Constants.CGeoAccu) && (APoint.ZCoor <= Math.Max(ASurface[I - 1].ZCoor, ASurface[I].ZCoor) + 0.1*Constants.CGeoAccu))) { result = I; break; } } } return result; } // ======================================================================================================================= // Date ID Modification // 2010-05-17 Bkm Created // ======================================================================================================================= public static double GetSignedDistanceOnBoundary(TPoints AStartPoint, TPoints AEndPoint, TPoints[] ASurface) { double result; if (AStartPoint.XCoor < AEndPoint.XCoor) { result = GetDistanceOnBoundary(AStartPoint, AEndPoint, ASurface); } else { result = -GetDistanceOnBoundary(AStartPoint, AEndPoint, ASurface); } return result; } // ======================================================================================================================= // Date ID Modification // 2008-08-13 Best Created // ======================================================================================================================= public static double GetDistanceOnBoundary(TPoints AStartPoint, TPoints AEndPoint, TPoints[] ASurface) { double result; int i; TPoints LStartPoint; // pointer TPoints LEndPoint; // pointer int LStartCurveNumber; int LEndCurveNumber; bool LContinue; result = 0; LEndCurveNumber = 0; LStartPoint = AStartPoint; LEndPoint = AEndPoint; if ((Math.Abs(AStartPoint.XCoor - AEndPoint.XCoor) < Constants.CGeoAccu)) { // Point lies on a vertical line result = (AEndPoint.ZCoor - AStartPoint.ZCoor); } else { // both points are not on a vetical line if ((AStartPoint.XCoor > AEndPoint.XCoor)) { LStartPoint = AEndPoint; LEndPoint = AStartPoint; } LStartCurveNumber = GetCurveNumberFromPoint(AStartPoint, ASurface); LContinue = (LStartCurveNumber > 0); if (LContinue) { LEndCurveNumber = GetCurveNumberFromPoint(AEndPoint, ASurface); LContinue = ((LEndCurveNumber >= LStartCurveNumber) && (LEndCurveNumber <= ASurface.Length)); } if (LContinue) { if ((LStartCurveNumber == LEndCurveNumber)) { result = Distance2D(AStartPoint.XCoor, AStartPoint.ZCoor, AEndPoint.XCoor, AEndPoint.ZCoor); } else { // The first surface result = Distance2D(LStartPoint.XCoor, LStartPoint.ZCoor, ASurface[LStartCurveNumber].XCoor, ASurface[LStartCurveNumber].ZCoor); i = LStartCurveNumber + 1; // The surfaces in between while ((i < LEndCurveNumber)) { result = result + Distance2D(ASurface[i].XCoor, ASurface[i].ZCoor, ASurface[i - 1].XCoor, ASurface[i - 1].ZCoor); i ++; } // The last surface result = result + Distance2D(ASurface[LEndCurveNumber - 1].XCoor, ASurface[LEndCurveNumber - 1].ZCoor, LEndPoint.XCoor, LEndPoint.ZCoor); } } } return result; } // ======================================================================================================================= // Date ID Modification // 2005-04-18 Best Created // 2010-04-18 Can handle negative values for distance, then search to the left of the starting point // ======================================================================================================================= public static void GetPointOnBoundaryAtDistance(ref TPoints ANewPoint, TPoints AOldPoint, double ADistance, TPoints[] ASurface) { int i; int LStartCurveNumber; bool LContinue; double LSom; double LDistance; double LRatio; double LXStart; double LZStart; int LNCurves; LSom = 0; LNCurves = ASurface.Length; LStartCurveNumber = GetCurveNumberFromPoint(AOldPoint, ASurface); LContinue = (LStartCurveNumber > 0); i = LStartCurveNumber; LXStart = AOldPoint.XCoor; LZStart = AOldPoint.ZCoor; if (ADistance >= 0) { while (LContinue && (i < LNCurves)) { LDistance = Distance2D(LXStart, LZStart, ASurface[i].XCoor, ASurface[i].ZCoor); if ((LDistance > ADistance - LSom)) { // punt ligt op curve LRatio = (ADistance - LSom)/LDistance; ANewPoint.XCoor = LXStart + LRatio*(ASurface[i].XCoor - LXStart); ANewPoint.ZCoor = LZStart + LRatio*(ASurface[i].ZCoor - LZStart); LContinue = false; } else { LSom = LSom + LDistance; LXStart = ASurface[i].XCoor; LZStart = ASurface[i].ZCoor; i ++; } } } else { ADistance = -ADistance; while (LContinue && (i > 0)) { LDistance = Distance2D(LXStart, LZStart, ASurface[i - 1].XCoor, ASurface[i - 1].ZCoor); if ((LDistance > ADistance - LSom)) { // punt ligt op curve LRatio = (ADistance - LSom)/LDistance; ANewPoint.XCoor = LXStart + LRatio*(ASurface[i - 1].XCoor - LXStart); ANewPoint.ZCoor = LZStart + LRatio*(ASurface[i - 1].ZCoor - LZStart); LContinue = false; } else { LSom = LSom + LDistance; LXStart = ASurface[i - 1].XCoor; LZStart = ASurface[i - 1].ZCoor; i -= 1; } } } } // ======================================================================================================================= // Date ID Modification // 2008-04-14 Best Created // ======================================================================================================================= public static void SortIntersectionPoints(List AZIntersect, List AXIntersect) { // sort intersection points for (int i = 0; i < AXIntersect.Count; i++) { for (int j = i; j < AXIntersect.Count; j++) { if (AXIntersect[i] > AXIntersect[j]) { SwapDoubles(AXIntersect, i, j); SwapDoubles(AZIntersect, i, j); } } } } // ======================================================================================================================= // Date ID Modification // 2008-04-01 Best Created // ======================================================================================================================= public static void RemovePointFromArray(int APointNo, ref double[] AX, ref double[] AY) { int LnDim; int i; LnDim = AX.Length; // Copy the values from higher points for (i = APointNo; i <= LnDim - 2; i ++) { AX[i] = AX[i + 1]; AY[i] = AY[i + 1]; } // Remove redundant point AX = new double[LnDim - 1]; AY = new double[LnDim - 1]; } // ======================================================================================================================= // Date ID Modification // 2008-04-01 Best Created // ======================================================================================================================= public static void RemoveDoubleIntersectionPoints(List AXIntersect, List AYIntersect) { var indices = new List(); // Put the indexes with double values in LDoubleCount for (int i = 1; i < AXIntersect.Count; i++) { if ((Distance2D(AXIntersect[i - 1], AYIntersect[i - 1], AXIntersect[i], AYIntersect[i]) < Constants.CSameSliceDist)) { indices.Add(i); } } // Remove the double points for (int i = indices.Count - 1; i >= 0; i--) { AXIntersect.RemoveAt(indices[i]); AYIntersect.RemoveAt(indices[i]); } } // ======================================================================================================================= // // Description : Checks to see if two given lines intersect, and // returns coordinates of intersections if they do. // // Name Type Function // ---- ---- -------- // Parameters - IN : X1,Y1 Double First point on first line // X2,Y2 Double Second point on first line // X3,Y3 Double First point on second line // X4,Y4 Double Second point on second line // Parameters - OUT : XS,YS Double Coordinates of intersection // if one was found. // // Pre-condition : None // // Post-condition : If Result of function is TRUE, then Xs, Ys // are the coordinates of the intersection point. // IF function result is FALSE, then Xs and Ys // are both 0.0. // // Comments : The function checks if the lines intersect // within the given begin and end points. If the // lines intersect further on (which will always // be the case unless both lines are parallel), // no intersection point is found. // // Example call : If IntersectLines(0,0,2,2, 0,2,2,0, Xs, Ys) then // etc. etc. // 2008-08-07 Best // ======================================================================================================================= public static bool IntersectLines(double x1, double y1, double x2, double y2, double x3, double y3, double x4, double y4, ref double xs, ref double ys) { bool result; double dx1; double dy1; double dx2; double dy2; double noem; double u; double v; // Method : Use parameter equations for both lines : // line 1 : X = x1 + u*(x2 - x1) = x1 + u*dx1 // Y = y1 + u*(y2 - y1) = y1 + u*dy1 and // line 2 : X = x3 + v*(x4 - x3) = x3 + v*dx2 // Y = y3 + v*(y4 - y3) = y3 + v*dy2 // Combining the equations for X and Y leaves 2 equations with // unknowns u and v. Solving for u, an intersection is found if // u can be sloved, and 0 <= u <= 1. Substituting u in equations // for line 1 gives the intersection point. If no value for u can // be calculated (due to denominator = 0), the lines are parrallel // dx2(y1-y3) - dy2(x1-x3) // u = ------------------------- // dx1*dy2 - dx2*dy1 // If 0 <= u <= 1, the intersection lies within line 1. A second // check needs to be made to ensure that the intersection point also // lies on line 2 by checking whether 0 <= v <= 1. // v = (x1 + u*dx1 - x3) / dx2 result = false; xs = 0; ys = 0; dx1 = x2 - x1; dy1 = y2 - y1; dx2 = x4 - x3; dy2 = y4 - y3; noem = dx1*dy2 - dy1*dx2; if (Math.Abs(noem) > 1E-10) { u = (dx2*(y1 - y3) - dy2*(x1 - x3))/noem; if ((u >= 0.0) && (u <= 1.0)) { v = (dx1*(y1 - y3) - dy1*(x1 - x3))/noem; if ((v >= 0.0) && (v <= 1.0)) { result = true; xs = x1 + u*dx1; ys = y1 + u*dy1; } } } return result; } // ======================================================================================================================= // // Procedure : Intersect_Circle_line(Xm, Ym, R, X1, Y1, X2, Y2, // ns, xs, ys) // // Declared in unit : BREKEN (Internal) // // Description : Voegt aan de array Snijp de X-koord. toe van de // snijpunten van de cirkel bepaald door Xm, Ym en // R met het lijnstuk bepaald door de punten X1,Y1 // en X2,Y2. // // Last Update 941114 : (Dek) Test for 0 <= U <= 1 changed in test for // -eps <= U <= 1+eps because in some special // cases an intersection was NOT found, while there // should be one. Occured when circle intersects // at precisely a geometry point. // 960411 : (Eg) Ook Y-coordinaat snijpunt toegevoegd. // Nieuwe naam en naar GDMATH // 970410 : (Best) Naar DUtils // 970930 : (Bka/Best) Results for ys1 and ys2 corrected to ys[1] // and ys[2] instead of xs[1] and xs[2] // // Name Type Function // ---- ---- -------- // Parameters - IN : Xm,Ym Real Circle midpoint // R Real Radius // X1,Y1 Real 1st point on a line // X2,Y2 Real 2nd point on a line // - OUT : Ns Integer Number of intersection // points // Xs Arr2 X-co. of intersection points. // Ys Arr2 Y-co. of intersection points. // // Pre-condition : None // // Post-condition : None // // Comments : An Accuracy of 1.0E-8 is used to determine the // the number of intersection points. // The line may consist only One point! // // Example call : Intersect_Circle_line(Xm,Ym,R,X1,Y1,X2,Y2,N,Xs,Ys) // ======================================================================================================================= public static void Intersect_Circle_line(double Xm, double Ym, double R, double x1, double x2, double y1, double y2, ref int Ns, ref double xs1, ref double ys1, ref double xs2, ref double ys2) { double ABC_Eps = 1E-8; double dx1; double dx2; double dy1; double dy2; double A; double B; double C; double D; double U; var xs = new double[2 + 1]; var ys = new double[2 + 1]; Ns = 0; // Oplossing door invullen parametrische vergelijking van lijn : // X = x1 + u*(x2 - x1) // Y = y1 + u*(y2 - y1) // in vergelijking van cirkel : // (X - xm)^2 + (Y - ym)^2 = r^2 // en vervolgens parameter u oplossen uit vierkantsverg. // - geen oplossing => geen snijpunt // - een of twee oplossingen => snijpunt als -� <= u <= 1+�. // Als (x1,y1) an (x2,y2) alleen een punt weergeven, doe niets. dx1 = x2 - x1; dx2 = x1 - Xm; dy1 = y2 - y1; dy2 = y1 - Ym; if ((Math.Abs(dx1) > ABC_Eps) || (Math.Abs(dy1) > ABC_Eps)) { // There is a real line to intersect with the circle. A = dx1*dx1 + dy1*dy1; B = 2*(dx1*dx2 + dy1*dy2); C = dx2*dx2 + dy2*dy2 - R*R; D = B*B - 4*A*C; if (D > ABC_Eps) { // 2 snijpunten U = (-B + Math.Sqrt(D))/(2*A); if ((U >= -ABC_Eps) && (U <= (1.0 + ABC_Eps))) { Ns ++; xs[Ns] = x1 + U*dx1; ys[Ns] = y1 + U*dy1; } U = (-B - Math.Sqrt(D))/(2*A); if ((U >= -ABC_Eps) && (U <= (1.0 + ABC_Eps))) { Ns ++; xs[Ns] = x1 + U*dx1; ys[Ns] = y1 + U*dy1; } } else if (Math.Abs(D) <= ABC_Eps) { // 1 raakpunt U = -B/(2*A); if ((U >= -ABC_Eps) && (U <= 1.0 + ABC_Eps)) { Ns ++; xs[Ns] = x1 + U*dx1; ys[Ns] = y1 + U*dy1; } } } // Line is one point; Check distance else if (Math.Abs(Math.Pow(dx2, 2) + Math.Pow(dy2, 2) - Math.Pow(R, 2)) < ABC_Eps) { Ns ++; xs[Ns] = x1; ys[Ns] = y1; } xs1 = xs[1]; ys1 = ys[1]; xs2 = xs[2]; ys2 = ys[2]; } // ======================================================================================================================= // Date ID Modification // 2008-04-10 Best Created // ======================================================================================================================= public static void SnijCirkelLijnstuk(double Xm, double Zm, double Rad, double Xl, double Zl, double Xr, double Zr, List AXIntersect, List AZIntersect) { double xs1 = 0; double xs2 = 0; double Zs1 = 0; double Zs2 = 0; int N = 0; if (((Xr < (Xm - Rad)) || (Xl > (Xm + Rad)) || ((Zl < (Zm - Rad)) && (Zr < (Zm - Rad))))) { N = 0; } else { Intersect_Circle_line(Xm, Zm, Rad, Xl, Xr, Zl, Zr, ref N, ref xs1, ref Zs1, ref xs2, ref Zs2); if (N > 0) { AXIntersect.Add(xs1); AZIntersect.Add(Zs1); } if (N > 1) { AXIntersect.Add(xs2); AZIntersect.Add(Zs2); } } } // ======================================================================================================================= // Date ID Modification // 2008-03-31 Best Copied from geolib // ======================================================================================================================= public static double Distance2D(double X1, double Y1, double X2, double Y2) { double result; double dx; double dy; dx = X1 - X2; dy = Y1 - Y2; result = Math.Sqrt(dx*dx + dy*dy); return result; } // ======================================================================================================================= // Date ID Modification // 2008-03-31 Best Copied from geolib // ======================================================================================================================= public static double DistanceToLine(double XS, double YS, double X1, double Y1, double X2, double Y2) { double result; double LAC; double LAB; double LBC; double LAClose; // afstand van A tot hoogtelijn op lijnstuk A - C // A, B, C is een driehoek LAB = Distance2D(X1, Y1, X2, Y2); LAC = Distance2D(X1, Y1, XS, YS); LBC = Distance2D(X2, Y2, XS, YS); if ((LAB != 0)) { // bereken hoogtelijn LAClose = (Math.Pow(LAC, 2) - Math.Pow(LBC, 2) + Math.Pow(LAB, 2))/(2*LAB); if (LAClose <= 0) { result = LAC; } else { if (LAClose >= LAB) { result = LBC; } else { result = Math.Sqrt(Math.Abs(Math.Pow(LAC, 2) - Math.Pow(LAClose, 2))); } } } else { result = LAC; } return result; } // ======================================================================================================================= // Date ID Modification // 2008-04-02 Best Created // ======================================================================================================================= public static double LinInpolY(double Xs, double Ys, double Xe, double Ye, double Xi) { double result; if (Math.Abs(Xe - Xs) < 1E-15) { result = (Ye + Ys)*0.5; } else { result = Ys + (Ye - Ys)*((Xi - Xs)/(Xe - Xs)); } return result; } // ======================================================================================================================= // Date ID Modification // 2008-04-02 Best Created // ======================================================================================================================= public static double ZTopAtSurface(double Ax, TPoints[] ASurface) { double result; int i; bool LFound; double x1; double Z1; double X2; double Z2; double LZMax; LZMax = -Constants.CMaxCoor; LFound = false; if (!(ASurface != null)) { LFound = false; } // check if xcoordinate is not outside surface reach x1 = ASurface[ASurface.GetLowerBound(0)].XCoor; X2 = ASurface[ASurface.GetUpperBound(0)].XCoor; // if points are somewhere on the surface if ((Ax >= x1) && (Ax <= X2)) { // If points are on a point or on a vertical line find them for (i = ASurface.GetLowerBound(0); i <= ASurface.GetUpperBound(0); i ++) { x1 = ASurface[i].XCoor; if ((Math.Abs(Ax - x1) < Constants.CGeoAccu)) { LFound = true; LZMax = Math.Max(LZMax, ASurface[i].ZCoor); } } // X is somwhere between the Points if (!LFound) { x1 = ASurface[ASurface.GetLowerBound(0)].XCoor; Z1 = ASurface[ASurface.GetLowerBound(0)].ZCoor; for (i = ASurface.GetLowerBound(0) + 1; i <= ASurface.GetUpperBound(0); i ++) { X2 = ASurface[i].XCoor; Z2 = ASurface[i].ZCoor; if ((Ax >= x1) && (Ax <= X2)) { LFound = true; LZMax = LinInpolY(x1, Z1, X2, Z2, Ax); break; } x1 = X2; Z1 = Z2; } } // Must must be somewhere Debug.Assert(LFound); } result = LZMax; return result; } // ======================================================================================================================= // Date ID Modification // 2008-04-02 Best Created // ======================================================================================================================= public static double ZBottomAtSurface(double Ax, TPoints[] ASurface) { double result; int i; bool LFound; double x1; double Z1; double X2; double Z2; double LZMin; LZMin = Constants.CMaxCoor; LFound = false; // check if xcoordinate is not outside surface reach x1 = ASurface[ASurface.GetLowerBound(0)].XCoor; X2 = ASurface[ASurface.GetUpperBound(0)].XCoor; // if points are somewhere on the surface if ((Ax >= x1) && (Ax <= X2)) { // If points are on a point or on a vertical line find them for (i = ASurface.GetLowerBound(0); i <= ASurface.GetUpperBound(0); i ++) { x1 = ASurface[i].XCoor; if ((Math.Abs(Ax - x1) < Constants.CGeoAccu)) { LFound = true; LZMin = Math.Min(LZMin, ASurface[i].ZCoor); } } // X is somwhere between the Points if (!LFound) { x1 = ASurface[ASurface.GetLowerBound(0)].XCoor; Z1 = ASurface[ASurface.GetLowerBound(0)].ZCoor; for (i = ASurface.GetLowerBound(0) + 1; i <= ASurface.GetUpperBound(0); i ++) { X2 = ASurface[i].XCoor; Z2 = ASurface[i].ZCoor; if ((Ax >= x1) && (Ax <= X2)) { LFound = true; LZMin = LinInpolY(x1, Z1, X2, Z2, Ax); break; } x1 = X2; Z1 = Z2; } } // Must must be somewhere Debug.Assert(LFound); } result = LZMin; return result; } // ======================================================================================================================= // Description Try to find a vertical line on a boundary // Date ID Modification // 2008-04-03 Best Created // ======================================================================================================================= public static void FindVertical(TPoints[] ASurface, ref bool IsVerticalFound, ref double AX, ref double AZ1, ref double AZ2) { double x1; double z1; double x2; double z2; int i; IsVerticalFound = false; for (i = ASurface.GetLowerBound(0); i < ASurface.GetUpperBound(0); i ++) { x1 = ASurface[i].XCoor; z1 = ASurface[i].ZCoor; x2 = ASurface[i + 1].XCoor; z2 = ASurface[i + 1].ZCoor; if ((Math.Abs(x1 - x2) < Constants.CGeoEps)) { if (IsVerticalFound) { // is it the same vertical? if ((Math.Abs(AX - x2) < Constants.CGeoEps)) { // yes, add the last z AZ2 = z2; } } else { // not yet found, this is the first initialise IsVerticalFound = true; AX = x1; AZ1 = z1; AZ2 = z2; } } } } // ======================================================================================================================= // Date ID Modification // 2008-04-10 Best Created // ======================================================================================================================= public static double ZAtXOnCircle(double Ax, double Axmid, double AZmid, double ARadius) { double result; double Dx; double Dy2; try { Dx = Axmid - Ax; Dy2 = ARadius*ARadius - Dx*Dx; if ((Dy2 <= 1.0E-6)) { result = AZmid; } else { result = AZmid - Math.Sqrt(Dy2); } } catch { result = AZmid; } return result; } // ======================================================================================================================= // Date ID Modification // 2008-04-10 Best Created // ======================================================================================================================= public static void SwapDoubles(ref double ADouble1, ref double ADouble2) { double LTmp; LTmp = ADouble1; ADouble1 = ADouble2; ADouble2 = LTmp; } public static void SwapDoubles(List list, int index1, int index2) { double temp = list[index1]; list[index1] = list[index2]; list[index2] = temp; } public static bool InPolygon(List Layerpoints, double AX, double AY) { var polygon = new List(); for (int i = 0; i < Layerpoints.Count; i++) { var Point = new Point2D(Layerpoints[i].XCoor, Layerpoints[i].ZCoor); polygon.Add(Point); } return ((Routines2D.CheckIfPointIsInPolygon(polygon, AX, AY) != PointInPolygon.OutsidePolygon)); } // ======================================================================================================================= // Date ID Modification // 2008-08-07 Best Created // ======================================================================================================================= public static bool InPolygon(PolyPoint[] APolynome, double AX, double AY) { bool result; const double eps = 1.0E-20; int i; bool LReady; double LDeltaX; double LDeltaY; double LAbsX; double LAbsY; double LAL1; double LAL2; double LAL3; double LSom; int LNPolyPoints; result = false; LNPolyPoints = APolynome.GetUpperBound(0) + 1; if ((LNPolyPoints > 2)) { LDeltaX = APolynome[0].X - AX; LDeltaY = APolynome[0].Y - AY; LAbsX = Math.Abs(LDeltaX); LAbsY = Math.Abs(LDeltaY); if ((LAbsX < eps) && (LAbsY < eps)) { result = true; } else { LAL1 = Math.Atan2(LDeltaY, LDeltaX); LSom = 0; LReady = false; i = 1; while ((i < LNPolyPoints) && !LReady) { LDeltaX = APolynome[i].X - AX; LDeltaY = APolynome[i].Y - AY; LAbsX = Math.Abs(LDeltaX); LAbsY = Math.Abs(LDeltaY); if ((LAbsX < eps) && (LAbsY < eps)) { result = true; LReady = true; } LAL2 = Math.Atan2(LDeltaY, LDeltaX); LAL3 = LAL2 - LAL1; // Normeer lAL3 tussen - PI en + PI if ((LAL3 < -Math.PI)) { LAL3 = LAL3 + 2.0*Math.PI; } if ((LAL3 > Math.PI)) { LAL3 = LAL3 - 2.0*Math.PI; } LSom = LSom + LAL3; LAL1 = LAL2; i ++; } if (!LReady) { if ((LSom > 1.9*Math.PI) || (LSom < -1.9*Math.PI)) { result = true; } else { result = false; } } } } return result; } // ======================================================================================================================= // Date ID Modification // 2009-04-09 Created // ======================================================================================================================= public static double DetermineArea(List Points) { double som = 0; if (Points.Count > 2) { for (int i = 1; i < Points.Count; i++) { som = som + (Points[i - 1].XCoor - Points[i].XCoor)*(Points[i - 1].ZCoor + Points[i].ZCoor)/2; } } return som; } // ======================================================================================================================= // Date ID Modification // 2009-04-09 Created // ======================================================================================================================= public static int FindLayerNumberFromPoint(TLayerRecord[] ASoilLayers, double AXCoor, double AZCoor) { int result = 0; var LXIntersect = new List(); var LZIntersect = new List(); var LNNumbers = new List(); result = -1; // Find the intersection points with all layers for (int i = 1; i <= ASoilLayers.Length; i++) { if (((ASoilLayers[i - 1].MinX <= AXCoor) && (ASoilLayers[i - 1].MaxX >= AXCoor) && (ASoilLayers[i - 1].MinZ <= AZCoor) && (ASoilLayers[i - 1].MaxZ >= AZCoor))) { // (* There is a posibility that point is in this layer} // see if there are intersections with the horizontal to the right of x // if this number is odd the point lies in the area , if the number is zero or even // the point lies outside the area *) double LXStart = AXCoor; double LZStart = AZCoor; double LXEnd = ASoilLayers[i - 1].MaxX + 10; double LZEnd = AZCoor; if (InPolygon(ASoilLayers[i - 1].Layerpoints, AXCoor, AZCoor)) { // Point Lies in Layer i LNNumbers.Add(i - 1); } } } // If nothing found, point lies outside geometry // If one point found that is the layer number if (LNNumbers.Count == 1) { result = LNNumbers[0]; } // if there are subareas, the point lies in the layer // with the smaller area if (LNNumbers.Count > 1) { double LMinArea = ASoilLayers[LNNumbers[0]].Area; result = LNNumbers[0]; for (int i = 1; i < LNNumbers.Count; i++) { if (LMinArea > ASoilLayers[LNNumbers[i]].Area) { LMinArea = ASoilLayers[LNNumbers[i]].Area; result = LNNumbers[i]; } } } return result; } /// /// Determines all x coordinates from a surface at a certain level /// /// /// /// public static List GetAllXAtSurface(TPoints[] SurfaceLine, double ZCoor) { int i; double LXStart; double LZStart; double LXEnd; double LZEnd; double lXs; var result = new List(); for (i = SurfaceLine.GetLowerBound(0) + 1; i <= SurfaceLine.GetUpperBound(0); i++) { if (((ZCoor >= SurfaceLine[i - 1].ZCoor) && (ZCoor <= SurfaceLine[i].ZCoor)) || ((ZCoor <= SurfaceLine[i - 1].ZCoor) && (ZCoor >= SurfaceLine[i].ZCoor))) { LXStart = SurfaceLine[i - 1].XCoor; LZStart = SurfaceLine[i - 1].ZCoor; LXEnd = SurfaceLine[i].XCoor; LZEnd = SurfaceLine[i].ZCoor; if (sameValue(LZStart, ZCoor) && sameValue(LZEnd, ZCoor)) { lXs = LXEnd; } else { lXs = LinInpolY(LZStart, LXStart, LZEnd, LXEnd, ZCoor); } bool LFound = false; for (int j = 0; j < result.Count; j++) { if (sameValue(result[j], lXs)) { LFound = true; break; } } if (!LFound) { result.Add(lXs); } } } return result; } // ======================================================================================================================= // Date ID Modification // 2009-04-09 Created // ======================================================================================================================= public static double GetZCoordinateSurfaceAtX(TPoints[] SurfaceLine, double AXCoor) { double result; int i; double LXStart; double LZStart; double LXEnd; double LZEnd; result = 0; for (i = SurfaceLine.GetLowerBound(0) + 1; i <= SurfaceLine.GetUpperBound(0); i ++) { if (((AXCoor >= SurfaceLine[i - 1].XCoor) && (AXCoor <= SurfaceLine[i].XCoor)) || ((AXCoor <= SurfaceLine[i - 1].XCoor) && (AXCoor >= SurfaceLine[i].XCoor))) { LXStart = SurfaceLine[i - 1].XCoor; LZStart = SurfaceLine[i - 1].ZCoor; LXEnd = SurfaceLine[i].XCoor; LZEnd = SurfaceLine[i].ZCoor; result = LinInpolY(LXStart, LZStart, LXEnd, LZEnd, AXCoor); break; } } return result; } // ======================================================================================================================= // Date ID Modification // 2009-04-15 Created // ======================================================================================================================= public static void GetnextpointOnSurface(TPoints[] SurfaceLine, int Asign, ref double AXCoor, ref double AZCoor) { int i; double LXStart; double LZStart; double LXEnd; double LZEnd; bool LFirst; bool LFound; bool LFinish = false; LFound = false; LXStart = AXCoor; LZStart = AZCoor; LXEnd = AXCoor; LZEnd = AZCoor; for (i = SurfaceLine.GetLowerBound(0) + 1; i <= SurfaceLine.GetUpperBound(0); i ++) { if (((AXCoor >= SurfaceLine[i - 1].XCoor) && (AXCoor < SurfaceLine[i].XCoor)) || ((AXCoor <= SurfaceLine[i - 1].XCoor) && (AXCoor > SurfaceLine[i].XCoor))) { LXStart = SurfaceLine[i - 1].XCoor; LZStart = SurfaceLine[i - 1].ZCoor; LXEnd = SurfaceLine[i].XCoor; LZEnd = SurfaceLine[i].ZCoor; // first check of one of the points is the start x and z coor} if ((Asign > 0)) { if (sameValue(LXStart, AXCoor) && sameValue(LZStart, AZCoor)) { AXCoor = LXEnd; AZCoor = LZEnd; LFinish = true; break; } } else { if (sameValue(LXEnd, AXCoor) && sameValue(LZEnd, AZCoor)) { AXCoor = LXStart; AZCoor = LZStart; LFinish = true; break; } } LFound = true; break; } } if ((LFound) && (!LFinish)) { if ((Asign < 0)) { LFirst = (LXStart <= LXEnd); } else { LFirst = (LXStart >= LXEnd); } if (LFirst) { AXCoor = LXStart; AZCoor = LZStart; } else { AXCoor = LXEnd; AZCoor = LZEnd; } } } // ======================================================================================================================= // Date ID Modification // 2009-04-09 Created // ======================================================================================================================= public static void IntersectLineSoileArea(TLayerRecord ASoilLayer, double XStart, double ZStart, double XEnd, double ZEnd, List XIntersect, List ZIntersect) { double LXS = 0; double LZS = 0; double Lx1 = 0; double LZ1 = 0; double LX2 = 0; double LZ2 = 0; // ASoilLayer is a closed polygon for (int i = 1; i < ASoilLayer.NumberOfLayerPoints; i ++) { Lx1 = ASoilLayer.Layerpoints[i - 1].XCoor; LZ1 = ASoilLayer.Layerpoints[i - 1].ZCoor; LX2 = ASoilLayer.Layerpoints[i].XCoor; LZ2 = ASoilLayer.Layerpoints[i].ZCoor; // If an intersection point is found add these to the array if (IntersectLines(Lx1, LZ1, LX2, LZ2, XStart, ZStart, XEnd, ZEnd, ref LXS, ref LZS)) { if ((!(sameValue(Lx1, LXS) && sameValue(LZ1, LZS))) || (!(sameValue(LX2, LXS) && sameValue(LZ2, LZS)))) { addIntersectionpointsIfNotAlreadyExcists(XIntersect, ZIntersect, LXS, LZS); } } } // the last line if a soil area is not close Lx1 = ASoilLayer.Layerpoints[ASoilLayer.NumberOfLayerPoints - 1].XCoor; LZ1 = ASoilLayer.Layerpoints[ASoilLayer.NumberOfLayerPoints - 1].ZCoor; LX2 = ASoilLayer.Layerpoints[0].XCoor; LZ2 = ASoilLayer.Layerpoints[0].ZCoor; // If an intersection point is found add these to the array if (IntersectLines(Lx1, LZ1, LX2, LZ2, XStart, ZStart, XEnd, ZEnd, ref LXS, ref LZS)) { addIntersectionpointsIfNotAlreadyExcists(XIntersect, ZIntersect, LXS, LZS); } } // ======================================================================================================================= // Date ID Modification // 2009-04-09 Created // ======================================================================================================================= public static void IntersectLineSoileAreas(TLayerRecord[] ASoilLayers, double XStart, double ZStart, double XEnd, double ZEnd, List XIntersect, List ZIntersect) { int LNlayers = 0; var LXIntersect = new List(); var LZIntersect = new List(); LNlayers = ASoilLayers.Length; // Fins the intersection points with all layers for (int i = 1; i <= LNlayers; i ++) { IntersectLineSoileArea(ASoilLayers[i - 1], XStart, ZStart, XEnd, ZEnd, LXIntersect, LZIntersect); if ((LXIntersect.Count > 0)) { for (int j = 0; j < LXIntersect.Count; j++) { addIntersectionpointsIfNotAlreadyExcists(XIntersect, ZIntersect, LXIntersect[j], LZIntersect[j]); } } } // Sort the intersection points SortIntersectionPoints(ZIntersect, XIntersect); } // ======================================================================================================================= // // Date ID Modification // 2010-04-15 Created // ======================================================================================================================= public static double FindIntersectionsWithPLLine(double AXCoor, TPoints[] APLLines) { double result; int i; result = -999.0; for (i = APLLines.GetLowerBound(0); i < APLLines.GetUpperBound(0); i ++) { if ((AXCoor >= APLLines[i].XCoor) && (AXCoor <= APLLines[i + 1].XCoor)) { result = LinInpolY(APLLines[i].XCoor, APLLines[i].ZCoor, APLLines[i + 1].XCoor, APLLines[i + 1].ZCoor, AXCoor); } } return result; } // ======================================================================================================================= // Description Find intersection from waternet lines with vertical. // Note that a waterline needs not go from left to right // Date ID Modification // 2010-04-15 Created // ======================================================================================================================= public static void FindIntersectionsWithWateline(double AXCoor, TPoints[] AWaterline, ref int ANIntersections, ref double[] ZCoors) { int i; double LMinX; double LMaxX; ANIntersections = 0; ZCoors = null; for (i = AWaterline.GetLowerBound(0); i < AWaterline.GetUpperBound(0); i ++) { LMinX = Math.Min(AWaterline[i].XCoor, AWaterline[i + 1].XCoor); LMaxX = Math.Max(AWaterline[i].XCoor, AWaterline[i + 1].XCoor); if ((AXCoor >= LMinX) && (AXCoor <= LMaxX)) { ANIntersections ++; ZCoors = new double[ANIntersections]; ZCoors[ANIntersections - 1] = LinInpolY(AWaterline[i].XCoor, AWaterline[i].ZCoor, AWaterline[i + 1].XCoor, AWaterline[i + 1].ZCoor, AXCoor); } } } private static bool sameValue(double aVal1, double aVal2) { return (Math.Abs(aVal1 - aVal2) < Constants.CAlmostZero); } private static void addIntersectionpointsIfNotAlreadyExcists(List XIntersect, List ZIntersect, double AXIntersect, Double AZIntersect) { bool lfound = false; for (int i = 0; i < XIntersect.Count; i++) { if (sameValue(XIntersect[i], AXIntersect) && sameValue(ZIntersect[i], AZIntersect)) { lfound = true; break; } } if (!lfound) { XIntersect.Add(AXIntersect); ZIntersect.Add(AZIntersect); } } } // end MStabDatafunctions }