# [Proj] Computing geodesics with Matlab; accuracy of spherical distances.

Charles Karney charles.karney at sri.com
Sun Dec 15 11:40:40 EST 2013

```About a year ago, I released a native Matlab implementation of the
geodesic algorithms in GeographicLib.  The standalone package is
available here:

http://www.mathworks.com/matlabcentral/fileexchange/39108

This works with Octave and Matlab.  The solution is vectorized; solving
1000000 inverse problems in parallel takes a few seconds on my home
machine.

This capability allow you to extend Tobler's 1964 paper

http://dx.doi.org/10.1111/j.0033-0124.1964.009_q.x

comparing spherical and ellipsoidal distances.  For example, a rather
accurate way of estimating the ellipsoidal distance is to transfer the
problem to a sphere of radius (2*a+b)/3 with the spherical latitude set
to

atan((1-f)^(3/2)*tan(phi))

(which is roughly the rectifying latitude).  The RMS relative error in
the distance is 0.05% and the max relative error is 0.11% (WGS84
ellipsoid).  The following short Matlab program illustrates
this result

ell=defaultellipsoid;                   % WGS84
m=1000000;                              % number of points
% sample points uniformly
lat1=asind(2*rand(m,1)-1); lon1=180*(2*rand(m,1)-1);
lat2=asind(2*rand(m,1)-1); lon2=180*(2*rand(m,1)-1);
% the ellipsoidal distance
sell=geoddistance(lat1,lon1,lat2,lon2,ell);
% the mean sphere
a=ell(1); f=ecc2flat(ell(2)); b=a*(1-f); sph=[(2*a+b)/3,0];
p=1.5;                                  % convert the latitude
mu1=atand((1-f)^p*tand(lat1)); mu2=atand((1-f)^p*tand(lat2));
% the spherical distance
ssph=geoddistance(mu1,lon1,mu2,lon2,sph);
err=(ssph-sell)./sell;                  % the relative error
100*[sqrt(mean(err(:).^2)),max(abs(err(:)))]

If the geographic latitude is used (set p=0 above), the errors are RMS =
0.18%, max = 0.56%.  (Note that Tobler gives too optimistic an estimate
of the error because he only uses 27 pairs of points in his tests,