function [ret] = bearing( LAT1, LONG1, LAT2, LONG2 )
  LAT1r=LAT1*pi/180;
  LONG1r=LONG1*pi/180;
  LAT2r=LAT2*pi/180;
  LONG2r=LONG2*pi/180;

  x=atan2( sin(LONG2r-LONG1r)*cos(LAT2r), cos(LAT1r)*sin(LAT2r)-sin(LAT1r)*cos(LAT2r)*cos(LONG1r-LONG2r))*180/pi;
  if x<0
    x+=360;
  endif
  ret=x;
endfunction
