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

Paul Ramsey pramsey at refractions.net
Thu Dec 19 17:37:11 PST 2002

```Norman,
A valiant try, but it did not seem to have any effect when I tested
it... could be my testing I suppose.
P.

Norman Vine wrote:
> 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));
> ....
> }
>
>
>
> _______________________________________________
> postgis-users mailing list
> postgis-users at postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-users

--
__
/
| Paul Ramsey
| Refractions Research
| Email: pramsey at refractions.net
| Phone: (250) 885-0632
\_

```