PostGIS  2.1.10dev-r@@SVN_REVISION@@
double distance_ellipse_calculation ( double  lat1,
double  long1,
double  lat2,
double  long2,
SPHEROID sphere 
)

Definition at line 233 of file lwgeom_spheroid.c.

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

Referenced by distance_ellipse().

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

Here is the call graph for this function:

Here is the caller graph for this function: