% using https://www.mathworks.com/matlabcentral/fileexchange/44848-geod-yet-another-geodetic-toolbox

lat = linspace(0, 90)';
ell = get_ellipsoid('wgs84');

M = get_radius_meridional(lat, ell);
N = get_radius_normal(lat, ell);

Rc = get_radius_gaussian(lat, ell);
Rm = get_radius_mean(lat, ell);

R = get_radius_geocentric(lat, ell);

%tmp = [N R Rc Rm M];
%lab = {'Transverse', 'Geocentric', 'Gaussian', 'Mean', 'Meridional'};
tmp = [N R Rc M];
lab = {'Transverse', 'Geocentric', 'Gaussian & Mean', 'Meridional'};

%%
figure
hold on
%plot(lat, tmp./1000)
plot(lat, tmp(:,1)./1000, '-r')
plot(lat, tmp(:,2)./1000, '-k')
plot(lat, tmp(:,3)./1000, '-g')
plot(lat, tmp(:,4)./1000, '-b')
xlim([0 90])
xlabel('Latitude (degrees)')
ylabel('Radii of curvature (km)')
h = legend(lab, 'Location','northwest');
  set(h, 'Box','off')
grid on
hline(ell.a/1000, ':k')
hline(ell.b/1000, '--k')
text(77, ell.a/1000, sprintf('equatorial\nradius'))
text(82, ell.b/1000, sprintf('polar\nradius'))