Home > matlab > vreckon.m

vreckon

PURPOSE ^

RECKON - Using the WGS-84 Earth ellipsoid, travel a given distance along

SYNOPSIS ^

function [lat2,lon2,a21] = vreckon(lat1,lon1,s,a12)

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Generated on Sun 28-Jan-2018 22:54:10 by m2html © 2005