function [d] = greatCircleDistance(lat1,lon1,lat2,lon2)
% 
% Calculates the great-circle distance between 2 lat/lon (decimal degrees)
% using the ‘haversine’ formula – that is, the shortest distance over the 
% earth’s surface – giving an ‘as-the-crow-flies’ distance between the
% points (ignoring elavation changes).
% 
% Haversine formula:
%   a = sin²(Δφ/2) + cos φ1 ⋅ cos φ2 ⋅ sin²(Δλ/2)
%   c = 2 ⋅ atan2( √a, √(1−a) )
%   d = R ⋅ c
% 
% where	φ is latitude, λ is longitude, R is earth’s radius (6,371km)
% angles need to be in radians to pass to trig functions
% 
% Ben Raanan, MBARI Dec 12, 2020

R = 6371e3; % earth’s radius in metres

lt1 = deg2rad(lat1);
lt2 = deg2rad(lat2);
dLat = deg2rad(lat2-lat1);
dLon = deg2rad(lon2-lon1);

a = sin(dLat/2).^2 + cos(lt1) .* cos(lt2) .* sin(dLon/2).^2;
c = 2 .* atan2(sqrt(a), sqrt(1-a));
d = R * c;
