function [s,lon,lat] = m_geodesic(lon1,lat1,lon2,lat2,Npoints,varargin)% M_GEODESIC - Compute points along geodesics of an ellipsoidal earth.%%       [S,LON,LAT] = m_geodesic(lon1,lat1,lon2,lat2,Npoints,spheroid)%% lat1 = GEODETIC latitude of first point (degrees)% lon1 = longitude of first point (degrees)% lat2, lon2 = second point (degrees)% Npoints = number of points along geodesic% spheroid = (Optional) spheroid, defaults to 'wgs84'% S = distance in meters  % LON,LAT = points on geodesics.%% Note that inputs can be the same size, or a mixture of scalars and matrices.% However, output curves will be on columns of a matrix (i.e. form of% inputs will not be preserved).%% See M_FDIST, M_IDIST, M_LLDIST% R. pawlowicz 10/Jan/2005%% Do the equivalent of meshgrid-type calls to make all% matrix sizes match. To do this we can have no more than% two different numbers in any particular dimension, and one% of these has to be a '1'.allsize=[size(lon1);size(lat1);size(lon2);size(lat2)];i1=ones(1,size(allsize,2));for k=1:size(allsize,2), rs=unique(allsize(:,k)); if length(rs)==2 & rs(1)==1,   j1=i1;j1(k)=rs(2);   if allsize(1,k)==1,lon1=repmat(lon1,j1); end;   if allsize(2,k)==1,lat1=repmat(lat1,j1); end;   if allsize(3,k)==1,lon2=repmat(lon2,j1); end;   if allsize(4,k)==1,lat2=repmat(lat2,j1); end; elseif length(rs)>2,  error('incompatible array sizes!');   end;end;% Get the distances and bearings[S,A12,A21]=m_idist(lon1,lat1,lon2,lat2,varargin{:});s=[0:Npoints-1]'/(Npoints-1)*S(:)';[lon,lat]=m_fdist(lon1(:)',lat1(:)',A12(:)',s,varargin{:});