Line data    Source code

       1              : /* contrib/earthdistance/earthdistance.c */
       2              : 
       3              : #include "postgres.h"
       4              : 
       5              : #include <math.h>
       6              : 
       7              : #include "utils/geo_decls.h"  /* for Point */
       8              : 
       9              : /* X/Open (XSI) requires <math.h> to provide M_PI, but core POSIX does not */
      10              : #ifndef M_PI
      11              : #define M_PI 3.14159265358979323846
      12              : #endif
      13              : 
      14            1 : PG_MODULE_MAGIC_EXT(
      15              :                     .name = "earthdistance",
      16              :                     .version = PG_VERSION
      17              : );
      18              : 
      19              : /* Earth's radius is in statute miles. */
      20              : static const double EARTH_RADIUS = 3958.747716;
      21              : static const double TWO_PI = 2.0 * M_PI;
      22              : 
      23              : 
      24              : /******************************************************
      25              :  *
      26              :  * degtorad - convert degrees to radians
      27              :  *
      28              :  * arg: double, angle in degrees
      29              :  *
      30              :  * returns: double, same angle in radians
      31              :  ******************************************************/
      32              : 
      33              : static double
      34           96 : degtorad(double degrees)
      35              : {
      36           96 :     return (degrees / 360.0) * TWO_PI;
      37              : }
      38              : 
      39              : /******************************************************
      40              :  *
      41              :  * geo_distance_internal - distance between points
      42              :  *
      43              :  * args:
      44              :  *   a pair of points - for each point,
      45              :  *     x-coordinate is longitude in degrees west of Greenwich
      46              :  *     y-coordinate is latitude in degrees above equator
      47              :  *
      48              :  * returns: double
      49              :  *   distance between the points in miles on earth's surface
      50              :  ******************************************************/
      51              : 
      52              : static double
      53           24 : geo_distance_internal(Point *pt1, Point *pt2)
      54              : {
      55              :     double      long1,
      56              :                 lat1,
      57              :                 long2,
      58              :                 lat2;
      59              :     double      longdiff;
      60              :     double      sino;
      61              : 
      62              :     /* convert degrees to radians */
      63              : 
      64           24 :     long1 = degtorad(pt1->x);
      65           24 :     lat1 = degtorad(pt1->y);
      66              : 
      67           24 :     long2 = degtorad(pt2->x);
      68           24 :     lat2 = degtorad(pt2->y);
      69              : 
      70              :     /* compute difference in longitudes - want < 180 degrees */
      71           24 :     longdiff = fabs(long1 - long2);
      72           24 :     if (longdiff > M_PI)
      73            0 :         longdiff = TWO_PI - longdiff;
      74              : 
      75           24 :     sino = sqrt(sin(fabs(lat1 - lat2) / 2.) * sin(fabs(lat1 - lat2) / 2.) +
      76           24 :                 cos(lat1) * cos(lat2) * sin(longdiff / 2.) * sin(longdiff / 2.));
      77           24 :     if (sino > 1.)
      78            0 :         sino = 1.;
      79              : 
      80           24 :     return 2. * EARTH_RADIUS * asin(sino);
      81              : }
      82              : 
      83              : 
      84              : /******************************************************
      85              :  *
      86              :  * geo_distance - distance between points
      87              :  *
      88              :  * args:
      89              :  *   a pair of points - for each point,
      90              :  *     x-coordinate is longitude in degrees west of Greenwich
      91              :  *     y-coordinate is latitude in degrees above equator
      92              :  *
      93              :  * returns: float8
      94              :  *   distance between the points in miles on earth's surface
      95              :  ******************************************************/
      96              : 
      97            2 : PG_FUNCTION_INFO_V1(geo_distance);
      98              : 
      99              : Datum
     100           24 : geo_distance(PG_FUNCTION_ARGS)
     101              : {
     102           24 :     Point      *pt1 = PG_GETARG_POINT_P(0);
     103           24 :     Point      *pt2 = PG_GETARG_POINT_P(1);
     104              :     float8      result;
     105              : 
     106           24 :     result = geo_distance_internal(pt1, pt2);
     107           24 :     PG_RETURN_FLOAT8(result);
     108              : }