function out = angdist(lat1,long1,lat2,long2);

% angdist.m
%
% calculates angular distance theta between points on a sphere specified 
% by latitude and longitude --
%
% syntax: theta = angdist(lat1,long1,lat2,long2);
%
% can have vector arguments
%
% uses formula:
%
% theta = arccos{cos(lat1)cos(lat2)[cos(long1)cos(long2) +...
%                   sin(long1)sin(long2)] + sin(lat1)sin(lat2)}
%
% which is devised by taking the dot product of two vectors and substituting
% polar coordinates. 
%
% Written by Greg Balco -- UW Cosmogenic Nuclide Lab
% balcs@u.washington.edu
% First version, March, 2001
% Checked October, 2005
% Part of the CRONUS-Earth online calculators: 
%      http://hess.ess.washington.edu/math
%
% Copyright 2001-2007, University of Washington
% All rights reserved
% Developed in part with funding from the National Science Foundation.
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License, version 2,
% as published by the Free Software Foundation (www.fsf.org).

lata = d2r(lat1);longa = d2r(long1);
latb = d2r(lat2);longb = d2r(long2);

theta = acos( (cos(lata).*cos(latb).*((cos(longa).*cos(longb)) + (sin(longa).*sin(longb)))) + (sin(lata).*sin(latb)));
out_1 = r2d(theta);

% make vector horizontal

if size(out_1,1) > size(out_1,2);
	out = out_1';
else;
	out = out_1;
end;

