function Chord_Dists_Mat=EarthChordDistances(llPoints)
% function Chord_Dists=EarthChordDistances(llPoints)
%Input is N by 2, each row a (long, lat) pair, -180<long<180; -90<lat<90.
%Output is a N by N matrix of chordal circle distances in KM (approximating the
%earth as a sphere), where the (i,j) entry is the distance between the ith
%and jth rows of the input N by 2 matric. So the diagonal is zero. 
%radius of the earth is taken as 6378.137



RR=6378.137; %radius of the earth in km
NN=length(llPoints(:,1));
%make a NN^2 by 4 matrix. Each row one of the NN^2 possible sets of
%two points:

Pts_paired_Vec=[kron(llPoints, ones(NN,1)), kron(ones(NN,1), llPoints)];

Half_Angles_AsVec=asin(sqrt(sin((Pts_paired_Vec(:,2)-Pts_paired_Vec(:,4))*pi/180/2).^2 + cos(Pts_paired_Vec(:,2)*pi/180).*cos(Pts_paired_Vec(:,4)*pi/180).*sin(abs(Pts_paired_Vec(:,1)-Pts_paired_Vec(:,3))*pi/180/2).^2));
Chords_as_vec=2*RR*sin(Half_Angles_AsVec);
Chord_Dists_Mat=reshape(Chords_as_vec, NN, NN);

