Index: applications/editors/josm/plugins/utilsplugin/src/UtilsPlugin/SimplifyWayAction.java
===================================================================
--- applications/editors/josm/plugins/utilsplugin/src/UtilsPlugin/SimplifyWayAction.java	(revision 5142)
+++ applications/editors/josm/plugins/utilsplugin/src/UtilsPlugin/SimplifyWayAction.java	(revision 5226)
@@ -103,8 +103,8 @@
 		for (int i = from+1; i < to; i++) {
 			Node n = wnew.nodes.get(i);
-			double xte = radtometers(linedist(
+			double xte = Math.abs(EARTH_RAD * xtd(
 				fromN.coor.lat(), fromN.coor.lon(),
-				n.coor.lat(), n.coor.lon(),
-				toN.coor.lat(), toN.coor.lon()));
+				toN.coor.lat(), toN.coor.lon(),
+				n.coor.lat(), n.coor.lon()));
 			if (xte > xtemax) {
 				xtemax = xte;
@@ -120,166 +120,24 @@
 	}
 
-	/* ---------------------------------------------------------------------- 
-	 * Everything below this comment was converted from C to Java by Frederik
-	 * Ramm. The original sources are the files grtcirc.c and smplrout.c from 
-	 * the gpsbabel source code (www.gpsbabel.org), which is under GPL. The
-	 * relevant code portions have been written by Robert Lipe.
-	 * 
-	 * Method names have been left unchanged where possible.
+	public static double EARTH_RAD = 6378137.0;
+
+	/* From Aviaton Formulary v1.3
+	 * http://williams.best.vwh.net/avform.htm
 	 */
-	
-	public static double EARTH_RAD = 6378137.0;
-	public static double radmiles = EARTH_RAD*100.0/2.54/12.0/5280.0;
-
-	public static double[] crossproduct(double[] v1, double[] v2) {
-		double[] rv = new double[3];
-		rv[0] = v1[1]*v2[2]-v2[1]*v1[2];
-		rv[1] = v1[2]*v2[0]-v2[2]*v1[0];
-		rv[2] = v1[0]*v2[1]-v1[1]*v2[0];
-		return rv;
+	public static double dist(double lat1, double lon1, double lat2, double lon2) {
+		return 2*Math.asin(Math.sqrt(Math.pow(Math.sin((lat1-lat2)/2), 2) + 
+			Math.cos(lat1)*Math.cos(lat2)*Math.pow(Math.sin((lon1-lon2)/2), 2)));
 	}
 
-	public static double dotproduct(double[] v1, double[] v2) {
-		return v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2];
+	public static double course(double lat1, double lon1, double lat2, double lon2) {
+		return Math.atan2(Math.sin(lon1-lon2)*Math.cos(lat2),
+		    Math.cos(lat1)*Math.sin(lat2)-Math.sin(lat1)*Math.cos(lat2)*Math.cos(lon1-lon2)) % (2*Math.PI);
 	}
 
-	public static double radtomiles(double rads) {
-		return (rads*radmiles);
-	}
-
-	public static double radtometers(double rads) {
-		return (rads * EARTH_RAD);
-	}
-	
-	public static double veclen(double[] vec) {
-		return Math.sqrt(vec[0]*vec[0]+vec[1]*vec[1]+vec[2]*vec[2]);
-	}
-
-	public static double gcdist(double lat1, double lon1, double lat2, double lon2) 
-	{
-		double res;
-		double sdlat, sdlon;
-
-		sdlat = Math.sin((lat1 - lat2) / 2.0);
-		sdlon = Math.sin((lon1 - lon2) / 2.0);
-
-		res = Math.sqrt(sdlat * sdlat + Math.cos(lat1) * Math.cos(lat2) * sdlon * sdlon);
-
-		if (res > 1.0) {
-			res = 1.0;
-		} else if (res < -1.0) {
-			res = -1.0;
-		}
-
-		res = Math.asin(res);
-		return 2.0 * res;
-	}
-
-	static double linedist(double lat1, double lon1, double lat2, double lon2, double lat3, double lon3) {
-
-		double dot;
-
-		/* degrees to radians */
-		lat1 = Math.toRadians(lat1);  lon1 = Math.toRadians(lon1);
-		lat2 = Math.toRadians(lat2);  lon2 = Math.toRadians(lon2);
-		lat3 = Math.toRadians(lat3);  lon3 = Math.toRadians(lon3);
-
-		/* polar to ECEF rectangular */
-		double[] v1 = new double[3];
-		double[] v2 = new double[3];
-		double[] v3 = new double[3];
-		v1[0] = Math.cos(lon1)*Math.cos(lat1); v1[1] = Math.sin(lat1); v1[2] = Math.sin(lon1)*Math.cos(lat1);
-		v2[0] = Math.cos(lon2)*Math.cos(lat2); v2[1] = Math.sin(lat2); v2[2] = Math.sin(lon2)*Math.cos(lat2);
-		v3[0] = Math.cos(lon3)*Math.cos(lat3); v3[1] = Math.sin(lat3); v3[2] = Math.sin(lon3)*Math.cos(lat3);
-
-		/* 'va' is the axis; the line that passes through the center of the earth
-		 * and is perpendicular to the great circle through point 1 and point 2 
-		 * It is computed by taking the cross product of the '1' and '2' vectors.*/
-		double[] va = crossproduct(v1, v2);
-		double la = veclen(va);
-
-		if (la != 0) {
-			va[0] /= la;
-			va[1] /= la;
-			va[2] /= la;
-
-			/* dot is the component of the length of '3' that is along the axis.
-			 * What's left is a non-normalized vector that lies in the plane of 
-			 * 1 and 2. */
-
-			dot = dotproduct(v3, va);
-
-			double[] vp = new double[3];
-			vp[0]=v3[0]-dot*va[0];
-			vp[1]=v3[1]-dot*va[1];
-			vp[2]=v3[2]-dot*va[2];
-
-			double lp = veclen(vp);
-
-			if (lp != 0) {
-
-				/* After this, 'p' is normalized */
-				vp[0] /= lp;
-				vp[1] /= lp;
-				vp[2] /= lp;
-
-				double[] cp1 = crossproduct(v1, vp);
-				double dp1 = dotproduct(cp1, va);
-
-				double[] cp2 = crossproduct(v2, vp);
-				double dp2 = dotproduct(cp2, va);
-
-				if ( dp1 >= 0 && dp2 >= 0 ) {
-					/* rather than call gcdist and all its sines and cosines and
-					 * worse, we can get the angle directly.  It's the arctangent
-					 * of the length of the component of vector 3 along the axis 
-					 * divided by the length of the component of vector 3 in the 
-					 * plane.  We already have both of those numbers. 
-					 * 
-					 * atan2 would be overkill because lp and Math.abs are both
-					 * known to be positive. */
-					return Math.atan(Math.abs(dot)/lp); 
-				}
-
-				/* otherwise, get the distance from the closest endpoint */
-				double c1 = dotproduct(v1, vp);
-				double c2 = dotproduct(v2, vp);
-				dp1 = Math.abs(dp1);
-				dp2 = Math.abs(dp2);
-
-				/* This is a hack.  d$n$ is proportional to the sine of the angle
-				 * between point $n$ and point p.  That preserves orderedness up
-				 * to an angle of 90 degrees.  c$n$ is proportional to the cosine
-				 * of the same angle; if the angle is over 90 degrees, c$n$ is
-				 * negative.  In that case, we flop the sine across the y=1 axis
-				 * so that the resulting value increases as the angle increases. 
-				 * 
-				 * This only works because all of the points are on a unit sphere. */
-
-				if (c1 < 0) {
-					dp1 = 2 - dp1;
-				}
-				if (c2 < 0) {
-					dp2 = 2 - dp2;
-				}
-
-				if (Math.abs(dp1) < Math.abs(dp2)) {
-					return gcdist(lat1,lon1,lat3,lon3);  
-				} else {
-					return gcdist(lat2,lon2,lat3,lon3);
-				}
-			} else {
-				/* lp is 0 when 3 is 90 degrees from the great circle */
-				return Math.PI/2;
-			}    
-		} else {
-			/* la is 0 when 1 and 2 are either the same point or 180 degrees apart */
-			dot = dotproduct(v1, v2);
-			if (dot >= 0) { 
-				return gcdist(lat1,lon1,lat3,lon3);
-			} else {
-				return 0;
-			}
-		}
+	public static double xtd(double lat1, double lon1, double lat2, double lon2, double lat3, double lon3) {
+		double dist_AD = dist(lat1, lon1, lat3, lon3);
+		double crs_AD = course(lat1, lon1, lat3, lon3);
+		double crs_AB = course(lat1, lon1, lat2, lon2);
+		return Math.asin(Math.sin(dist_AD)*Math.sin(crs_AD-crs_AB));
 	}
 }