in lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/Plane.java [1749:2278]
public void recordBounds(
final PlanetModel planetModel, final XYZBounds boundsInfo, final Membership... bounds) {
// Basic plan is to do three intersections of the plane and the planet.
// For min/max x, we intersect a vertical plane such that y = 0.
// For min/max y, we intersect a vertical plane such that x = 0.
// For min/max z, we intersect a vertical plane that is chosen to go through the high point of
// the arc. For clarity, load local variables with good names
final double A = this.x;
final double B = this.y;
final double C = this.z;
// Do Z. This can be done simply because it is symmetrical.
if (!boundsInfo.isSmallestMinZ(planetModel) || !boundsInfo.isLargestMaxZ(planetModel)) {
// System.err.println(" computing Z bound");
// Compute Z bounds for this arc
// With ellipsoids, we really have only one viable way to do this computation.
// Specifically, we compute an appropriate vertical plane, based on the current plane's x-y
// orientation, and then intersect it with this one and with the ellipsoid. This gives us
// zero, one, or two points to use as bounds.
// There is one special case: horizontal circles. These require TWO vertical planes: one for
// the x, and one for the y, and we use all four resulting points in the bounds computation.
if ((Math.abs(A) >= MINIMUM_RESOLUTION || Math.abs(B) >= MINIMUM_RESOLUTION)) {
// NOT a degenerate case
// System.err.println(" not degenerate");
final Plane normalizedZPlane = constructNormalizedZPlane(A, B);
final GeoPoint[] points =
findIntersections(planetModel, normalizedZPlane, bounds, NO_BOUNDS);
for (final GeoPoint point : points) {
assert planetModel.pointOnSurface(point);
// System.err.println(" Point = "+point+";
// this.evaluate(point)="+this.evaluate(point)+";
// normalizedZPlane.evaluate(point)="+normalizedZPlane.evaluate(point));
addPoint(boundsInfo, bounds, point);
}
} else {
// Since a==b==0, any plane including the Z axis suffices.
// System.err.println(" Perpendicular to z");
GeoPoint[] points = findIntersections(planetModel, normalYPlane, NO_BOUNDS, NO_BOUNDS);
if (points.length == 0) {
points = findIntersections(planetModel, normalXPlane, NO_BOUNDS, NO_BOUNDS);
}
if (points.length == 0) {
boundsInfo.addZValue(new GeoPoint(0.0, 0.0, -this.z));
} else {
boundsInfo.addZValue(points[0]);
}
}
}
// First, compute common subexpressions
final double k =
1.0
/ ((x * x + y * y) * planetModel.xyScaling * planetModel.xyScaling
+ z * z * planetModel.zScaling * planetModel.zScaling);
final double abSquared = planetModel.xyScaling * planetModel.xyScaling;
final double cSquared = planetModel.zScaling * planetModel.zScaling;
final double ASquared = A * A;
final double BSquared = B * B;
final double CSquared = C * C;
final double r = 2.0 * D * k;
final double rSquared = r * r;
if (!boundsInfo.isSmallestMinX(planetModel) || !boundsInfo.isLargestMaxX(planetModel)) {
// For min/max x, we need to use lagrange multipliers.
//
// For this, we need grad(F(x,y,z)) = (dF/dx, dF/dy, dF/dz).
//
// Minimize and maximize f(x,y,z) = x, with respect to g(x,y,z) = Ax + By + Cz - D and
// h(x,y,z) = x^2/xyScaling^2 + y^2/xyScaling^2 + z^2/zScaling^2 - 1
//
// grad(f(x,y,z)) = (1,0,0)
// grad(g(x,y,z)) = (A,B,C)
// grad(h(x,y,z)) = (2x/xyScaling^2,2y/xyScaling^2,2z/zScaling^2)
//
// Equations we need to simultaneously solve:
//
// grad(f(x,y,z)) = l * grad(g(x,y,z)) + m * grad(h(x,y,z))
// g(x,y,z) = 0
// h(x,y,z) = 0
//
// Equations:
// 1 = l*A + m*2x/xyScaling^2
// 0 = l*B + m*2y/xyScaling^2
// 0 = l*C + m*2z/zScaling^2
// Ax + By + Cz + D = 0
// x^2/xyScaling^2 + y^2/xyScaling^2 + z^2/zScaling^2 - 1 = 0
//
// Solve for x,y,z in terms of (l, m):
//
// x = ((1 - l*A) * xyScaling^2 ) / (2 * m)
// y = (-l*B * xyScaling^2) / ( 2 * m)
// z = (-l*C * zScaling^2)/ (2 * m)
//
// Two equations, two unknowns:
//
// A * (((1 - l*A) * xyScaling^2 ) / (2 * m)) + B * ((-l*B * xyScaling^2) / ( 2 * m)) + C *
// ((-l*C * zScaling^2)/ (2 * m)) + D = 0
//
// and
//
// (((1 - l*A) * xyScaling^2 ) / (2 * m))^2/xyScaling^2 + ((-l*B * xyScaling^2) / ( 2 *
// m))^2/xyScaling^2 + ((-l*C * zScaling^2)/ (2 * m))^2/zScaling^2 - 1 = 0
//
// Simple: solve for l and m, then find x from it.
//
// (a) Use first equation to find l in terms of m.
//
// A * (((1 - l*A) * xyScaling^2 ) / (2 * m)) + B * ((-l*B * xyScaling^2) / ( 2 * m)) + C *
// ((-l*C * zScaling^2)/ (2 * m)) + D = 0
// A * ((1 - l*A) * xyScaling^2 ) + B * (-l*B * xyScaling^2) + C * (-l*C * zScaling^2) + D * 2
// * m = 0
// A * xyScaling^2 - l*A^2* xyScaling^2 - B^2 * l * xyScaling^2 - C^2 * l * zScaling^2 + D * 2
// * m = 0
// - l *(A^2* xyScaling^2 + B^2 * xyScaling^2 + C^2 * zScaling^2) + (A * xyScaling^2 + D * 2 *
// m) = 0
// l = (A * xyScaling^2 + D * 2 * m) / (A^2* xyScaling^2 + B^2 * xyScaling^2 + C^2 *
// zScaling^2)
// l = A * xyScaling^2 / (A^2* xyScaling^2 + B^2 * xyScaling^2 + C^2 * zScaling^2) + m * 2 * D
// / (A^2* xyScaling^2 + B^2 * xyScaling^2 + C^2 * zScaling^2)
//
// For convenience:
//
// k = 1.0 / (A^2* xyScaling^2 + B^2 * xyScaling^2 + C^2 * zScaling^2)
//
// Then:
//
// l = A * xyScaling^2 * k + m * 2 * D * k
// l = k * (A*xyScaling^2 + m*2*D)
//
// For further convenience:
//
// q = A*xyScaling^2*k
// r = 2*D*k
//
// l = (r*m + q)
// l^2 = (r^2 * m^2 + 2*r*m*q + q^2)
//
// (b) Simplify the second equation before substitution
//
// (((1 - l*A) * xyScaling^2 ) / (2 * m))^2/xyScaling^2 + ((-l*B * xyScaling^2) / ( 2 *
// m))^2/xyScaling^2 + ((-l*C * zScaling^2)/ (2 * m))^2/zScaling^2 - 1 = 0
// ((1 - l*A) * xyScaling^2 )^2/xyScaling^2 + (-l*B * xyScaling^2)^2/xyScaling^2 + (-l*C *
// zScaling^2)^2/zScaling^2 = 4 * m^2
// (1 - l*A)^2 * xyScaling^2 + (-l*B)^2 * xyScaling^2 + (-l*C)^2 * zScaling^2 = 4 * m^2
// (1 - 2*l*A + l^2*A^2) * xyScaling^2 + l^2*B^2 * xyScaling^2 + l^2*C^2 * zScaling^2 = 4 *
// m^2
// xyScaling^2 - 2*A*xyScaling^2*l + A^2*xyScaling^2*l^2 + B^2*xyScaling^2*l^2 +
// C^2*zScaling^2*l^2 - 4*m^2 = 0
//
// (zScaling) Substitute for l, l^2
//
// xyScaling^2 - 2*A*xyScaling^2*(r*m + q) + A^2*xyScaling^2*(r^2 * m^2 + 2*r*m*q + q^2) +
// B^2*xyScaling^2*(r^2 * m^2 + 2*r*m*q + q^2) + C^2*zScaling^2*(r^2 * m^2 + 2*r*m*q + q^2) -
// 4*m^2 = 0
// xyScaling^2 - 2*A*xyScaling^2*r*m - 2*A*xyScaling^2*q + A^2*xyScaling^2*r^2*m^2 +
// 2*A^2*xyScaling^2*r*q*m +
// A^2*xyScaling^2*q^2 + B^2*xyScaling^2*r^2*m^2 + 2*B^2*xyScaling^2*r*q*m +
// B^2*xyScaling^2*q^2 + C^2*zScaling^2*r^2*m^2 + 2*C^2*zScaling^2*r*q*m + C^2*zScaling^2*q^2
// - 4*m^2 = 0
//
// (d) Group
//
// m^2 * [A^2*xyScaling^2*r^2 + B^2*xyScaling^2*r^2 + C^2*zScaling^2*r^2 - 4] +
// m * [- 2*A*xyScaling^2*r + 2*A^2*xyScaling^2*r*q + 2*B^2*xyScaling^2*r*q +
// 2*C^2*zScaling^2*r*q] +
// [xyScaling^2 - 2*A*xyScaling^2*q + A^2*xyScaling^2*q^2 + B^2*xyScaling^2*q^2 +
// C^2*zScaling^2*q^2] = 0
// Useful subexpressions for this bound
final double q = A * abSquared * k;
final double qSquared = q * q;
// Quadratic equation
final double a =
ASquared * abSquared * rSquared
+ BSquared * abSquared * rSquared
+ CSquared * cSquared * rSquared
- 4.0;
final double b =
-2.0 * A * abSquared * r
+ 2.0 * ASquared * abSquared * r * q
+ 2.0 * BSquared * abSquared * r * q
+ 2.0 * CSquared * cSquared * r * q;
final double c =
abSquared
- 2.0 * A * abSquared * q
+ ASquared * abSquared * qSquared
+ BSquared * abSquared * qSquared
+ CSquared * cSquared * qSquared;
if (Math.abs(a) >= MINIMUM_RESOLUTION_SQUARED) {
final double sqrtTerm = b * b - 4.0 * a * c;
if (Math.abs(sqrtTerm) < MINIMUM_RESOLUTION_SQUARED) {
// One solution
final double m = -b / (2.0 * a);
// Valid?
if (Math.abs(m) >= MINIMUM_RESOLUTION) {
final double l = r * m + q;
// x = ((1 - l*A) * xyScaling^2 ) / (2 * m)
// y = (-l*B * xyScaling^2) / ( 2 * m)
// z = (-l*C * zScaling^2)/ (2 * m)
final double denom0 = 0.5 / m;
final GeoPoint thePoint =
new GeoPoint(
(1.0 - l * A) * abSquared * denom0,
-l * B * abSquared * denom0,
-l * C * cSquared * denom0);
// Math is not quite accurate enough for this
// assert planetModel.pointOnSurface(thePoint): "Point: "+thePoint+";
// Planetmodel="+planetModel+"; A="+A+" B="+B+" C="+C+" D="+D+" planetfcn="+
// (thePoint.x*thePoint.x*planetModel.inverseXYScaling*planetModel.inverseXYScaling +
// thePoint.y*thePoint.y*planetModel.inverseXYScaling*planetModel.inverseXYScaling +
// thePoint.z*thePoint.z*planetModel.inverseZScaling*planetModel.inverseZScaling);
// assert evaluateIsZero(thePoint): "Evaluation of point: "+evaluate(thePoint);
addPoint(boundsInfo, bounds, thePoint);
} else {
// This is a plane of the form A=n B=0 C=0. We can set a bound only by noting the D
// value.
boundsInfo.addXValue(-D / A);
}
} else if (sqrtTerm > 0.0) {
// Two solutions
final double sqrtResult = Math.sqrt(sqrtTerm);
final double commonDenom = 0.5 / a;
final double m1 = (-b + sqrtResult) * commonDenom;
assert Math.abs(a * m1 * m1 + b * m1 + c) < MINIMUM_RESOLUTION;
final double m2 = (-b - sqrtResult) * commonDenom;
assert Math.abs(a * m2 * m2 + b * m2 + c) < MINIMUM_RESOLUTION;
if (Math.abs(m1) >= MINIMUM_RESOLUTION || Math.abs(m2) >= MINIMUM_RESOLUTION) {
final double l1 = r * m1 + q;
final double l2 = r * m2 + q;
// x = ((1 - l*A) * xyScaling^2 ) / (2 * m)
// y = (-l*B * xyScaling^2) / ( 2 * m)
// z = (-l*C * zScaling^2)/ (2 * m)
final double denom1 = 0.5 / m1;
final double denom2 = 0.5 / m2;
final GeoPoint thePoint1 =
new GeoPoint(
(1.0 - l1 * A) * abSquared * denom1,
-l1 * B * abSquared * denom1,
-l1 * C * cSquared * denom1);
final GeoPoint thePoint2 =
new GeoPoint(
(1.0 - l2 * A) * abSquared * denom2,
-l2 * B * abSquared * denom2,
-l2 * C * cSquared * denom2);
// Math is not quite accurate enough for this
// assert planetModel.pointOnSurface(thePoint1): "Point1: "+thePoint1+";
// Planetmodel="+planetModel+"; A="+A+" B="+B+" C="+C+" D="+D+" planetfcn="+
// (thePoint1.x*thePoint1.x*planetModel.inverseXYScaling*planetModel.inverseXYScaling +
// thePoint1.y*thePoint1.y*planetModel.inverseXYScaling*planetModel.inverseXYScaling +
// thePoint1.z*thePoint1.z*planetModel.inverseZScaling*planetModel.inverseZScaling);
// assert planetModel.pointOnSurface(thePoint2): "Point1: "+thePoint2+";
// Planetmodel="+planetModel+"; A="+A+" B="+B+" C="+C+" D="+D+" planetfcn="+
// (thePoint2.x*thePoint2.x*planetModel.inverseXYScaling*planetModel.inverseXYScaling +
// thePoint2.y*thePoint2.y*planetModel.inverseXYScaling*planetModel.inverseXYScaling +
// thePoint2.z*thePoint2.z*planetModel.inverseZScaling*planetModel.inverseZScaling);
// assert evaluateIsZero(thePoint1): "Evaluation of point1: "+evaluate(thePoint1);
// assert evaluateIsZero(thePoint2): "Evaluation of point2: "+evaluate(thePoint2);
addPoint(boundsInfo, bounds, thePoint1);
addPoint(boundsInfo, bounds, thePoint2);
} else {
// This is a plane of the form A=n B=0 C=0. We can set a bound only by noting the D
// value.
boundsInfo.addXValue(-D / A);
}
} else {
// No solutions
}
} else if (Math.abs(b) > MINIMUM_RESOLUTION_SQUARED) {
// a = 0, so m = - zScaling / b
final double m = -c / b;
final double l = r * m + q;
// x = ((1 - l*A) * xyScaling^2 ) / (2 * m)
// y = (-l*B * xyScaling^2) / ( 2 * m)
// z = (-l*C * zScaling^2)/ (2 * m)
final double denom0 = 0.5 / m;
final GeoPoint thePoint =
new GeoPoint(
(1.0 - l * A) * abSquared * denom0,
-l * B * abSquared * denom0,
-l * C * cSquared * denom0);
// Math is not quite accurate enough for this
// assert planetModel.pointOnSurface(thePoint): "Point: "+thePoint+";
// Planetmodel="+planetModel+"; A="+A+" B="+B+" C="+C+" D="+D+" planetfcn="+
// (thePoint.x*thePoint.x*planetModel.inverseXYScaling*planetModel.inverseXYScaling +
// thePoint.y*thePoint.y*planetModel.inverseXYScaling*planetModel.inverseXYScaling +
// thePoint.z*thePoint.z*planetModel.inverseZScaling*planetModel.inverseZScaling);
// assert evaluateIsZero(thePoint): "Evaluation of point: "+evaluate(thePoint);
addPoint(boundsInfo, bounds, thePoint);
} else {
// Something went very wrong; a = b = 0
}
}
// Do Y
if (!boundsInfo.isSmallestMinY(planetModel) || !boundsInfo.isLargestMaxY(planetModel)) {
// For min/max x, we need to use lagrange multipliers.
//
// For this, we need grad(F(x,y,z)) = (dF/dx, dF/dy, dF/dz).
//
// Minimize and maximize f(x,y,z) = y, with respect to g(x,y,z) = Ax + By + Cz - D and
// h(x,y,z) = x^2/xyScaling^2 + y^2/xyScaling^2 + z^2/zScaling^2 - 1
//
// grad(f(x,y,z)) = (0,1,0)
// grad(g(x,y,z)) = (A,B,C)
// grad(h(x,y,z)) = (2x/xyScaling^2,2y/xyScaling^2,2z/zScaling^2)
//
// Equations we need to simultaneously solve:
//
// grad(f(x,y,z)) = l * grad(g(x,y,z)) + m * grad(h(x,y,z))
// g(x,y,z) = 0
// h(x,y,z) = 0
//
// Equations:
// 0 = l*A + m*2x/xyScaling^2
// 1 = l*B + m*2y/xyScaling^2
// 0 = l*C + m*2z/zScaling^2
// Ax + By + Cz + D = 0
// x^2/xyScaling^2 + y^2/xyScaling^2 + z^2/zScaling^2 - 1 = 0
//
// Solve for x,y,z in terms of (l, m):
//
// x = (-l*A * xyScaling^2 ) / (2 * m)
// y = ((1 - l*B) * xyScaling^2) / ( 2 * m)
// z = (-l*C * zScaling^2)/ (2 * m)
//
// Two equations, two unknowns:
//
// A * ((-l*A * xyScaling^2 ) / (2 * m)) + B * (((1 - l*B) * xyScaling^2) / ( 2 * m)) + C *
// ((-l*C * zScaling^2)/ (2 * m)) + D = 0
//
// and
//
// ((-l*A * xyScaling^2 ) / (2 * m))^2/xyScaling^2 + (((1 - l*B) * xyScaling^2) / ( 2 *
// m))^2/xyScaling^2 + ((-l*C * zScaling^2)/ (2 * m))^2/zScaling^2 - 1 = 0
//
// Simple: solve for l and m, then find y from it.
//
// (a) Use first equation to find l in terms of m.
//
// A * ((-l*A * xyScaling^2 ) / (2 * m)) + B * (((1 - l*B) * xyScaling^2) / ( 2 * m)) + C *
// ((-l*C * zScaling^2)/ (2 * m)) + D = 0
// A * (-l*A * xyScaling^2 ) + B * ((1-l*B) * xyScaling^2) + C * (-l*C * zScaling^2) + D * 2 *
// m = 0
// -A^2*l*xyScaling^2 + B*xyScaling^2 - l*B^2*xyScaling^2 - C^2*l*zScaling^2 + D*2*m = 0
// - l *(A^2* xyScaling^2 + B^2 * xyScaling^2 + C^2 * zScaling^2) + (B * xyScaling^2 + D * 2 *
// m) = 0
// l = (B * xyScaling^2 + D * 2 * m) / (A^2* xyScaling^2 + B^2 * xyScaling^2 + C^2 *
// zScaling^2)
// l = B * xyScaling^2 / (A^2* xyScaling^2 + B^2 * xyScaling^2 + C^2 * zScaling^2) + m * 2 * D
// / (A^2* xyScaling^2 + B^2 * xyScaling^2 + C^2 * zScaling^2)
//
// For convenience:
//
// k = 1.0 / (A^2* xyScaling^2 + B^2 * xyScaling^2 + C^2 * zScaling^2)
//
// Then:
//
// l = B * xyScaling^2 * k + m * 2 * D * k
// l = k * (B*xyScaling^2 + m*2*D)
//
// For further convenience:
//
// q = B*xyScaling^2*k
// r = 2*D*k
//
// l = (r*m + q)
// l^2 = (r^2 * m^2 + 2*r*m*q + q^2)
//
// (b) Simplify the second equation before substitution
//
// ((-l*A * xyScaling^2 ) / (2 * m))^2/xyScaling^2 + (((1 - l*B) * xyScaling^2) / ( 2 *
// m))^2/xyScaling^2 + ((-l*C * zScaling^2)/ (2 * m))^2/zScaling^2 - 1 = 0
// (-l*A * xyScaling^2 )^2/xyScaling^2 + ((1 - l*B) * xyScaling^2)^2/xyScaling^2 + (-l*C *
// zScaling^2)^2/zScaling^2 = 4 * m^2
// (-l*A)^2 * xyScaling^2 + (1 - l*B)^2 * xyScaling^2 + (-l*C)^2 * zScaling^2 = 4 * m^2
// l^2*A^2 * xyScaling^2 + (1 - 2*l*B + l^2*B^2) * xyScaling^2 + l^2*C^2 * zScaling^2 = 4 *
// m^2
// A^2*xyScaling^2*l^2 + xyScaling^2 - 2*B*xyScaling^2*l + B^2*xyScaling^2*l^2 +
// C^2*zScaling^2*l^2 - 4*m^2 = 0
//
// (zScaling) Substitute for l, l^2
//
// A^2*xyScaling^2*(r^2 * m^2 + 2*r*m*q + q^2) + xyScaling^2 - 2*B*xyScaling^2*(r*m + q) +
// B^2*xyScaling^2*(r^2 * m^2 + 2*r*m*q + q^2) + C^2*zScaling^2*(r^2 * m^2 + 2*r*m*q + q^2) -
// 4*m^2 = 0
// A^2*xyScaling^2*r^2*m^2 + 2*A^2*xyScaling^2*r*q*m + A^2*xyScaling^2*q^2 + xyScaling^2 -
// 2*B*xyScaling^2*r*m - 2*B*xyScaling^2*q + B^2*xyScaling^2*r^2*m^2 +
// 2*B^2*xyScaling^2*r*q*m + B^2*xyScaling^2*q^2 + C^2*zScaling^2*r^2*m^2 +
// 2*C^2*zScaling^2*r*q*m + C^2*zScaling^2*q^2 - 4*m^2 = 0
//
// (d) Group
//
// m^2 * [A^2*xyScaling^2*r^2 + B^2*xyScaling^2*r^2 + C^2*zScaling^2*r^2 - 4] +
// m * [2*A^2*xyScaling^2*r*q - 2*B*xyScaling^2*r + 2*B^2*xyScaling^2*r*q +
// 2*C^2*zScaling^2*r*q] +
// [A^2*xyScaling^2*q^2 + xyScaling^2 - 2*B*xyScaling^2*q + B^2*xyScaling^2*q^2 +
// C^2*zScaling^2*q^2] = 0
// System.err.println(" computing Y bound");
// Useful subexpressions for this bound
final double q = B * abSquared * k;
final double qSquared = q * q;
// Quadratic equation
final double a =
ASquared * abSquared * rSquared
+ BSquared * abSquared * rSquared
+ CSquared * cSquared * rSquared
- 4.0;
final double b =
2.0 * ASquared * abSquared * r * q
- 2.0 * B * abSquared * r
+ 2.0 * BSquared * abSquared * r * q
+ 2.0 * CSquared * cSquared * r * q;
final double c =
ASquared * abSquared * qSquared
+ abSquared
- 2.0 * B * abSquared * q
+ BSquared * abSquared * qSquared
+ CSquared * cSquared * qSquared;
if (Math.abs(a) >= MINIMUM_RESOLUTION_SQUARED) {
final double sqrtTerm = b * b - 4.0 * a * c;
if (Math.abs(sqrtTerm) < MINIMUM_RESOLUTION_SQUARED) {
// One solution
final double m = -b / (2.0 * a);
// Valid?
if (Math.abs(m) >= MINIMUM_RESOLUTION) {
final double l = r * m + q;
// x = (-l*A * xyScaling^2 ) / (2 * m)
// y = ((1.0-l*B) * xyScaling^2) / ( 2 * m)
// z = (-l*C * zScaling^2)/ (2 * m)
final double denom0 = 0.5 / m;
final GeoPoint thePoint =
new GeoPoint(
-l * A * abSquared * denom0,
(1.0 - l * B) * abSquared * denom0,
-l * C * cSquared * denom0);
// Math is not quite accurate enough for this
// assert planetModel.pointOnSurface(thePoint): "Point: "+thePoint+";
// Planetmodel="+planetModel+"; A="+A+" B="+B+" C="+C+" D="+D+" planetfcn="+
// (thePoint1.x*thePoint.x*planetModel.inverseXYScaling*planetModel.inverseXYScaling +
// thePoint.y*thePoint.y*planetModel.inverseXYScaling*planetModel.inverseXYScaling +
// thePoint.z*thePoint.z*planetModel.inverseZScaling*planetModel.inverseZScaling);
// assert evaluateIsZero(thePoint): "Evaluation of point: "+evaluate(thePoint);
addPoint(boundsInfo, bounds, thePoint);
} else {
// This is a plane of the form A=0 B=n C=0. We can set a bound only by noting the D
// value.
boundsInfo.addYValue(-D / B);
}
} else if (sqrtTerm > 0.0) {
// Two solutions
final double sqrtResult = Math.sqrt(sqrtTerm);
final double commonDenom = 0.5 / a;
final double m1 = (-b + sqrtResult) * commonDenom;
assert Math.abs(a * m1 * m1 + b * m1 + c) < MINIMUM_RESOLUTION;
final double m2 = (-b - sqrtResult) * commonDenom;
assert Math.abs(a * m2 * m2 + b * m2 + c) < MINIMUM_RESOLUTION;
if (Math.abs(m1) >= MINIMUM_RESOLUTION || Math.abs(m2) >= MINIMUM_RESOLUTION) {
final double l1 = r * m1 + q;
final double l2 = r * m2 + q;
// x = (-l*A * xyScaling^2 ) / (2 * m)
// y = ((1.0-l*B) * xyScaling^2) / ( 2 * m)
// z = (-l*C * zScaling^2)/ (2 * m)
final double denom1 = 0.5 / m1;
final double denom2 = 0.5 / m2;
final GeoPoint thePoint1 =
new GeoPoint(
-l1 * A * abSquared * denom1,
(1.0 - l1 * B) * abSquared * denom1,
-l1 * C * cSquared * denom1);
final GeoPoint thePoint2 =
new GeoPoint(
-l2 * A * abSquared * denom2,
(1.0 - l2 * B) * abSquared * denom2,
-l2 * C * cSquared * denom2);
// Math is not quite accurate enough for this
// assert planetModel.pointOnSurface(thePoint1): "Point1: "+thePoint1+";
// Planetmodel="+planetModel+"; A="+A+" B="+B+" C="+C+" D="+D+" planetfcn="+
// (thePoint1.x*thePoint1.x*planetModel.inverseXYScaling*planetModel.inverseXYScaling +
// thePoint1.y*thePoint1.y*planetModel.inverseXYScaling*planetModel.inverseXYScaling +
// thePoint1.z*thePoint1.z*planetModel.inverseZScaling*planetModel.inverseZScaling);
// assert planetModel.pointOnSurface(thePoint2): "Point2: "+thePoint2+";
// Planetmodel="+planetModel+"; A="+A+" B="+B+" C="+C+" D="+D+" planetfcn="+
// (thePoint2.x*thePoint2.x*planetModel.inverseXYScaling*planetModel.inverseXYScaling +
// thePoint2.y*thePoint2.y*planetModel.inverseXYScaling*planetModel.inverseXYScaling +
// thePoint2.z*thePoint2.z*planetModel.inverseZScaling*planetModel.inverseZScaling);
// assert evaluateIsZero(thePoint1): "Evaluation of point1: "+evaluate(thePoint1);
// assert evaluateIsZero(thePoint2): "Evaluation of point2: "+evaluate(thePoint2);
addPoint(boundsInfo, bounds, thePoint1);
addPoint(boundsInfo, bounds, thePoint2);
} else {
// This is a plane of the form A=0 B=n C=0. We can set a bound only by noting the D
// value.
boundsInfo.addYValue(-D / B);
}
} else {
// No solutions
}
} else if (Math.abs(b) > MINIMUM_RESOLUTION_SQUARED) {
// a = 0, so m = - zScaling / b
final double m = -c / b;
final double l = r * m + q;
// x = ( -l*A * xyScaling^2 ) / (2 * m)
// y = ((1-l*B) * xyScaling^2) / ( 2 * m)
// z = (-l*C * zScaling^2)/ (2 * m)
final double denom0 = 0.5 / m;
final GeoPoint thePoint =
new GeoPoint(
-l * A * abSquared * denom0,
(1.0 - l * B) * abSquared * denom0,
-l * C * cSquared * denom0);
// Math is not quite accurate enough for this
// assert planetModel.pointOnSurface(thePoint): "Point: "+thePoint+";
// Planetmodel="+planetModel+"; A="+A+" B="+B+" C="+C+" D="+D+" planetfcn="+
// (thePoint.x*thePoint.x*planetModel.inverseXYScaling*planetModel.inverseXYScaling +
// thePoint.y*thePoint.y*planetModel.inverseXYScaling*planetModel.inverseXYScaling +
// thePoint.z*thePoint.z*planetModel.inverseZScaling*planetModel.inverseZScaling);
// assert evaluateIsZero(thePoint): "Evaluation of point: "+evaluate(thePoint);
addPoint(boundsInfo, bounds, thePoint);
} else {
// Something went very wrong; a = b = 0
}
}
}