code: public static double cosineFormulaRad(double lat1, double lon1, double lat2, double lon2) {
      if(lat1 == lat2 && lon1 == lon2) {
        return 0.;
      }
      final DoubleWrapper tmp = new DoubleWrapper(); // To return cosine
      final double slat1 = sinAndCos(lat1, tmp), clat1 = tmp.value;
      final double slat2 = sinAndCos(lat2, tmp), clat2 = tmp.value;
      final double a = slat1 * slat2 + clat1 * clat2 * cos(lon2 - lon1);
      return a < .9999_9999_9999_999 ? acos(a) : 0;
    }
docs: Compute the approximate great-circle distance of two points using the
  Spherical law of cosines.
  <p>
  Complexity: 6 trigonometric functions. Note that acos
  is rather expensive apparently - roughly atan + sqrt.
  <p>
  Reference:
  <p>
  R. W. Sinnott,<br>
  Virtues of the Haversine<br>
  Sky and Telescope 68(2)
  
  @param lat1 Latitude of first point in degree
  @param lon1 Longitude of first point in degree
  @param lat2 Latitude of second point in degree
  @param lon2 Longitude of second point in degree
  @return Distance on unit sphere