drawComplex.R

# drawComplex.R

# gridパッケージが必要
if (!require('grid')) {
  install.packages('grid')
  stopifnot(require('grid'))
}

# noDefaultとmakeFunctionが必要
stopifnot(exists("noDefault") && is.function(noDefault))
stopifnot(exists("makeFunction") && is.function(makeFunction))

# 複素関数の値（wの値）を成分とする複素数マトリックスを作成し、関連情報と共に返す
makeXyzw <-
  function (
    # 複素関数の式、コールオブジェクトを渡す
    expr,
    # 複素変数の名前、文字列で指定する
    varname = "z",
    # 複素変数のガウス平面（x-y-座標）における描画対象域
    xlim = c(-1, 1), ylim = c(-1, 1),
    # 描画対象域のx,y方向を幾つの少区間に分割するか（分割数）
    xdiv = 20, ydiv = 20)
{
  x0 <- xlim[1]
  x1 <- xlim[2]
  y0 <- ylim[1]
  y1 <- ylim[2]
  xStep <- (x1 - x0)/xdiv
  yStep <- (y1 - y0)/ydiv

  xBreaks <- seq(from = x0, to = x1, by = xStep)
  yBreaks <- seq(from = y0, to = y1, by = yStep)

  # 関数の仮引数リストの作成
  fplist <- list(noDefault()) # デフォルト値なし1変数
  names(fplist) <- c(varname)  # 引数名（複素独立変数名）を指定
  # 関数の作成
  fun <- makeFunction(fplist, expr) 

  wMatrix <- matrix(NA, xdiv + 1, ydiv + 1)
  for (i in 1:(xdiv + 1)) {
    for (j in 1:(ydiv + 1)) {
      wMatrix[i, j] <- fun(xBreaks[i] + yBreaks[j]*1i)
    }
  }
  # 参考のために独立変数zのマトリックスも添える
  zMatrix <- outer(xBreaks, yBreaks, function(x, y) x + y*1i)
  list(
      xlim = xlim,
      ylim = ylim,
      xStep = xStep,
      yStep = yStep,
      x = xBreaks,
      y = yBreaks,
      z = zMatrix,
      w = wMatrix
      )
}  

# 複素関数を描く
# 描画方法："map", "re", "im", "mod", "arg"
# map : 写像図、x方向／y方向の直線群に対する像曲線群を描く
# re : 値の実部を描く
# im : 値の虚部を描く
# mod : 値の絶対値を描く
# arg : 値の偏角を描く
drawComplex <-
  function (
    # 複素関数の式、引数としてzの式を書く
    expr,
    # 描画方法、シンボルでもよい
    way = "map",
    # 複素変数のガウス平面（x-y-座標）における描画対象域、原点中心正方形の辺長の半分
    half = 1,
    # 複素変数のガウス平面（x-y-座標）における描画対象域、x-y-共通
    lim = c(-1*half, half),
    # 複素変数のガウス平面（x-y-座標）における描画対象域、x-y-個別
    xlim = lim, ylim = lim,
    # 描画対象域のx,y方向を幾つの少区間に分割するか（分割数）、x-y-共通
    div = 20,
    # 描画対象域のx,y方向を幾つの少区間に分割するか（分割数）、x-y-個別
    xdiv = div, ydiv = div,
    # perspのオプション
    theta = 30, phi = 30, expand = 0.5, col = "skyblue",
    # 関数の情報を戻り値として出力するか
    output.info = FALSE
    )
{
  # 引数の式を評価せずに取り出す
  expr <- substitute(expr)
  # way引数の取得とチェック
  way <- substitute(way)
  stopifnot(is.name(way) || is.character(way))
  way <- as.character(way)
  stopifnot(any(c("map", "re", "im", "mod", "arg") == way))
  
  # 関数値のマトリックスを作成する
  xyzw  <- makeXyzw(expr,
                    varname = "z",
                    xlim = xlim, ylim = ylim,
                    xdiv = xdiv, ydiv = ydiv
                    )
  # 基本的なパラメータをセット
  w <- xyzw$w
  ulim <- range(Re(xyzw$w))
  vlim <- range(Im(xyzw$w))
  rangeWidth <- ulim[2] - ulim[1]
  rangeHeight <- vlim[2] - vlim[1]
  rangeSize <- max(rangeWidth, rangeHeight)
  u0 <- ulim[1]
  v0 <- vlim[1]
  
  if (way == "map") {
    # 写像図による描画
    # 使用する変数: w, u0, v0, rangeSize, rangeWidth, rangeHeight
    uOffset <- (rangeSize - rangeWidth)/(2*rangeSize)
    vOffset <- (rangeSize - rangeHeight)/(2*rangeSize)
    
    grid.newpage()
    grid.rect(width = rangeWidth/rangeSize, height = rangeHeight/rangeSize)
    # x方向直線の像曲線
    for (j in 1:(ydiv + 1)) {
      u <- Re(w[ ,j])
      v <- Im(w[ ,j])
      grid.lines((u - u0)/rangeSize + uOffset, (v - v0)/rangeSize + vOffset , gp = gpar(col="red"))
    }
    # y方向直線の像曲線
    for (i in 1:(xdiv + 1)) {
      u <- Re(w[i, ])
      v <- Im(w[i, ])
      grid.lines((u - u0)/rangeSize + uOffset, (v - v0)/rangeSize + vOffset, gp = gpar(col="blue"))
    }
  } else {
    # perspによる3Dグラフ描画
    # 使用する変数: way, x, y, w, xlim, ylim
    funName <- paste(toupper(substr(way, 1, 1)), substr(way, 2, 10), sep="")
    fun <- eval(as.name(funName))

    # perspの引数リストを作成する
    args <- list(
      x = xyzw$x,
      y = xyzw$y,
      z = fun(xyzw$w),
      xlab = "x",
      ylab = "y",
      zlab = paste(funName, "(", deparse(expr), ")", sep=""),
      xlim = xlim,
      ylim = ylim,
      theta = theta,
      phi = phi,
      expand = expand,
      col = col
      )
    # persp関数の呼び出し
    do.call("persp", args)
  }
  # 情報の出力（オプショナル）
  if (output.info) {
    return(
      list (
        xlim = xlim,
        ylim = ylim,
        ulim = ulim,
        vlim = vlim,
        rangeRect = list(
          left = u0,
          bottom = v0,
          width = rangeWidth, 
          height = rangeHeight
          )
        )
      )
  }
}