function D = dec_deg_2_miles_matrix(L)

%Takes locations in decimal-degree coordinates and calculates pairwise
%distances in miles. 

% INPUTS: L = (xi,yi) = (lon(i),lat(i)) decimal-degree coords.

%*******************************************************************
% NOTE: For the US, xi is NEGATIVE and yi is POSITIVE. Be sure they
%       are in the right order.
%*******************************************************************

% OUTPUT: D = nx(n-1)/2 vector of distinct pairwise distances

% DIST_VEC computes the vector of distinct pairwise distances
%          for a given location vector, loc.

%dbstop if warning

n = length(L(:,1)) ; % number of locations

D = zeros(n,n);

%Compute pairwise distances

for i = 1:n
    
    for j =1:n
        
       lon1 = L(i,1); 
       lat1 = L(i,2); 
       
      lon2 = L(j,1); 
      lat2 = L(j,2); 

      if (lon1 == lon2)&(lat1 == lat2)
          
          D(i,j) = D(i,j);
          
      else          

          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);
          D(i,j) = dist * 60 * 1.1515;

      end

   end
     
end