applied patricia tozer's patch (distance function was taking acos of something

just slightly outside [-1,1]).

git-svn-id: http://svn.osgeo.org/postgis/trunk@439 b70326c6-7e19-0410-871a-916f4a2858ee

diff --git a/postgis_proj.c b/postgis_proj.c
index d76f3db8c6a61fa991ca1c55c44abfb83428a9e7..e3cf680456f5cba8a2664dd0ec22b1e33115e785 100644 (file)
--- a/postgis_proj.c
+++ b/postgis_proj.c
@@ -8,9 +8,13 @@
  *
  * This is free software; you can redistribute and/or modify it under
  * the terms of hte GNU General Public Licence. See the COPYING file.
- * 
+ *
  **********************************************************************
  * $Log$
+ * Revision 1.9  2004/02/04 02:53:20  dblasby
+ * applied patricia tozer's patch (distance function was taking acos of something
+ * just slightly outside [-1,1]).
+ *
  * Revision 1.8  2003/12/04 18:58:35  dblasby
  * changed david skae to skea
  *
@@ -50,9 +54,9 @@
 
 
 
-//use the WKT definition of an ellipsoid 
-// ie. SPHEROID["name",A,rf] or SPHEROID("name",A,rf) 
-//       SPHEROID["GRS_1980",6378137,298.257222101] 
+//use the WKT definition of an ellipsoid
+// ie. SPHEROID["name",A,rf] or SPHEROID("name",A,rf)
+//       SPHEROID["GRS_1980",6378137,298.257222101]
 // wkt says you can use "(" or "["
 
 PG_FUNCTION_INFO_V1(ellipsoid_in);
@@ -66,11 +70,11 @@ Datum ellipsoid_in(PG_FUNCTION_ARGS)
 
        memset(sphere,0, sizeof(SPHEROID));
 
-       if (strstr(str,"SPHEROID") !=  str ) 
+       if (strstr(str,"SPHEROID") !=  str )
        {
                 elog(ERROR,"SPHEROID parser - doesnt start with SPHEROID");
                 pfree(sphere);
-                PG_RETURN_NULL();      
+                PG_RETURN_NULL();
        }
 
        nitems = sscanf(str,"SPHEROID[\"%19[^\"]\",%lf,%lf]",sphere->name,&sphere->a,&rf);
@@ -91,7 +95,7 @@ Datum ellipsoid_in(PG_FUNCTION_ARGS)
        sphere->e = sqrt(sphere->e_sq);
 
        PG_RETURN_POINTER(sphere);
-       
+
 }
 
 
@@ -99,7 +103,7 @@ Datum ellipsoid_in(PG_FUNCTION_ARGS)
 PG_FUNCTION_INFO_V1(ellipsoid_out);
 Datum ellipsoid_out(PG_FUNCTION_ARGS)
 {
-       SPHEROID  *sphere = (SPHEROID *) PG_GETARG_POINTER(0);  
+       SPHEROID  *sphere = (SPHEROID *) PG_GETARG_POINTER(0);
        char    *result;
 
        result = palloc(MAX_DIGS_DOUBLE + MAX_DIGS_DOUBLE + 20 + 9 + 2);
@@ -123,12 +127,12 @@ double deltaLongitude(double azimuth, double sigma, double tsm,SPHEROID *sphere)
        das = cos(azimuth)*cos(azimuth);
        C = sphere->f/16.0 * das * (4.0 + sphere->f * (4.0 - 3.0 * das));
        // compute the difference in longitude
-       
+
        ctsm = cos(tsm);
        DL = ctsm + C * cos(sigma) * (-1.0 + 2.0 * ctsm*ctsm);
        DL = sigma + C * sin(sigma) * DL;
        return (1.0 - C) * sphere->f * sin(azimuth) * DL;
-}  
+}
 
 
 //support function for distance calc
@@ -142,7 +146,7 @@ double mu2(double azimuth,SPHEROID *sphere)
 
        e2 = sqrt(sphere->a*sphere->a-sphere->b*sphere->b)/sphere->b;
        return cos(azimuth)*cos(azimuth) * e2*e2;
-}  
+}
 
 
 //support function for distance calc
@@ -153,7 +157,7 @@ double mu2(double azimuth,SPHEROID *sphere)
 double bigA(double u2)
 {
        return 1.0 + u2/256.0 * (64.0 + u2 * (-12.0 + 5.0 * u2));
-}  
+}
 
 
 //support function for distance calc
@@ -164,7 +168,7 @@ double bigA(double u2)
 double bigB(double u2)
 {
        return u2/512.0 * (128.0 + u2 * (-64.0 + 37.0 * u2));
-}  
+}
 
 
 //given 2 lat/longs and ellipse, find the distance
@@ -184,7 +188,9 @@ double      distance_ellipse(double lat1, double long1,
        double u2,A,B;
        double dsigma;
 
-       int iterations;  
+       double TEMP;
+
+       int iterations;
 
        L1 = atan((1.0 - sphere->f ) * tan( lat1) );
        L2 = atan((1.0 - sphere->f ) * tan( lat2) );
@@ -202,7 +208,21 @@ double     distance_ellipse(double lat1, double long1,
                cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosdl1;
                sigma = acos(cosSigma);
                azimuthEQ = asin((cosU1 * cosU2 * sindl1)/sin(sigma));
-               tsm = acos(cosSigma - (2.0 * sinU1 * sinU2)/(cos(azimuthEQ)*cos(azimuthEQ)));
+
+                       // patch from patrica tozer to handle minor mathematical stability problem
+               TEMP = cosSigma - (2.0 * sinU1 * sinU2)/(cos(azimuthEQ)*cos(azimuthEQ));
+               if(TEMP > 1)
+               {
+                          TEMP = 1;
+               }
+               else if(TEMP < -1)
+               {
+                          TEMP = -1;
+               }
+        tsm = acos(TEMP);
+
+
+               //tsm = acos(cosSigma - (2.0 * sinU1 * sinU2)/(cos(azimuthEQ)*cos(azimuthEQ)));
                dl2 = deltaLongitude(azimuthEQ, sigma, tsm,sphere);
                dl3 = dl1 - (dl + dl2);
                dl1 = dl + dl2;
@@ -211,12 +231,12 @@ double    distance_ellipse(double lat1, double long1,
                iterations++;
        } while ( (iterations<999) && (fabs(dl3) > 1.0e-032));
 
-          // compute expansions A and B 
+          // compute expansions A and B
        u2 = mu2(azimuthEQ,sphere);
        A = bigA(u2);
        B = bigB(u2);
 
-       // compute length of geodesic 
+       // compute length of geodesic
        dsigma = B * sin(sigma) * (cos(tsm) + (B*cosSigma*(-1.0 + 2.0 * (cos(tsm)*cos(tsm))))/4.0);
        return sphere->b * (A * (sigma - dsigma));
 }
@@ -237,7 +257,7 @@ double length2d_ellipse_linestring(LINE3D   *line, SPHEROID         *sphere)
        for (i=1; i<line->npoints;i++)
        {
                to = &line->points[i];
-       
+
                dist +=  distance_ellipse(frm->y*M_PI/180.0 , frm->x*M_PI/180.0,
                                                to->y*M_PI/180.0 , to->x*M_PI/180.0,
                                                sphere);
@@ -262,7 +282,7 @@ double length3d_ellipse_linestring(LINE3D   *line, SPHEROID         *sphere)
        for (i=1; i<line->npoints;i++)
        {
                to = &line->points[i];
-               
+
                dist_ellipse =  distance_ellipse(frm->y*M_PI/180.0 , frm->x*M_PI/180.0,
                                                to->y*M_PI/180.0 , to->x*M_PI/180.0,
                                                sphere);
@@ -281,7 +301,7 @@ double length3d_ellipse_linestring(LINE3D   *line, SPHEROID         *sphere)
 // find the "length of a geometry"
 // length2d(point) = 0
 // length2d(line) = length of line
-// length2d(polygon) = 0  
+// length2d(polygon) = 0
 // uses ellipsoidal math to find the distance
 //// x's are longitude, and y's are latitude - both in decimal degrees
 
@@ -289,8 +309,8 @@ PG_FUNCTION_INFO_V1(length_ellipsoid);
 Datum length_ellipsoid(PG_FUNCTION_ARGS)
 {
                GEOMETRY                      *geom = (GEOMETRY *)  PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
-               SPHEROID                *sphere = (SPHEROID *) PG_GETARG_POINTER(1);    
-       
+               SPHEROID                *sphere = (SPHEROID *) PG_GETARG_POINTER(1);
+
        int32                           *offsets1;
        char                            *o1;
        int32                           type1,j;
@@ -304,7 +324,7 @@ Datum length_ellipsoid(PG_FUNCTION_ARGS)
 
        for (j=0; j< geom->nobjs; j++)          //for each object in geom1
        {
-               o1 = (char *) geom +offsets1[j] ;  
+               o1 = (char *) geom +offsets1[j] ;
                type1=  geom->objType[j];
                if (type1 == LINETYPE)  //LINESTRING
                {
@@ -322,7 +342,7 @@ Datum length_ellipsoid(PG_FUNCTION_ARGS)
 // find the "length of a geometry"
 // length3d(point) = 0
 // length3d(line) = length of line
-// length3d(polygon) = 0  
+// length3d(polygon) = 0
 // uses ellipsoidal math to find the distance on the XY plane, then
 //  uses simple Pythagoras theorm to find the 3d distance on each segment
 // x's are longitude, and y's are latitude - both in decimal degrees
@@ -331,8 +351,8 @@ PG_FUNCTION_INFO_V1(length3d_ellipsoid);
 Datum length3d_ellipsoid(PG_FUNCTION_ARGS)
 {
                GEOMETRY                      *geom = (GEOMETRY *)  PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
-               SPHEROID                *sphere = (SPHEROID *) PG_GETARG_POINTER(1);    
-       
+               SPHEROID                *sphere = (SPHEROID *) PG_GETARG_POINTER(1);
+
        int32                           *offsets1;
        char                            *o1;
        int32                           type1,j;
@@ -346,7 +366,7 @@ Datum length3d_ellipsoid(PG_FUNCTION_ARGS)
 
        for (j=0; j< geom->nobjs; j++)          //for each object in geom1
        {
-               o1 = (char *) geom +offsets1[j] ;  
+               o1 = (char *) geom +offsets1[j] ;
                type1=  geom->objType[j];
                if (type1 == LINETYPE)  //LINESTRING
                {
@@ -362,7 +382,7 @@ Datum length3d_ellipsoid(PG_FUNCTION_ARGS)
 PG_FUNCTION_INFO_V1(distance_ellipsoid);
 Datum distance_ellipsoid(PG_FUNCTION_ARGS)
 {
-               SPHEROID                *sphere = (SPHEROID *) PG_GETARG_POINTER(2);    
+               SPHEROID                *sphere = (SPHEROID *) PG_GETARG_POINTER(2);
                GEOMETRY                      *geom1 = (GEOMETRY *)  PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
                GEOMETRY                      *geom2 = (GEOMETRY *)  PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
 
@@ -377,7 +397,7 @@ Datum distance_ellipsoid(PG_FUNCTION_ARGS)
                elog(ERROR,"optimistic_overlap:Operation on two GEOMETRIES with different SRIDs\n");
                PG_RETURN_NULL();
        }
-       
+
        if (geom1->type != POINTTYPE)
        {
                elog(ERROR,"optimistic_overlap: first arg isnt a point\n");