use Vincenty's formula in spoint_dist(): this vindicates a fixed ./expected/circle_extended.out

Author: Markus Nullmeier

expected/circle.out

@@ -165,9 +165,9 @@ SELECT 180.0*dist('<( 23h 30m 00s , 0d 0m 0s), 1.0d>'::scircle,'<( 1h 0m 00s , 0
 (1 row)
 
 SELECT 180.0*dist('<( 0h 40m 00s , 0d 0m 0s), 1.0d>'::scircle,'<( 0h 50m 00s , 0d 0m 0s),1.0d>'::scircle)/pi();
-     ?column?      
--------------------
- 0.499999999999953
+ ?column? 
+----------
+      0.5
 (1 row)
 
 SELECT 180.0*dist('<( 0h 40m 00s , 0d 0m 0s), 1.5d>'::scircle,'<( 0h 50m 00s , 0d 0m 0s),1.5d>'::scircle)/pi();
@@ -407,9 +407,9 @@ SELECT 180.0*('<( 23h 30m 00s , 0d 0m 0s), 1.0d>'::scircle<->'<( 1h 0m 00s , 0d
 (1 row)
 
 SELECT 180.0*('<( 0h 40m 00s , 0d 0m 0s), 1.0d>'::scircle<->'<( 0h 50m 00s , 0d 0m 0s),1.0d>'::scircle)/pi();
-     ?column?      
--------------------
- 0.499999999999953
+ ?column? 
+----------
+      0.5
 (1 row)
 
 SELECT 180.0*('<( 0h 40m 00s , 0d 0m 0s), 1.5d>'::scircle<->'<( 0h 50m 00s , 0d 0m 0s),1.5d>'::scircle)/pi();

expected/circle_extended.out

@@ -305,6 +305,6 @@ select count(sc) from scircle_data where '<(0d,0d),2.1d>'::scircle ~ sc;
 select count(spoint_data.sp) from spoint_data,scircle_data where spoint_data.sp @ scircle_data.sc;
  count 
 -------
- 76682
+ 78842
 (1 row)
 

expected/points.out

@@ -261,7 +261,7 @@ SELECT spoint(0.0 , 2.141592653589793116);
 SELECT dist('( 0h 2m 30s , 0d 0m 0s)'::spoint,'( 0h 0m 30s , 0d 0m 0s)'::spoint);
         dist         
 ---------------------
- 0.00872664625996925
+ 0.00872664625997165
 (1 row)
 
 SELECT dist('( 0h 2m 30s , 0d 0m 0s)'::spoint,'( 0h 2m 30s , 10d 0m 0s)'::spoint);
@@ -609,7 +609,7 @@ SELECT '( 0h 2m 30s , 0d 0m 0s)'::spoint<>'( 12h 2m 30s , 45d 0m 0s)'::spoint;
 SELECT '( 0h 2m 30s , 0d 0m 0s)'::spoint<->'( 0h 0m 30s , 0d 0m 0s)'::spoint;
       ?column?       
 ---------------------
- 0.00872664625996925
+ 0.00872664625997165
 (1 row)
 
 SELECT '( 0h 2m 30s , 0d 0m 0s)'::spoint<->'( 0h 2m 30s , 10d 0m 0s)'::spoint;

point.c

@@ -158,27 +158,22 @@ spherepoint_from_long_lat(PG_FUNCTION_ARGS)
 	PG_RETURN_POINTER(p);
 }
 
+static double
+norm2(double a, double b)
+{
+	return sqrt(a * a + b * b);
+}
+
 float8
 spoint_dist(const SPoint *p1, const SPoint *p2)
 {
 	float8		dl = p1->lng - p2->lng;
-	float8		f = ((sin(p1->lat) * sin(p2->lat) +
-					  cos(p1->lat) * cos(p2->lat) * cos(dl)));
-
-	if (FPeq(f, 1.0))
-	{
-		/* for small distances */
-		Vector3D	v1,	v2,	v3;
-
-		spoint_vector3d(&v1, p1);
-		spoint_vector3d(&v2, p2);
-		vector3d_cross(&v3, &v1, &v2);
-		f = vector3d_length(&v3);
-	}
-	else
-	{
-		f = acos(f);
-	}
+	/* use Vincenty's formula for the inverse geodesic problem on the sphere */
+	float8 f = atan2(norm2(cos(p2->lat) * sin(dl),
+							 cos(p1->lat) * sin(p2->lat)
+						   - sin(p1->lat) * cos(p2->lat) * cos(dl)),
+					   sin(p1->lat) * sin(p2->lat)
+					 + cos(p1->lat) * cos(p2->lat) * cos(dl));
 	if (FPzero(f))
 	{
 		return 0.0;