dms_to_rad <- function(d, m, s) (d + m / 60 + s / 3600) * pi / 180




great_circle_distance <- function(lat1, long1, lat2, long2) {
   a <- sin(0.5 * (lat2 - lat1))
   b <- sin(0.5 * (long2 - long1))
   12742 * asin(sqrt(a * a + cos(lat1) * cos(lat2) * b * b))
}





great_circle_distance(
   dms_to_rad(36,  7, 28.10), dms_to_rad( 86, 40, 41.50),   
   dms_to_rad(33, 56, 32.98), dms_to_rad(118, 24, 29.05))