function [ang]=bearing(lat1,lon1,lat2,lon2)
% Input lat lon of two adjacent points in degrees
%   lat1=-30 ; % lat1 in degrees
%   lon1=50  ; % lon1 in degrees
%   lat2=-40 ; % lat2 in degrees
%   lon2=50  ; % lon2 in degrees
%% local pi, d2r, r2d, lat1r, lon1r, lat2r, lon2r, ang_tmp, lat1_tmp, lon1_tmp
% calculates the bearing (angle of movement over the sphere) between two points on a sphere
% output angle in degrees
% D.R. Thatcher 15/07/2014
% Magdalena Carranza: coded in MATLAB from bearing.ncl, downloaded from
% Colin Zarzycki's Github:
% https://github.com/zarzycki/ExTraTrack/tree/master/functions
% March 26, 2019: returns NaN if storm is not moving (otherwise it would
% return 0!)

% % Define constants
%   pi = 4.*atan(1.);
%   d2r = pi/180;
%   r2d = (1/d2r);
% 
% % convert lat/lon points from degrees to radians
%   lat1r = lat1*d2r;
%   lon1r = lon1*d2r;
%   lat2r = lat2*d2r;
%   lon2r = lon2*d2r;
%   
% % compute angle of motion based on two lat/lon points
% ang_tmp = r2d*atan2(sin((lon2r-lon1r))*cos(lat2r),... 
%                   cos(lat1r)*sin(lat2r)-sin(lat1r)*cos(lat2r)*cos(lat2r-lat1r)); % atan2(dy,dx)
% % convert angles from (-180 to +180) to (0 to 360)
% ang=wrapTo360(ang_tmp)  % gives angle in met convention
%                                       % 0= towars North; 
%                                       % 90 = moving East;
%                                       % 180 = moving South;
%                                       % 270= moving West;

% To test what the code is doing, using distance.m function instead:  
R=earthRadius('km');
Y=distance(lat1,lon1,lat2,lon1,[R 0]); if lat2<lat1; Y=-Y;end 
X=distance(lat1,lon1,lat1,lon2,[R 0]); if lon2<lon1; X=-X;end
ang_tmp = rad2deg(atan2(X,Y));
ang= wrapTo360(ang_tmp);

% if Storm is still, return NaN:
% (otherwise, it would return 0)
if lat1==lat2 && lon1==lon2
    ang=NaN;
end

% gives angle in met convention
                                      % 0= towars North; 
                                      % 90 = moving East;
                                      % 180 = moving South;
                                      % 270= moving West;
end