PostGIS  2.1.10dev-r@@SVN_REVISION@@
 double spheroid_distance ( const GEOGRAPHIC_POINT * a, const GEOGRAPHIC_POINT * b, const SPHEROID * spheroid )

Computes the shortest distance along the surface of the spheroid between two points.

Based on Vincenty's formula for the geodetic inverse problem as described in "Geocentric Datum of Australia Technical Manual", Chapter 4. Tested against: http://mascot.gdbc.gov.bc.ca/mascot/util1a.html and http://www.ga.gov.au/nmd/geodesy/datums/vincenty_inverse.jsp

Parameters
 a - location of first point. b - location of second point. s - spheroid to calculate on
Returns
spheroidal distance between a and b in spheroid units.

Definition at line 59 of file lwspheroid.c.

60 {
61  double lambda = (b->lon - a->lon);
62  double f = spheroid->f;
63  double omf = 1 - spheroid->f;
64  double u1, u2;
65  double cos_u1, cos_u2;
66  double sin_u1, sin_u2;
67  double big_a, big_b, delta_sigma;
68  double alpha, sin_alpha, cos_alphasq, c;
69  double sigma, sin_sigma, cos_sigma, cos2_sigma_m, sqrsin_sigma, last_lambda, omega;
70  double cos_lambda, sin_lambda;
71  double distance;
72  int i = 0;
73
74  /* Same point => zero distance */
75  if ( geographic_point_equals(a, b) )
76  {
77  return 0.0;
78  }
79
80  u1 = atan(omf * tan(a->lat));
81  cos_u1 = cos(u1);
82  sin_u1 = sin(u1);
83  u2 = atan(omf * tan(b->lat));
84  cos_u2 = cos(u2);
85  sin_u2 = sin(u2);
86
87  omega = lambda;
88  do
89  {
90  cos_lambda = cos(lambda);
91  sin_lambda = sin(lambda);
92  sqrsin_sigma = POW2(cos_u2 * sin_lambda) +
93  POW2((cos_u1 * sin_u2 - sin_u1 * cos_u2 * cos_lambda));
94  sin_sigma = sqrt(sqrsin_sigma);
95  cos_sigma = sin_u1 * sin_u2 + cos_u1 * cos_u2 * cos_lambda;
96  sigma = atan2(sin_sigma, cos_sigma);
97  sin_alpha = cos_u1 * cos_u2 * sin_lambda / sin(sigma);
98
99  /* Numerical stability issue, ensure asin is not NaN */
100  if ( sin_alpha > 1.0 )
101  alpha = M_PI_2;
102  else if ( sin_alpha < -1.0 )
103  alpha = -1.0 * M_PI_2;
104  else
105  alpha = asin(sin_alpha);
106
107  cos_alphasq = POW2(cos(alpha));
108  cos2_sigma_m = cos(sigma) - (2.0 * sin_u1 * sin_u2 / cos_alphasq);
109
110  /* Numerical stability issue, cos2 is in range */
111  if ( cos2_sigma_m > 1.0 )
112  cos2_sigma_m = 1.0;
113  if ( cos2_sigma_m < -1.0 )
114  cos2_sigma_m = -1.0;
115
116  c = (f / 16.0) * cos_alphasq * (4.0 + f * (4.0 - 3.0 * cos_alphasq));
117  last_lambda = lambda;
118  lambda = omega + (1.0 - c) * f * sin(alpha) * (sigma + c * sin(sigma) *
119  (cos2_sigma_m + c * cos(sigma) * (-1.0 + 2.0 * POW2(cos2_sigma_m))));
120  i++;
121  }
122  while ( (i < 999) && (lambda != 0.0) && (fabs((last_lambda - lambda)/lambda) > 1.0e-9) );
123
124  u2 = spheroid_mu2(alpha, spheroid);
125  big_a = spheroid_big_a(u2);
126  big_b = spheroid_big_b(u2);
127  delta_sigma = big_b * sin_sigma * (cos2_sigma_m + (big_b / 4.0) * (cos_sigma * (-1.0 + 2.0 * POW2(cos2_sigma_m)) -
128  (big_b / 6.0) * cos2_sigma_m * (-3.0 + 4.0 * sqrsin_sigma) * (-3.0 + 4.0 * POW2(cos2_sigma_m))));
129
130  distance = spheroid->b * big_a * (sigma - delta_sigma);
131
132  /* Algorithm failure, distance == NaN, fallback to sphere */
133  if ( distance != distance )
134  {
135  lwerror("spheroid_distance returned NaN: (%.20g %.20g) (%.20g %.20g) a = %.20g b = %.20g",a->lat, a->lon, b->lat, b->lon, spheroid->a, spheroid->b);
136  return spheroid->radius * sphere_distance(a, b);
137  }
138
139  return distance;
140 }
static double spheroid_big_b(double u2)
Definition: lwspheroid.c:40
double sphere_distance(const GEOGRAPHIC_POINT *s, const GEOGRAPHIC_POINT *e)
Given two points on a unit sphere, calculate their distance apart in radians.
Definition: lwgeodetic.c:897
static double spheroid_big_a(double u2)
Definition: lwspheroid.c:35
double b
Definition: liblwgeom.h:270
Definition: liblwgeom.h:274
static double spheroid_mu2(double alpha, const SPHEROID *s)
Definition: lwspheroid.c:29
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:67
double f
Definition: liblwgeom.h:271
#define POW2(x)
Definition: lwgeodetic.h:28
Datum distance(PG_FUNCTION_ARGS)
int geographic_point_equals(const GEOGRAPHIC_POINT *g1, const GEOGRAPHIC_POINT *g2)
Definition: lwgeodetic.c:147
double a
Definition: liblwgeom.h:269

Here is the call graph for this function:

Here is the caller graph for this function: