LUCENE-8076: It is possible to create planet models that are not normalized and this messes up surface distance computations. Therefore, normalize these computations. Committed on behalf of Ignacio Vera.

lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/PlanetModel.java

@@ -60,6 +60,11 @@ 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.
   
@@ -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);

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(){