PostGIS  2.5.0dev-r@@SVN_REVISION@@
double spheroid_direction ( const GEOGRAPHIC_POINT r,
const GEOGRAPHIC_POINT s,
const SPHEROID spheroid 
)

Computes the direction of the geodesic joining two points on the spheroid.

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
r- location of first point
s- location of second point
Returns
azimuth of line joining r and s

Definition at line 282 of file lwspheroid.c.

References SPHEROID::f, GEOGRAPHIC_POINT::lat, GEOGRAPHIC_POINT::lon, and POW2.

Referenced by lwgeom_azumith_spheroid(), and ptarray_area_spheroid().

283 {
284  int i = 0;
285  double lambda = s->lon - r->lon;
286  double omf = 1 - spheroid->f;
287  double u1 = atan(omf * tan(r->lat));
288  double cos_u1 = cos(u1);
289  double sin_u1 = sin(u1);
290  double u2 = atan(omf * tan(s->lat));
291  double cos_u2 = cos(u2);
292  double sin_u2 = sin(u2);
293 
294  double omega = lambda;
295  double alpha, sigma, sin_sigma, cos_sigma, cos2_sigma_m, sqr_sin_sigma, last_lambda;
296  double sin_alpha, cos_alphasq, C, alphaFD;
297  do
298  {
299  sqr_sin_sigma = POW2(cos_u2 * sin(lambda)) +
300  POW2((cos_u1 * sin_u2 - sin_u1 * cos_u2 * cos(lambda)));
301  sin_sigma = sqrt(sqr_sin_sigma);
302  cos_sigma = sin_u1 * sin_u2 + cos_u1 * cos_u2 * cos(lambda);
303  sigma = atan2(sin_sigma, cos_sigma);
304  sin_alpha = cos_u1 * cos_u2 * sin(lambda) / sin(sigma);
305 
306  /* Numerical stability issue, ensure asin is not NaN */
307  if ( sin_alpha > 1.0 )
308  alpha = M_PI_2;
309  else if ( sin_alpha < -1.0 )
310  alpha = -1.0 * M_PI_2;
311  else
312  alpha = asin(sin_alpha);
313 
314  cos_alphasq = POW2(cos(alpha));
315  cos2_sigma_m = cos(sigma) - (2.0 * sin_u1 * sin_u2 / cos_alphasq);
316 
317  /* Numerical stability issue, cos2 is in range */
318  if ( cos2_sigma_m > 1.0 )
319  cos2_sigma_m = 1.0;
320  if ( cos2_sigma_m < -1.0 )
321  cos2_sigma_m = -1.0;
322 
323  C = (spheroid->f / 16.0) * cos_alphasq * (4.0 + spheroid->f * (4.0 - 3.0 * cos_alphasq));
324  last_lambda = lambda;
325  lambda = omega + (1.0 - C) * spheroid->f * sin(alpha) * (sigma + C * sin(sigma) *
326  (cos2_sigma_m + C * cos(sigma) * (-1.0 + 2.0 * POW2(cos2_sigma_m))));
327  i++;
328  }
329  while ( (i < 999) && (lambda != 0) && (fabs((last_lambda - lambda) / lambda) > 1.0e-9) );
330 
331  alphaFD = atan2((cos_u2 * sin(lambda)),
332  (cos_u1 * sin_u2 - sin_u1 * cos_u2 * cos(lambda)));
333  if (alphaFD < 0.0)
334  {
335  alphaFD = alphaFD + 2.0 * M_PI;
336  }
337  if (alphaFD > 2.0 * M_PI)
338  {
339  alphaFD = alphaFD - 2.0 * M_PI;
340  }
341  return alphaFD;
342 }
double f
Definition: liblwgeom.h:314
#define POW2(x)
Definition: lwgeodetic.h:42

Here is the caller graph for this function: