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