PostGIS  2.1.10dev-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 155 of file lwspheroid.c.

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

Referenced by lwgeom_azumith_spheroid(), and ptarray_area_spheroid().

156 {
157  int i = 0;
158  double lambda = s->lon - r->lon;
159  double omf = 1 - spheroid->f;
160  double u1 = atan(omf * tan(r->lat));
161  double cos_u1 = cos(u1);
162  double sin_u1 = sin(u1);
163  double u2 = atan(omf * tan(s->lat));
164  double cos_u2 = cos(u2);
165  double sin_u2 = sin(u2);
166
167  double omega = lambda;
168  double alpha, sigma, sin_sigma, cos_sigma, cos2_sigma_m, sqr_sin_sigma, last_lambda;
169  double sin_alpha, cos_alphasq, C, alphaFD;
170  do
171  {
172  sqr_sin_sigma = POW2(cos_u2 * sin(lambda)) +
173  POW2((cos_u1 * sin_u2 - sin_u1 * cos_u2 * cos(lambda)));
174  sin_sigma = sqrt(sqr_sin_sigma);
175  cos_sigma = sin_u1 * sin_u2 + cos_u1 * cos_u2 * cos(lambda);
176  sigma = atan2(sin_sigma, cos_sigma);
177  sin_alpha = cos_u1 * cos_u2 * sin(lambda) / sin(sigma);
178
179  /* Numerical stability issue, ensure asin is not NaN */
180  if ( sin_alpha > 1.0 )
181  alpha = M_PI_2;
182  else if ( sin_alpha < -1.0 )
183  alpha = -1.0 * M_PI_2;
184  else
185  alpha = asin(sin_alpha);
186
187  cos_alphasq = POW2(cos(alpha));
188  cos2_sigma_m = cos(sigma) - (2.0 * sin_u1 * sin_u2 / cos_alphasq);
189
190  /* Numerical stability issue, cos2 is in range */
191  if ( cos2_sigma_m > 1.0 )
192  cos2_sigma_m = 1.0;
193  if ( cos2_sigma_m < -1.0 )
194  cos2_sigma_m = -1.0;
195
196  C = (spheroid->f / 16.0) * cos_alphasq * (4.0 + spheroid->f * (4.0 - 3.0 * cos_alphasq));
197  last_lambda = lambda;
198  lambda = omega + (1.0 - C) * spheroid->f * sin(alpha) * (sigma + C * sin(sigma) *
199  (cos2_sigma_m + C * cos(sigma) * (-1.0 + 2.0 * POW2(cos2_sigma_m))));
200  i++;
201  }
202  while ( (i < 999) && (lambda != 0) && (fabs((last_lambda - lambda) / lambda) > 1.0e-9) );
203
204  alphaFD = atan2((cos_u2 * sin(lambda)),
205  (cos_u1 * sin_u2 - sin_u1 * cos_u2 * cos(lambda)));
206  if (alphaFD < 0.0)
207  {
208  alphaFD = alphaFD + 2.0 * M_PI;
209  }
210  if (alphaFD > 2.0 * M_PI)
211  {
212  alphaFD = alphaFD - 2.0 * M_PI;
213  }
214  return alphaFD;
215 }
double f
Definition: liblwgeom.h:271
#define POW2(x)
Definition: lwgeodetic.h:28

Here is the caller graph for this function: