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. (See Erastothenes Project Page ). 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');