PostGIS  2.2.7dev-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 268 of file lwspheroid.c.

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

Referenced by lwgeom_azumith_spheroid(), and ptarray_area_spheroid().

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

Here is the caller graph for this function: