[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[dennou-ruby:001827] Re: GPhys::EP_Flux



塚原様:

「流れ」ベクトル用のメソッドを作りました。ggraph.rb の DCLExtモ
ジュールに入れてあります(CVSのみ)。ugvect の代りに使うものですが、
前に話した通り、factor で全体の長さを微調整できます。

また、xintv, yintv で指定する間引き機能をつけました。
ugvect を呼んでますので、矩形座標専用です(ugvect)。
さらに言えば、itr=1 でないと、ちゃんと流れになりません。
将来的には ugvect を使わないで地図投影に対応するものも欲しいです
ね。

factor=1.0は x, y 成分の長さの最大が平均グリッド間隔になるよう
にという ugvect 的仕様です。流れを表すため、ベクトルの縦横比を、
(U 座標の縦横比/V座標の縦横比)で調整します。速度場やフラック
スの表示に好適。

    def ug_flow_vect( fx, fy, factor=1.0, xintv=1, yintv=1)
      raise ArgumentError,"Expect 2D arrays" if fx.rank != 2 || fy.rank != 2
      raise ArgumentError,"fx.shape != fy.shape" if fx.shape != fy.shape
      raise ArgumentError,"xintv must be a positive integer" if xintv < 0
      raise ArgumentError,"yintv must be a positive integer" if yintv < 0
      nx, ny = fx.shape
      if xintv >= 2
	idx = NArray.int(nx/xintv).indgen!*xintv  # [0,xintv,2*xintv,..]
	fx = fx[idx, true]
	fy = fy[idx, true]
      end
      if yintv >= 2
	idx = NArray.int(ny/yintv).indgen!*yintv  # [0,yintv,2*yintv,..]
	fx = fx[true, idx]
	fy = fy[true, idx]
      end
      nx, ny = fx.shape  # again, because of xintv & yintv 
      vx0,vx1,vy0,vy1 = DCL.sgqvpt
      ux0,ux1,uy0,uy1 = DCL.sgqwnd
      dvx = (vx1-vx0)/nx
      dvy = (vy1-vy0)/ny
      ax = (vx1-vx0)/(ux1-ux0)   # factor to convert from U to V coordinate
      ay = (vy1-vy0)/(uy1-uy0)   # factor to convert from U to V coordinate
      fxmx = fx.abs.max
      fymx = fy.abs.max
      raise "fx has no data or all zero" if fxmx == 0
      raise "fy has no data or all zero" if fymx == 0
      cn = [ dvx/(ax*fxmx),  dvy/(ay*fymx) ].min  # normarization constant
      lnrmal = DCL.ugpget('LNRMAL')  # save
      DCL.ugpset('LNRMAL', false)
      DCL.ugvect( factor*cn*ax*fx, factor*cn*ay*fy )
      DCL.ugpset('LNRMAL', lnrmal)    # recover
    end

以下、テストプログラムです。U座標の縦横比にかかわらず、流れが
追えるのがわかると思います。

---
require "numru/ggraph"
include NumRu

nx = 20
ny = 10
ux = 10.0
uy = 1.0
fx = ux*NArray.sfloat(nx,ny).indgen!
fy = uy*NArray.sfloat(nx,ny).indgen!

DCL.gropn(1)

DCL.grfrm
DCL.grswnd(0,ux,0,uy)
DCL.grsvpt(0.2,0.8,0.2,0.8)
DCL.grstrn(1)
DCL.grstrf
DCL.usdaxs
DCLExt.ug_flow_vect(fx, fy, 2.0)

DCL.grfrm
DCL.grswnd(0,ux,0,uy)
DCL.grsvpt(0.2,0.8,0.2,0.8)
DCL.grstrn(1)
DCL.grstrf
DCL.usdaxs
DCLExt.ug_flow_vect(fx, fy, 1.0, 2, 1)

DCL.grcls
---

堀之内