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

[dennou-ruby:002412] Re: integrate



堀之内です.

> 竹広です. 
> 
> GPhys のマニュアルみて, なんとなくできそうだけど
> 自信ないのでしつもんしちゃいます. 
> 
> GPhys では重みが定義されている座標に関する積分をちゃんと計算してくれる
> メソッドが用意されてますでしょうか? integrate てのがそれっぽいけど
> 
> 具体例としては, Gaussian 格子点の緯度座標を
> 
>   lat		座標
>   lat_weight	Gaussian 荷重
> 
> とした場合, です. 

何もせず普通に GPhys で NetCDF 開いた場合,格子点はセルでなく点
を代表することになってます(コンベンション対応が追い付いてない).
その場合に使われるデフォルトの積分法は台形公式です.これが,「セ
ル型」の軸の場合,セル代表点(普通は中心)な格子に沿った積分は,
セルの幅を掛けて足していく,重み方式になります.Axis クラスの機
能ですので,axis.rb に該当部分があります.概念的には,
http://ruby.gfd-dennou.org/products/gphys/gphys-intro-J/ の7枚目
のスライドをみてください.セルの幅は格子の両端の距離です.セル型
の Axis の場合は,格子の境界の位置情報は必須です.

ガウス緯度の場合は,普通はセルの境界データを作らないでしょうから,
ちょっと工夫が要ります.解決方法は以下の2パターン考えられますが,
今回のお勧めは2番目です.

  パターン1: 緯度軸セル型に適合するようにする.
     例えば,座標変数は度でなく,sin緯度にし,ガウス重みから
     ±1 の区間でのセル境界の位置を決め,緯度軸をセル型にする.

  パターン2:デフォルトの積分法に頼らない.
     実は台形公式だろうがなんだろうが,Axis クラスでは公式に
     合わせた積分重みを導出したうえで積和をとるという方式に
     してます.で,デフォルトの重み計算は @integ_weight という内
     部変数が定義されていないときにのみ発動されます.従って,
     予め @integ_weight を定義すればやりたい放題です.今の場合は
     ガウス重みを @integ_weight にすればいいんです.

この機能は GPhys 作成の初期に作ったんですが,今まで自分では 
@integ_weight を指定するケースは使ったことがないので,今回見てみ
たらバグを発見しました.一部タイポで @integ_wieghtとなってるとこ
ろがあります.また, @integ_weight を設定するメソッドを作ってな
かったので,合わせて修正して CVS コミットします.なお,本来
重みは単位付であるべきなのに,NArray になっちゃってます.
VArray 用にしないと単位がきちんと更新されない.これも変えないと...

なお,明後日から夏休みをとります.明日も電脳 ruby メールに反応す
る時間がとれるかわかりません.悪しからず.