I have included fixes to declare some floating point constants as double

bmomjian · bmomjian · commit 24bebf0b725e · 2002-09-20T03:47:22.000Z
instead of int, change the calculation method to use the haversine
formula
which is more accurrate for short distances, added a grant to public for
geo_distance and added a regression test.

I will resubmit the earth distance stuff based on cube after 7.3 is
released.

Bruno Wolff III
diff --git a/contrib/earthdistance/Makefile b/contrib/earthdistance/Makefile
@@ -1,4 +1,4 @@
-# $Header: /cvsroot/pgsql/contrib/earthdistance/Makefile,v 1.11 2001/09/06 10:49:29 petere Exp $
+# $Header: /cvsroot/pgsql/contrib/earthdistance/Makefile,v 1.12 2002/09/20 03:47:22 momjian Exp $
 
 subdir = contrib/earthdistance
 top_builddir = ../..
@@ -7,5 +7,6 @@ include $(top_builddir)/src/Makefile.global
 MODULES = earthdistance
 DATA_built = earthdistance.sql
 DOCS = README.earthdistance
+REGRESS = earthdistance
 
 include $(top_srcdir)/contrib/contrib-global.mk
diff --git a/contrib/earthdistance/README.earthdistance b/contrib/earthdistance/README.earthdistance
@@ -1,3 +1,13 @@
+---------------------------------------------------------------------
+I corrected a bug in the geo_distance code where two double constants
+were declared as int. I changed the distance function to use the
+haversine formula which is more accurate for small distances.
+I added a regression test to the package. I added a grant statement
+to give execute access for geo_distance to public.
+
+Bruno Wolff III
+September 2002
+---------------------------------------------------------------------
 Date: Wed, 1 Apr 1998 15:19:32 -0600 (CST)
 From: Hal Snyder <hal@vailsys.com>
 To: vmehr@ctp.com
diff --git a/contrib/earthdistance/earthdistance.c b/contrib/earthdistance/earthdistance.c
@@ -6,8 +6,8 @@
 
 
 /* Earth's radius is in statute miles. */
-const int	EARTH_RADIUS = 3958.747716;
-const int	TWO_PI = 2.0 * M_PI;
+const double	EARTH_RADIUS = 3958.747716;
+const double	TWO_PI = 2.0 * M_PI;
 
 double	   *geo_distance(Point *pt1, Point *pt2);
 
@@ -50,6 +50,7 @@ geo_distance(Point *pt1, Point *pt2)
 				long2,
 				lat2;
 	double		longdiff;
+        double		sino;
 	double	   *resultp = palloc(sizeof(double));
 
 	/* convert degrees to radians */
@@ -65,8 +66,10 @@ geo_distance(Point *pt1, Point *pt2)
 	if (longdiff > M_PI)
 		longdiff = TWO_PI - longdiff;
 
-	*resultp = EARTH_RADIUS * acos
-		(sin(lat1) * sin(lat2) + cos(lat1) * cos(lat2) * cos(longdiff));
+        sino = sqrt(sin(fabs(lat1-lat2)/2.)*sin(fabs(lat1-lat2)/2.) +
+                cos(lat1) * cos(lat2) * sin(longdiff/2.)*sin(longdiff/2.));
+        if (sino > 1.) sino = 1.;
+        *resultp = 2. * EARTH_RADIUS * asin(sino);
 
 	return resultp;
 }
diff --git a/contrib/earthdistance/earthdistance.sql.in b/contrib/earthdistance/earthdistance.sql.in
@@ -1,24 +1,24 @@
+begin;
 
 --------------- geo_distance
 
-DROP FUNCTION geo_distance (point, point);
-CREATE FUNCTION geo_distance (point, point) RETURNS float8
-  AS 'MODULE_PATHNAME' LANGUAGE 'c'
-  WITH (isstrict);
-
-SELECT geo_distance ('(1,2)'::point, '(3,4)'::point);
+CREATE OR REPLACE FUNCTION geo_distance (point, point) RETURNS float8
+  LANGUAGE 'c' IMMUTABLE STRICT AS 'MODULE_PATHNAME';
 
 --------------- geo_distance as operator <@>
 
-DROP OPERATOR <@> (point, point);
 CREATE OPERATOR <@> (
   leftarg = point,
   rightarg = point,
   procedure = geo_distance,
   commutator = <@>
 );
 
--- ( 87.6, 41.8) is in Chicago
--- (106.7, 35.1) is in Albuquerque
--- The cities are about 1100 miles apart
-SELECT '(87.6,41.8)'::point <@> '(106.7,35.1)'::point;
+--
+-- By default this function is made executable by anyone. To restrict
+-- access by default, comment out the following grant command.
+--
+
+grant execute on function geo_distance(point, point) to public;
+
+commit;
