# distanza fra due punti definiti dai vettori posizione
dist <- function(p1, p2) {
   dp = p1 - p2
   sqrt( as.numeric( t(dp) %*% dp ) )
}

# forza SUL punto materiale di massa m1 situato in p1 
# dovuta al punto materiale di massa m2 situato in p2

forza <- function(m1, p1, m2, p2) {
   G = 6.67e-11 # N m^2 / kg^2
   
   return( -G*m1*m2 / dist(p1,p2)^3 * (p1-p2) )
}