[dennou-ruby:002124] Re: ruby資 源共有化の加速

[Shigenori OTSUKA <otsuka@xxxxxxxxxxxxxxxxxx>]

大塚です。

> 堀之内です。
> 
> > ありがとうございます。
> > では報告させていただきたいと思います。
> 
> よろしくお願いします。急ぎませんが、プログラムに記載する標題を私
> まで送ってください。
>  
> > #もうひとつ、流線を描くルーチンもありました。
> 
> 流線関数のコンターでなく、（発散成分があってもいい）流線なんです
> よね？ それは素晴らしい。以前から DCL にほしいと思ってましたが、
> 私としては実現されるんだったら DCL の中でも外でも構いません。
> （どうせ Ruby 経由でしか使わないので）
> 
> もしよければ、MLにソースを流すか、どこかに置いて貰えませんか。
> 

標題についてはちょっと考えさせてください。
そのうち送ります。

流線の描画ルーチンですが、とりあえず今の物を添付します。
いいかげんな書式ですがお許しください。

2種類のサブルーチンがありますが、最初のもの(自動描画)を
使うとよいと思います。

線の密度や矢羽の間隔がオプションで変えられない状態なので、
それだけでも改良したいと思っています。

================================
京都大学理学研究科
気象学研究室 M1
大塚成徳 Shigenori OTSUKA
email: otsuka@xxxxxxxxxxxxxxxxxx

#!/usr/bin/env ruby

=begin
等間隔グリッドデータに対する流線の描画ルーチン
2004.9.9 - 2004.9.11
by Shigenori OTSUKA

def streamline(u,v,xx,yy,dl,index)
  u : x方向成分格子データ
  v : y方向成分格子データ
  xx: x座標格子データ
  yy: y座標格子データ
  dl: 格子間隔
  index: 線の種類

=end

require "narray"
require "numru/dcl"
include NumRu
include NMath

def streamline(u,v,xx,yy,dl,index)
  #等間隔等方グリッドデータに対する流線のDCL描画ルーチン
  #積分開始点の適正間隔の計算
  istart = 20
  if u.shape[0]>u.shape[1] then
    #nstart = (u.shape[0] / 40).to_i
    nstart = (u.shape[0] / istart).to_i
  else
    #nstart = (u.shape[1] / 40).to_i
    nstart = (u.shape[1] / istart).to_i
  end
  nstart2 = nstart**2  
  tarrow = 60.0 * 60.0 * 48.0 #積分48時間おきに矢羽を描く 
  xy = []
  if u.shape[0]>u.shape[1] then
    step = u.shape[0]/100
  else
    step = u.shape[1]/100
  end
  step = 1 if step < 1
  NArray.sfloat((u.shape[0]/step).to_i).indgen!(1,step).each{|i|
    NArray.sfloat((u.shape[1]/step).to_i).indgen!(1,step).each{|j|
      skipflag = false
      if xy.length != 0 then
        a = NArray.to_na(xy)
        ax = a[0,true]
        ay = a[1,true]
        if (((ax - i)**2 + (ay - j)**2).lt(nstart2)).max >= 1 then 
          skipflag = true
        end
      end
      skipflag = true if i >= (u.shape[0]-1)
      skipflag = true if j >= (u.shape[1]-1)
      #p [i,j]
      if !skipflag then
        sx = i
        sy = j
        #print "ok\n" 
        dtime = [0.1,-0.1]
        xy_temp = []
        dtime.each{ |dt|
          x = sx
          y = sy
          if dt > 0 then
            t = 0.0
          else
            t = tarrow
          end
          loopflag = false
          for n in 0..10000
            xy_temp.push [x,y]
            flag = false
            u0 = interp_2D(u,x,y)
            v0 = interp_2D(v,x,y)
            while ((u0*dt/dl).abs>1 or (v0*dt/dl).abs>1)
              dt = dt * 0.9
              #print "dt has benn changed to #{dt}\n"
            end
            while ((u0*dt/dl).abs<0.5 and (v0*dt/dl).abs<0.5)
              dt = dt * 1.1
              #print "dt has benn changed to #{dt}\n"
            end
            while (x+u0*dt/dl<0.0 or x+u0*dt/dl>((u.shape[0]-1).to_f) or y+v0*dt/dl<0.0 or y+v0*dt/dl>((u.shape[1]-1).to_f))
              dt = dt * 0.9
              flag = true
              #print "dt has benn changed to #{dt} to draw within the window\n"
            end
            dx = u0 * dt / dl
            dy = v0 * dt / dl
            t += dt.abs

            xx0 = []
            yy0 = []
            [[x,y],[x+dx,y+dy]].each{|label|
              xx0.push interp_2D(xx,label[0],label[1])
              yy0.push interp_2D(yy,label[0],label[1])
            }
            if t > tarrow then
              if dt >= 0 then
                DCL.sglset("LPROP",false)
                DCL.sglazu(xx0[0],yy0[0],xx0[1],yy0[1],1,index)
              else
                DCL.sglset("LPROP",false)
                DCL.sglazu(xx0[1],yy0[1],xx0[0],yy0[0],1,index)
              end
              t = t - tarrow
            else
              DCL.sglnzu(xx0[0],yy0[0],xx0[1],yy0[1],index)
            end
            x = x + dx
            y = y + dy
            if xy.length != 0 then
              a = NArray.to_na(xy)
              ax = a[0,true]
              ay = a[1,true]
              if (((ax - x)**2 + (ay - y)**2).lt(nstart.to_f/2.0)).max >= 1 then 
              #if (((ax - x)**2 + (ay - y)**2).lt(1)).max >= 1 then 
                flag = true
              end
            end
            if loopflag then
              if xy_temp.length < 10 then
                ntempend = xy_temp.length / 2
              else
                ntempend = xy_temp.length-10
              end
              #for ntemp in 0..ntempend
              #  if((xy_temp[ntemp][0]-x)**2 + (xy_temp[ntemp][1]-y)**2) < 1 then
              #    flag = true
              #  end
              #end
              if xy_temp[0..ntempend].length != 0 then
                a = NArray.to_na(xy_temp[0..ntempend])
                ax = a[0,true]
                ay = a[1,true]
                #if (((ax - x)**2 + (ay - y)**2).lt(nstart.to_f/4.0)).max >= 1 then 
                if (((ax - x)**2 + (ay - y)**2).lt(1)).max >= 1 then 
                  flag = true
                end
              end
            end
            loopflag = true if ((sx-x)**2+(sy-y)**2) > 2
            break if flag
            break if x<=0.0
            break if x>=(u.shape[0]-1).to_f
            break if y<=0.0
            break if y>=(u.shape[1]-1).to_f
            break if u0**2 + v0**2 < 0.05#0.5
          end
        }
        xy_temp.each{|temp|
          xy.push temp
        }
        xy_temp = nil
        GC.start
      end 
    }
  }
  xy = nil
  GC.start
end

def streamline_manual(u,v,xx,yy,points,index)
  #等間隔等方グリッドデータに対する流線のDCL描画ルーチン
  #開始点手動設定

  tarrow = 5.0
  xy = []
  points.each{|pnts|
    sx = pnts[0]
    sy = pnts[1]
    #print "ok\n" 
    dtime = [0.1,-0.1]
    xy_temp = []
    dtime.each{ |dt|
      x = sx
      y = sy
      if dt > 0 then
        t = 0.0
      else
        t = tarrow
      end
      loopflag = false
      for n in 0..10000
        xy_temp.push [x,y]
        flag = false
        u0 = interp_2D(u,x,y)
        v0 = interp_2D(v,x,y)
        while ((u0*dt).abs>1 or (v0*dt).abs>1)
          dt = dt * 0.9
          #print "dt has benn changed to #{dt}\n"
        end
        while ((u0*dt).abs<0.5 and (v0*dt).abs<0.5)
          dt = dt * 1.1
          #print "dt has benn changed to #{dt}\n"
        end
        while (x+u0*dt<0.0 or x+u0*dt>((u.shape[0]-1).to_f) or y+v0*dt<0.0 or y+v0*dt>((u.shape[1]-1).to_f))
          dt = dt * 0.9
          flag = true
          #print "dt has benn changed to #{dt} to draw within the window\n"
        end
        dx = u0 * dt
        dy = v0 * dt
        t += dt.abs
        
        xx0 = []
        yy0 = []
        [[x,y],[x+dx,y+dy]].each{|label|
          xx0.push interp_2D(xx,label[0],label[1])
          yy0.push interp_2D(yy,label[0],label[1])
        }
        if t > tarrow then
          if dt >= 0 then
            DCL.sglset("LPROP",false)
            #DCL.sglazu(x,y,x+dx,y+dy,1,1)
            DCL.sglazu(xx0[0],yy0[0],xx0[1],yy0[1],1,index)
          else
            DCL.sglset("LPROP",false)
            #DCL.sglazu(x+dx,y+dy,x,y,1,1)
            DCL.sglazu(xx0[1],yy0[1],xx0[0],yy0[0],1,index)
          end
          t = t - tarrow
        else
        #DCL.sglnzu(x,y,x+dx,y+dy,1)
          DCL.sglnzu(xx0[0],yy0[0],xx0[1],yy0[1],index)
        end
        x = x + dx
        y = y + dy
        if xy.length != 0 then
          a = NArray.to_na(xy)
          ax = a[0,true]
          ay = a[1,true]
          if (((ax - x)**2 + (ay - y)**2).lt(1)).max >= 1 then 
            flag = true
          end
        end
        #xy.each{|xy1|
        #  if xy1[0]-x<1 and xy1[1]-y<1 then
        #    if((xy1[0]-x)**2 + (xy1[1]-y)**2) < 1 then
        #      flag = true
        #    end
        #  end
        #}
        if loopflag then
          if xy_temp.length < 10 then
            ntempend = xy_temp.length / 2
          else
            ntempend = xy_temp.length-10
          end
          #for ntemp in 0..ntempend
          #  if((xy_temp[ntemp][0]-x)**2 + (xy_temp[ntemp][1]-y)**2) < 1 then
          #    flag = true
          #  end
          #end
          if xy_temp[0..ntempend].length != 0 then
            a = NArray.to_na(xy_temp[0..ntempend])
            ax = a[0,true]
            ay = a[1,true]
            if (((ax - x)**2 + (ay - y)**2).lt(1)).max >= 1 then 
              flag = true
            end
          end
        end
        loopflag = true if ((sx-x)**2+(sy-y)**2) > 10
        break if flag
        break if x<=0.0
        break if x>=(u.shape[0]-1).to_f
        break if y<=0.0
        break if y>=(u.shape[1]-1).to_f
        break if u0**2 + v0**2 < 0.5
      end
    }
    xy_temp.each{|temp|
      xy.push temp
    }
  }
end

def streamline_core(u,v,xx,yy,i,j)

end

def interp_2D(data,x,y)
  # linear interpolation
  data0 = 0
  if x.ceil == x.floor then
    if y.ceil == y.floor then
      data0 = data[x.to_i,y.to_i]
    else
      data0 = data[x.to_i,y.ceil]*(y-y.floor) + data[x.to_i,y.floor]*(y.ceil-y)
    end
  elsif y.ceil == y.floor then
    data0 = data[x.ceil,y.to_i]*(x-x.floor) + data[x.floor,y.to_i]*(x.ceil-x)
  else
    #print "x:real,y:real #{x.floor}<x<#{x.ceil},#{y.floor}<y<#{y.ceil}\n"
    data0 = (data[x.ceil,y.ceil]*(x-x.floor) + data[x.floor,y.ceil]*(x.ceil-x))*(y-y.floor)
    data0 +=(data[x.ceil,y.floor]*(x-x.floor) + data[x.floor,y.floor]*(x.ceil-x))*(y.ceil-y)
  end
  return data0
end

#test----------------------------------------------
#
#u = NArray.sfloat(100,100)
#v = NArray.sfloat(100,100)
#ベクトル場の設定
#for i in 0...100
#  for j in 0...100
#    u[i,j] = -10.0 * Math.sin(Math.atan2(i-50,j-50))
#    v[i,j] = 10.0 * Math.cos(Math.atan2(i-50,j-50))
#  end
#end
#
#
#ベクトル場の描画
#DCL.gropn(1)
#DCL.grfrm
#DCL.grswnd(0.0,99.0,0.0,99.0)
#DCL.grsvpt(0.2,0.8,0.2,0.8)
#DCL.grstrf
#DCL.usdaxs
#for i in 0...(u.shape[0])
#  for j in 0...(v.shape[1])
#    if i%5==0 and j%5==0 then
#      DCL.sglau(i,j,i+u[i,j]/10.0,j+v[i,j]/10.0)
#    end
#  end
#end
#
#流線の描画
#DCL.grfrm
#DCL.grswnd(0.0,99.0,0.0,99.0)
#DCL.grsvpt(0.2,0.8,0.2,0.8)
#DCL.grstrf
#DCL.usdaxs
#streamline(u,v)
#DCL.grcls