Branch data     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                 :             : }