[postgis-users] Selecting geocoded postal codes within a given distance ( miles for example )

[Norman Vine]

Paul Ramsey writes:
> 
> You are exercising the function in a way which breaks it :) The 
> workaround is to use two points which are *very*nearly* horizontal, but 
> not quite ((0,0) (1,0.000001)). I am pretty sure I found this bug once 
> as well -- tbe the geodetic length code is kinda complex are rarely 
> used, so I do not think we ever fixed it. Anyone with a fix will be 
> greeted with open arms.

I believe the problem(s) is a divide by zero in this case
< following off the top of my head needs checking >

//postgis_proj.c
double distance_ellipse(...)
{
......
 do {
  double sinSigma;  // new
 
  cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosdl1;
  sigma = acos(cosSigma);

  // new check for zero
  sinSigma = sin(sigma)
  if ( sinSigma < EPS )
     azimuthEQ = 0.0;
  else
     azimuthEQ = asin((cosU1 * cosU2 * sindl1)/sinSigma);
  cosAzimuthEQ = cos(azimuthEQ);
  // new check for zero
  if ( cosAzimuthEQ < EPS )
      tsm = 0.0;
  else
      tsm = acos(cosSigma - (2.0 * sinU1 * sinU2)/(cosAzimuthEQ*cosAzimuthEQ));
  dl2 = deltaLongitude(azimuthEQ, sigma, tsm,sphere);
  dl3 = dl1 - (dl + dl2);
  dl1 = dl + dl2;
  cosdl1 = cos(dl1);
  sindl1 = sin(dl1);
  iterations++;
 } while ( (iterations<999) && (fabs(dl3) > 1.0e-032));
....
}