PostGIS  2.5.1dev-r@@SVN_REVISION@@

## ◆ spheroid_distance()

 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 186 of file lwspheroid.c.

187 {
188  double lambda = (b->lon - a->lon);
189  double f = spheroid->f;
190  double omf = 1 - spheroid->f;
191  double u1, u2;
192  double cos_u1, cos_u2;
193  double sin_u1, sin_u2;
194  double big_a, big_b, delta_sigma;
195  double alpha, sin_alpha, cos_alphasq, c;
196  double sigma, sin_sigma, cos_sigma, cos2_sigma_m, sqrsin_sigma, last_lambda, omega;
197  double cos_lambda, sin_lambda;
198  double distance;
199  int i = 0;
200
201  /* Same point => zero distance */
202  if ( geographic_point_equals(a, b) )
203  {
204  return 0.0;
205  }
206
207  u1 = atan(omf * tan(a->lat));
208  cos_u1 = cos(u1);
209  sin_u1 = sin(u1);
210  u2 = atan(omf * tan(b->lat));
211  cos_u2 = cos(u2);
212  sin_u2 = sin(u2);
213
214  omega = lambda;
215  do
216  {
217  cos_lambda = cos(lambda);
218  sin_lambda = sin(lambda);
219  sqrsin_sigma = POW2(cos_u2 * sin_lambda) +
220  POW2((cos_u1 * sin_u2 - sin_u1 * cos_u2 * cos_lambda));
221  sin_sigma = sqrt(sqrsin_sigma);
222  cos_sigma = sin_u1 * sin_u2 + cos_u1 * cos_u2 * cos_lambda;
223  sigma = atan2(sin_sigma, cos_sigma);
224  sin_alpha = cos_u1 * cos_u2 * sin_lambda / sin(sigma);
225
226  /* Numerical stability issue, ensure asin is not NaN */
227  if ( sin_alpha > 1.0 )
228  alpha = M_PI_2;
229  else if ( sin_alpha < -1.0 )
230  alpha = -1.0 * M_PI_2;
231  else
232  alpha = asin(sin_alpha);
233
234  cos_alphasq = POW2(cos(alpha));
235  cos2_sigma_m = cos(sigma) - (2.0 * sin_u1 * sin_u2 / cos_alphasq);
236
237  /* Numerical stability issue, cos2 is in range */
238  if ( cos2_sigma_m > 1.0 )
239  cos2_sigma_m = 1.0;
240  if ( cos2_sigma_m < -1.0 )
241  cos2_sigma_m = -1.0;
242
243  c = (f / 16.0) * cos_alphasq * (4.0 + f * (4.0 - 3.0 * cos_alphasq));
244  last_lambda = lambda;
245  lambda = omega + (1.0 - c) * f * sin(alpha) * (sigma + c * sin(sigma) *
246  (cos2_sigma_m + c * cos(sigma) * (-1.0 + 2.0 * POW2(cos2_sigma_m))));
247  i++;
248  }
249  while ( (i < 999) && (lambda != 0.0) && (fabs((last_lambda - lambda)/lambda) > 1.0e-9) );
250
251  u2 = spheroid_mu2(alpha, spheroid);
252  big_a = spheroid_big_a(u2);
253  big_b = spheroid_big_b(u2);
254  delta_sigma = big_b * sin_sigma * (cos2_sigma_m + (big_b / 4.0) * (cos_sigma * (-1.0 + 2.0 * POW2(cos2_sigma_m)) -
255  (big_b / 6.0) * cos2_sigma_m * (-3.0 + 4.0 * sqrsin_sigma) * (-3.0 + 4.0 * POW2(cos2_sigma_m))));
256
257  distance = spheroid->b * big_a * (sigma - delta_sigma);
258
259  /* Algorithm failure, distance == NaN, fallback to sphere */
260  if ( distance != distance )
261  {
262  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);
263  return spheroid->radius * sphere_distance(a, b);
264  }
265
266  return distance;
267 }
static double spheroid_big_b(double u2)
Definition: lwspheroid.c:60
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:917
static double spheroid_big_a(double u2)
Definition: lwspheroid.c:55
double b
Definition: liblwgeom.h:316
Definition: liblwgeom.h:320
static double spheroid_mu2(double alpha, const SPHEROID *s)
Definition: lwspheroid.c:49
double f
Definition: liblwgeom.h:317
#define POW2(x)
Definition: lwgeodetic.h:42
Datum distance(PG_FUNCTION_ARGS)
int geographic_point_equals(const GEOGRAPHIC_POINT *g1, const GEOGRAPHIC_POINT *g2)
Definition: lwgeodetic.c:161
double a
Definition: liblwgeom.h:315
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:190
Here is the call graph for this function:
Here is the caller graph for this function: