Calculation of the Earth's size - Noon project

The following program was written to calculate the circumference of the Earth from data collected during the spring noon project. This program has an advantage over the spreadsheet approach since the values for every site can be calculated at once. It is then possible to compare sites and ask questions about the data. Any feedback regarding this program would be greatly appreciated. 

% This program acccepts data from the Noon project and calculates
% the circumference of Earth.

% Data should be in the form of a matrix with the first column being
% the latitude of the site, with negative signs for southern
% latitudes, the second column being the longitude, with negative signs
% for western longitudes and the third column being the length of a shadow
% cast from a 1-meter standard.

clear

% Read the file
% Change these 2 lines to match your data file!!!
load('noon.dat');

% Assign columns of file to vectors
n=size(noon,1);
lat=noon(:,1);
long=noon(:,2);
shad=noon(:,3);

% Calculate circumference comparing each site (i) to all 
% other sites (j)
for i=1:n;
        for j=1:n;
                if lat(j)<1; 
                
% Calculate angle from shadow length
% assumes 1 meter gnomon, 57.3 converts radians to degrees 
                        ang(j)=90+atan(shad(j)/100)*57.3;
                else
                        ang(j)=90-atan(shad(j)/100)*57.3;
                end
% if at same latitude, skip calculations and put zero in                
                if lat(i)==lat(j);
                        ns(i,j)=0;
                        cir(i,j)=0;
                        
%If two different latitudes, but same shadow, skip calcs
                elseif shad(i)==shad(j);
                        ns(i,j)=0;
                        cir(i,j)=0;
                else
                
% Calculate north-south latitude, 111.133 km per degree of latitude
                        ns(i,j)=abs(lat(i)-lat(j))*111.133;

%Calculate circumference of earth
                        cir(i,j)=abs(ns(i,j)*360/(ang(j)-ang(i)));
                        
% Calculate percent error based on 40008 as true circumference
                        per(i,j)=(cir(i,j)-40008)*100/40008;
                end
        end
end

% Calculate averages and standard deviation for nonzero circumferences

for i=1:n
        circ2=nonzeros(cir(:,i));
        avg(i)=mean(circ2);     
        stdev(i)=std(circ2);
end

% Printout results
disp('     Summary of Circumference Results');
disp('Lat       Long     Avg.Cir   St.Dev  Perc.Err');
result(1,:)=reshape(lat,1,n);
result(2,:)=reshape(long,1,n);
result(3,:)=avg;
result(4,:)=stdev;
perct=(avg-40008)*100/40008;
result(5,:)=perct;
fid = fopen('noonresult','w');
fprintf(fid,'%4.2f  %4.2f  %5.0f  %3.0f  %3.0f\n',result);
fclose(fid);
fprintf('%4.2f    %4.2f    %5.0f      %3.0f      %3.0f\n',result);

disp('Average Circumference from all site averages');
avgall=sum(avg)/n;
disp(avgall);
disp ('Overall Percent error');
perctall=(avgall-40008)*100/40008;
disp (perctall);

% Plot the Average Circumference plotted by each site
% vs. the latitude of that site

plot(lat,avg);
title('Average Circumference by Site');
xlabel('Latitude');
ylabel('Avg. Circumference');