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           2 : 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         192 : degtorad(double degrees)
      35             : {
      36         192 :     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          48 : 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          48 :     long1 = degtorad(pt1->x);
      65          48 :     lat1 = degtorad(pt1->y);
      66             : 
      67          48 :     long2 = degtorad(pt2->x);
      68          48 :     lat2 = degtorad(pt2->y);
      69             : 
      70             :     /* compute difference in longitudes - want < 180 degrees */
      71          48 :     longdiff = fabs(long1 - long2);
      72          48 :     if (longdiff > M_PI)
      73           0 :         longdiff = TWO_PI - longdiff;
      74             : 
      75          48 :     sino = sqrt(sin(fabs(lat1 - lat2) / 2.) * sin(fabs(lat1 - lat2) / 2.) +
      76          48 :                 cos(lat1) * cos(lat2) * sin(longdiff / 2.) * sin(longdiff / 2.));
      77          48 :     if (sino > 1.)
      78           0 :         sino = 1.;
      79             : 
      80          48 :     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           4 : PG_FUNCTION_INFO_V1(geo_distance);
      98             : 
      99             : Datum
     100          48 : geo_distance(PG_FUNCTION_ARGS)
     101             : {
     102          48 :     Point      *pt1 = PG_GETARG_POINT_P(0);
     103          48 :     Point      *pt2 = PG_GETARG_POINT_P(1);
     104             :     float8      result;
     105             : 
     106          48 :     result = geo_distance_internal(pt1, pt2);
     107          48 :     PG_RETURN_FLOAT8(result);
     108             : }