RECKON - Using the WGS-84 Earth ellipsoid, travel a given distance along a given azimuth starting at a given initial point, and return the endpoint within a few millimeters of accuracy, using Vincenty's algorithm. USAGE: [lat2,lon2] = vreckon(lat1, lon1, s, a12) VARIABLES: lat1 = inital latitude (degrees) lon1 = initial longitude (degrees) s = distance (meters) a12 = intial azimuth (degrees) lat2, lon2 = second point (degrees) a21 = reverse azimuth (degrees), at final point facing back toward the intial point Original algorithm source: T. Vincenty, "Direct and Inverse Solutions of Geodesics on the Ellipsoid with Application of Nested Equations", Survey Review, vol. 23, no. 176, April 1975, pp 88-93. Available at: http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf Notes: (1) The Vincenty reckoning algorithm was transcribed verbatim into JavaScript by Chris Veness. It was modified and translated to Matlab by Michael Kleder. Mr. Veness's website is: http://www.movable-type.co.uk/scripts/latlong-vincenty-direct.html (2) Error correcting code, polar error corrections, WGS84 ellipsoid parameters, testing, and comments by Michael Kleder. (3) By convention, when starting at a pole, the longitude of the initial point (otherwise meaningless) determines the longitude line along which to traverse, and hence the longitude of the final point. (4) The convention noted in (3) above creates a discrepancy with VDIST when the the intial or final point is at a pole. In the VDIST function, when traversing from a pole, the azimuth is 0 when heading away from the south pole and 180 when heading away from the north pole. In contrast, this VRECKON function uses the azimuth as noted in (3) above when traversing away form a pole. (5) In testing, where the traversal subtends no more than 178 degrees, this function correctly inverts the VDIST function to within 0.2 millimeters of distance, 5e-10 degrees of forward azimuth, and 5e-10 degrees of reverse azimuth. Precision reduces as test points approach antipodal because the precision of VDIST is reduced for nearly antipodal points. (A warning is given by VDIST.) (6) Tested but no warranty. Use at your own risk. (7) Ver 1.0, Michael Kleder, November 2007
0001 function [lat2,lon2,a21] = vreckon(lat1,lon1,s,a12) 0002 % RECKON - Using the WGS-84 Earth ellipsoid, travel a given distance along 0003 % a given azimuth starting at a given initial point, and return 0004 % the endpoint within a few millimeters of accuracy, using 0005 % Vincenty's algorithm. 0006 % 0007 % USAGE: 0008 % [lat2,lon2] = vreckon(lat1, lon1, s, a12) 0009 % 0010 % VARIABLES: 0011 % lat1 = inital latitude (degrees) 0012 % lon1 = initial longitude (degrees) 0013 % s = distance (meters) 0014 % a12 = intial azimuth (degrees) 0015 % lat2, lon2 = second point (degrees) 0016 % a21 = reverse azimuth (degrees), at final point facing back toward the 0017 % intial point 0018 % 0019 % Original algorithm source: 0020 % T. Vincenty, "Direct and Inverse Solutions of Geodesics on the Ellipsoid 0021 % with Application of Nested Equations", Survey Review, vol. 23, no. 176, 0022 % April 1975, pp 88-93. 0023 % Available at: http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf 0024 % 0025 % Notes: 0026 % (1) The Vincenty reckoning algorithm was transcribed verbatim into 0027 % JavaScript by Chris Veness. It was modified and translated to Matlab 0028 % by Michael Kleder. Mr. Veness's website is: 0029 % http://www.movable-type.co.uk/scripts/latlong-vincenty-direct.html 0030 % (2) Error correcting code, polar error corrections, WGS84 ellipsoid 0031 % parameters, testing, and comments by Michael Kleder. 0032 % (3) By convention, when starting at a pole, the longitude of the initial 0033 % point (otherwise meaningless) determines the longitude line along 0034 % which to traverse, and hence the longitude of the final point. 0035 % (4) The convention noted in (3) above creates a discrepancy with VDIST 0036 % when the the intial or final point is at a pole. In the VDIST 0037 % function, when traversing from a pole, the azimuth is 0 when 0038 % heading away from the south pole and 180 when heading away from the 0039 % north pole. In contrast, this VRECKON function uses the azimuth as 0040 % noted in (3) above when traversing away form a pole. 0041 % (5) In testing, where the traversal subtends no more than 178 degrees, 0042 % this function correctly inverts the VDIST function to within 0.2 0043 % millimeters of distance, 5e-10 degrees of forward azimuth, 0044 % and 5e-10 degrees of reverse azimuth. Precision reduces as test 0045 % points approach antipodal because the precision of VDIST is reduced 0046 % for nearly antipodal points. (A warning is given by VDIST.) 0047 % (6) Tested but no warranty. Use at your own risk. 0048 % (7) Ver 1.0, Michael Kleder, November 2007 0049 0050 % Input check: 0051 if abs(lat1)>90 0052 error('Input latitude must be between -90 and 90 degrees, inclusive.') 0053 end 0054 a = 6378137; % semimajor axis 0055 b = 6356752.31424518; % semiminor axis 0056 f = 1/298.257223563; % flattening coefficient 0057 lat1 = lat1 * .1745329251994329577e-1; % intial latitude in radians 0058 lon1 = lon1 * .1745329251994329577e-1; % intial longitude in radians 0059 % correct for errors at exact poles by adjusting 0.6 millimeters: 0060 kidx = abs(pi/2-abs(lat1)) < 1e-10; 0061 if any(kidx) 0062 lat1(kidx) = sign(lat1(kidx))*(pi/2-(1e-10)); 0063 end 0064 alpha1 = a12 * .1745329251994329577e-1; % inital azimuth in radians 0065 sinAlpha1 = sin(alpha1); 0066 cosAlpha1 = cos(alpha1); 0067 tanU1 = (1-f) * tan(lat1); 0068 cosU1 = 1 / sqrt(1 + tanU1*tanU1); 0069 sinU1 = tanU1*cosU1; 0070 sigma1 = atan2(tanU1, cosAlpha1); 0071 sinAlpha = cosU1 * sinAlpha1; 0072 cosSqAlpha = 1 - sinAlpha*sinAlpha; 0073 uSq = cosSqAlpha * (a*a - b*b) / (b*b); 0074 A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq))); 0075 B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq))); 0076 sigma = s / (b*A); 0077 sigmaP = 2*pi; 0078 while (abs(sigma-sigmaP) > 1e-12) 0079 cos2SigmaM = cos(2*sigma1 + sigma); 0080 sinSigma = sin(sigma); 0081 cosSigma = cos(sigma); 0082 deltaSigma = B*sinSigma*(cos2SigmaM+B/4*(cosSigma*(-1+... 0083 2*cos2SigmaM*cos2SigmaM)-... 0084 B/6*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+... 0085 4*cos2SigmaM*cos2SigmaM))); 0086 sigmaP = sigma; 0087 sigma = s / (b*A) + deltaSigma; 0088 end 0089 tmp = sinU1*sinSigma - cosU1*cosSigma*cosAlpha1; 0090 lat2 = atan2(sinU1*cosSigma + cosU1*sinSigma*cosAlpha1,... 0091 (1-f)*sqrt(sinAlpha*sinAlpha + tmp*tmp)); 0092 lambda = atan2(sinSigma*sinAlpha1, cosU1*cosSigma - ... 0093 sinU1*sinSigma*cosAlpha1); 0094 C = f/16*cosSqAlpha*(4+f*(4-3*cosSqAlpha)); 0095 L = lambda - (1-C) * f * sinAlpha * (sigma + C*sinSigma*(cos2SigmaM+... 0096 C*cosSigma*(-1+2*cos2SigmaM*cos2SigmaM))); 0097 lon2 = lon1 + L; 0098 % output degrees 0099 lat2 = rad2deg(lat2); 0100 lon2 = rad2deg(lon2); 0101 lon2 = mod(lon2,360); % follow [0,360] convention 0102 if nargout > 2 0103 a21 = atan2(sinAlpha, -tmp); 0104 a21 = 180 + rad2deg(a21); % note direction reversal 0105 a21=mod(a21,360); 0106 end 0107 0108 end % function