function theta = angdistr(lat1,long1,lat2,long2)

% angdistr.m
%
% calculates angular distance theta between points on a sphere specified 
% by latitude and longitude --
%
% syntax: out = angdistr(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@bgc.org
% First version of 'angdist.m', March, 2001
% Checked October, 2005
% Revised Oct 2015 as 'angdistr.m' to function in radians not degrees. 
% Part of the CRONUS-Earth online calculators: 
%      http://hess.ess.washington.edu/
%
% Copyright 2001-2007, University of Washington
% 2007- Greg Balco, Berkeley Geochronology Center
% 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).

theta = acos( (cos(lat1).*cos(lat2).*((cos(long1).*cos(long2)) + (sin(long1).*sin(long2)))) + (sin(lat1).*sin(lat2)));

if size(theta,1) > size(theta,2); 
    % Is vertical, transpose
    theta = theta';
end;
