2010-09-20 22:08:53 +02:00
|
|
|
/* contrib/earthdistance/earthdistance.c */
|
2006-03-11 05:38:42 +01:00
|
|
|
|
2001-02-10 03:31:31 +01:00
|
|
|
#include "postgres.h"
|
|
|
|
|
1998-06-16 05:55:15 +02:00
|
|
|
#include <math.h>
|
|
|
|
|
2008-04-20 03:05:52 +02:00
|
|
|
#include "utils/geo_decls.h" /* for Point */
|
2001-02-10 03:31:31 +01:00
|
|
|
|
2006-10-19 22:08:03 +02:00
|
|
|
#ifndef M_PI
|
|
|
|
#define M_PI 3.14159265358979323846
|
|
|
|
#endif
|
|
|
|
|
2006-05-31 00:12:16 +02:00
|
|
|
PG_MODULE_MAGIC;
|
|
|
|
|
1998-06-16 05:55:15 +02:00
|
|
|
/* Earth's radius is in statute miles. */
|
2008-04-20 03:05:52 +02:00
|
|
|
static const double EARTH_RADIUS = 3958.747716;
|
|
|
|
static const double TWO_PI = 2.0 * M_PI;
|
2000-06-15 20:55:34 +02:00
|
|
|
|
1998-06-16 05:55:15 +02:00
|
|
|
|
|
|
|
/******************************************************
|
|
|
|
*
|
|
|
|
* degtorad - convert degrees to radians
|
|
|
|
*
|
1999-05-25 18:15:34 +02:00
|
|
|
* arg: double, angle in degrees
|
1998-06-16 05:55:15 +02:00
|
|
|
*
|
1999-05-25 18:15:34 +02:00
|
|
|
* returns: double, same angle in radians
|
1998-06-16 05:55:15 +02:00
|
|
|
******************************************************/
|
|
|
|
|
|
|
|
static double
|
1999-05-25 18:15:34 +02:00
|
|
|
degtorad(double degrees)
|
|
|
|
{
|
1998-06-16 05:55:15 +02:00
|
|
|
return (degrees / 360.0) * TWO_PI;
|
|
|
|
}
|
|
|
|
|
|
|
|
/******************************************************
|
|
|
|
*
|
2008-04-21 03:11:43 +02:00
|
|
|
* geo_distance_internal - distance between points
|
1998-06-16 05:55:15 +02:00
|
|
|
*
|
|
|
|
* args:
|
1999-05-25 18:15:34 +02:00
|
|
|
* a pair of points - for each point,
|
|
|
|
* x-coordinate is longitude in degrees west of Greenwich
|
|
|
|
* y-coordinate is latitude in degrees above equator
|
1998-06-16 05:55:15 +02:00
|
|
|
*
|
2008-04-21 03:11:43 +02:00
|
|
|
* returns: double
|
1999-05-25 18:15:34 +02:00
|
|
|
* distance between the points in miles on earth's surface
|
1998-06-16 05:55:15 +02:00
|
|
|
******************************************************/
|
|
|
|
|
2008-04-21 03:11:43 +02:00
|
|
|
static double
|
|
|
|
geo_distance_internal(Point *pt1, Point *pt2)
|
1999-05-25 18:15:34 +02:00
|
|
|
{
|
|
|
|
double long1,
|
|
|
|
lat1,
|
|
|
|
long2,
|
|
|
|
lat2;
|
|
|
|
double longdiff;
|
2002-11-08 21:20:22 +01:00
|
|
|
double sino;
|
1998-06-16 05:55:15 +02:00
|
|
|
|
|
|
|
/* convert degrees to radians */
|
|
|
|
|
1999-05-25 18:15:34 +02:00
|
|
|
long1 = degtorad(pt1->x);
|
|
|
|
lat1 = degtorad(pt1->y);
|
1998-06-16 05:55:15 +02:00
|
|
|
|
1999-05-25 18:15:34 +02:00
|
|
|
long2 = degtorad(pt2->x);
|
|
|
|
lat2 = degtorad(pt2->y);
|
1998-06-16 05:55:15 +02:00
|
|
|
|
|
|
|
/* compute difference in longitudes - want < 180 degrees */
|
1999-05-25 18:15:34 +02:00
|
|
|
longdiff = fabs(long1 - long2);
|
1998-06-16 05:55:15 +02:00
|
|
|
if (longdiff > M_PI)
|
|
|
|
longdiff = TWO_PI - longdiff;
|
|
|
|
|
2003-08-04 02:43:34 +02:00
|
|
|
sino = sqrt(sin(fabs(lat1 - lat2) / 2.) * sin(fabs(lat1 - lat2) / 2.) +
|
Phase 3 of pgindent updates.
Don't move parenthesized lines to the left, even if that means they
flow past the right margin.
By default, BSD indent lines up statement continuation lines that are
within parentheses so that they start just to the right of the preceding
left parenthesis. However, traditionally, if that resulted in the
continuation line extending to the right of the desired right margin,
then indent would push it left just far enough to not overrun the margin,
if it could do so without making the continuation line start to the left of
the current statement indent. That makes for a weird mix of indentations
unless one has been completely rigid about never violating the 80-column
limit.
This behavior has been pretty universally panned by Postgres developers.
Hence, disable it with indent's new -lpl switch, so that parenthesized
lines are always lined up with the preceding left paren.
This patch is much less interesting than the first round of indent
changes, but also bulkier, so I thought it best to separate the effects.
Discussion: https://postgr.es/m/E1dAmxK-0006EE-1r@gemulon.postgresql.org
Discussion: https://postgr.es/m/30527.1495162840@sss.pgh.pa.us
2017-06-21 21:35:54 +02:00
|
|
|
cos(lat1) * cos(lat2) * sin(longdiff / 2.) * sin(longdiff / 2.));
|
2003-08-04 02:43:34 +02:00
|
|
|
if (sino > 1.)
|
|
|
|
sino = 1.;
|
1998-06-16 05:55:15 +02:00
|
|
|
|
2008-04-21 03:11:43 +02:00
|
|
|
return 2. * EARTH_RADIUS * asin(sino);
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
/******************************************************
|
|
|
|
*
|
|
|
|
* geo_distance - distance between points
|
|
|
|
*
|
|
|
|
* args:
|
|
|
|
* a pair of points - for each point,
|
|
|
|
* x-coordinate is longitude in degrees west of Greenwich
|
|
|
|
* y-coordinate is latitude in degrees above equator
|
|
|
|
*
|
|
|
|
* returns: float8
|
|
|
|
* distance between the points in miles on earth's surface
|
|
|
|
******************************************************/
|
|
|
|
|
|
|
|
PG_FUNCTION_INFO_V1(geo_distance);
|
|
|
|
|
|
|
|
Datum
|
|
|
|
geo_distance(PG_FUNCTION_ARGS)
|
|
|
|
{
|
|
|
|
Point *pt1 = PG_GETARG_POINT_P(0);
|
|
|
|
Point *pt2 = PG_GETARG_POINT_P(1);
|
|
|
|
float8 result;
|
|
|
|
|
|
|
|
result = geo_distance_internal(pt1, pt2);
|
2008-04-20 03:05:52 +02:00
|
|
|
PG_RETURN_FLOAT8(result);
|
1998-06-16 05:55:15 +02:00
|
|
|
}
|