[dennou-ruby:002611] [Ruby_VTK] drawing stream line & contour surface using Ruby VTK

[Youhei SASAKI <uwabami@xxxxxxxxxxxxxxxxxxxx>]

西澤様, 竹広様
Cc: dennou-ruby

佐々木です.

Ruby_VTK + easy_vtk を用いて流線と等値面を描いてみました.

# 西澤さん, ありがとうございます.
# スゴイ, 素晴しいです. 思わず興奮して鼻血出しそうでした.

とりあえず, できた絵とスクリプト添付します.

絵は「渦度の南北成分の等値面 + 流線(赤道面からスタート)」です. 等値面
の色使いは, とりあえずツッコマないで下さい

# 緑にして, 透明度上げたら増えるワカメみたいと石渡さんに言われた, とか.
# 赤だと内臓みたいだ, と林さんに言われたとか. むー.

西澤様, easy_vtk について幾つか質問したいのですが,

 1. 流線に色をつけたり, 流線の太さを変えたりする事は可能ですか?

 2. 等値面を描く際に, 特定の領域だけ描かせる事は可能ですか?
    
    例えば, 流線はきちんと描いてみたけど, 境界層がジャマなんで
    渦度の等値面は境界層抜いて描いてみたい, とか.

 3. 流線描く際のイテレータについて,
    まだ streamline の使い方がよくわかって無いのかもしれませんけど.

    Rdoc かけた easy_vtk によれば

     streamline( point, time, ds, options={})

    なので, point で指定されている座標からスルスルっと流線を描いて
    るのですよね. この point の指定について繰り返し処理させる, と.
    points について, なんらかの制限はありますか?

    例えば, 計算分解能が低く, 指定した point が存在しない場合は, どう
    なるのでしょう. 流線が描けないだけでしょうか? 

---
佐々木 洋平 (Sasaki Youhei)
北海道大学 大学院理学研究科 地球惑星科学専攻 
地球惑星流体科学講座 地球流体力学研究室 
e-mailto: Youhei SASAKI<uwabami@xxxxxxxxxxxxxxxxxxxx> 

=begin
$Lastupdate: 2006/05/10 18:41:10 $ Youhei SASAKI
流線(磁力線)お絵描きのテストクリプト
=end

# library included
require 'easy_vtk'
require 'numru/gphys'
include NumRu
include NMath

# File I/O
fname = 'case1-T21N16-11.nc'
x_vname = 'vlon'
y_vname = 'vlat'
z_vname = 'vrad'
x_gphys = GPhys::IO.open(fname,x_vname)
x_gphys = x_gphys.cut('t'=>11)
x_gphys = x_gphys.cyclic_ext(0, 360)
y_gphys = GPhys::IO.open(fname,y_vname)
y_gphys = y_gphys.cut('t'=>11)
y_gphys = y_gphys.cyclic_ext(0, 360)
z_gphys = GPhys::IO.open(fname,z_vname)
z_gphys = z_gphys.cut('t'=>11)
z_gphys = z_gphys.cyclic_ext(0, 360)

## 座標の設定
# ラジアンに変換
lon = x_gphys.coord(0).val*PI/180
lat = y_gphys.coord(1).val*PI/180
rad = z_gphys.coord(2).val

# 球→カルテジャン変換
nlon,nlat,nrad = x_gphys.shape
lon.reshape!( nlon, 1, 1 )
lat.reshape!( 1, nlat, 1 )
rad.reshape!( 1, 1, nrad )
x = rad*cos(lat)*cos(lon)
y = rad*cos(lat)*sin(lon)
z = NArray.sfloat(nlon,1,1).fill(1)*rad*sin(lat)

x_gphys = x_gphys.val
y_gphys = y_gphys.val
z_gphys = z_gphys.val
x_vec = - x_gphys*cos(lon) + ( y_gphys*sin(lat) + z_gphys*cos(lat) )*cos(lon)
y_vec = x_gphys*sin(lon) + ( y_gphys*sin(lat) + z_gphys*cos(lat) )*sin(lon)
z_vec = y_gphys*cos(lat) + z_gphys*sin(lat)

x.reshape!( nlon*nlat*nrad )
y.reshape!( nlon*nlat*nrad )
z.reshape!( nlon*nlat*nrad )
vec = NArray.sfloat( nlon*nlat*nrad, 3 )
vec[true,0] = x_vec.reshape!( nlon*nlat*nrad )
vec[true,1] = y_vec.reshape!( nlon*nlat*nrad )
vec[true,2] = z_vec.reshape!( nlon*nlat*nrad )

# 描画開始
EasyVtk::init
EasyVtk::set_size( 800, 600 ) # ウィンドウサイズ
EasyVtk::set_background( 1, 1, 1 ) # 背景色
# 格子点 および データ の設定
EasyVtk::set_grid( x, y, z, [nlon,nlat,nrad] )
EasyVtk::set_data( vec )

# 流線
npoints = 10
mpoints = 10
npoints.times{|i|
  r = rad[i*nrad/npoints]
  mpoints.times{|j|
    lamda = j*(PI/10)
    phi  = 0
    px = r*cos(lamda)*cos(phi)
    py = r*sin(lamda)*cos(phi)
    pz = r*sin(phi)
    EasyVtk::streamline( [px,py,pz], 10, 0.01 )
  }
}

## 等値面描画
# 渦度の等値面も描いてみる.
s_vname = 'zlat'
s_gphys = GPhys::IO.open(fname,s_vname)
s_gphys = s_gphys.cut('t'=>11)
s_gphys = s_gphys.cyclic_ext(0,360)
s_gphys = s_gphys.val.reshape!(nlon*nlat*nrad).to_va
EasyVtk::set_data(s_gphys)
EasyVtk::surface(-30, 'opacity'=>0.3, 'color'=>[0,0,0.9])

## 球面
# 外核
#EasyVtk::sphere( rad[0], 'resolution'=>21 , 'opacity'=>0.1 )
# 内核
EasyVtk::sphere( rad[-1], 'resolution'=>21, 'color'=>[0.5,0.5,0.5], 'opacity'=>0.9 )

# GUIスタート
EasyVtk::gui_start
#EasyVtk::save_file("velocity.png")