code: public static final double radBearing(double lat1, double lon1, double lat2, double lon2)
      {
          checkLatitude(lat1);
          checkLongitude(lon1);
          checkLatitude(lat2);
          checkLongitude(lon2);
          double dep = departure(lat1, lat2);
          double aa = dep*(lon2-lon1);
          double bb = lat2-lat1;
          double dd = Math.atan2(aa, bb);
          if (dd < 0)
          {
              dd += 2*Math.PI;
          }
          return dd;
      }
docs: Returns bearing from (lat1, lon1) to (lat2, lon2) in radians
  @param lat1
  @param lon1
  @param lat2
  @param lon2
  @return