function [d,alpha]=dist(lon1,lat1,lon2,lat2)
% function [d,alpha]=dist(lon1,lat1,lon2,lat2)
%
% distance [km] between points (lon1,lat1) and (lon2,lat2)
% and angle 'alpha' [degrees] of connecting line. longitude is
% east-longitude or the negative of west-longitude.
% spherical approx. (note by cw - this is a significant error)
%        alpha           (2) relative to (1)
% -------------------------------------------------
%          0                  eastward
%         90                  northward
%        180                  westward
%        270                  southward
% -------------------------------------------------
%
% r. schlitzer  (5/18/89)
fac=111.194929; alpha=0; dlon=lon2-lon1; dlat=lat2-lat1;
ddlon=fac*cos(2*pi/360*(lat1+lat2)/2)*dlon; ddlat=fac*dlat;
if ddlon==0 & dlat>0,     alpha=90;
elseif ddlon==0 & dlat<0, alpha=-90;
else,                     alpha=360/(2*pi)*atan(ddlat/ddlon); end
if ddlon<0, alpha=alpha+180; end, if alpha<0, alpha=alpha+360; end
d=sqrt(ddlat^2+ddlon^2); 
