Как нарисовать заказное векторное поле в базе R?

Question

Я хочу сгенерировать график, который бы отображал «векторное поле» на графике, используя базу R.

В этом разделе скрипта будет сгенерирован график с разными кругами (изменяя только их радиус, чтобы сделать база векторного поля).

r = 100 # Set the maximum radius to make an empty plot
x = function(ang) r* cos(ang) # To draw circles
y = function(ang) r* sin(ang) # To draw circles
nb = seq(from = 0,to = (2*pi),length.out = 100) # To have a sequence to draw the circles
plot(x(nb),y(nb), # Empty plot 
     asp = 1, 
     type = "n", 
     bg = "black", 
     col = "black", 
     pch =21, main = "Circle",
     ylab = "Y values",
     xlab ="X values")
abline(h=0,v=0) # Draw axes
for (i in seq(0,100,by = 5)) { # Draw a series of circles
  r = i
  points(x(nb),y(nb), 
         type = "l", 
         lwd = 1.0,
         lty = 3)
}

# DRAWING TE VECTORS ----------------------
by = 10 # Define a "resolution" to see better the circles (This value will be smaller to be more precise)

changex = seq(0,100, by =by) # For each circle draw a radius with this sequence 

current = -1 # This is to "flip" the orientation of the vectors
mag = current* seq(100,0, by = current*by)
arrows(x0 = changex, y0 = 0, # Draw the vectors 
       x1 = changex, y1 = mag,
       code = 2,
       length=0.1,
       angle=40)

Остальная часть кода пытается печатать векторы при изменении угла на графике:

xycircle <- function(ang,r) { # function to draw position on the circle
  x = r*cos(ang)
  y = r*sin(ang)
  return(list(x,y))
}

pilist = c(#0,1/4*pi,#1/2*pi, # List of PI values to go around the circle
           #pi, #3/4*pi,
           #3/2*pi,
           2*pi)
for (pip in 1:length(pilist)) { # Going around the circle

ang = pilist[pip] # extract 1 angle value to draw 

abline(a=0,b=tan(ang), lty = 3, lwd = 3) # Get a line that will show the angle selected 
r = seq(0,100, by = by) # List of radius
mag = current* seq(-100,-0, by = by) # Magnitude of the vectors 

for (i in 1:length(r)) {  # Draw vectors when the angle changes  
  arrows(x0 = xycircle(ang,r[i])[[1]], # Base position of the vector (tangent to the circle)
         y0 = xycircle(ang,r[i])[[2]],
         x1 = cos(atan2(r[i],mag[i])-ang)*sqrt(r[i]^2+mag[i]^2), # Position of the tip of the vector (x)
         y1 = sin(atan2(r[i],mag[i])-ang)*sqrt(r[i]^2+mag[i]^2), # Position of the tip of the vector (y)
         code = 2, # Change the arrow head 
         length = 0.1,
         angle = 40)
}

}

Как видите, когда я совершаю полный оборот, векторы не совпадают с исходными векторами (они должны ...).

Но когда я начинаю поворачиваться по кругу (скажем, 1/4 * пи), это прекрасно.

Как можно было бы заставить векторы «поворачиваться» вокруг круга (в зависимости от углов), чтобы они вращались вокруг круга так, чтобы векторы всегда были перпендикулярны круги (как на последнем графике, но все углы).

Answer 1

Наконец-то я получил его на работу

# parameter list ----------------------------------------------------------
by = 10
current = -1
invert.speed = TRUE
circleefrom = 10
circleeto = 28
# See 1 circle from 3 (2 for 4, 3 for 5...)
resolution = 15 # Which is also the number of arrows
resolution.circles = 100
maxcurrent = 10
mincurrent = 0


# Define functions --------------------------------------------------------

# Make a list of radius to draw various circles 
r = seq(from = circleefrom, to = circleeto,
        length.out = resolution)

# Make a function that will allow to draw circles and extract the values from the x,y position
xycircle <- function(ang,r) {
  x = r*cos(ang)
  y = r*sin(ang)
  return(list(x,y))
}

# Samples enough point for the circles (100 at least)
nb = seq(from = 0, # starts at 0 
         to = (2*pi), # Does a FULL revolution around the circle 
         length.out = resolution.circles)


# Plot circles ------------------------------------------------------------
# Create empty canva 
plot(x = xycircle(nb,max(r))[[1]],
     y = xycircle(nb,max(r))[[2]],
     asp = 1, 
     type = "n", 
     bg = "black", 
     col = "black", 
     pch =21, main = "Rheotaxis experiment",
     ylab = "Y values",
     xlab ="X values")

# Draw the axes 
abline(h = 0,
       v = 0)

# This will draw the circles 
dbcicle = NULL
for (i in seq(from = circleefrom, to = circleeto,
              length.out = resolution)) {
  points(x = xycircle(nb,i)[[1]],
         y = xycircle(nb,i)[[2]],
         col = "blue",
         type = "l", 
         lwd = 1.0,
         lty = 3)
tmp = as.data.frame.list(x = c(xycircle(nb,i),i), col.names = c("x","y","radius"))
dbcicle = rbind(dbcicle,tmp)
}

# Add a black center to the design 
plotrix::draw.circle(0, 0, radius =  circleefrom, 
                     nv = 1000, 
                     border = NULL,
                     col = "black", lty = 1, lwd = 1)
# Add circle at the contour 
plotrix::draw.circle(0, 0, radius =  circleeto, 
                     nv = 1000, 
                     border = NULL,
                     col = NA, lty = 1, lwd = 1)


# Current specification ---------------------------------------------------
# Create the simulated current 
if(invert.speed) {
  mag = current * seq(from = mincurrent,
                      to = maxcurrent, 
                      length.out = resolution + 1)

} else {
mag = current * seq(from = maxcurrent,
                    to = mincurrent, 
                    length.out = resolution + 1)
}
# mag = current * rep(16,11)


# Draw the vecotrs of current  --------------------------------------------


# Get different angle values 
pilist = seq(from = 1/2*pi,
             to = 2*pi,
             by = 1/2*pi)
pilist = seq(from = 0,
             to = 2*pi,
             length.out = resolution)

dbcicle2 = NULL
for (pip in 1:length(pilist)) {
  ang = pilist[pip]
  abline(a=0,b=tan(ang), lty = 3, lwd = 3)
  for (i in 1:length(r)) {
    arrows(x0 = xycircle(ang,r[i])[[1]], 
           y0 = xycircle(ang,r[i])[[2]],
           x1 = xycircle(c(ang-atan2(mag[i],r[i])),sqrt(r[i]^2+mag[i]^2))[[1]], 
           y1 = xycircle(c(ang-atan2(mag[i],r[i])),sqrt(r[i]^2+mag[i]^2))[[2]],
           col = "blue",
           code = 2, 
           length = 0.1,
           angle = 40)
  tmp = as.data.frame.list(x = c(xycircle(ang,r[i]),r[i],-1*mag[i],ang+pi/2), 
                           col.names = c("x","y","radius","magnitude","ang"))
  dbcicle2 = rbind(dbcicle2,tmp)
  }
}