function distance = great_circle(x,y)

%GREAT_CIRCLE calculates great-circel distance in miles between two points
%with coordinates in decimal degrees.

%INPUTS: x = (x1,x2) = [lat1,lon1] in decimal degrees
%        y = (y1,y2) = [lat2,lon2] in decimal degrees

%%EXAMPLE: From: San Francisco (x = [37.68,-122.308]) 
%            To: New York (y = [40.670,-73.944])
% 
% dist = great_circle(x,y)
% 
% dist = 2565.9 (miles)


lat1 = x(1);
lon1 = x(2);
lat2 = y(1);
lon2 = y(2);


theta = lon1 - lon2;
dist = sin(deg2rad(lat1)) * sin(deg2rad(lat2)) + ...
    cos(deg2rad(lat1)) * cos(deg2rad(lat2)) * cos(deg2rad(theta));
dist = acos(dist);
dist = rad2deg(dist);
distance = dist * 60 * 1.1515;