ベクトル場を矢印で可視化したい

昔書いたコードを掘り起こした。どこを参考に書いたのかは忘れた。以下は二次元ベクトル場計算の例(調和振動子のポテンシャァル想定)。

gradient.harmonic.potential <- function(coordinate, x0, y0, k)
{
  #Evaluate : -grad(Potential)
  x <- coordinate[1]
  y <- coordinate[2]
  -c(vx=sum(k*(x-x0)), vy=sum(k*(y-y0)))
}

上のベクトル場を２つの調和振動子のポテンシャルに対して適当なグリッドポイントで計算＆可視化。
赤丸の箇所(X0, Y0)が二つの調和振動子の各々の均衡点(力が0になる点)。一方、実際にplotされているベクトル場はその合力。

#Setting for potential position
K <-  c(1, 1)*(0.1)
X0 <- c(25, 75)
Y0 <- c(25, 50)
#grid setting for plot
grid <- expand.grid(x=seq(1, 100, 10), y=seq(1, 100, 10))
#evaluate vector filed of potential and plot
vector.filed <- cbind(grid, t(apply(grid, 1, function(points){gradient.harmonic.potential(points, X0, Y0, K)})))
p <- ggplot(vector.filed, aes(x=x, y=y))  
print(p + geom_segment(aes(xend=x+vx, yend=y+vy), arrow = arrow(length = unit(0.5,"cm"))) + 
  geom_point(data=data.frame(x=X0,y=Y0), colour="red", size = 8.0) 
)