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