%GREATCIRCLE  Finds the shortest distance between two points on a sphere.
%
%    Given the inputs of latitude and longitude (in degrees) for points 
%    A and B, this function returns a vector with the angle separating these 
%    two points (in radians) and the distance separating the two points 
%    (in kilometres).
%
%    To call this function, enter
%
%    >> greatcircle(A,B)
%
%    where A and B are 2-dimensional vectors conataining information 
%    as [lat long].
%
%    Follows p135 of {Lay and Wallace (1995), Modern Global Seismology}.
%
%       So a ham sandwich walks into a bar. The bartender says, 
%       "Sorry, we don't serve food here."
%

function [distance] = greatcircle(X1,X2)

A = abs(deg2rad(X2(2)-X1(2)));
b = pi/2 - deg2rad(X1(1));
c = pi/2 - deg2rad(X2(1));

a = acos(cos(b) * cos(c) + sin(b) * sin(c) * cos(A));

Delta = 6371 * a;

distance = [a Delta];