PostGIS  2.3.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 359 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().

360 {
361  double omf = 1 - spheroid->f;
362  double tan_u1 = omf * tan(r->lat);
363  double u1 = atan(tan_u1);
364  double sigma, last_sigma, delta_sigma, two_sigma_m;
365  double sigma1, sin_alpha, alpha, cos_alphasq;
366  double u2, A, B;
367  double lat2, lambda, lambda2, C, omega;
368  int i = 0;
369 
370  if (azimuth < 0.0)
371  {
372  azimuth = azimuth + M_PI * 2.0;
373  }
374  if (azimuth > (M_PI * 2.0))
375  {
376  azimuth = azimuth - M_PI * 2.0;
377  }
378 
379  sigma1 = atan2(tan_u1, cos(azimuth));
380  sin_alpha = cos(u1) * sin(azimuth);
381  alpha = asin(sin_alpha);
382  cos_alphasq = 1.0 - POW2(sin_alpha);
383 
384  u2 = spheroid_mu2(alpha, spheroid);
385  A = spheroid_big_a(u2);
386  B = spheroid_big_b(u2);
387 
388  sigma = (distance / (spheroid->b * A));
389  do
390  {
391  two_sigma_m = 2.0 * sigma1 + sigma;
392  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))))));
393  last_sigma = sigma;
394  sigma = (distance / (spheroid->b * A)) + delta_sigma;
395  i++;
396  }
397  while (i < 999 && fabs((last_sigma - sigma) / sigma) > 1.0e-9);
398 
399  lat2 = atan2((sin(u1) * cos(sigma) + cos(u1) * sin(sigma) *
400  cos(azimuth)), (omf * sqrt(POW2(sin_alpha) +
401  POW2(sin(u1) * sin(sigma) - cos(u1) * cos(sigma) *
402  cos(azimuth)))));
403  lambda = atan2((sin(sigma) * sin(azimuth)), (cos(u1) * cos(sigma) -
404  sin(u1) * sin(sigma) * cos(azimuth)));
405  C = (spheroid->f / 16.0) * cos_alphasq * (4.0 + spheroid->f * (4.0 - 3.0 * cos_alphasq));
406  omega = lambda - (1.0 - C) * spheroid->f * sin_alpha * (sigma + C * sin(sigma) *
407  (cos(two_sigma_m) + C * cos(sigma) * (-1.0 + 2.0 * POW2(cos(two_sigma_m)))));
408  lambda2 = r->lon + omega;
409  g->lat = lat2;
410  g->lon = lambda2;
411  return LW_SUCCESS;
412 }
static double spheroid_big_b(double u2)
Definition: lwspheroid.c:60
static double spheroid_big_a(double u2)
Definition: lwspheroid.c:55
double b
Definition: liblwgeom.h:313
#define LW_SUCCESS
Definition: liblwgeom.h:79
static double spheroid_mu2(double alpha, const SPHEROID *s)
Definition: lwspheroid.c:49
double f
Definition: liblwgeom.h:314
#define POW2(x)
Definition: lwgeodetic.h:42
Datum distance(PG_FUNCTION_ARGS)

Here is the call graph for this function:

Here is the caller graph for this function: