kungfoo / geohash-java

Implementation of GeoHashes in java. We try to be/stay compliant to the spec, as far as possible.
Other
982 stars 310 forks source link

VincentGeodesy.moveInDirection - issue when close to the international date line #9

Closed ChrisWooldridge closed 12 years ago

ChrisWooldridge commented 13 years ago

If I use:

WGS84Point point = new WGS84Point(0.0, 179.999);
point = VincentyGeodesy.moveInDirection(point , 90.0, 1000);

I get

Exception in thread "main" java.lang.IllegalArgumentException: The supplied coordinates (-3.4624406578435726E-18,180.00798315284126) are out of range. at ch.hsr.geohash.WGS84Point.(WGS84Point.java:26) at ch.hsr.geohash.util.VincentyGeodesy.moveInDirection(VincentyGeodesy.java:81) at WeatherLocationTest.(WeatherLocationTest.java:52) at WeatherLocationTest.main(WeatherLocationTest.java:28)

I believe a quick fix is:

public static WGS84Point moveInDirection(WGS84Point point, double bearingInDegrees, double distanceInMeters) {

if (bearingInDegrees < 0 || bearingInDegrees > 360) {
  throw new IllegalArgumentException("direction must be in (0,360)");
}

double a = 6378137, b = 6356752.3142, f = 1 / 298.257223563; // WGS-84
// ellipsiod
double alpha1 = bearingInDegrees * degToRad;
double sinAlpha1 = Math.sin(alpha1), cosAlpha1 = Math.cos(alpha1);

double tanU1 = (1 - f) * Math.tan(point.getLatitude() * degToRad);
double cosU1 = 1 / Math.sqrt((1 + tanU1 * tanU1)), sinU1 = tanU1 * cosU1;
double sigma1 = Math.atan2(tanU1, cosAlpha1);
double sinAlpha = cosU1 * sinAlpha1;
double cosSqAlpha = 1 - sinAlpha * sinAlpha;
double uSq = cosSqAlpha * (a * a - b * b) / (b * b);
double A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)));
double B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)));

double sinSigma = 0, cosSigma = 0, cos2SigmaM = 0;
double sigma = distanceInMeters / (b * A), sigmaP = 2 * Math.PI;
while (Math.abs(sigma - sigmaP) > 1e-12) {
  cos2SigmaM = Math.cos(2 * sigma1 + sigma);
  sinSigma = Math.sin(sigma);
  cosSigma = Math.cos(sigma);
  double deltaSigma = B
  * sinSigma
  * (cos2SigmaM + B
          / 4
          * (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM) - B / 6 * cos2SigmaM
                  * (-3 + 4 * sinSigma * sinSigma) * (-3 + 4 * cos2SigmaM * cos2SigmaM)));
  sigmaP = sigma;
  sigma = distanceInMeters / (b * A) + deltaSigma;
}

double tmp = sinU1 * sinSigma - cosU1 * cosSigma * cosAlpha1;
double lat2 = Math.atan2(sinU1 * cosSigma + cosU1 * sinSigma * cosAlpha1, (1 - f)
                         * Math.sqrt(sinAlpha * sinAlpha + tmp * tmp));
double lambda = Math.atan2(sinSigma * sinAlpha1, cosU1 * cosSigma - sinU1 * sinSigma * cosAlpha1);
double C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha));
double L = lambda - (1 - C) * f * sinAlpha
* (sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)));

double newLat = lat2 / degToRad;
double newLon = point.getLongitude() + L / degToRad;

newLon = (newLon > 180.0 ? 360.0 - newLon : newLon);
newLon = (newLon < -180.0 ? 360.0 + newLon : newLon);

return new WGS84Point(newLat, newLon);

}

kungfoo commented 13 years ago

Did you add the two lines at the end? Can you either create a clone and publish the patch there or create a patch for the fix? It would be a lot easier to see the changes there. This is in fact very simple with git:

$ git diff

Will yield the changes in patch format.

kungfoo commented 12 years ago

This has in fact been fixed.