File Exchange

image thumbnail

VDISTINV: Find the endpoint of a geodesic on the ellipsoidal earth

version 1.0.0.0 (4.22 KB) by Michael Kleder
Obtains a final location on ellipsoidal earth, given a start point, an azimuth and a distance.

2 Downloads

Updated 24 Apr 2006

No License

This function solves the "forward geodesic problem," which is to compute the endpoint of a geodesic (shortest-distance) path on the ellipsoidal earth, given the start point, a path length, and a starting azimuth.

In response to a user request, this function numerically inverts the previously published VDIST function, which is contained within. (VDIST solves the "inverse geodesic problem" of computing the distance and azimuth between two known points, using Vincenty's 1975 algorithm.) Although numerically inverting VDIST is not a very efficient way to solve the forward problem, it provides a quick solution and also serves as an example of how to invert a function numerically in Matlab.

OTHER NOTES: The distance between essentially antipodal points on an ellipsoid is very sensitive to small deviations in azimuth, so such points should be avoided. (A warning is given.) For other cases, precision is set to about one part in 10^12. This function has been tested but no warranty is provided; use at your own risk. Written by Michael Kleder, April 2006.

Cite As

Michael Kleder (2020). VDISTINV: Find the endpoint of a geodesic on the ellipsoidal earth (https://www.mathworks.com/matlabcentral/fileexchange/10821-vdistinv-find-the-endpoint-of-a-geodesic-on-the-ellipsoidal-earth), MATLAB Central File Exchange. Retrieved .

Comments and Ratings (4)

The code below is much faster and did not have any errors while the VDISTINV did not seem to always work.

felipe diniz

John Lansberry

It would be interesting to compare this to the NOAA code that I translated to an m-file and vectorized.

function [lat2_out, lon2_out, baz_out] = wgs84invdist(varargin)
%------------------------------------------------------------------------------
% function [lat2, lon2, baz] = wgs84invdist(lat1, lon1, faz, s [, deg_in, deg_out])
%
% This code is part of the geodetic tool kit, which can be downloaded
% from the National Oceanic and Atmospheric Administration (NOAA)
% National Geodetic Survey (ngs) website at:
% http://www.ngs.noaa.gov/index.shtml.
%
% Input Description
%
% lat1 initial latitude (rad or deg)
% lon1 initial longitude (rad or deg)
% faz forward azimuth (rad or deg)
% s surface distance along ellipsoid (m)
% deg_in selects whether input angular quantities are in degrees or
% radians (deg = 1 for degrees, radians is default)
% deg_out selects whether output angular quantities are in degrees or
% radians (deg = 1 for degrees, radians is default)
%
% Output Description
%
% lat2 final latitude (rad or deg)
% lon2 final longitude (rad or deg)
% baz backward azimuth (rad or deg)
%------------------------------------------------------------------------------

%------------------------------------------------------------------------------
% Original source code comments:
%
% Solution of the geodetic direct problem after T.Vincenty
% modified Rainsford's method with Helmert's elliptical terms
% effective in any azimuth and at any distance short of antipodal
% standpoint/forepoint must not be the geographic pole
%
% a is the semi-major axis of the reference ellipsoid
% f is the flattening of the reference ellipsoid
% latitudes and longitudes in radians positive north and east
% forward azimuths at both points returned in radians from north
%
% programmed for cdc-6600 by LCDR L.Pfeifer NGS Rockville MD 18FEB75
% modified for system 360 by John G Gergen NGS Rockville MD 7507
%------------------------------------------------------------------------------

%-------------------------------------------------------------------------------
% check input arguments
%-------------------------------------------------------------------------------
error(nargchk(4, 6, nargin));

%-------------------------------------------------------------------------------
% check for input errors (e.g., inputs have incorrect number of rows or columns)
%-------------------------------------------------------------------------------
rows = zeros(1,nargin);
cols = zeros(1,nargin);

for ii=1:nargin,
[rows(ii), cols(ii)] = size(varargin{ii});
end

maxrows = max(rows);

% make sure each input has either 1 row or n rows
inputerr = any(rows ~= 1 & rows ~= maxrows);

if inputerr,
error('inputs must have either 1 row or n rows')
end

% all inputs must have 1 column
inputerr = any(cols ~= 1);

if inputerr
error('all inputs must have one column');
end

%-------------------------------------------------------------------------------
% expand inputs as necessary
%-------------------------------------------------------------------------------
expand = any(rows == 1) & (maxrows > 1);

if expand,
for ii=1:nargin,
if rows(ii) == 1,
varargin{ii} = repmat(varargin{ii}, maxrows, 1);
end
end
end

%------------------------------------------------------------------------------
% extract inputs from varargin
%------------------------------------------------------------------------------
lat1 = varargin{1};
lon1 = varargin{2};
faz = varargin{3};
s = varargin{4};

deg_in = zeros(maxrows, 1);
deg_out = zeros(maxrows, 1);

if nargin >= 5,
deg_in = varargin{5};
if nargin == 6,
deg_out = varargin{6};
end
end

%------------------------------------------------------------------------------
% WGS-84 defining parameters.
%------------------------------------------------------------------------------
a = 6378137.0;
f = 1.0 / 298.257223563;

%------------------------------------------------------------------------------
% Miscellaneous parameters.
%------------------------------------------------------------------------------
rad2deg = 180 / pi;
deg2rad = pi / 180;

eps = 0.5e-13;

zero = 0.0; one = 1.0; two = 2.0;three = 3.0; four = 4.0; six = 6.0;
three_eighths = 3.0 / 8.0;sixteen = 16.0;twopi = 2.0 * pi;

%------------------------------------------------------------------------------
% Input conversions.
%------------------------------------------------------------------------------
lat1(deg_in == 1) = lat1(deg_in == 1) * deg2rad;
lon1(deg_in == 1) = lon1(deg_in == 1) * deg2rad;
faz(deg_in == 1) = faz(deg_in == 1) * deg2rad;

%------------------------------------------------------------------------------
% find cases where s = 0 and remove them from computations
%------------------------------------------------------------------------------
szero = s == 0;

lat2_out = zeros(size(lat1));
lon2_out = zeros(size(lat1));
baz_out = zeros(size(lat1));

lat2_out(szero) = lat1(szero);
lon2_out(szero) = lon1(szero);

lat1 = lat1(~szero);
lon1 = lon1(~szero);
faz = faz(~szero);

%------------------------------------------------------------------------------
% Main routine.
%------------------------------------------------------------------------------
r = one - f;
tu = r * sin(lat1) ./ cos(lat1);
sf = sin(faz);
cf = cos(faz);

baz = zeros(size(faz));

if (cf ~= zero),
baz = atan2(tu, cf) * two;
end

cu = one ./ sqrt(tu .* tu + one);
su = tu .* cu;
sa = cu .* sf;
c2a = -sa .* sa + one;
x = sqrt((one / r / r - one) * c2a + one) + one;
x = (x - two) ./ x;
c = one - x;
c = (x .* x / four + one) ./ c;
d = (three_eighths * x .* x - one) .* x;
tu = s ./ a ./ c / r;
y = tu;

repeat = 1;

while repeat == 1,

sy = sin(y);
cy = cos(y);
cz = cos(baz + y);
e = cz .* cz * two - one;
c = y;
x = e .* cy;
y = e + e - one;
y = (((sy .* sy * four - three) .* y .* cz .* d / six + x) .* d / four - cz) .* sy .* d + tu;

if all(abs(y - c) <= eps),
break
end

end

baz = cu .* cy .* cf - su .* sy;
c = r * sqrt(sa .* sa + baz .* baz);
d = su .* cy + cu .* sy .* cf;
lat2 = atan2(d, c);
c = cu .* cy - su .* sy .* cf;
x = atan2(sy .* sf, c);
c = ((-three * c2a + four) .* f + four) .* c2a .* f / sixteen;
d = ((e .* cy .* c + cz) .* sy .* c + y) .* sa;
lon2 = lon1 + x - (one - c) .* d .* f;
baz = atan2(sa, baz) + pi;

%------------------------------------------------------------------------------
% Map to outputs.
%------------------------------------------------------------------------------
lat2_out(~szero) = lat2;
lon2_out(~szero) = lon2;
baz_out(~szero) = baz;

ii = lon2_out < -pi;
lon2_out(ii) = lon2_out(ii) + twopi;

ii = lon2_out > pi;
lon2_out(ii) = lon2_out(ii) - twopi;

%------------------------------------------------------------------------------
% Output conversions.
%------------------------------------------------------------------------------
lat2_out(deg_out == 1) = lat2_out(deg_out == 1) * rad2deg;
lon2_out(deg_out == 1) = lon2_out(deg_out == 1) * rad2deg;
baz_out(deg_out == 1) = baz_out(deg_out == 1) * rad2deg;

if nargout == 1,
lat2_out = [lat2_out lon2_out baz_out];
end

Barbara Kobzik

Very helpful function, numerically inverts VDIST. Due to numerical optimization, it is a bit slower (according to author 2ms for VDIST, 250ms for VDISTINV), but accuracy very good. Author was extremely helpful when contacted for further questions. Thank you very much!

MATLAB Release Compatibility
Created with R14SP3
Compatible with any release
Platform Compatibility
Windows macOS Linux