Drawing quiver plots with R

Our extreme solvers site solvers/rplotpage/quiver
allows one to draw a quiver plot or arrows plot over a displayed image of a function in two variables of the form $$fxy= f(x,y)$$. It is based on R routines by Ripley and Hand and illustrated in the addictedtor.org gallery from which our qpy code is based.

Here is an example for the function (or expression) $$fxy= \sin(x) + \cos(y)$$. Here we have the settings xlo=ylo = -2.0, xhi=yhi=2.0 and xby=yby = 0.2. The contour color is gray and the pallette for the image is terrain.colors. Here is the output image.

Note that we have shown a test code to generate quiver plots previously.

Generating quiver plots using R under Python

Here is Python code to generate an R script which in turn does the actual image generation.

"""
quiver-test.py
"""

def quiverfunc(fxy, x0, x1, xby, y0, y1, yby, pallette="terrain.colors", contourcolor="gray"):
   S = """
par.uin <- function() 
  # determine scale of inches/userunits in x and y
  # from http://tolstoy.newcastle.edu.au/R/help/01c/2714.html
  # Brian Ripley Tue 20 Nov 2001 - 20:13:52 EST
 {
    u <- par("usr") 
    p <- par("pin") 
    c(p[1]/(u[2] - u[1]), p[2]/(u[4] - u[3]))
 }

quiver2 <- function(expr,
                     x,
                     y,
                     nlevels=20, 
                     length=0.05, 
                     ...){

    z <- expand.grid(x,y) 
    xx  <- x
    x   <- z[,1]
    yy  <- y
    y   <- z[,2]

    fxy <- eval(expr) 
    grad_x <- eval(D(expr, "x")) 
    grad_y <- eval(D(expr, "y")) 

  dim(fxy) <- c(length(xx), length(yy)) 
  dim(grad_x) <- dim(fxy) 
  dim(grad_y) <- dim(fxy) 

  maxlen <- min(diff(xx), diff(yy)) * .9 
  grad_x <- grad_x / max(grad_x) * maxlen 
  grad_y <- grad_y / max(grad_y) * maxlen 

  filled.contour(xx, yy, fxy, nlevels=nlevels, 
    plot.axes = { 
      contour(xx, yy, fxy, add=T, col="gray", 
              nlevels=nlevels, drawlabels=FALSE) 

      arrows(x0  = x, 
             x1  = x + grad_x,
             y0  = y,
             y1  = y + grad_y,
             length = length*min(par.uin())) 

      axis(1) 
      axis(2) 
    },
    ...)
}

f <- expression( #expr) 
x <- seq(#x0, #x1,by= #xby) 
y <- seq(#y0, #y1,by= #yby) 
par(mar=c(3,3,3,3)) 
quiver2(f,x,y, color.palette=#pallette) 
graphics.off()
"""
   S = S.replace("#expr", fxy)
   S = S.replace("#x0", x0)
   S = S.replace("#x1", x1)
   S = S.replace("#xby",xby)
   S = S.replace("#y0", y0)
   S = S.replace("#y1", y1)
   S = S.replace("#yby",yby)
   S = S.replace("#pallette",pallette)
   S = S.replace("#contourcol", contourcolor)
   return S 

Rcode = 'png("temp.png", 6 *72, 6 *72)'
Rcode +=  quiverfunc("(3*x^2 + y) * exp(-x^2-y^2)", "-2", "2", "0.2", "-2", "2", "0.2", pallette="terrain.colors", contourcolor="gray")
print Rcode

Save the file to "test-quiver2.py". Then to run it, type python quiver2.py > out.R. Then issue R < out.R --no-save. The generated image is in temp.png. Finally, display it using the imagemagick tool: display temp.png.