using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.Linq;
using Deltares.Geometry;
using Deltares.Geotechnics;
using Deltares.Geotechnics.Soils;
// =======================================================================================================================
// Class name:
//
// Description:
//
// Copyright (c) 2005-2009 Deltares
//
// Date ID Modification
// 2008-03-31 Best Created
// =======================================================================================================================
namespace Deltares.Stability.Calculation2
{
public struct PolyPoint
{
public double X;
public double Y;
}
// end PolyPoint
}
namespace Deltares.Stability.Calculation2
{
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, GeometryPointString ASurface)
{
int result;
int I;
result = -1;
for (I = 1; I < ASurface.Points.Count; I ++)
{
if ((APoint.XCoor >= ASurface.Points[I - 1].X - 0.1*Constants.CGeoAccu) && (APoint.XCoor <= ASurface.Points[I].X + 0.1*Constants.CGeoAccu))
{
if (((APoint.ZCoor >= Math.Min(ASurface.Points[I - 1].Z, ASurface.Points[I].Z) - 0.1*Constants.CGeoAccu) && (APoint.ZCoor <= Math.Max(ASurface.Points[I - 1].Z, ASurface.Points[I].Z) + 0.1*Constants.CGeoAccu)))
{
result = I;
break;
}
}
}
return result;
}
public static double GetSignedDistanceOnBoundary(GeometryPoint startPoint, GeometryPoint endPoint, GeometryPointString surfaceLine)
{
var start = new TPoints(startPoint.X, startPoint.Z);
var end = new TPoints(endPoint.X, endPoint.Z);
return GetSignedDistanceOnBoundary(start, end, surfaceLine);
}
///
/// get signed distance on boundary
///
///
///
///
///
public static double GetSignedDistanceOnBoundary(TPoints startPoint, TPoints endPoint, GeometryPointString surfaceLine)
{
double result;
result = startPoint.XCoor < endPoint.XCoor
? GetDistanceOnBoundary(startPoint, endPoint, surfaceLine)
: -GetDistanceOnBoundary(startPoint, endPoint, surfaceLine);
return result;
}
// =======================================================================================================================
// Date ID Modification
// 2008-08-13 Best Created
// =======================================================================================================================
public static double GetDistanceOnBoundary(TPoints AStartPoint, TPoints AEndPoint, GeometryPointString 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(LStartPoint, ASurface);
LContinue = (LStartCurveNumber > 0);
if (LContinue)
{
LEndCurveNumber = GetCurveNumberFromPoint(LEndPoint, ASurface);
LContinue = ((LEndCurveNumber >= LStartCurveNumber) && (LEndCurveNumber <= ASurface.Points.Count));
}
if (LContinue)
{
if ((LStartCurveNumber == LEndCurveNumber))
{
result = Distance2D(LStartPoint.XCoor, LStartPoint.ZCoor, LEndPoint.XCoor, LEndPoint.ZCoor);
}
else
{
// The first surface
result = Distance2D(LStartPoint.XCoor, LStartPoint.ZCoor, ASurface.Points[LStartCurveNumber].X, ASurface.Points[LStartCurveNumber].Z);
i = LStartCurveNumber + 1;
// The surfaces in between
while ((i < LEndCurveNumber))
{
result = result + Distance2D(ASurface.Points[i].X, ASurface.Points[i].Z, ASurface.Points[i - 1].X, ASurface.Points[i - 1].Z);
i ++;
}
// The last surface
result = result + Distance2D(ASurface.Points[LEndCurveNumber - 1].X, ASurface.Points[LEndCurveNumber - 1].Z, LEndPoint.XCoor, LEndPoint.ZCoor);
}
}
}
return result;
}
public static void GetPointOnBoundaryAtDistance(ref GeometryPoint ANewPoint, GeometryPoint AOldPoint, double ADistance, GeometryPointString ASurface)
{
var newPoint = new TPoints(ANewPoint.X, ANewPoint.Z);
var oldPoint = new TPoints(AOldPoint.X, AOldPoint.Z);
GetPointOnBoundaryAtDistance(ref newPoint, oldPoint, ADistance, ASurface);
ANewPoint.X = newPoint.XCoor;
ANewPoint.Z = newPoint.ZCoor;
}
// =======================================================================================================================
// 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, GeometryPointString ASurface)
{
int i;
int LStartCurveNumber;
bool LContinue;
double LSom;
double LDistance;
double LRatio;
double LXStart;
double LZStart;
int LNCurves;
LSom = 0;
LNCurves = ASurface.Points.Count;
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.Points[i].X, ASurface.Points[i].Z);
if ((LDistance > ADistance - LSom))
{
// punt ligt op curve
LRatio = (ADistance - LSom)/LDistance;
ANewPoint.XCoor = LXStart + LRatio*(ASurface.Points[i].X - LXStart);
ANewPoint.ZCoor = LZStart + LRatio*(ASurface.Points[i].Z - LZStart);
LContinue = false;
}
else
{
LSom = LSom + LDistance;
LXStart = ASurface.Points[i].X;
LZStart = ASurface.Points[i].Z;
i ++;
}
}
}
else
{
ADistance = -ADistance;
while (LContinue && (i > 0))
{
LDistance = Distance2D(LXStart, LZStart, ASurface.Points[i - 1].X, ASurface.Points[i - 1].Z);
if ((LDistance > ADistance - LSom))
{
// punt ligt op curve
LRatio = (ADistance - LSom)/LDistance;
ANewPoint.XCoor = LXStart + LRatio*(ASurface.Points[i - 1].X - LXStart);
ANewPoint.ZCoor = LZStart + LRatio*(ASurface.Points[i - 1].Z - LZStart);
LContinue = false;
}
else
{
LSom = LSom + LDistance;
LXStart = ASurface.Points[i - 1].X;
LZStart = ASurface.Points[i - 1].Z;
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];
}
///
/// get intersection points with circle and line containing 2 points
///
///
///
///
public static GeometryPointString GetIntersectionPointsCircleAndline(Circle circle, GeometryPointString pointString)
{
double ABC_Eps = 1E-8;
double A;
double B;
double C;
double D;
double U;
var xs = new double[2 + 1];
var zs = new double[2 + 1];
int 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.
double dx1 = pointString.Points[1].X - pointString.Points[0].X;
double dx2 = pointString.Points[0].X - circle.Xcentre;
double dz1 = pointString.Points[1].Z - pointString.Points[0].Z;
double dz2 = pointString.Points[0].Z - circle.Zcentre;
if ((Math.Abs(dx1) > ABC_Eps) || (Math.Abs(dz1) > ABC_Eps))
{
// There is a real line to intersect with the circle.
A = dx1*dx1 + dz1*dz1;
B = 2*(dx1*dx2 + dz1*dz2);
C = dx2*dx2 + dz2*dz2 - circle.Radius*circle.Radius;
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] = pointString.Points[0].X + U*dx1;
zs[Ns] = pointString.Points[0].Z + U*dz1;
}
U = (-B - Math.Sqrt(D))/(2*A);
if ((U >= -ABC_Eps) && (U <= (1.0 + ABC_Eps)))
{
Ns++;
xs[Ns] = pointString.Points[0].X + U*dx1;
zs[Ns] = pointString.Points[0].Z + U*dz1;
}
}
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] = pointString.Points[0].X + U*dx1;
zs[Ns] = pointString.Points[0].Z + U*dz1;
}
}
}
// Line is one point; Check distance
else if (Math.Abs(Math.Pow(dx2, 2) + Math.Pow(dz2, 2) - Math.Pow(circle.Radius, 2)) < ABC_Eps)
{
Ns++;
xs[Ns] = pointString.Points[0].X;
zs[Ns] = pointString.Points[0].Z;
}
var gString = new GeometryPointString();
gString.Points.Add(new GeometryPoint(xs[1], 0, zs[1]));
gString.Points.Add(new GeometryPoint(xs[2], 0, zs[2]));
return gString;
}
// =======================================================================================================================
// 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;
}
///
/// get top or bottom at the surfaceline
///
///
///
///
///
public static double GetZTopOrBottomAtSurface(double Ax, GeometryPointString surfaceLine, bool top)
{
double result;
int i;
bool LFound;
double x1;
double Z1;
double X2;
double Z2;
double maxMin = -Constants.CMaxCoor;
LFound = false;
if (surfaceLine == null)
{
return maxMin;
}
// check if xcoordinate is not outside surface reach
x1 = surfaceLine.Points.First().X;
X2 = surfaceLine.Points.Last().X;
// 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 = 0; i < surfaceLine.Points.Count; i++)
{
x1 = surfaceLine.Points[i].X;
if ((Math.Abs(Ax - x1) < Constants.CGeoAccu))
{
LFound = true;
if (top)
{
maxMin = Math.Max(maxMin, surfaceLine.Points[i].Z);
}
else
{
maxMin = Math.Min(maxMin, surfaceLine.Points[i].Z);
}
}
}
// X is somwhere between the Points
if (!LFound)
{
x1 = surfaceLine.Points.First().X;
Z1 = surfaceLine.Points.First().Z;
for (i = 1; i < surfaceLine.Points.Count; i++)
{
X2 = surfaceLine.Points[i].X;
Z2 = surfaceLine.Points[i].Z;
if ((Ax >= x1) && (Ax <= X2))
{
LFound = true;
maxMin = LinInpolY(x1, Z1, X2, Z2, Ax);
break;
}
x1 = X2;
Z1 = Z2;
}
}
// must be somewhere
Debug.Assert(LFound);
}
result = maxMin;
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
// 2013 Dijk, L. updated
// =======================================================================================================================
public static Tuple FindVertical(GeometryPointString surfaceLine)
{
bool isVerticalFound = false;
double AX = 0;
double AZ1 = 0;
double AZ2 = 0;
int i;
for (i = 0; i < surfaceLine.Points.Count - 1; i++)
{
double x1 = surfaceLine.Points[i].X;
double z1 = surfaceLine.Points[i].Z;
double x2 = surfaceLine.Points[i + 1].X;
double z2 = surfaceLine.Points[i + 1].Z;
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;
}
}
}
return new Tuple(isVerticalFound, AX, AZ1, AZ2);
}
///
/// find vertical line
///
///
///
public static GeometryPointString FindVerticalLine(GeometryPointString surfaceLine)
{
var vertical = new GeometryPointString();
int i;
for (i = 0; i < surfaceLine.Points.Count - 1; i++)
{
double x1 = surfaceLine.Points[i].X;
double z1 = surfaceLine.Points[i].Z;
double x2 = surfaceLine.Points[i + 1].X;
double z2 = surfaceLine.Points[i + 1].Z;
if ((Math.Abs(x1 - x2) < Constants.CGeoEps))
{
vertical.Points.Add(new GeometryPoint(x1, 0, z1));
vertical.Points.Add(new GeometryPoint(x1, 0, z2));
break;
}
}
return vertical;
}
// =======================================================================================================================
// 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;
}
// =======================================================================================================================
// 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(IList 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
// =======================================================================================================================
//@@new: Layers
public static int FindLayerNumberFromPoint(List 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.Count; i++)
{
var minX = ASoilLayers[i - 1].GeometrySurface.OuterLoop.GetMinX();
var maxX = ASoilLayers[i - 1].GeometrySurface.OuterLoop.GetMaxX();
var minZ = ASoilLayers[i - 1].GeometrySurface.OuterLoop.GetMinZ();
var maxZ = ASoilLayers[i - 1].GeometrySurface.OuterLoop.GetMaxZ();
if (((minX <= AXCoor) && (maxX >= AXCoor) && (minZ <= AZCoor) && (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 = maxX + 10;
double LZEnd = AZCoor;
IntersectLineSoileArea(ASoilLayers[i - 1], LXStart, LZStart, LXEnd, LZEnd, LXIntersect, LZIntersect);
//@ Unsupported function or procedure: 'Odd'
if (LXIntersect.Count%2 == 1)
{
// 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)
{
//@@new: Layers
double LMinArea = Math.Abs(DetermineArea(ASoilLayers[LNNumbers[0]].GeometrySurface.OuterLoop.Points));
result = LNNumbers[0];
for (int i = 1; i < LNNumbers.Count; i++)
{
if (LMinArea > Math.Abs(DetermineArea(ASoilLayers[LNNumbers[i]].GeometrySurface.OuterLoop.Points)))
{
LMinArea = Math.Abs(DetermineArea(ASoilLayers[LNNumbers[i]].GeometrySurface.OuterLoop.Points));
result = LNNumbers[i];
}
}
}
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;
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;
LFound = true;
break;
}
}
if (LFound)
{
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
// =======================================================================================================================
//@@new: Layers
public static void IntersectLineSoileArea(SoilLayer2D ASoilLayer, double XStart, double ZStart, double XEnd, double ZEnd, List XIntersect, List ZIntersect)
{
double LXS = 0;
double LZS = 0;
// ASoilLayer is a closed polygon
for (int i = 1; i < ASoilLayer.GeometrySurface.OuterLoop.Points.Count; i ++)
{
double Lx1 = ASoilLayer.GeometrySurface.OuterLoop.Points[i - 1].X;
double LZ1 = ASoilLayer.GeometrySurface.OuterLoop.Points[i - 1].Z;
double LX2 = ASoilLayer.GeometrySurface.OuterLoop.Points[i].X;
double LZ2 = ASoilLayer.GeometrySurface.OuterLoop.Points[i].Z;
// 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))
{
XIntersect.Add(LXS);
ZIntersect.Add(LZS);
}
}
}
// =======================================================================================================================
// Date ID Modification
// 2009-04-09 Created
// =======================================================================================================================
//@@new: Layers
public static void IntersectLineSoileAreas(List 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.Count;
// 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++)
{
XIntersect.Add(LXIntersect[j]);
ZIntersect.Add(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);
}
}
}
///
/// get unique intersection points of circle with surfaceline
///
///
///
///
///
///
public static List GetUniqueIntersectionPoints(double xCenterPoint, double zCenterPoint, double radius, GeometryPointString surfaceLine)
{
// find intersection points
var intersectionPoints = new List();
for (int i = 1; i < surfaceLine.Points.Count; i++)
{
var circle = new Circle()
{
Xcentre = xCenterPoint,
Zcentre = zCenterPoint,
Radius = radius,
};
var gString = new GeometryPointString();
gString.Points.Add(new GeometryPoint(surfaceLine.Points[i - 1].X, 0, surfaceLine.Points[i - 1].Z));
gString.Points.Add(new GeometryPoint(surfaceLine.Points[i].X, 0, surfaceLine.Points[i].Z));
var result = GetIntersectionPointLineSegmentSlipcircle(circle, surfaceLine);
foreach (var geometryPoint in result)
{
intersectionPoints.Add(geometryPoint);
}
}
return intersectionPoints;
}
///
/// get intersection points of the line segment and the slip circle
///
///
///
///
internal static List GetIntersectionPointLineSegmentSlipcircle(Circle circle, GeometryPointString gString)
{
var xCenterPoint = circle.Xcentre;
var zCenterPoint = circle.Zcentre;
var radius = circle.Radius;
double lxStart = gString.Points[0].X;
double lzStart = gString.Points[0].Z;
double lxEnd = gString.Points[1].X;
double lzEnd = gString.Points[1].Z;
var gpoints = new List();
if (((lxEnd < (xCenterPoint - radius)) || (lxStart > (xCenterPoint + radius)) || ((lzStart < (zCenterPoint - radius)) && (lzEnd < (zCenterPoint - radius)))))
{
return gpoints;
}
else
{
var intersections = GetIntersectionPointsCircleAndline(circle, gString);
if (intersections.Points.Count > 1)
{
gpoints.Add(intersections.Points[1]);
gpoints.Add(intersections.Points[0]);
}
else if (intersections.Points.Count > 0)
{
gpoints.Add(intersections.Points[0]);
}
return gpoints;
}
}
private static double DetermineArea(IList geometryPoints)
{
double som = 0;
if (geometryPoints.Count > 2)
{
for (int i = 1; i < geometryPoints.Count; i++)
{
som = som + (geometryPoints[i - 1].X - geometryPoints[i].X)*(geometryPoints[i - 1].Z + geometryPoints[i].Z)/2;
}
}
return som;
}
}
// end MStabDatafunctions
}