diff --git a/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/PlanetModel.java b/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/PlanetModel.java index 37297d3172c..2aabfc1b824 100644 --- a/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/PlanetModel.java +++ b/lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/PlanetModel.java @@ -60,7 +60,12 @@ public class PlanetModel implements SerializableObject { public final double flattening; /** The square ratio */ public final double squareRatio; - + /** The scale of the planet */ + public final double scale; + /** The inverse of scale */ + public final double inverseScale; + + // We do NOT include radius, because all computations in geo3d are in radians, not meters. // Compute north and south pole for planet model, since these are commonly used. @@ -99,6 +104,8 @@ public class PlanetModel implements SerializableObject { this.MAX_X_POLE = new GeoPoint(ab, 1.0, 0.0, 0.0, 0.0, 0.0); this.MIN_Y_POLE = new GeoPoint(ab, 0.0, -1.0, 0.0, 0.0, -Math.PI * 0.5); this.MAX_Y_POLE = new GeoPoint(ab, 0.0, 1.0, 0.0, 0.0, Math.PI * 0.5); + this.scale = (2.0 * ab + c)/3.0; + this.inverseScale = 1.0 / scale; this.minimumPoleDistance = Math.min(surfaceDistance(NORTH_POLE, SOUTH_POLE), surfaceDistance(MIN_X_POLE, MAX_X_POLE)); } @@ -315,14 +322,13 @@ public class PlanetModel implements SerializableObject { lambda = L + (1.0 - C) * flattening * sinAlpha * (sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1.0 + 2.0 * cos2SigmaM *cos2SigmaM))); } while (Math.abs(lambda-lambdaP) >= Vector.MINIMUM_RESOLUTION && ++iterLimit < 100); - final double uSq = cosSqAlpha * this.squareRatio; final double A = 1.0 + uSq / 16384.0 * (4096.0 + uSq * (-768.0 + uSq * (320.0 - 175.0 * uSq))); final double B = uSq / 1024.0 * (256.0 + uSq * (-128.0 + uSq * (74.0 - 47.0 * uSq))); final double deltaSigma = B * sinSigma * (cos2SigmaM + B / 4.0 * (cosSigma * (-1.0 + 2.0 * cos2SigmaM * cos2SigmaM)- B / 6.0 * cos2SigmaM * (-3.0 + 4.0 * sinSigma * sinSigma) * (-3.0 + 4.0 * cos2SigmaM * cos2SigmaM))); - return c * A * (sigma - deltaSigma); + return c * inverseScale * A * (sigma - deltaSigma); } /** Compute new point given original point, a bearing direction, and an adjusted angle (as would be computed by @@ -360,7 +366,7 @@ public class PlanetModel implements SerializableObject { double cosσ; double Δσ; - double σ = dist / (c * A); + double σ = dist / (c * inverseScale * A); double σʹ; double iterations = 0; do { @@ -370,9 +376,8 @@ public class PlanetModel implements SerializableObject { Δσ = B * sinσ * (cos2σM + B / 4.0 * (cosσ * (-1.0 + 2.0 * cos2σM * cos2σM) - B / 6.0 * cos2σM * (-3.0 + 4.0 * sinσ * sinσ) * (-3.0 + 4.0 * cos2σM * cos2σM))); σʹ = σ; - σ = dist / (c * A) + Δσ; + σ = dist / (c * inverseScale * A) + Δσ; } while (Math.abs(σ - σʹ) >= Vector.MINIMUM_RESOLUTION && ++iterations < 100); - double x = sinU1 * sinσ - cosU1 * cosσ * cosα1; double φ2 = Math.atan2(sinU1 * cosσ + cosU1 * sinσ * cosα1, (1.0 - flattening) * Math.sqrt(sinα * sinα + x * x)); double λ = Math.atan2(sinσ * sinα1, cosU1 * cosσ - sinU1 * sinσ * cosα1); diff --git a/lucene/spatial3d/src/test/org/apache/lucene/spatial3d/geom/GeoExactCircleTest.java b/lucene/spatial3d/src/test/org/apache/lucene/spatial3d/geom/GeoExactCircleTest.java index 859cdedbced..2f758f1169a 100644 --- a/lucene/spatial3d/src/test/org/apache/lucene/spatial3d/geom/GeoExactCircleTest.java +++ b/lucene/spatial3d/src/test/org/apache/lucene/spatial3d/geom/GeoExactCircleTest.java @@ -52,6 +52,31 @@ public class GeoExactCircleTest extends RandomGeo3dShapeGenerator{ assertTrue(c.isWithin(gp)); } + @Test + public void testSurfacePointOnBearingScale(){ + double ab = 1.6; + double c = 0.7; + PlanetModel p1 = PlanetModel.WGS84; + PlanetModel p2 = new PlanetModel(0.5 * PlanetModel.WGS84.ab, 0.5 * PlanetModel.WGS84.c ); + GeoPoint point1P1 = new GeoPoint(p1, 0, 0); + GeoPoint point2P1 = new GeoPoint(p1, 1, 1); + GeoPoint point1P2 = new GeoPoint(p2, point1P1.getLatitude(), point1P1.getLongitude()); + GeoPoint point2P2 = new GeoPoint(p2, point2P1.getLatitude(), point2P1.getLongitude()); + + double dist = 0.2* Math.PI; + double bearing = 0.2 * Math.PI; + + GeoPoint new1 = p1.surfacePointOnBearing(point2P1, dist, bearing); + GeoPoint new2 = p2.surfacePointOnBearing(point2P2, dist, bearing); + + assertEquals(new1.getLatitude(), new2.getLatitude(), 1e-12); + assertEquals(new1.getLongitude(), new2.getLongitude(), 1e-12); + //This is true if surfaceDistance return results always in radians + double d1 = p1.surfaceDistance(point1P1, point2P1); + double d2 = p2.surfaceDistance(point1P2, point2P2); + assertEquals(d1, d2, 1e-12); + } + @Test @Repeat(iterations = 100) public void RandomPointBearingWGS84Test(){