Age Owner Branch data 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 : :
164 tgl@sss.pgh.pa.us 14 :CBC 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
9601 bruce@momjian.us 34 : 96 : degtorad(double degrees)
35 : : {
9944 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
6347 tgl@sss.pgh.pa.us 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 : :
9601 bruce@momjian.us 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);
9944 72 [ - + ]: 24 : if (longdiff > M_PI)
9944 bruce@momjian.us 73 :UBC 0 : longdiff = TWO_PI - longdiff;
74 : :
8069 bruce@momjian.us 75 :CBC 24 : sino = sqrt(sin(fabs(lat1 - lat2) / 2.) * sin(fabs(lat1 - lat2) / 2.) +
2999 tgl@sss.pgh.pa.us 76 : 24 : cos(lat1) * cos(lat2) * sin(longdiff / 2.) * sin(longdiff / 2.));
8069 bruce@momjian.us 77 [ - + ]: 24 : if (sino > 1.)
8069 bruce@momjian.us 78 :UBC 0 : sino = 1.;
79 : :
6347 tgl@sss.pgh.pa.us 80 :CBC 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);
6348 107 : 24 : PG_RETURN_FLOAT8(result);
108 : : }
|
