% eda11_02
% distances on a sphere as an example
% of small number approximations
clear all;

RTOD = 180/pi;
DTOR = pi/180;

lat1 = 30*DTOR;
lon1 = 30*DTOR;
N = 101;
lat2min = 30*DTOR;
lat2max = 40*DTOR;
lat2 = lat2min + (lat2max-lat2min)*[1:N-1]'/(N-1);
lon2 = lat2;
latbar = (lat1 + lat2)/2;
Dlon = lon2 - lon1;
Dlat = lat2 - lat1;

cr = sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2).*cos(Dlon);
r = acos(cr);

coslatbar2 = cos(latbar).^2;
Dlat2 = (Dlat).^2;
Dlon2 = (Dlon).^2;
t1 = Dlat2 + (Dlon).^2 .* coslatbar2;
r1 = sqrt( t1 );

t2 = Dlat2.*Dlon2.*(0.25-(coslatbar2/6));
r2 = sqrt( t1-t2 );

figure(1);
clf;
subplot(2,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [30, 40, 0, 15 ] );
xlabel('latitude (deg)');
ylabel('r (deg)');
plot( RTOD*lat2, RTOD*r, '-', 'LineWidth', 5, 'Color', [0.8,0.8,0.8] );
plot( RTOD*lat2, RTOD*r1, 'k:', 'LineWidth', 2 );
plot( RTOD*lat2, RTOD*r2, 'k-', 'LineWidth', 2 );

subplot(2,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [30, 40, 0, 0.02 ] );
xlabel('latitude (deg)');
ylabel('error (deg)');
plot( RTOD*lat2, RTOD*(r1-r), 'k:', 'LineWidth', 2 );
plot( RTOD*lat2, RTOD*(r2-r), 'k-', 'LineWidth', 2 );
