|

# earthdistance.c

Go to the documentation of this file.
```00001 /* contrib/earthdistance/earthdistance.c */
00002
00003 #include "postgres.h"
00004
00005 #include <math.h>
00006
00007 #include "utils/geo_decls.h"    /* for Point */
00008
00009 #ifndef M_PI
00010 #define M_PI 3.14159265358979323846
00011 #endif
00012
00013
00014 PG_MODULE_MAGIC;
00015
00016 /* Earth's radius is in statute miles. */
00017 static const double EARTH_RADIUS = 3958.747716;
00018 static const double TWO_PI = 2.0 * M_PI;
00019
00020
00021 /******************************************************
00022  *
00023  * degtorad - convert degrees to radians
00024  *
00025  * arg: double, angle in degrees
00026  *
00027  * returns: double, same angle in radians
00028  ******************************************************/
00029
00030 static double
00032 {
00033     return (degrees / 360.0) * TWO_PI;
00034 }
00035
00036 /******************************************************
00037  *
00038  * geo_distance_internal - distance between points
00039  *
00040  * args:
00041  *   a pair of points - for each point,
00042  *     x-coordinate is longitude in degrees west of Greenwich
00043  *     y-coordinate is latitude in degrees above equator
00044  *
00045  * returns: double
00046  *   distance between the points in miles on earth's surface
00047  ******************************************************/
00048
00049 static double
00050 geo_distance_internal(Point *pt1, Point *pt2)
00051 {
00052     double      long1,
00053                 lat1,
00054                 long2,
00055                 lat2;
00056     double      longdiff;
00057     double      sino;
00058
00059     /* convert degrees to radians */
00060
00061     long1 = degtorad(pt1->x);
00062     lat1 = degtorad(pt1->y);
00063
00064     long2 = degtorad(pt2->x);
00065     lat2 = degtorad(pt2->y);
00066
00067     /* compute difference in longitudes - want < 180 degrees */
00068     longdiff = fabs(long1 - long2);
00069     if (longdiff > M_PI)
00070         longdiff = TWO_PI - longdiff;
00071
00072     sino = sqrt(sin(fabs(lat1 - lat2) / 2.) * sin(fabs(lat1 - lat2) / 2.) +
00073             cos(lat1) * cos(lat2) * sin(longdiff / 2.) * sin(longdiff / 2.));
00074     if (sino > 1.)
00075         sino = 1.;
00076
00077     return 2. * EARTH_RADIUS * asin(sino);
00078 }
00079
00080
00081 /******************************************************
00082  *
00083  * geo_distance - distance between points
00084  *
00085  * args:
00086  *   a pair of points - for each point,
00087  *     x-coordinate is longitude in degrees west of Greenwich
00088  *     y-coordinate is latitude in degrees above equator
00089  *
00090  * returns: float8
00091  *   distance between the points in miles on earth's surface
00092  *
00093  * If float8 is passed-by-value, the oldstyle version-0 calling convention
00094  * is unportable, so we use version-1.  However, if it's passed-by-reference,
00095  * continue to use oldstyle.  This is just because we'd like earthdistance
00096  * to serve as a canary for any unintentional breakage of version-0 functions
00097  * with float8 results.
00098  ******************************************************/
00099
00100 #ifdef USE_FLOAT8_BYVAL
00101
00102 Datum       geo_distance(PG_FUNCTION_ARGS);
00103
00104 PG_FUNCTION_INFO_V1(geo_distance);
00105
00106 Datum
00107 geo_distance(PG_FUNCTION_ARGS)
00108 {
00109     Point      *pt1 = PG_GETARG_POINT_P(0);
00110     Point      *pt2 = PG_GETARG_POINT_P(1);
00111     float8      result;
00112
00113     result = geo_distance_internal(pt1, pt2);
00114     PG_RETURN_FLOAT8(result);
00115 }
00116 #else                           /* !USE_FLOAT8_BYVAL */
00117
00118 double     *geo_distance(Point *pt1, Point *pt2);
00119
00120 double *
00121 geo_distance(Point *pt1, Point *pt2)
00122 {
00123     double     *resultp = palloc(sizeof(double));
00124
00125     *resultp = geo_distance_internal(pt1, pt2);
00126     return resultp;
00127 }
00128
00129 #endif   /* USE_FLOAT8_BYVAL */
```