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

Norman Vine nhv at cape.com
Thu Dec 19 15:53:23 PST 2002

```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));
....
}

```