PostGIS  2.2.7dev-r@@SVN_REVISION@@
int spheroid_project ( const GEOGRAPHIC_POINT r,
const SPHEROID spheroid,
double  distance,
double  azimuth,
GEOGRAPHIC_POINT g 
)

Given a location, an azimuth and a distance, computes the location of the projected point.

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

Parameters
r- location of first point.
distance- distance in meters.
azimuth- azimuth in radians.
Returns
s - location of projected point.

Definition at line 345 of file lwspheroid.c.

References SPHEROID::b, SPHEROID::f, GEOGRAPHIC_POINT::lat, GEOGRAPHIC_POINT::lon, LW_SUCCESS, POW2, spheroid_big_a(), spheroid_big_b(), and spheroid_mu2().

Referenced by lwgeom_project_spheroid(), and ptarray_area_spheroid().

346 {
347  double omf = 1 - spheroid->f;
348  double tan_u1 = omf * tan(r->lat);
349  double u1 = atan(tan_u1);
350  double sigma, last_sigma, delta_sigma, two_sigma_m;
351  double sigma1, sin_alpha, alpha, cos_alphasq;
352  double u2, A, B;
353  double lat2, lambda, lambda2, C, omega;
354  int i = 0;
355 
356  if (azimuth < 0.0)
357  {
358  azimuth = azimuth + M_PI * 2.0;
359  }
360  if (azimuth > (M_PI * 2.0))
361  {
362  azimuth = azimuth - M_PI * 2.0;
363  }
364 
365  sigma1 = atan2(tan_u1, cos(azimuth));
366  sin_alpha = cos(u1) * sin(azimuth);
367  alpha = asin(sin_alpha);
368  cos_alphasq = 1.0 - POW2(sin_alpha);
369 
370  u2 = spheroid_mu2(alpha, spheroid);
371  A = spheroid_big_a(u2);
372  B = spheroid_big_b(u2);
373 
374  sigma = (distance / (spheroid->b * A));
375  do
376  {
377  two_sigma_m = 2.0 * sigma1 + sigma;
378  delta_sigma = B * sin(sigma) * (cos(two_sigma_m) + (B / 4.0) * (cos(sigma) * (-1.0 + 2.0 * POW2(cos(two_sigma_m)) - (B / 6.0) * cos(two_sigma_m) * (-3.0 + 4.0 * POW2(sin(sigma))) * (-3.0 + 4.0 * POW2(cos(two_sigma_m))))));
379  last_sigma = sigma;
380  sigma = (distance / (spheroid->b * A)) + delta_sigma;
381  i++;
382  }
383  while (i < 999 && fabs((last_sigma - sigma) / sigma) > 1.0e-9);
384 
385  lat2 = atan2((sin(u1) * cos(sigma) + cos(u1) * sin(sigma) *
386  cos(azimuth)), (omf * sqrt(POW2(sin_alpha) +
387  POW2(sin(u1) * sin(sigma) - cos(u1) * cos(sigma) *
388  cos(azimuth)))));
389  lambda = atan2((sin(sigma) * sin(azimuth)), (cos(u1) * cos(sigma) -
390  sin(u1) * sin(sigma) * cos(azimuth)));
391  C = (spheroid->f / 16.0) * cos_alphasq * (4.0 + spheroid->f * (4.0 - 3.0 * cos_alphasq));
392  omega = lambda - (1.0 - C) * spheroid->f * sin_alpha * (sigma + C * sin(sigma) *
393  (cos(two_sigma_m) + C * cos(sigma) * (-1.0 + 2.0 * POW2(cos(two_sigma_m)))));
394  lambda2 = r->lon + omega;
395  g->lat = lat2;
396  g->lon = lambda2;
397  return LW_SUCCESS;
398 }
static double spheroid_big_b(double u2)
Definition: lwspheroid.c:46
static double spheroid_big_a(double u2)
Definition: lwspheroid.c:41
double b
Definition: liblwgeom.h:298
#define LW_SUCCESS
Definition: liblwgeom.h:65
static double spheroid_mu2(double alpha, const SPHEROID *s)
Definition: lwspheroid.c:35
double f
Definition: liblwgeom.h:299
#define POW2(x)
Definition: lwgeodetic.h:27
Datum distance(PG_FUNCTION_ARGS)

Here is the call graph for this function:

Here is the caller graph for this function: