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 : :
6529 tgl@sss.pgh.pa.us 14 :CBC 1 : 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
9091 bruce@momjian.us 31 : 96 : degtorad(double degrees)
32 : : {
9434 33 : 96 : 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
5837 tgl@sss.pgh.pa.us 50 : 24 : 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 : :
9091 bruce@momjian.us 61 : 24 : long1 = degtorad(pt1->x);
62 : 24 : lat1 = degtorad(pt1->y);
63 : :
64 : 24 : long2 = degtorad(pt2->x);
65 : 24 : lat2 = degtorad(pt2->y);
66 : :
67 : : /* compute difference in longitudes - want < 180 degrees */
68 : 24 : longdiff = fabs(long1 - long2);
9434 69 [ - + ]: 24 : if (longdiff > M_PI)
9434 bruce@momjian.us 70 :UBC 0 : longdiff = TWO_PI - longdiff;
71 : :
7559 bruce@momjian.us 72 :CBC 24 : sino = sqrt(sin(fabs(lat1 - lat2) / 2.) * sin(fabs(lat1 - lat2) / 2.) +
2489 tgl@sss.pgh.pa.us 73 : 24 : cos(lat1) * cos(lat2) * sin(longdiff / 2.) * sin(longdiff / 2.));
7559 bruce@momjian.us 74 [ - + ]: 24 : if (sino > 1.)
7559 bruce@momjian.us 75 :UBC 0 : sino = 1.;
76 : :
5837 tgl@sss.pgh.pa.us 77 :CBC 24 : 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 : 2 : PG_FUNCTION_INFO_V1(geo_distance);
95 : :
96 : : Datum
97 : 24 : geo_distance(PG_FUNCTION_ARGS)
98 : : {
99 : 24 : Point *pt1 = PG_GETARG_POINT_P(0);
100 : 24 : Point *pt2 = PG_GETARG_POINT_P(1);
101 : : float8 result;
102 : :
103 : 24 : result = geo_distance_internal(pt1, pt2);
5838 104 : 24 : PG_RETURN_FLOAT8(result);
105 : : }
|
