/** * Compute Distance and Bearings Between Two Points on The Earth * @author Robert John Morton * @version 20 August 1998 revised 30 March 2010 */ /* Solution to the geodetic inverse problem after T. Vincenty. Modified Rains- ford's Method with Helmert's elliptical terms. Effective in any azimuth and at any distance short of antipodal. Neither the observer nor the observed can be at a geographic pole. Latitudes and longitudes in radians, positive north and east. Distance between observer and observed is given in metres. Forward azimuths at observer and observed returned in radians clockwise from north. Originally programmed for CDC6600 by Lcdr L Pfeifer NGS Rock- ville MD 18FEB75. Modified for IBM 360 by John G Gergen NGS Rockville MD JUL75 */ class dandb{ private static final double π = 3.14159265358979323846264338327950288, // convergence limit for distance accuracy better than 1mm δ = 0.5e-13, // polar limit (1 minute < right-angle) P = 1.57050543858623089763516774317834; private static dandb db; // reference to current instance of this class private double α1, // forward azimuth α2, // back-azimuth s, // ellipsoidal distance observer to observed a, // equitorial radius of the Earth in metres f, // flattening factor r; // equitorial radius * r = polar radius private int ic; /* Iteration count for vincenty loop: 3 to 6 iterations are normally sufficient to achieve convergence. */ // set which geo-ellipsoid to use (above) dandb(double a, double f){ this.a = a; // equitorial radius of the Earth in metres this.f = f; // flattening factor r = 1 - f; // equitorial radius * r = polar radius db = this; // reference to current instance of this class } public void setaf(double a, double f){ this.a = a; // set new value for semi-major axis this.f = f; // set new value for flattening r = 1 - f; // equitorial radius * r = polar radius } public double getDist(){return s;} // get ellipsoidal distance in metres // return bearing of observed from observer public double getForwardAzimuth(){return α1;} // return bearing of observer from observed public double getBackAzimuth(){return α2;} public int getIterations(){return ic;} // to check for non-convergence /* Code fragment for calling the Vincenty dandb class: dandb db = dandb.getCurrent(); int x = db.vincenty(φ1, λ1, φ2, λ2); // compute distance & bearings if(x == 0) { // if the dandb computation is valid if((s = db.getDist()) != 0) { // provided distance isn't zero s /= 1000; // if you need distance in kilometres, or... s /= 6366197.723857773; // if you need distance in mean Earth-radians α1 = db.getForwardAzimuth(); // both azumuths are valid // else the points are coincident so dandb is invalid so leave α1, α2 as they were last time. α2 = db.getBackAzimuth(); } } Note that if s = 0, then the observer and the observed are coincident and therefore the azimuths are meaningless. In this case, for flight control applications, leave both azimuths at their previous values. If the computation returns a non-zero error code, look at the premature returns in the following method to see their meanings. */ public int vincenty( // VINCENTY INVERSE METHOD double φ1, // latitude of observer double λ1, // longitude of observer double φ2, // latitude of observed double λ2 // longitude of observed ) { if(φ1 > P) return 1; // observer's latitude > 89°59'00"N else if(φ1 < -P) return 2; // observer's latitude > 89°59'00"S else if(φ2 > P) return 3; // observed's latitude > 89°59'00"N else if(φ2 < -P) return 4; // observed's latitude > 89°59'00"S else if(λ1 > π || λ1 < -π) // observer's longitude > 180° E or W return 5; else if(λ2 > π || λ2 < -π) // observed's longitude > 180° E or W return 6; double Δλ = λ2 - λ1, // longitude difference A = r * Math.tan(φ1), // tan(reduced latitude of observer) B = r * Math.tan(φ2), // tan(reduced latitude of observed) C = 1 / Math.sqrt(1 + A * A), // cos(reduced latitude of observer) D = C * A, // sin(reduced latitude of the observer) E = 1 / Math.sqrt(1 + B * B), // cos(reduced latitude of observed) θ = Δλ, // prime θ before entering loop χ, // previous value of θ for each iteration F, G, H, I, J, K, ψ, e, y, c; s = C * E; // prime ellipsoidal distance (observer to observed) α2 = s * B; // prime the back-azimuth α1 = α2 * A; // prime the forward azimuth ic = 0; // zero the iteration counter /* ITERATE THE APPROPRIATE EQUATIONS Until there is no further signi- ficant change in θ or until 100 iterations, whichever occurs first. */ do{ H = Math.sin(θ); I = Math.cos(θ); A = E * H; B = α2 - D * E * I; F = Math.sqrt(A * A + B * B); G = s * I + α1; y = Math.atan2(F, G); J = s * H / F; K = 1 - J * J; ψ = α1 + α1; if(K > 0) ψ = G - ψ / K; e = ψ * ψ * 2 - 1; c = ((4 - 3 * K) * f + 4) * K * f / 16; χ = θ; // save old θ before computing new value θ = (1 - c) * ((e * G * c + ψ) * F * c + y) * J * f + Δλ } while( Math.abs(χ - θ) > δ && ++ic < 100 ); if(ic >= 100) // if θ didn't converge sufficiently within 100 iterations return 7; // return code 7 = non-convergence // distance calculation θ = Math.sqrt((1 / (r * r) - 1) * K + 1) + 1; θ = (θ - 2) / θ; c = (θ * θ / 4 + 1) / (1 - θ); θ = (0.375 * θ * θ - 1) * θ; s = ((((F * F * 4 - 3) * (1 - e - e) * ψ * θ / 6 - e * G) * θ / 4 + ψ) * F * θ + y) * c * a * r; if(Double.isNaN(s)) { // if distance calculation indeterminate if(Math.abs(Δλ) > 1) // points are antipodal (invalid) return 8; // return code 8 = antipodal s = 0; // the points must be coincident } else { // else distance is determinable, so // compute forward and back azimuths α1 = Math.atan2(A, B); α2 = Math.atan2(C * H, α2 * I - D * E) + π; } return 0; // return code 0 = success } }