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

[Takeshi Horinouchi <horinout@xxxxxxxxxxxxxxxxxx>]

塚原様

堀之内です。

> > 「流れ」ベクトル用のメソッドを作りました。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
-----------------