PostGIS  2.5.0dev-r@@SVN_REVISION@@

◆ distance_ellipse_calculation()

double distance_ellipse_calculation ( double  lat1,
double  long1,
double  lat2,
double  long2,
SPHEROID sphere 
)

Definition at line 246 of file lwgeom_spheroid.c.

References SPHEROID::b, bigA(), bigB(), deltaLongitude(), SPHEROID::f, LWGEOM_length2d_ellipsoid(), mu2(), and PG_FUNCTION_INFO_V1().

Referenced by distance_ellipse().

248 {
249  /*
250  * Code is taken from David Skea
251  * Geographic Data BC, Province of British Columbia, Canada.
252  * Thanks to GDBC and David Skea for allowing this to be
253  * put in PostGIS.
254  */
255 
256  double L1,L2,sinU1,sinU2,cosU1,cosU2;
257  double dl,dl1,dl2,dl3,cosdl1,sindl1;
258  double cosSigma,sigma,azimuthEQ,tsm;
259  double u2,A,B;
260  double dsigma;
261 
262  double TEMP;
263 
264  int iterations;
265 
266 
267  L1 = atan((1.0 - sphere->f ) * tan( lat1) );
268  L2 = atan((1.0 - sphere->f ) * tan( lat2) );
269  sinU1 = sin(L1);
270  sinU2 = sin(L2);
271  cosU1 = cos(L1);
272  cosU2 = cos(L2);
273 
274  dl = long2- long1;
275  dl1 = dl;
276  cosdl1 = cos(dl);
277  sindl1 = sin(dl);
278  iterations = 0;
279  do
280  {
281  cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosdl1;
282  sigma = acos(cosSigma);
283  azimuthEQ = asin((cosU1 * cosU2 * sindl1)/sin(sigma));
284 
285  /*
286  * Patch from Patrica Tozer to handle minor
287  * mathematical stability problem
288  */
289  TEMP = cosSigma - (2.0 * sinU1 * sinU2)/
290  (cos(azimuthEQ)*cos(azimuthEQ));
291  if (TEMP > 1)
292  {
293  TEMP = 1;
294  }
295  else if (TEMP < -1)
296  {
297  TEMP = -1;
298  }
299  tsm = acos(TEMP);
300 
301 
302  /* (old code?)
303  tsm = acos(cosSigma - (2.0 * sinU1 * sinU2)/(cos(azimuthEQ)*cos(azimuthEQ)));
304  */
305 
306  dl2 = deltaLongitude(azimuthEQ, sigma, tsm,sphere);
307  dl3 = dl1 - (dl + dl2);
308  dl1 = dl + dl2;
309  cosdl1 = cos(dl1);
310  sindl1 = sin(dl1);
311  iterations++;
312  }
313  while ( (iterations<999) && (fabs(dl3) > 1.0e-032));
314 
315  /* compute expansions A and B */
316  u2 = mu2(azimuthEQ,sphere);
317  A = bigA(u2);
318  B = bigB(u2);
319 
320  /* compute length of geodesic */
321  dsigma = B * sin(sigma) * (cos(tsm) +
322  (B*cosSigma*(-1.0 + 2.0 * (cos(tsm)*cos(tsm))))/4.0);
323  return sphere->b * (A * (sigma - dsigma));
324 }
double b
Definition: liblwgeom.h:316
double bigB(double u2)
double mu2(double azimuth, SPHEROID *sphere)
double f
Definition: liblwgeom.h:317
double bigA(double u2)
double deltaLongitude(double azimuth, double sigma, double tsm, SPHEROID *sphere)
uint32_t e
Definition: geobuf.h:57
Here is the call graph for this function:
Here is the caller graph for this function: