/
vecfield2d.R
100 lines (97 loc) · 3.65 KB
/
vecfield2d.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
#' Vector Fields.
#'
#' Draw a vector field associated to a system of at least two ODE.
#'
#' @param coords a matrix with two columns or more that is optionally used to
#' alternatively define the coordinates of the vector field.
#' @param FUN the function that describes the dynamical system (see details).
#' @param args the parameters of the dynamical system (see details).
#' @param ndim number of dimension of the system. If `NULL` the values is
#' based on `coords` and `slice`
#' @param slices a vector of 2 elements providing the dimensions to be
#' displayed, (default set to c(1,2)).
#' @param fixed the values used for non drawn dimension, if `NULL` the values
#' will be set to 0.
#' @param cex.x the magnification coefficient to be used for lengths of
#' vectors along the x axis.
#' @param cex.y the magnification coefficient to be used for lengths of
#' vectors along the y axis.
#' @param log a logical. If `TRUE`, the lengths of arrows are log-transformed.
#' @param add a logical. If `TRUE`, the vector field is added on the current plot.
#' @param ... additional arguments to be passed to [arrows()].
#'
#' @keywords empty plot
#'
#' @export
#'
#' @details
#' The `FUN` function must be a function of at least two arguments. The
#' first argument must contain the dynamical variables as a vector and the
#' second arguments must contain all the other parameters that shapes the
#' dynamical system. When some dimensions are missing, the order of `fixed`
#' is the one in `FUN` once the drawn dimension are withdrawn.
#'
#' @examples
#' systLin <- function(X, beta){
#' Y <- matrix(0,ncol=2)
#' Y[1L] <- beta[1,1]*X[1L]+beta[1,2]*X[2L]
#' Y[2L] <- beta[2,1]*X[1L]+beta[2,2]*X[2L]
#' return(Y)
#' }
#' seqx <- seq(-2,2,0.31)
#' seqy <- seq(-2,2,0.31)
#' beta1 <- matrix(c(0,-1,1,0),2)
#' # Plot 1:
#' vecfield2d(coords=expand.grid(seqx, seqy), FUN=systLin,
#' args=list(beta=beta1))
#' # Plot 2:
#' graphics::par(mar=c(2,2,2,2))
#' vecfield2d(coords=expand.grid(seqx, seqy), FUN=systLin,
#' args=list(beta=beta1), cex.x=0.35, cex.arr=0.25,
#' border=NA,cex.hh=1, cex.shr=0.6, col=8)
#' graphics::abline(v=0,h=0)
vecfield2d <- function(coords, FUN, args = NULL, ndim = NULL,
slices = c(1, 2), fixed = NULL, cex.x = 0.25, cex.y = cex.x,
log = FALSE, add = FALSE, ...) {
##-- Format checking
grid <- as.matrix(coords)
stopifnot(ncol(grid) > 1)
stopifnot(length(slices) == 2)
##-- Generating a grid based on coords / slices / fixed
if (is.null(ndim))
ndim <- max(ncol(grid), length(fixed) + length(slices), slices)
gridin <- matrix(0, ncol = ndim, nrow = nrow(grid))
#
if (ncol(grid) < ndim) {
gridin[, slices] <- grid
if (is.null(fixed))
fixed <- rep(0, ndim - 2)
grid_val <- matrix(rep(fixed, each = nrow(grid)), ncol = length(fixed))
gridin[, -slices] <- grid_val
} else gridin <- grid
##---
gridout <- gridin * 0
fun_names <- names(formals(FUN))
for (i in seq_len(nrow(gridin))) {
args[[fun_names[1L]]] <- gridin[i, ]
gridout[i, ] <- do.call(FUN, args)
}
if (!add)
plot0(range(gridin[, slices[1L]]), range(gridin[, slices[2L]]))
## ----
for (i in seq_len(nrow(gridout))) {
dstx <- cex.x * gridout[i, 1L]
dsty <- cex.y * gridout[i, 2L]
if (log) {
dstx <- log(abs(dstx) + 1)
if (gridout[i, 1L] < 0)
dstx <- -dstx
dsty <- log(abs(dsty) + 1)
if (gridout[i, 2L] < 0)
dsty <- -dsty
}
arrows2(gridin[i, 1L], gridin[i, 2L], gridin[i, 1L] + dstx,
gridin[i, 2L] + dsty, ...)
}
invisible(NULL)
}