LUCENE-6544: Geo3D: Regularize path & polygon construction, add PlanetModel.surfaceDistance(), cache lat & lon in GeoPoint, add thread-safety where missing -- Geo3dShape.

git-svn-id: https://svn.apache.org/repos/asf/lucene/dev/trunk@1686285 13f79535-47bb-0310-9956-ffa450edef68

lucene/CHANGES.txt

@@ -95,6 +95,11 @@ New Features
 * LUCENE-6577: Give earlier and better error message for invalid CRC.
   (Robert Muir)
 
+* LUCENE-6544: Geo3D: (1) Regularize path & polygon construction, (2) add
+  PlanetModel.surfaceDistance() (ellipsoidal calculation), (3) cache lat & lon
+  in GeoPoint, (4) add thread-safety where missing -- Geo3dShape. (Karl Wright,
+  David Smiley)
+
 API Changes
 
 * LUCENE-6508: Simplify Lock api, there is now just 

lucene/spatial/src/java/org/apache/lucene/spatial/spatial4j/Geo3dShape.java

@@ -18,6 +18,7 @@ package org.apache.lucene.spatial.spatial4j;
  */
 
 import com.spatial4j.core.context.SpatialContext;
+import com.spatial4j.core.distance.DistanceUtils;
 import com.spatial4j.core.shape.Point;
 import com.spatial4j.core.shape.Rectangle;
 import com.spatial4j.core.shape.Shape;
@@ -31,20 +32,21 @@ import org.apache.lucene.spatial.spatial4j.geo3d.GeoShape;
 import org.apache.lucene.spatial.spatial4j.geo3d.PlanetModel;
 
 /**
- * A 3D planar geometry based Spatial4j Shape implementation.
+ * A Spatial4j Shape wrapping a {@link GeoShape} ("Geo3D") -- a 3D planar geometry based Spatial4j Shape implementation.
+ * Geo3D implements shapes on the surface of a sphere or ellipsoid.
  *
  * @lucene.experimental
  */
 public class Geo3dShape implements Shape {
+  /** The required size of this adjustment depends on the actual planetary model chosen.
+   * This value is big enough to account for WGS84. */
+  protected static final double ROUNDOFF_ADJUSTMENT = 0.05;
 
   public final SpatialContext ctx;
   public final GeoShape shape;
   public final PlanetModel planetModel;
 
-  private Rectangle boundingBox = null;
+  private volatile Rectangle boundingBox = null; // lazy initialized
-
-  public final static double RADIANS_PER_DEGREE = Math.PI / 180.0;
-  public final static double DEGREES_PER_RADIAN = 1.0 / RADIANS_PER_DEGREE;
 
   public Geo3dShape(final GeoShape shape, final SpatialContext ctx) {
     this(PlanetModel.SPHERE, shape, ctx);
@@ -72,10 +74,10 @@ public class Geo3dShape implements Shape {
   protected SpatialRelation relate(Rectangle r) {
     // Construct the right kind of GeoArea first
     GeoArea geoArea = GeoAreaFactory.makeGeoArea(planetModel,
-        r.getMaxY() * RADIANS_PER_DEGREE,
+        r.getMaxY() * DistanceUtils.DEGREES_TO_RADIANS,
-        r.getMinY() * RADIANS_PER_DEGREE,
+        r.getMinY() * DistanceUtils.DEGREES_TO_RADIANS,
-        r.getMinX() * RADIANS_PER_DEGREE,
+        r.getMinX() * DistanceUtils.DEGREES_TO_RADIANS,
-        r.getMaxX() * RADIANS_PER_DEGREE);
+        r.getMaxX() * DistanceUtils.DEGREES_TO_RADIANS);
     int relationship = geoArea.getRelationship(shape);
     if (relationship == GeoArea.WITHIN)
       return SpatialRelation.WITHIN;
@@ -91,7 +93,7 @@ public class Geo3dShape implements Shape {
 
   protected SpatialRelation relate(Point p) {
     // Create a GeoPoint
-    GeoPoint point = new GeoPoint(planetModel, p.getY()*RADIANS_PER_DEGREE, p.getX()*RADIANS_PER_DEGREE);
+    GeoPoint point = new GeoPoint(planetModel, p.getY()* DistanceUtils.DEGREES_TO_RADIANS, p.getX()* DistanceUtils.DEGREES_TO_RADIANS);
     if (shape.isWithin(point)) {
       // Point within shape
       return SpatialRelation.CONTAINS;
@@ -99,13 +101,12 @@ public class Geo3dShape implements Shape {
     return SpatialRelation.DISJOINT;
   }
 
-  // The required size of this adjustment depends on the actual planetary model chosen.
+
-  // This value is big enough to account for WGS84.
-  protected final double ROUNDOFF_ADJUSTMENT = 0.05;
   
   @Override
   public Rectangle getBoundingBox() {
-    if (boundingBox == null) {
+    Rectangle bbox = this.boundingBox;//volatile read once
+    if (bbox == null) {
       Bounds bounds = shape.getBounds(null);
       double leftLon;
       double rightLon;
@@ -113,24 +114,25 @@ public class Geo3dShape implements Shape {
         leftLon = -180.0;
         rightLon = 180.0;
       } else {
-        leftLon = bounds.getLeftLongitude().doubleValue() * DEGREES_PER_RADIAN;
+        leftLon = bounds.getLeftLongitude().doubleValue() * DistanceUtils.RADIANS_TO_DEGREES;
-        rightLon = bounds.getRightLongitude().doubleValue() * DEGREES_PER_RADIAN;
+        rightLon = bounds.getRightLongitude().doubleValue() * DistanceUtils.RADIANS_TO_DEGREES;
       }
       double minLat;
       if (bounds.checkNoBottomLatitudeBound()) {
         minLat = -90.0;
       } else {
-        minLat = bounds.getMinLatitude().doubleValue() * DEGREES_PER_RADIAN;
+        minLat = bounds.getMinLatitude().doubleValue() * DistanceUtils.RADIANS_TO_DEGREES;
       }
       double maxLat;
       if (bounds.checkNoTopLatitudeBound()) {
         maxLat = 90.0;
       } else {
-        maxLat = bounds.getMaxLatitude().doubleValue() * DEGREES_PER_RADIAN;
+        maxLat = bounds.getMaxLatitude().doubleValue() * DistanceUtils.RADIANS_TO_DEGREES;
       }
-      boundingBox = new RectangleImpl(leftLon, rightLon, minLat, maxLat, ctx).getBuffered(ROUNDOFF_ADJUSTMENT, ctx);
+      bbox = new RectangleImpl(leftLon, rightLon, minLat, maxLat, ctx).getBuffered(ROUNDOFF_ADJUSTMENT, ctx);
+      this.boundingBox = bbox;
     }
-    return boundingBox;
+    return bbox;
   }
 
   @Override
@@ -140,17 +142,17 @@ public class Geo3dShape implements Shape {
 
   @Override
   public double getArea(SpatialContext ctx) {
-    throw new RuntimeException("Unimplemented");
+    throw new UnsupportedOperationException();
   }
 
   @Override
   public Point getCenter() {
-    throw new RuntimeException("Unimplemented");
+    throw new UnsupportedOperationException();
   }
 
   @Override
   public Shape getBuffered(double distance, SpatialContext ctx) {
-    return this;
+    throw new UnsupportedOperationException();
   }
 
   @Override
@@ -160,7 +162,7 @@ public class Geo3dShape implements Shape {
 
   @Override
   public String toString() {
-    return "Geo3dShape{planetmodel=" + planetModel+", shape="+shape + '}';
+    return "Geo3dShape{planetmodel=" + planetModel + ", shape=" + shape + '}';
   }
 
   @Override
@@ -168,7 +170,7 @@ public class Geo3dShape implements Shape {
     if (!(other instanceof Geo3dShape))
       return false;
     Geo3dShape tr = (Geo3dShape)other;
-    return tr.planetModel.equals(planetModel) && tr.shape.equals(shape);
+    return tr.ctx.equals(ctx) && tr.planetModel.equals(planetModel) && tr.shape.equals(shape);
   }
 
   @Override

lucene/spatial/src/java/org/apache/lucene/spatial/spatial4j/geo3d/GeoCircle.java

@@ -39,31 +39,40 @@ public class GeoCircle extends GeoBaseExtendedShape implements GeoDistanceShape,
       throw new IllegalArgumentException("Cutoff angle out of bounds");
     final double cosAngle = Math.cos(cutoffAngle);
     this.center = new GeoPoint(planetModel, lat, lon);
-    final double magnitude = center.magnitude();
     // In an ellipsoidal world, cutoff distances make no sense, unfortunately.  Only membership
     // can be used to make in/out determination.
     this.cutoffAngle = cutoffAngle;
-    // The plane's normal vector needs to be normalized, since we compute D on that basis
+    // Compute two points on the circle, with the right angle from the center.  We'll use these
-    this.circlePlane = new SidedPlane(center, center.normalize(), -cosAngle * magnitude);
+    // to obtain the perpendicular plane to the circle.
-
+    double upperLat = lat + cutoffAngle;
+    double upperLon = lon;
+    if (upperLat > Math.PI * 0.5) {
+      upperLon += Math.PI;
+      if (upperLon > Math.PI)
+        upperLon -= 2.0 * Math.PI;
+      upperLat = Math.PI - upperLat;
+    }
+    double lowerLat = lat - cutoffAngle;
+    double lowerLon = lon;
+    if (lowerLat < -Math.PI * 0.5) {
+      lowerLon += Math.PI;
+      if (lowerLon > Math.PI)
+        lowerLon -= 2.0 * Math.PI;
+      lowerLat = -Math.PI - lowerLat;
+    }
+    final GeoPoint upperPoint = new GeoPoint(planetModel, upperLat, upperLon);
+    final GeoPoint lowerPoint = new GeoPoint(planetModel, lowerLat, lowerLon);
+    // Construct normal plane
+    final Plane normalPlane = new Plane(upperPoint, center);
+    // Construct a sided plane that goes through the two points and whose normal is in the normalPlane.
+    this.circlePlane = SidedPlane.constructNormalizedPerpendicularSidedPlane(center, normalPlane, upperPoint, lowerPoint);
+    if (circlePlane == null)
+      throw new RuntimeException("Couldn't construct circle plane.  Cutoff angle = "+cutoffAngle+"; upperPoint = "+upperPoint+"; lowerPoint = "+lowerPoint);
     // Compute a point on the circle boundary.
     if (cutoffAngle == Math.PI)
       this.edgePoints = new GeoPoint[0];
     else {
-      // We already have circle plane, which is the definitive determination of the edge of the "circle".
+      this.edgePoints = new GeoPoint[]{upperPoint};
-      // Next, compute vertical plane going through origin and the center point (C = 0, D = 0).
-      Plane verticalPlane = Plane.constructNormalizedVerticalPlane(this.center.x, this.center.y);
-      if (verticalPlane == null) {
-        verticalPlane = new Plane(1.0,0.0);
-      }
-      // Finally, use Plane.findIntersections() to find the intersection points.
-      final GeoPoint edgePoint = this.circlePlane.getSampleIntersectionPoint(planetModel, verticalPlane);
-      if (edgePoint == null) {
-        throw new RuntimeException("Could not find edge point for circle at lat="+lat+" lon="+lon+" cutoffAngle="+cutoffAngle+" planetModel="+planetModel);
-      }
-      //if (Math.abs(circlePlane.evaluate(edgePoint)) > 1e-10)
-      //    throw new RuntimeException("Computed an edge point that does not satisfy circlePlane equation! "+circlePlane.evaluate(edgePoint));
-      this.edgePoints = new GeoPoint[]{edgePoint};
     }
   }
 

lucene/spatial/src/java/org/apache/lucene/spatial/spatial4j/geo3d/GeoConvexPolygon.java

@@ -34,13 +34,14 @@ public class GeoConvexPolygon extends GeoBaseExtendedShape implements GeoMembers
   protected final BitSet isInternalEdges;
 
   protected SidedPlane[] edges = null;
-  protected boolean[] internalEdges = null;
   protected GeoPoint[][] notableEdgePoints = null;
 
   protected GeoPoint[] edgePoints = null;
 
   protected double fullDistance = 0.0;
 
+  protected boolean isDone = false;
+  
   /**
    * Create a convex polygon from a list of points.  The first point must be on the
    * external edge.
@@ -48,19 +49,20 @@ public class GeoConvexPolygon extends GeoBaseExtendedShape implements GeoMembers
   public GeoConvexPolygon(final PlanetModel planetModel, final List<GeoPoint> pointList) {
     super(planetModel);
     this.points = pointList;
-    this.isInternalEdges = null;
+    this.isInternalEdges = new BitSet();
-    donePoints(false);
+    done(false);
   }
 
   /**
    * Create a convex polygon from a list of points, keeping track of which boundaries
    * are internal.  This is used when creating a polygon as a building block for another shape.
    */
-  public GeoConvexPolygon(final PlanetModel planetModel, final List<GeoPoint> pointList, final BitSet internalEdgeFlags, final boolean returnEdgeInternal) {
+  public GeoConvexPolygon(final PlanetModel planetModel, final List<GeoPoint> pointList, final BitSet internalEdgeFlags,
+                          final boolean returnEdgeInternal) {
     super(planetModel);
     this.points = pointList;
     this.isInternalEdges = internalEdgeFlags;
-    donePoints(returnEdgeInternal);
+    done(returnEdgeInternal);
   }
 
   /**
@@ -69,16 +71,9 @@ public class GeoConvexPolygon extends GeoBaseExtendedShape implements GeoMembers
    */
   public GeoConvexPolygon(final PlanetModel planetModel, final double startLatitude, final double startLongitude) {
     super(planetModel);
-    points = new ArrayList<GeoPoint>();
+    points = new ArrayList<>();
     isInternalEdges = new BitSet();
-    // Argument checking
+    points.add(new GeoPoint(planetModel, startLatitude, startLongitude));
-    if (startLatitude > Math.PI * 0.5 || startLatitude < -Math.PI * 0.5)
-      throw new IllegalArgumentException("Latitude out of range");
-    if (startLongitude < -Math.PI || startLongitude > Math.PI)
-      throw new IllegalArgumentException("Longitude out of range");
-
-    final GeoPoint p = new GeoPoint(planetModel, startLatitude, startLongitude);
-    points.add(p);
   }
 
   /**
@@ -87,36 +82,39 @@ public class GeoConvexPolygon extends GeoBaseExtendedShape implements GeoMembers
    *
    * @param latitude       is the latitude of the next point.
    * @param longitude      is the longitude of the next point.
-   * @param isInternalEdge is true if the edge just added should be considered "internal", and not
+   * @param isInternalEdge is true if the edge just added with this point should be considered "internal", and not
    *                       intersected as part of the intersects() operation.
    */
   public void addPoint(final double latitude, final double longitude, final boolean isInternalEdge) {
-    // Argument checking
+    if (isDone)
-    if (latitude > Math.PI * 0.5 || latitude < -Math.PI * 0.5)
+      throw new IllegalStateException("Can't call addPoint() if done() already called");
-      throw new IllegalArgumentException("Latitude out of range");
+    if (isInternalEdge)
-    if (longitude < -Math.PI || longitude > Math.PI)
+      isInternalEdges.set(points.size() - 1);
-      throw new IllegalArgumentException("Longitude out of range");
+    points.add(new GeoPoint(planetModel, latitude, longitude));
-
-    final GeoPoint p = new GeoPoint(planetModel, latitude, longitude);
-    isInternalEdges.set(points.size(), isInternalEdge);
-    points.add(p);
   }
 
   /**
    * Finish the polygon, by connecting the last added point with the starting point.
    */
-  public void donePoints(final boolean isInternalReturnEdge) {
+  public void done(final boolean isInternalReturnEdge) {
+    if (isDone)
+      throw new IllegalStateException("Can't call done() more than once");
     // If fewer than 3 points, can't do it.
     if (points.size() < 3)
       throw new IllegalArgumentException("Polygon needs at least three points.");
+
+    if (isInternalReturnEdge)
+      isInternalEdges.set(points.size() - 1);
+
+    isDone = true;
+    
     // Time to construct the planes.  If the polygon is truly convex, then any adjacent point
     // to a segment can provide an interior measurement.
     edges = new SidedPlane[points.size()];
     notableEdgePoints = new GeoPoint[points.size()][];
-    internalEdges = new boolean[points.size()];
+
     for (int i = 0; i < points.size(); i++) {
       final GeoPoint start = points.get(i);
-      final boolean isInternalEdge = (isInternalEdges != null ? (i == isInternalEdges.size() ? isInternalReturnEdge : isInternalEdges.get(i)) : false);
       final GeoPoint end = points.get(legalIndex(i + 1));
       final double distance = start.arcDistance(end);
       if (distance > fullDistance)
@@ -126,7 +124,6 @@ public class GeoConvexPolygon extends GeoBaseExtendedShape implements GeoMembers
       //System.out.println("Created edge "+sp+" using start="+start+" end="+end+" check="+check);
       edges[i] = sp;
       notableEdgePoints[i] = new GeoPoint[]{start, end};
-      internalEdges[i] = isInternalEdge;
     }
     createCenterPoint();
   }
@@ -183,7 +180,7 @@ public class GeoConvexPolygon extends GeoBaseExtendedShape implements GeoMembers
     for (int edgeIndex = 0; edgeIndex < edges.length; edgeIndex++) {
       final SidedPlane edge = edges[edgeIndex];
       final GeoPoint[] points = this.notableEdgePoints[edgeIndex];
-      if (!internalEdges[edgeIndex]) {
+      if (!isInternalEdges.get(edgeIndex)) {
         //System.err.println(" non-internal edge "+edge);
         // Edges flagged as 'internal only' are excluded from the matching
         // Construct boundaries
@@ -250,14 +247,9 @@ public class GeoConvexPolygon extends GeoBaseExtendedShape implements GeoMembers
     GeoConvexPolygon other = (GeoConvexPolygon) o;
     if (!super.equals(other))
       return false;
-    if (other.points.size() != points.size())
+    if (!other.isInternalEdges.equals(isInternalEdges))
       return false;
-
+    return (other.points.equals(points));
-    for (int i = 0; i < points.size(); i++) {
-      if (!other.points.get(i).equals(points.get(i)))
-        return false;
-    }
-    return true;
   }
 
   @Override

lucene/spatial/src/java/org/apache/lucene/spatial/spatial4j/geo3d/GeoPath.java

@@ -18,6 +18,7 @@ package org.apache.lucene.spatial.spatial4j.geo3d;
  */
 
 import java.util.ArrayList;
+import java.util.Collections;
 import java.util.List;
 
 /**
@@ -31,15 +32,24 @@ import java.util.List;
 public class GeoPath extends GeoBaseExtendedShape implements GeoDistanceShape {
   
   public final double cutoffAngle;
-  public final double sinAngle;
+
-  public final double cosAngle;
+  public final double sinAngle; // sine of cutoffAngle
+  public final double cosAngle; // cosine of cutoffAngle
 
   public final List<GeoPoint> points = new ArrayList<GeoPoint>();
   
-  public final List<SegmentEndpoint> endPoints = new ArrayList<SegmentEndpoint>();
+  public List<SegmentEndpoint> endPoints;
-  public final List<PathSegment> segments = new ArrayList<PathSegment>();
+  public List<PathSegment> segments;
 
-  public GeoPoint[] edgePoints = null;
+  public GeoPoint[] edgePoints;
+
+  public boolean isDone = false;
+  
+  public GeoPath(final PlanetModel planetModel, final double maxCutoffAngle, final GeoPoint[] pathPoints) {
+    this(planetModel, maxCutoffAngle);
+    Collections.addAll(points, pathPoints);
+    done();
+  }
   
   public GeoPath(final PlanetModel planetModel, final double maxCutoffAngle) {
     super(planetModel);
@@ -51,16 +61,20 @@ public class GeoPath extends GeoBaseExtendedShape implements GeoDistanceShape {
   }
 
   public void addPoint(double lat, double lon) {
-    if (lat < -Math.PI * 0.5 || lat > Math.PI * 0.5)
+    if (isDone)
-      throw new IllegalArgumentException("Latitude out of range");
+      throw new IllegalStateException("Can't call addPoint() if done() already called");
-    if (lon < -Math.PI || lon > Math.PI)
-      throw new IllegalArgumentException("Longitude out of range");
     points.add(new GeoPoint(planetModel, lat, lon));
   }
   
   public void done() {
+    if (isDone)
+      throw new IllegalStateException("Can't call done() twice");
     if (points.size() == 0)
       throw new IllegalArgumentException("Path must have at least one point");
+    isDone = true;
+
+    endPoints = new ArrayList<>(points.size());
+    segments = new ArrayList<>(points.size());
     // Compute an offset to use for all segments.  This will be based on the minimum magnitude of
     // the entire ellipsoid.
     final double cutoffOffset = this.sinAngle * planetModel.getMinimumMagnitude();
@@ -80,21 +94,32 @@ public class GeoPath extends GeoBaseExtendedShape implements GeoDistanceShape {
     
     if (segments.size() == 0) {
       // Simple circle
-      final SegmentEndpoint onlyEndpoint = new SegmentEndpoint(points.get(0), cutoffOffset);
+      double lat = points.get(0).getLatitude();
+      double lon = points.get(0).getLongitude();
+      // Compute two points on the circle, with the right angle from the center.  We'll use these
+      // to obtain the perpendicular plane to the circle.
+      double upperLat = lat + cutoffAngle;
+      double upperLon = lon;
+      if (upperLat > Math.PI * 0.5) {
+        upperLon += Math.PI;
+        if (upperLon > Math.PI)
+          upperLon -= 2.0 * Math.PI;
+        upperLat = Math.PI - upperLat;
+      }
+      double lowerLat = lat - cutoffAngle;
+      double lowerLon = lon;
+      if (lowerLat < -Math.PI * 0.5) {
+        lowerLon += Math.PI;
+        if (lowerLon > Math.PI)
+          lowerLon -= 2.0 * Math.PI;
+        lowerLat = -Math.PI - lowerLat;
+      }
+      final GeoPoint upperPoint = new GeoPoint(planetModel, upperLat, upperLon);
+      final GeoPoint lowerPoint = new GeoPoint(planetModel, lowerLat, lowerLon);
+
+      final SegmentEndpoint onlyEndpoint = new SegmentEndpoint(points.get(0), upperPoint, lowerPoint);
       endPoints.add(onlyEndpoint);
-      // Find an edgepoint
+      this.edgePoints = new GeoPoint[]{upperPoint};
-      // We already have circle plane, which is the definitive determination of the edge of the "circle".
-      // Next, compute vertical plane going through origin and the center point (C = 0, D = 0).
-      Plane verticalPlane = Plane.constructNormalizedVerticalPlane(onlyEndpoint.point.x, onlyEndpoint.point.y);
-      if (verticalPlane == null) {
-        verticalPlane = new Plane(1.0,0.0);
-      }
-      // Finally, use Plane.findIntersections() to find the intersection points.
-      final GeoPoint edgePoint = onlyEndpoint.circlePlane.getSampleIntersectionPoint(planetModel, verticalPlane);
-      if (edgePoint == null) {
-        throw new RuntimeException("Could not find edge point for path endpoint="+onlyEndpoint.point+" cutoffOffset="+cutoffOffset+" planetModel="+planetModel);
-      }
-      this.edgePoints = new GeoPoint[]{edgePoint};
       return;
     }
     
@@ -386,17 +411,9 @@ public class GeoPath extends GeoBaseExtendedShape implements GeoDistanceShape {
     GeoPath p = (GeoPath) o;
     if (!super.equals(p))
       return false;
-    if (endPoints.size() != p.endPoints.size())
-      return false;
     if (cutoffAngle != p.cutoffAngle)
       return false;
-    for (int i = 0; i < endPoints.size(); i++) {
+    return points.equals(p.points);
-      SegmentEndpoint point = endPoints.get(i);
-      SegmentEndpoint point2 = p.endPoints.get(i);
-      if (!point.equals(point2))
-        return false;
-    }
-    return true;
   }
 
   @Override
@@ -404,7 +421,7 @@ public class GeoPath extends GeoBaseExtendedShape implements GeoDistanceShape {
     int result = super.hashCode();
     long temp = Double.doubleToLongBits(cutoffAngle);
     result = 31 * result + (int) (temp ^ (temp >>> 32));
-    result = 31 * result + endPoints.hashCode();
+    result = 31 * result + points.hashCode();
     return result;
   }
 
@@ -443,11 +460,12 @@ public class GeoPath extends GeoBaseExtendedShape implements GeoDistanceShape {
     /** Constructor for case (1).
      * Generate a simple circle cutoff plane.
      */
-    public SegmentEndpoint(final GeoPoint point, final double cutoffOffset) {
+    public SegmentEndpoint(final GeoPoint point, final GeoPoint upperPoint, final GeoPoint lowerPoint) {
       this.point = point;
-      final double magnitude = point.magnitude();
+      // Construct normal plane
-      // Normalize vector to make D value correct
+      final Plane normalPlane = new Plane(upperPoint, point);
-      this.circlePlane = new SidedPlane(point, point.normalize(), -Math.sqrt(magnitude * magnitude - cutoffOffset * cutoffOffset));
+      // Construct a sided plane that goes through the two points and whose normal is in the normalPlane.
+      this.circlePlane = SidedPlane.constructNormalizedPerpendicularSidedPlane(point, normalPlane, upperPoint, lowerPoint);
       this.cutoffPlanes = new Membership[0];
       this.notablePoints = new GeoPoint[0];
     }

lucene/spatial/src/java/org/apache/lucene/spatial/spatial4j/geo3d/GeoPoint.java

@@ -24,9 +24,34 @@ package org.apache.lucene.spatial.spatial4j.geo3d;
  */
 public class GeoPoint extends Vector {
   
+  // By making lazily-evaluated variables be "volatile", we guarantee atomicity when they
+  // are updated.  This is necessary if we are using these classes in a multi-thread fashion,
+  // because we don't try to synchronize for the lazy computation.
+  
   /** This is the lazily-evaluated magnitude.  Some constructors include it, but others don't, and
-   * we try not to create extra computation by always computing it. */
+   * we try not to create extra computation by always computing it.  Does not need to be
-  protected double magnitude = Double.NEGATIVE_INFINITY;
+   * synchronized for thread safety, because depends wholly on immutable variables of this class. */
+  protected volatile double magnitude = Double.NEGATIVE_INFINITY;
+  /** Lazily-evaluated latitude.  Does not need to be
+   * synchronized for thread safety, because depends wholly on immutable variables of this class.  */
+  protected volatile double latitude = Double.NEGATIVE_INFINITY;
+  /** Lazily-evaluated longitude.   Does not need to be
+   * synchronized for thread safety, because depends wholly on immutable variables of this class.  */
+  protected volatile double longitude = Double.NEGATIVE_INFINITY;
+
+  /** Construct a GeoPoint from the trig functions of a lat and lon pair.
+   * @param planetModel is the planetModel to put the point on.
+   * @param sinLat is the sin of the latitude.
+   * @param sinLon is the sin of the longitude.
+   * @param cosLat is the cos of the latitude.
+   * @param cosLon is the cos of the longitude.
+   * @param lat is the latitude.
+   * @param lon is the longitude.
+   */
+  public GeoPoint(final PlanetModel planetModel, final double sinLat, final double sinLon, final double cosLat, final double cosLon, final double lat, final double lon) {
+    this(computeDesiredEllipsoidMagnitude(planetModel, cosLat * cosLon, cosLat * sinLon, sinLat),
+      cosLat * cosLon, cosLat * sinLon, sinLat, lat, lon);
+  }
   
   /** Construct a GeoPoint from the trig functions of a lat and lon pair.
    * @param planetModel is the planetModel to put the point on.
@@ -46,7 +71,26 @@ public class GeoPoint extends Vector {
    * @param lon is the longitude.
    */
   public GeoPoint(final PlanetModel planetModel, final double lat, final double lon) {
-    this(planetModel, Math.sin(lat), Math.sin(lon), Math.cos(lat), Math.cos(lon));
+    this(planetModel, Math.sin(lat), Math.sin(lon), Math.cos(lat), Math.cos(lon), lat, lon);
+  }
+  
+  /** Construct a GeoPoint from a unit (x,y,z) vector and a magnitude.
+   * @param magnitude is the desired magnitude, provided to put the point on the ellipsoid.
+   * @param x is the unit x value.
+   * @param y is the unit y value.
+   * @param z is the unit z value.
+   * @param lat is the latitude.
+   * @param lon is the longitude.
+   */
+  public GeoPoint(final double magnitude, final double x, final double y, final double z, double lat, double lon) {
+    super(x * magnitude, y * magnitude, z * magnitude);
+    this.magnitude = magnitude;
+    if (lat > Math.PI * 0.5 || lat < -Math.PI * 0.5)
+      throw new IllegalArgumentException("Latitude out of range");
+    if (lon < -Math.PI || lon > Math.PI)
+      throw new IllegalArgumentException("Longitude out of range");
+    this.latitude = lat;
+    this.longitude = lon;
   }
 
   /** Construct a GeoPoint from a unit (x,y,z) vector and a magnitude.
@@ -71,6 +115,8 @@ public class GeoPoint extends Vector {
   }
 
   /** Compute an arc distance between two points.
+   * Note: this is an angular distance, and not a surface distance, and is therefore independent of planet model.
+   * For surface distance, see {@link org.apache.lucene.spatial.spatial4j.geo3d.PlanetModel#surfaceDistance(GeoPoint, GeoPoint)}
    * @param v is the second point.
    * @return the angle, in radians, between the two points.
    */
@@ -79,19 +125,27 @@ public class GeoPoint extends Vector {
   }
 
   /** Compute the latitude for the point.
-   *@return the latitude.
+   * @return the latitude.
    */
   public double getLatitude() {
-    return Math.asin(z / magnitude() );
+    double lat = this.latitude;//volatile-read once
+    if (lat == Double.NEGATIVE_INFINITY)
+      this.latitude = lat = Math.asin(z / magnitude());
+    return lat;
   }
   
   /** Compute the longitude for the point.
    * @return the longitude value.  Uses 0.0 if there is no computable longitude.
    */
   public double getLongitude() {
-    if (Math.abs(x) < MINIMUM_RESOLUTION && Math.abs(y) < MINIMUM_RESOLUTION)
+    double lon = this.longitude;//volatile-read once
-      return 0.0;
+    if (lon == Double.NEGATIVE_INFINITY) {
-    return Math.atan2(y,x);
+      if (Math.abs(x) < MINIMUM_RESOLUTION && Math.abs(y) < MINIMUM_RESOLUTION)
+        this.longitude = lon = 0.0;
+      else
+        this.longitude = lon = Math.atan2(y,x);
+    }
+    return lon;
   }
   
   /** Compute the linear magnitude of the point.
@@ -99,9 +153,10 @@ public class GeoPoint extends Vector {
    */
   @Override
   public double magnitude() {
-    if (this.magnitude == Double.NEGATIVE_INFINITY) {
+    double mag = this.magnitude;//volatile-read once
-      this.magnitude = super.magnitude();
+    if (mag == Double.NEGATIVE_INFINITY) {
+      this.magnitude = mag = super.magnitude();
     }
-    return magnitude;
+    return mag;
   }
 }

lucene/spatial/src/java/org/apache/lucene/spatial/spatial4j/geo3d/PlanetModel.java

@@ -29,7 +29,7 @@ public class PlanetModel {
   /** Mean radius */
   public static final double WGS84_MEAN = 6371009.0;
   /** Polar radius */
-  public static final double WGS84_POLAR = 6356752.3;
+  public static final double WGS84_POLAR = 6356752.314245;
   /** Equatorial radius */
   public static final double WGS84_EQUATORIAL = 6378137.0;
   /** Planet model corresponding to WGS84 */
@@ -45,6 +45,8 @@ public class PlanetModel {
   public final double inverseC;
   public final double inverseAbSquared;
   public final double inverseCSquared;
+  public final double flattening;
+  public final double squareRatio;
   // 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.
@@ -56,10 +58,12 @@ public class PlanetModel {
     this.c = c;
     this.inverseAb = 1.0 / ab;
     this.inverseC = 1.0 / c;
+    this.flattening = (ab - c) * inverseAb;
+    this.squareRatio = (ab * ab - c * c) / (c * c);
     this.inverseAbSquared = inverseAb * inverseAb;
     this.inverseCSquared = inverseC * inverseC;
-    this.NORTH_POLE = new GeoPoint(c, 0.0, 0.0, 1.0);
+    this.NORTH_POLE = new GeoPoint(c, 0.0, 0.0, 1.0, Math.PI * 0.5, 0.0);
-    this.SOUTH_POLE = new GeoPoint(c, 0.0, 0.0, -1.0);
+    this.SOUTH_POLE = new GeoPoint(c, 0.0, 0.0, -1.0, -Math.PI * 0.5, 0.0);
   }
   
   /** Find the minimum magnitude of all points on the ellipsoid.
@@ -74,6 +78,67 @@ public class PlanetModel {
     return Math.max(this.ab, this.c);
   }
   
+  /** Compute surface distance between two points.
+   * @param p1 is the first point.
+   * @param p2 is the second point.
+   * @return the adjusted angle, when multiplied by the mean earth radius, yields a surface distance.  This will differ
+   * from GeoPoint.arcDistance() only when the planet model is not a sphere. @see {@link org.apache.lucene.spatial.spatial4j.geo3d.GeoPoint#arcDistance(GeoPoint)}
+   */
+  public double surfaceDistance(final GeoPoint p1, final GeoPoint p2) {
+    final double latA = p1.getLatitude();
+    final double lonA = p1.getLongitude();
+    final double latB = p2.getLatitude();
+    final double lonB = p2.getLongitude();
+
+    final double L = lonB - lonA;
+    final double oF = 1.0 - this.flattening;
+    final double U1 = Math.atan(oF * Math.tan(latA));
+    final double U2 = Math.atan(oF * Math.tan(latB));
+    final double sU1 = Math.sin(U1);
+    final double cU1 = Math.cos(U1);
+    final double sU2 = Math.sin(U2);
+    final double cU2 = Math.cos(U2);
+
+    double sigma, sinSigma, cosSigma;
+    double cos2Alpha, cos2SigmaM;
+    
+    double lambda = L;
+    double iters = 100;
+
+    do {
+      final double sinLambda = Math.sin(lambda);
+      final double cosLambda = Math.cos(lambda);
+      sinSigma = Math.sqrt((cU2 * sinLambda) * (cU2 * sinLambda) + (cU1 * sU2 - sU1 * cU2 * cosLambda)
+          * (cU1 * sU2 - sU1 * cU2 * cosLambda));
+      if (Math.abs(sinSigma) < Vector.MINIMUM_RESOLUTION)
+        return 0.0;
+
+      cosSigma = sU1 * sU2 + cU1 * cU2 * cosLambda;
+      sigma = Math.atan2(sinSigma, cosSigma);
+      final double sinAlpha = cU1 * cU2 * sinLambda / sinSigma;
+      cos2Alpha = 1.0 - sinAlpha * sinAlpha;
+      cos2SigmaM = cosSigma - 2.0 * sU1 * sU2 / cos2Alpha;
+
+      final double c = this.flattening * 0.625 * cos2Alpha * (4.0 + this.flattening * (4.0 - 3.0 * cos2Alpha));
+      final double lambdaP = lambda;
+      lambda = L + (1.0 - c) * this.flattening * sinAlpha * (sigma + c * sinSigma * (cos2SigmaM + c * cosSigma *
+          (-1.0 + 2.0 * cos2SigmaM * cos2SigmaM)));
+      if (Math.abs(lambda - lambdaP) < Vector.MINIMUM_RESOLUTION)
+        break;
+    } while (--iters > 0);
+
+    if (iters == 0)
+      return 0.0;
+
+    final double uSq = cos2Alpha * this.squareRatio;
+    final double A = 1.0 + uSq * 0.00006103515625 * (4096.0 + uSq * (-768.0 + uSq * (320.0 - 175.0 * uSq)));
+    final double B = uSq * 0.0009765625 * (256.0 + uSq * (-128.0 + uSq * (74.0 - 47.0 * uSq)));
+    final double deltaSigma = B * sinSigma * (cos2SigmaM + B * 0.25 * (cosSigma * (-1.0 + 2.0 * cos2SigmaM * cos2SigmaM) - B * 0.16666666666666666666667 * cos2SigmaM
+            * (-3.0 + 4.0 * sinSigma * sinSigma) * (-3.0 + 4.0 * cos2SigmaM * cos2SigmaM)));
+
+    return this.c * A * (sigma - deltaSigma);
+  }
+
   @Override
   public boolean equals(final Object o) {
     if (!(o instanceof PlanetModel))

lucene/spatial/src/java/org/apache/lucene/spatial/spatial4j/geo3d/SidedPlane.java

@@ -91,6 +91,8 @@ public class SidedPlane extends Plane implements Membership {
     final Vector normalVector, final Vector point1, final Vector point2) {
     final Vector pointsVector = new Vector(point1.x - point2.x, point1.y - point2.y, point1.z - point2.z);
     final Vector newNormalVector = new Vector(normalVector, pointsVector).normalize();
+    if (newNormalVector == null)
+      return null;
     // To construct the plane, we now just need D, which is simply the negative of the evaluation of the circle normal vector at one of the points.
     return new SidedPlane(insidePoint, newNormalVector, -newNormalVector.dotProduct(point1));
   }

lucene/spatial/src/test/org/apache/lucene/spatial/spatial4j/Geo3dShapeRectRelationTestCase.java

@@ -24,6 +24,7 @@ import java.util.Random;
 import com.carrotsearch.randomizedtesting.RandomizedContext;
 import com.spatial4j.core.context.SpatialContext;
 import com.spatial4j.core.distance.DistanceUtils;
+import com.spatial4j.core.shape.Circle;
 import com.spatial4j.core.shape.Point;
 import org.apache.lucene.spatial.spatial4j.geo3d.Bounds;
 import org.apache.lucene.spatial.spatial4j.geo3d.GeoBBox;
@@ -173,17 +174,16 @@ public abstract class Geo3dShapeRectRelationTestCase extends RandomizedShapeTest
       protected Geo3dShape generateRandomShape(Point nearP) {
         final Point centerPoint = randomPoint();
         final int maxDistance = random().nextInt(160) + 20;
+        final Circle pointZone = ctx.makeCircle(centerPoint, maxDistance);
         final int vertexCount = random().nextInt(3) + 3;
         while (true) {
           final List<GeoPoint> geoPoints = new ArrayList<>();
           while (geoPoints.size() < vertexCount) {
-            final Point point = randomPoint();
+            final Point point = randomPointIn(pointZone);
-            if (ctx.getDistCalc().distance(point,centerPoint) > maxDistance)
-              continue;
             final GeoPoint gPt = new GeoPoint(planetModel, point.getY() * DEGREES_TO_RADIANS, point.getX() * DEGREES_TO_RADIANS);
             geoPoints.add(gPt);
           }
-          final int convexPointIndex = random().nextInt(vertexCount);       //If we get this wrong, hopefully we get IllegalArgumentException
+          final int convexPointIndex = random().nextInt(vertexCount); //If we get this wrong, hopefully we get IllegalArgumentException
           try {
             final GeoShape shape = GeoPolygonFactory.makeGeoPolygon(planetModel, geoPoints, convexPointIndex);
             return new Geo3dShape(planetModel, shape, ctx);
@@ -217,18 +217,15 @@ public abstract class Geo3dShapeRectRelationTestCase extends RandomizedShapeTest
       protected Geo3dShape generateRandomShape(Point nearP) {
         final Point centerPoint = randomPoint();
         final int maxDistance = random().nextInt(160) + 20;
+        final Circle pointZone = ctx.makeCircle(centerPoint, maxDistance);
         final int pointCount = random().nextInt(5) + 1;
         final double width = (random().nextInt(89)+1) * DEGREES_TO_RADIANS;
         while (true) {
           try {
             final GeoPath path = new GeoPath(planetModel, width);
-            int i = 0;
+            while (path.points.size() < pointCount) {
-            while (i < pointCount) {
+              final Point nextPoint = randomPointIn(pointZone);
-              final Point nextPoint = randomPoint();
-              if (ctx.getDistCalc().distance(nextPoint,centerPoint) > maxDistance)
-                continue;
               path.addPoint(nextPoint.getY() * DEGREES_TO_RADIANS, nextPoint.getX() * DEGREES_TO_RADIANS);
-              i++;
             }
             path.done();
             return new Geo3dShape(planetModel, path, ctx);

lucene/spatial/src/test/org/apache/lucene/spatial/spatial4j/geo3d/GeoBBoxTest.java

@@ -53,6 +53,7 @@ public class GeoBBoxTest {
   public void testBBoxPointWithin() {
     GeoBBox box;
     GeoPoint gp;
+
     // Standard normal Rect box, not crossing dateline
     box = GeoBBoxFactory.makeGeoBBox(PlanetModel.SPHERE, 0.0, -Math.PI * 0.25, -1.0, 1.0);
     gp = new GeoPoint(PlanetModel.SPHERE, -0.1, 0.0);
@@ -65,6 +66,7 @@ public class GeoBBoxTest {
     assertFalse(box.isWithin(gp));
     gp = new GeoPoint(PlanetModel.SPHERE, -0.1, -1.1);
     assertFalse(box.isWithin(gp));
+
     // Standard normal Rect box, crossing dateline
     box = GeoBBoxFactory.makeGeoBBox(PlanetModel.SPHERE, 0.0, -Math.PI * 0.25, Math.PI - 1.0, -Math.PI + 1.0);
     gp = new GeoPoint(PlanetModel.SPHERE, -0.1, -Math.PI);
@@ -75,8 +77,9 @@ public class GeoBBoxTest {
     assertFalse(box.isWithin(gp));
     gp = new GeoPoint(PlanetModel.SPHERE, -0.1, -Math.PI + 1.1);
     assertFalse(box.isWithin(gp));
-    gp = new GeoPoint(PlanetModel.SPHERE, -0.1, -Math.PI - 1.1);
+    //bad lon: gp = new GeoPoint(PlanetModel.SPHERE, -0.1, -Math.PI - 1.1);
-    assertFalse(box.isWithin(gp));
+    //assertFalse(box.isWithin(gp));
+
     // Latitude zone rectangle
     box = GeoBBoxFactory.makeGeoBBox(PlanetModel.SPHERE, 0.0, -Math.PI * 0.25, -Math.PI, Math.PI);
     gp = new GeoPoint(PlanetModel.SPHERE, -0.1, -Math.PI);
@@ -87,8 +90,9 @@ public class GeoBBoxTest {
     assertFalse(box.isWithin(gp));
     gp = new GeoPoint(PlanetModel.SPHERE, -0.1, -Math.PI + 1.1);
     assertTrue(box.isWithin(gp));
-    gp = new GeoPoint(PlanetModel.SPHERE, -0.1, -Math.PI - 1.1);
+    //bad lon: gp = new GeoPoint(PlanetModel.SPHERE, -0.1, -Math.PI - 1.1);
-    assertTrue(box.isWithin(gp));
+    //assertTrue(box.isWithin(gp));
+
     // World
     box = GeoBBoxFactory.makeGeoBBox(PlanetModel.SPHERE, Math.PI * 0.5, -Math.PI * 0.5, -Math.PI, Math.PI);
     gp = new GeoPoint(PlanetModel.SPHERE, -0.1, -Math.PI);
@@ -99,8 +103,8 @@ public class GeoBBoxTest {
     assertTrue(box.isWithin(gp));
     gp = new GeoPoint(PlanetModel.SPHERE, -0.1, -Math.PI + 1.1);
     assertTrue(box.isWithin(gp));
-    gp = new GeoPoint(PlanetModel.SPHERE, -0.1, -Math.PI - 1.1);
+    //bad lat: gp = new GeoPoint(PlanetModel.SPHERE, -0.1, -Math.PI - 1.1);
-    assertTrue(box.isWithin(gp));
+    //assertTrue(box.isWithin(gp));
 
   }
 

lucene/spatial/src/test/org/apache/lucene/spatial/spatial4j/geo3d/GeoConvexPolygonTest.java

@@ -34,7 +34,7 @@ public class GeoConvexPolygonTest {
     c.addPoint(0.0, -0.6, false);
     c.addPoint(0.1, -0.5, false);
     c.addPoint(0.0, -0.4, false);
-    c.donePoints(false);
+    c.done(false);
     // Sample some points within
     gp = new GeoPoint(PlanetModel.SPHERE, 0.0, -0.5);
     assertTrue(c.isWithin(gp));
@@ -73,7 +73,7 @@ public class GeoConvexPolygonTest {
     c.addPoint(0.0, -0.6, false);
     c.addPoint(0.1, -0.5, false);
     c.addPoint(0.0, -0.4, false);
-    c.donePoints(false);
+    c.done(false);
 
     b = c.getBounds(null);
     assertFalse(b.checkNoLongitudeBound());

lucene/spatial/src/test/org/apache/lucene/spatial/spatial4j/geo3d/GeoModelTest.java

@@ -57,10 +57,10 @@ public class GeoModelTest {
     assertFalse(bounds.checkNoTopLatitudeBound());
     assertFalse(bounds.checkNoLongitudeBound());
     assertFalse(bounds.checkNoBottomLatitudeBound());
-    assertEquals(1.20, bounds.getMaxLatitude(), 0.01);
+    assertEquals(Math.PI * 0.25 + 0.01, bounds.getMaxLatitude(), 0.00001);
-    assertEquals(Math.PI * 0.25 - 0.01, bounds.getMinLatitude(), 0.01);
+    assertEquals(Math.PI * 0.25 - 0.01, bounds.getMinLatitude(), 0.00001);
-    assertEquals(-0.36, bounds.getLeftLongitude(), 0.01);
+    assertEquals(-0.0125, bounds.getLeftLongitude(), 0.0001);
-    assertEquals(0.36, bounds.getRightLongitude(), 0.01);
+    assertEquals(0.0125, bounds.getRightLongitude(), 0.0001);
 
     circle = new GeoCircle(scaledModel, Math.PI * 0.125, 0.0, 0.01);
     assertTrue(circle.isWithin(point2));
@@ -70,11 +70,11 @@ public class GeoModelTest {
     assertFalse(bounds.checkNoLongitudeBound());
     assertFalse(bounds.checkNoTopLatitudeBound());
     assertFalse(bounds.checkNoBottomLatitudeBound());
-    // Asymmetric, as expected
+    // Symmetric, as expected
-    assertEquals(Math.PI * 0.125 - 0.01, bounds.getMinLatitude(), 0.01);
+    assertEquals(Math.PI * 0.125 - 0.01, bounds.getMinLatitude(), 0.00001);
-    assertEquals(0.74, bounds.getMaxLatitude(), 0.01);
+    assertEquals(Math.PI * 0.125 + 0.01, bounds.getMaxLatitude(), 0.00001);
-    assertEquals(-0.18, bounds.getLeftLongitude(), 0.01);
+    assertEquals(-0.0089, bounds.getLeftLongitude(), 0.0001);
-    assertEquals(0.18, bounds.getRightLongitude(), 0.01);
+    assertEquals(0.0089, bounds.getRightLongitude(), 0.0001);
 
   }
 

lucene/spatial/src/test/org/apache/lucene/spatial/spatial4j/geo3d/GeoPointTest.java

@@ -30,10 +30,10 @@ public class GeoPointTest extends LuceneTestCase {
 
   @Test
   public void testConversion() {
-    testPointRoundTrip(PlanetModel.SPHERE, 90, 0, 1e-12);
+    testPointRoundTrip(PlanetModel.SPHERE, 90 * DistanceUtils.DEGREES_TO_RADIANS, 0, 1e-6);
-    testPointRoundTrip(PlanetModel.SPHERE, -90, 0, 1e-12);
+    testPointRoundTrip(PlanetModel.SPHERE, -90 * DistanceUtils.DEGREES_TO_RADIANS, 0, 1e-6);
-    testPointRoundTrip(PlanetModel.WGS84, 90, 0, 1e-12);
+    testPointRoundTrip(PlanetModel.WGS84, 90 * DistanceUtils.DEGREES_TO_RADIANS, 0, 1e-6);
-    testPointRoundTrip(PlanetModel.WGS84, -90, 0, 1e-12);
+    testPointRoundTrip(PlanetModel.WGS84, -90 * DistanceUtils.DEGREES_TO_RADIANS, 0, 1e-6);
 
     final int times = atLeast(100);
     for (int i = 0; i < times; i++) {
@@ -46,11 +46,35 @@ public class GeoPointTest extends LuceneTestCase {
 
   protected void testPointRoundTrip(PlanetModel planetModel, double pLat, double pLon, double epsilon) {
     final GeoPoint p1 = new GeoPoint(planetModel, pLat, pLon);
-    final GeoPoint p2 = new GeoPoint(planetModel, p1.getLatitude(), p1.getLongitude());
+    // In order to force the reverse conversion, we have to construct a geopoint from just x,y,z
-    double dist = p1.arcDistance(p2);
+    final GeoPoint p2 = new GeoPoint(p1.x, p1.y, p1.z);
+    // Now, construct the final point based on getLatitude() and getLongitude()
+    final GeoPoint p3 = new GeoPoint(planetModel, p2.getLatitude(), p2.getLongitude());
+    double dist = p1.arcDistance(p3);
     assertEquals(0, dist, epsilon);
   }
 
+  @Test
+  public void testSurfaceDistance() {
+    final int times = atLeast(100);
+    for (int i = 0; i < times; i++) {
+      final double p1Lat = (randomFloat() * 180.0 - 90.0) * DistanceUtils.DEGREES_TO_RADIANS;
+      final double p1Lon = (randomFloat() * 360.0 - 180.0) * DistanceUtils.DEGREES_TO_RADIANS;
+      final double p2Lat = (randomFloat() * 180.0 - 90.0) * DistanceUtils.DEGREES_TO_RADIANS;
+      final double p2Lon = (randomFloat() * 360.0 - 180.0) * DistanceUtils.DEGREES_TO_RADIANS;
+      final GeoPoint p1 = new GeoPoint(PlanetModel.SPHERE, p1Lat, p1Lon);
+      final GeoPoint p2 = new GeoPoint(PlanetModel.SPHERE, p2Lat, p2Lon);
+      final double arcDistance = p1.arcDistance(p2);
+      // Compute ellipsoid distance; it should agree for a sphere
+      final double surfaceDistance = PlanetModel.SPHERE.surfaceDistance(p1,p2);
+      assertEquals(arcDistance, surfaceDistance, 1e-6);
+    }
+  }
+
+  @Test(expected = IllegalArgumentException.class)
+  public void testBadLatLon() {
+    new GeoPoint(PlanetModel.SPHERE, 50.0, 32.2);
+  }
 }