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

Referenced by lwgeom_project_spheroid(), and ptarray_area_spheroid().

233 {
234  double omf = 1 - spheroid->f;
235  double tan_u1 = omf * tan(r->lat);
236  double u1 = atan(tan_u1);
237  double sigma, last_sigma, delta_sigma, two_sigma_m;
238  double sigma1, sin_alpha, alpha, cos_alphasq;
239  double u2, A, B;
240  double lat2, lambda, lambda2, C, omega;
241  int i = 0;
242
243  if (azimuth < 0.0)
244  {
245  azimuth = azimuth + M_PI * 2.0;
246  }
247  if (azimuth > (PI * 2.0))
248  {
249  azimuth = azimuth - M_PI * 2.0;
250  }
251
252  sigma1 = atan2(tan_u1, cos(azimuth));
253  sin_alpha = cos(u1) * sin(azimuth);
254  alpha = asin(sin_alpha);
255  cos_alphasq = 1.0 - POW2(sin_alpha);
256
257  u2 = spheroid_mu2(alpha, spheroid);
258  A = spheroid_big_a(u2);
259  B = spheroid_big_b(u2);
260
261  sigma = (distance / (spheroid->b * A));
262  do
263  {
264  two_sigma_m = 2.0 * sigma1 + sigma;
265  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))))));
266  last_sigma = sigma;
267  sigma = (distance / (spheroid->b * A)) + delta_sigma;
268  i++;
269  }
270  while (i < 999 && fabs((last_sigma - sigma) / sigma) > 1.0e-9);
271
272  lat2 = atan2((sin(u1) * cos(sigma) + cos(u1) * sin(sigma) *
273  cos(azimuth)), (omf * sqrt(POW2(sin_alpha) +
274  POW2(sin(u1) * sin(sigma) - cos(u1) * cos(sigma) *
275  cos(azimuth)))));
276  lambda = atan2((sin(sigma) * sin(azimuth)), (cos(u1) * cos(sigma) -
277  sin(u1) * sin(sigma) * cos(azimuth)));
278  C = (spheroid->f / 16.0) * cos_alphasq * (4.0 + spheroid->f * (4.0 - 3.0 * cos_alphasq));
279  omega = lambda - (1.0 - C) * spheroid->f * sin_alpha * (sigma + C * sin(sigma) *
280  (cos(two_sigma_m) + C * cos(sigma) * (-1.0 + 2.0 * POW2(cos(two_sigma_m)))));
281  lambda2 = r->lon + omega;
282  g->lat = lat2;
283  g->lon = lambda2;
284  return LW_SUCCESS;
285 }
static double spheroid_big_b(double u2)
Definition: lwspheroid.c:40
static double spheroid_big_a(double u2)
Definition: lwspheroid.c:35
double b
Definition: liblwgeom.h:270
#define LW_SUCCESS
Definition: liblwgeom.h:55
static double spheroid_mu2(double alpha, const SPHEROID *s)
Definition: lwspheroid.c:29
double f
Definition: liblwgeom.h:271
#define POW2(x)
Definition: lwgeodetic.h:28
Datum distance(PG_FUNCTION_ARGS)
#define PI
PI.

Here is the call graph for this function:

Here is the caller graph for this function: