Code covered by the BSD License

GREATCIRCLE and LOXODROME: "As the crow flies" and rhumb line path, distance and bearing

François Beauducel (view profile)

30 Jul 2012 (Updated )

Computes path coordinates, distance and bearing on Earth along great circle or rhumb line.

[lat,lon,dist,bear]=loxodrome(varargin)
```function [lat,lon,dist,bear]=loxodrome(varargin)
%LOXODROME Rhumb line path and distance.
%
%	[LAT,LON]=LOXODROME(LAT1,LON1,LAT2,LON2) computes the line between two
%	points defined by their spherical coordinates latitude and longitude
%	(in decimal degrees), crossing all meridians of longitude at the same
%	angle, i.e. a path derived from an initial and constant bearing.
%	LAT and LON contain vectors of coordinates of the way points.
%
%	[LAT,LON,DISTANCE,BEARING]=LOXODROME(...) returns also vector of
%	distances (in km) along the path way, and constant bearing angle (in
%	degrees from North, a scalar).
%
%	[...]=LOXODROME(...,N) uses N intermediate points (default is 100).
%
%	Example:
%
%	  contour(0:359,-89:90,topo,[0,0],'k')
%	  [lat,lon,dis,bear] = loxodrome(48.8,2.3,35.7,139.7);
%	  hold on, plot(lon,lat,'r','linewidth',2), hold off
%	  title(sprintf('Paris to Tokyo = %g km - bear = %g N',dis(end),bear))
%
%
%
%	Author: Francois Beauducel <beauducel@ipgp.fr>
%	Created: 2012-10-30
%	Updated: 2012-11-08
%
%	References:
%	  http://en.wikipedia.org/wiki/Rhumb_line
%
%	Acknowledgments: Jorge David Taramona Perea, Jacques Vernin.

%	Development history:
%	  [2012-11-08] correction of the bearing angle calculation, with
%	     consequence to distance values. Detection by Jacques Vernin.
%
%
%	Redistribution and use in source and binary forms, with or without
%	modification, are permitted provided that the following conditions are
%	met:
%
%	   * Redistributions of source code must retain the above copyright
%	     notice, this list of conditions and the following disclaimer.
%	   * Redistributions in binary form must reproduce the above copyright
%	     notice, this list of conditions and the following disclaimer in
%	     the documentation and/or other materials provided with the distribution
%
%	THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
%	AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
%	IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
%	ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
%	LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
%	CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
%	SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
%	INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
%	CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
%	ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
%	POSSIBILITY OF SUCH DAMAGE.

if nargin < 4
error('Number of input arguments not correct.')
else
lat1 = varargin{1};
lon1 = varargin{2};
lat2 = varargin{3};
lon2 = varargin{4};
end

if nargin < 5
n = 100;
else
n = varargin{5};
end

if any(~cellfun(@isnumeric,varargin)) | any(~cellfun(@isscalar,varargin))
error('All input arguments must be numeric scalars.')
end

% length of one degree of latitude (in km)
degkm = 6371*pi/180;

% defines linear vectors of latitudes and longitudes
lat = arcgdd(linspace(gdd(lat1),gdd(lat2),n));
lon = linspace(lon1,lon2,n);

if nargout > 2
% computes the bearing angle
bear = atan2d(lon2 - lon1,arcgdd(lat2) - arcgdd(lat1));

% computes the distance along the path
if lat1 ~= lat2
dist = degkm*(lat - lat1)./cosd(bear);
else
dist = degkm*(lon - lon1).*cosd(lat1);
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = gdd(x)
%GDD Gudermann function in degree.

y =  atand(sinh(x*pi/180));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = arcgdd(x)
%ARCGDD Gudermann inverse function in degree.

%y =  180/pi * log(tand(45 + x/2));
y = atanh(sind(x))*180/pi;```