loxodrome.m

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:
%
%	  load topo
%	  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))
%
%
%	See also GREATCIRCLE.
%
%	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.
%
%	Copyright (c) 2012, François Beauducel, covered by BSD License.
%	All rights reserved.
%
%	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;