いちからはじめる GMT その35 - psxyz

2010年 12月 1日 掲載

GMT には，3次元の描画をする psxyz というコマンドがある． 今回は，これを使ってみよう．

まずはサンプルデータ作りから． Lorenz (1963) を参考にして，Lorenz Attractor の軌跡データを作る． 以下のような awkスクリプト lorenz.awk を書いてみた．

BEGIN{

  sigma = 10
  b = 8 / 3
  gamma = 28

  dt = 0.01

  x = 0
  y = 1
  z = 0

  for(i = 0; i <= 6000; i++){
    print x, y, z

    x_n1 = x + fx(x, y, z) * dt
    y_n1 = y + fy(x, y, z) * dt
    z_n1 = z + fz(x, y, z) * dt

    x_n2 = x_n1 + fx(x_n1, y_n1, z_n1) * dt
    y_n2 = y_n1 + fy(x_n1, y_n1, z_n1) * dt
    z_n2 = z_n1 + fz(x_n1, y_n1, z_n1) * dt

    x = (x + x_n2) / 2
    y = (y + y_n2) / 2
    z = (z + z_n2) / 2
  }

}

function fx(x, y, z){
  return - sigma * x + sigma * y
}

function fy(x, y, z){
  return - x * z + gamma * x - y
}

function fz(x, y, z){
  return x * y - b * z
}

これを実行すると，6000ステップの時間積分データが出力される．

$ awk -f lorenz.awk > lorenz.txt

まずはこれを，psxyコマンドで，ふつうに 2次元の図に描いてみよう．

$ psxy lorenz.txt -Jx0.2 -R-20/20/-30/30 -Ba5 > lorenz.eps

さあ，次は 3次元の描画だ． とりあえず，psxyコマンドを psxyzコマンドで置き換えてみようか．

$ psxyz lorenz.txt -Jx0.2 -R-20/20/-30/30 -Ba5 > lorenz.eps

psxyコマンドと何も変わらない．

psxyzコマンドで 3次元の描画をするには，さらに

1. Z座標方向の縮尺を指定する -Jzオプション
2. Z座標方向の描画範囲（これは -Rオプションに追加）
3. 視点の方向を指定する -Eオプション

の 3つの指定を追加する必要がある．

では実際にためしてみよう． Z座標方向の縮尺は，-Jx と同じで -Jz0.2，描画範囲はデータを見て，0/50 とする． 視点の方向はいろいろためして，-E150/30（方位 150度，仰角 30度）とした． この視点の方向の指定については，次回くわしく説明する．

$ psxyz lorenz.txt -Jx0.2 -Jz0.2 -R-20/20/-30/30/0/50  -Ba5 -E150/30 > lorenz.eps

おお，かっこいい．

参考文献：

• Lorenz, E. N., 1963: Deterministic Nonperiodic Flow. J. atmos. Sci., 20, 130-141

[ KDOホームにもどる ] [ いちからはじめる GMT インデックスにもどる ]

(C) 2010 KDO