[dennou-ruby:002268] Re: Ruby-ISPACK 作成中

[Yuji Kitamura <kitamura@xxxxxxxxxxxxxxxxxx>]

北村です．

> 堀之内さんに言ったら作業変数を隠蔽するような上位のラッパーもほ 
> しい
> っていうことですのでそれもおいおいできればいいかなと思います。

自分が利用しているN2PACKの一部だけですが，ラップしてみまし 
た．
作業変数の隠蔽と，グリッドデータ，スペクトルデータが1次元 
配列でなくても
動くことを目標にしています．以下に実装例とテストプログラムを添付 
します．

サンプルでは，スペクトル変換のためのオブジェクトを作成したら，そ 
のオブ
ジェクトにデータを引数として与えれば，諸々の演算を行えるようにし 
てあり
ます．テストプログラムでは，ガウシアンを与えてスペクトル空間で 
x微分して
います．それなりにspmodelの雰囲気に近づいてはいると思いま 
す．
ちなみにx微分，y微分はN2PACKのラップではなく， 
自前で用意したものです．

どういう設計にするのがいいのかがよく分からないので，もっと適切な 
方法が
ありましたら教えていただければ幸いです．

require "narray"
require "ispack_raw"

class N2Pack
  def initialize(im, jm, nx, ny)
    @im, @jm, @nx, @ny = im, jm, nx, ny
    @itj = NArray.int(5)
    @iti = NArray.int(5)
    @tj = NArray.float(2*jm)
    @ti = NArray.float(2*im)
    @wkfft = NArray.float(im*jm)
    @itj, @tj, @iti, @ti = ISPACK::n2init(@jm, @im)
  end

  def forward(grid)
    shape = grid.shape
    grid.reshape!(grid.total)
    obj = ISPACK::n2g2sa(@ny ,@nx, @jm, @im, grid,
                 @wkfft, @itj, @tj, @iti, @ti)
    grid.reshape!(*shape)
    return obj
  end

  def backward(spectrum)
    shape = spectrum.shape
    spectrum.reshape!(spectrum.total)
    obj = ISPACK::n2s2ga(@ny ,@nx, @jm, @im, spectrum,
             @wkfft, @itj, @tj, @iti, @ti)
    spectrum.reshape!(*shape)
    return obj
  end

  def dx(spectrum)
    tmp = spectrum.reshape(2*@ny+1, 2*@nx+1)
    factor = -(NArray.float(1,tmp.shape[1]).indgen! - (tmp.shape[1]  
- 1)/2)
    tmp = tmp[-1..0, -1..0]
    return (tmp*factor).reshape!(tmp.total)
  end

  def dy(spectrum)
    tmp = spectrum.reshape(2*@ny+1, 2*@nx+1)
    factor = -(NArray.float(tmp.shape[0],1).indgen! - (tmp.shape[0]  
- 1)/2)
    tmp = tmp[-1..0, -1..0]
    return (tmp*factor).reshape!(tmp.total)
  end

end

=begin
###
### test program
###

require "numru/dcl"

include NumRu
include NMath

## グラフの描画
def draw(grid)
  DCL.grfrm
  DCL.grswnd(0.0, 1.0, 0.0, 1.0)
  DCL.uspfit
  DCL.grstrf

  DCL.ussttl('X', '', 'Y', '')
  DCL.uelset('ltone', true)
  DCL.uetone(grid)
  DCL.usdaxs
  DCL.udcntr(grid)
end

## 初期設定
i = 64
j = 64
nx = 31
ny = 31
x = NArray.float(1,i).indgen!/i
y = NArray.float(j,1).indgen!/j
grid = exp(-100.0*((x-0.5)*(x-0.5) + (y-0.5)*(y-0.5)))

DCL.gropn(1)

draw(grid)
## Fourier  変換 + 演算
n2pack = N2Pack.new(i, j, nx, ny)
grid = n2pack.backward(n2pack.dx(n2pack.forward(grid))).
       reshape!(j, i).transpose(1,0)

draw(grid)

DCL.grcls
=end

----
Yuji Kitamura, Ph.D. <kitamura@xxxxxxxxxxxxxxxxxx>
Division of Earth and Planetary Sciences,
Graduate School of Science, Kyoto University
Tel: +81-75-753-3933/Fax: +81-75-753-3715