/** * 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 Rainsford'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 returned in metres. Forward azimuths at observer and observed returned in radians clockwise from north. Originally programmed for CDC6600 by Lcdr L Pfeifer NGS Rockville 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 */ // CONSTRUCTOR: setS 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 } // get ellipsoidal distance in metres public double getDist(){return s;} //return bearing of observed from observer public double getForwardAzimuth(){return α1;} //return bearing of observer from observed public double getBackAzimuth(){return α2;} //to check for non-convergence public int getIterations(){return ic;} /* 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 α2 = db.getBackAzimuth(); //else the points are coincident } //so leave α1, α2 as they were last time } //dandb invalid so leave all as they were 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) // if observer's latitude > 89°59'00"N return 1; // return error code 1 else if(φ1 < -P) // if observer's latitude > 89°59'00"S return 2; // return error code 2 else if(φ2 > P) // if observed's latitude > 89°59'00"N return 3; // return error code 3 else if(φ2 < -P) // if observed's latitude > 89°59'00"S return 4; // return error code 4 else if(λ1 > π || λ1 < -π) // if observer's longitude > 180° E or W return 5; // return error code 5 else if(λ2 > π || λ2 < -π) // if observed's longitude > 180° E or W return 6; // return error code 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 χ, // holds 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 do { // do the following 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; θ = (1 - c) * ((e * G * c + ψ) * F * c + y) * J * f + Δλ χ = θ; // save old value of θ before computing its new value /* while there is still a further significant change in θ and while 100 iterations have not yet been completed */ } while(Math.abs(χ - θ) > δ && ++ic < 100); if(ic >= 100) // if θ didn't converge sufficiently within 100 iterations return 7; // return error code 7 // 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) // if the points are antipodal return 8; // return error code 8 s = 0; // otherwise the points must be coincident } else { // else the distance is determinate α1 = Math.atan2(A,B); // so compute forward and back azimuths α2 = Math.atan2(C * H,α2 * I - D * E) + π; } return 0; // success } }