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