[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[dennou-ruby:001838] Re: GPhys::EP_Flux
塚原様
堀之内です。
> > 「流れ」ベクトル用のメソッドを作りました。ggraph.rb の DCLExtモ
> > ジュールに入れてあります(CVSのみ)。ugvect の代りに使うものですが、
> > 前に話した通り、factor で全体の長さを微調整できます。
>
> 早速作ってくださったのですね. ありがとうございます.
お送りしたやつはどうもユニットベクトルがちゃんとかけないですね。
そもそも ugvect のユニットベクトル表示機能に頼っている限りは、ま
ず無理。ということで、できるだけ UG パックから独立させる方向に動
いてます。とりあえず、現状は以下のようになってます。名前の接頭辞
ug_ はもはやはずしました。ugpackのパラメターを併用しないんで。
-----------------
module NumRu
module DCLExt # 追加分のみ
# < flow vector package >
def __truncate(float, order=2)
# truncate (round) a floating number with the number digits
# specified by "order".
# e.g., if order=3, -0.012345 => -0.0123; 6.6666 => 6.67
exponent = 10**(-Math::log10(float.abs).floor+order-1)
(float * exponent).round.to_f/exponent
end
def unit_vect( vxuxratio, vyuyratio, # (V cood length)/(actual length)
uxunit=nil, uyunit=nil, # Unit vect len in U coord(float).
# If non nil, override v[xy]unit
vxunit=0.05, vyunit=0.05, # Unit vect len in V coord.
# Used only when u[xy]unit==nil
vxuoff=0.03, vyuoff=0.0, # specify v[xy]uloc by offset
vxuloc=nil, vyuloc=nil, # starting position of unit vect
annotate = true, # whether to show unit vect len
rsizet=nil,
index=3 # line index of unit vect
)
if uxunit && uyunit
vxunit = vxuxratio * uxunit
vyunit = vyuyratio * uyunit
else
uxunit = vxunit / vxuxratio
uyunit = vyunit / vyuyratio
end
uxunit = __truncate( (uxusv=uxunit) ) # round down to a good number
uyunit = __truncate( (uyusv=uyunit) ) # round down to a good number
vxunit = vxunit * (uxunit/uxusv)
vyunit = vyunit * (uyunit/uyusv)
if !(vxuloc && vyuloc)
vx0,vx1,vy0,vy1 = DCL.sgqvpt
vxuloc = vx1 + vxuoff
vyuloc = vy0 + vyuoff
end
DCL.sglazv( vxuloc, vyuloc, vxuloc+vxunit, vyuloc, 1, index )
DCL.sglazv( vxuloc, vyuloc, vxuloc, vyuloc+vyunit, 1, index )
if annotate
msg = "XUNIT=#{sprintf("%.2g",uxunit)} YUNIT=#{sprintf("%.2g",uyunit)}"
rsizet = DCL.uzpget('rsizel1') if !rsizet
before = uz_set_params({'rsizec1'=>rsizet})
DCL.uxsttl('b',msg,0.0)
uz_set_params(before)
end
end
def 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
vxuxratio = factor*cn*ax
vyuyratio = factor*cn*ay
before = ug_set_params( {'LNRMAL'=>false, 'LMSG'=>false,
'XFACT1'=>1.0, 'YFACT1'=>1.0} )
DCL.ugvect( vxuxratio*fx, vyuyratio*fy )
ug_set_params( before )
unit_vect_info = [ vxuxratio, vyuyratio, fx.max, fy.max ]
return unit_vect_info
end
end
end
-----------------
テスト
-----------------
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
uvectinfo = DCLExt.flow_vect(fx, fy, 3.0)
DCLExt.unit_vect( *uvectinfo ) # 最も長い矢印にほぼ合わせて表示
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
uvectinfo = DCLExt.flow_vect(fx, fy, 3.0)
DCLExt.unit_vect( *uvectinfo[0..1] ) # V座標での長さを約 0.05(default) に
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
uvectinfo = DCLExt.flow_vect(fx, fy, 1.0, 2, 1)
DCLExt.unit_vect( *uvectinfo[0..1] ) # V座標での長さを約 0.05(default) に
DCL.grcls
-----------------