to.radians<-function(degrees){
  degrees * pi / 180
}


#-----------------------------------------------------------------------------
# Vectorised havershine (lat1, long1) is a single point
# lats and longs vectors of other places
# returns a vector of distances from (lat1, long1) to each of the other places specified
#-----------------------------------------------------------------------------
vhaversine <- function(lat1, long1, lats, longs) {
  radius <- 6378      # radius of Earth in kilometers
  term1 <- sin((lats  - lat1)/2) ^ 2
  term2 <- cos(lat1) * cos(lats) * sin((longs - long1)/2) ^ 2
  the.terms <- term1 + term2
  delta.sigma <- 2 * atan2(sqrt(the.terms), sqrt(1-the.terms))
  distance <- radius * delta.sigma
  return(distance)
}

#-----------------------------------------------------------------------------
# My version calling vectorised havershine
#-----------------------------------------------------------------------------
just.R.mike.vectorised <- function(dframe){
  numrows <- nrow(dframe)

  lats  <- to.radians(dframe$x)
  longs <- to.radians(dframe$y)

  acc <- 0
  for (ii in 1:(numrows-1)) {
      acc <- acc + sum(vhaversine(lats[ii], longs[ii], lats[(ii+1):numrows], longs[(ii+1):numrows]))
  }

  acc / (numrows*(numrows-1)/2)
}


#-----------------------------------------------------------------------------
# Create data and do timings of the two methods
#-----------------------------------------------------------------------------
df <- read.csv("optim/Airport_Codes_mapped_to_Latitude_Longitude_in_the_United_States.csv", as.is=TRUE)
colnames(df) <- c('locationID', 'x', 'y')


system.time({print(just.R.mike.vectorised(df)) })

#> [1] 1869.744
#>   user  system elapsed
#>  10.021   3.000  13.036