TLA  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 GIC           1 : PG_MODULE_MAGIC;
      15 ECB             : 
      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 GIC          96 : degtorad(double degrees)
      32 ECB             : {
      33 GIC          96 :     return (degrees / 360.0) * TWO_PI;
      34 ECB             : }
      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 GIC          24 : geo_distance_internal(Point *pt1, Point *pt2)
      51 ECB             : {
      52                 :     double      long1,
      53                 :                 lat1,
      54                 :                 long2,
      55                 :                 lat2;
      56                 :     double      longdiff;
      57                 :     double      sino;
      58                 : 
      59                 :     /* convert degrees to radians */
      60                 : 
      61 GIC          24 :     long1 = degtorad(pt1->x);
      62 CBC          24 :     lat1 = degtorad(pt1->y);
      63 ECB             : 
      64 GIC          24 :     long2 = degtorad(pt2->x);
      65 CBC          24 :     lat2 = degtorad(pt2->y);
      66 ECB             : 
      67                 :     /* compute difference in longitudes - want < 180 degrees */
      68 GIC          24 :     longdiff = fabs(long1 - long2);
      69 CBC          24 :     if (longdiff > M_PI)
      70 LBC           0 :         longdiff = TWO_PI - longdiff;
      71 EUB             : 
      72 GIC          24 :     sino = sqrt(sin(fabs(lat1 - lat2) / 2.) * sin(fabs(lat1 - lat2) / 2.) +
      73 CBC          24 :                 cos(lat1) * cos(lat2) * sin(longdiff / 2.) * sin(longdiff / 2.));
      74              24 :     if (sino > 1.)
      75 LBC           0 :         sino = 1.;
      76 EUB             : 
      77 GIC          24 :     return 2. * EARTH_RADIUS * asin(sino);
      78 ECB             : }
      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 GIC           2 : PG_FUNCTION_INFO_V1(geo_distance);
      95 ECB             : 
      96                 : Datum
      97 GIC          24 : geo_distance(PG_FUNCTION_ARGS)
      98 ECB             : {
      99 GIC          24 :     Point      *pt1 = PG_GETARG_POINT_P(0);
     100 CBC          24 :     Point      *pt2 = PG_GETARG_POINT_P(1);
     101 ECB             :     float8      result;
     102                 : 
     103 GIC          24 :     result = geo_distance_internal(pt1, pt2);
     104 CBC          24 :     PG_RETURN_FLOAT8(result);
     105 ECB             : }