Introduction

# Haversine formula (https://en.wikipedia.org/wiki/Haversine_formula): 

my_dist <- function(long1, lat1, long2, lat2) {
  rad <- pi/180
  a1 <- lat1*rad
  a2 <- long1*rad
  b1 <- lat2*rad
  b2 <- long2*rad
  dlon <- b2 - a2
  dlat <- b1 - a1
  a <- (sin(dlat/2))^2 + cos(a1)*cos(b1)*(sin(dlon/2))^2
  c <- 2*atan2(sqrt(a), sqrt(1 - a))
  R <- 6378137
  d <- R*c
  return(d)
}

# Test hàm (https://andrew.hedges.name/experiments/haversine/): 
my_dist(53.36575, 7.348687, 53.36507, 7.348940)

## [1] 80.18433

# Use geosphere package: 

library(geosphere)
distHaversine(c(53.36575, 7.348687), c(53.36507, 7.348940))

## [1] 80.18433