```function [dist,lons,lats] = ll2dist(long,lat,N)
% ll2dist (originally M_LLDIST from the m_map package)
% Spherical earth distance between points in long/lat coordinates.
%   RANGE=M_LLDIST(LONG,LAT) gives the distance in kilometers between
%   successive points in the vectors LONG and LAT, computed
%   using the Haversine formula on a spherical earth of radius
%   6378.137km. Distances are probably good to better than 1% of the
%   "true" distance on the ellipsoidal earth
%
%   [RANGE,LONGS,LATS]=M_LLDIST(LONG,LAT,N) computes the N-point geodesics
%   between successive points. Each geodesic is returned on its
%   own row of length N+1.
%
% Rich Pawlowicz (rich@ocgy.ubc.ca) 6/Nov/00
% This software is provided "as is" without warranty of any kind. But
% it's mine, so you can't sell it.
%
% 30/Dec/2005 - added n-point geodesic computations, based on an algorithm
%               coded by Jeff Barton at Johns Hopkins APL in an m-file
%               I looked at at mathworks.com.

pi180=pi/180;

m=length(long)-1;

long1=reshape(long(1:end-1),m,1)*pi180;
long2=reshape(long(2:end)  ,m,1)*pi180;
lat1= reshape(lat(1:end-1) ,m,1)*pi180;
lat2= reshape(lat(2:end)   ,m,1)*pi180;

dlon = long2 - long1;
dlat = lat2 - lat1;
a = (sin(dlat/2)).^2 + cos(lat1) .* cos(lat2) .* (sin(dlon/2)).^2;
angles = 2 * atan2( sqrt(a), sqrt(1-a) );

if nargin==3 & nargout>1,   % Compute geodesics.

% Cartesian unit vectors in rows of v1,v2
v1=[cos(long1).*cos(lat1)   sin(long1).*cos(lat1)   sin(lat1) ];
v2=[cos(long2).*cos(lat2)   sin(long2).*cos(lat2)   sin(lat2) ];

% We want to get a unit vector tangent to the great circle.
n1=cross(v1,v2,2);
t1=cross(n1,v1,2);
t1=t1./repmat(sqrt(sum(t1.^2,2)),1,3);

lons=zeros(m,N+1);
lats=zeros(m,N+1);
for k=1:m,