[dennou-ruby:001815] Re: GPhys::EP_flux

[Takeshi Horinouchi <horinout@xxxxxxxxxxxxxxxxxx>]

堀之内です。

＃ 石岡さんdennou-ruby ML 入ってましたっけ？

塚原君は作業で頑張ってますので、私が代りにメールします。
作戦変更レポート。

● SHTLIB (in DCL/MATH2にある球面調和関数展開ライブラリー) について

等間隔データが対象ということで、等間隔化して提供されている客観解
析データの処理にはもってこいということで使う話でした。RubyDCLに
含まれてますので、すぐ使えるのも魅力。ただ、RubyDCLのSHTLIBイン
ターフェースにはちょっとバグがあるのを塚原君が発見し、直してくれ
ました。次期リリースには取り込まれます。（もしすぐ使いたいという
声があればすぐにリリースしますが、今のところ SHTLIB 以外のアップ
デートはありませんので。）

[[ここで石岡さんに質問です。初期化ルーチン shtint では、

   JM >= (MM+1)/2, IM >= MM+1 でなければならない

と書いてありますが、やってみると IM, JM がそれよりある程度小さく
てもうまく行くように見えます。ただし小さすぎると実際に変換→逆変
換で戻ってきませんし、下手すると segmentation fault がでます。実
際のところ、どうなんでしょう。特に、MM+1 の 「+1」は必要ですか。
とすると、例えば東西に10分割した場合に、波数が 5 まででなく 4 ま
でになってしまいます。ちなみに他のルーチンのドキュメントに少々間
違いを発見しました。そのうち塚原君から報告があると思います]]
(石岡さんへの質問おわり)

さて、SHTLIBはどうやら三角切断専用なようです。東西波数を南北より
も小さくとれないので、EP flux で出てくるような帯状平均物理量、す
なわち波数ゼロ成分だけに関する微分計算に使うには不向です。むりや
り東西方向に同じものをコピーして（つまり ary(ny) --> any(nx,ny)
のように広げちゃう）使うというのも考えたのですが、あまりに無駄が
多い(∵ nx=~2*ny)ので、やめようということになりました。

あと、下手な関数を与えるとギブス現象がでます。参考まで。ちゃんと
考察してませんが、極での特異性があるべき状態からはずれると起るよ
うに思います。よって、球面geometry的にリーズナブルな関数を与え、
かつ球面geometry的にリーズナブルな計算をしてる限りは問題にならな
いような気がします。ということで話はつきそうなもんですが、一応以
下に具体例をあげます。ちゃんと考えないとトラブるという実感が共有
して頂けると思います。

テストその１：  g(λ,φ) = cosφ (帯状一様)

shts2g で isw=1 とすると ∂g/∂φ が計算できますが、極で１になる
べきところがゼロになり、そのぶん相当シビアなギブスの角がでました。
一方、shtg2s で isw=1 とすると(cosφ)^-1∂(g cosφ)/∂φ が計算
できますが、この場合は問題ありません。実際、EP flux 計算で帯状平
均量に対して適用するのは、問題のなかったほうの後者です。

テストその２：  g(λ,φ) = sinλ cosφ

この場合は ∂g/∂φ はうまく行き、(cosφ)^-1∂(g cosφ)/∂φ は
激しいギブスの角がでました。

● で、微分はどうするの

球面というgeometryに特化したモジュールを作らないことになったので、
汎用の微分用モジュールをつくる。モジュール名は NumRu::Deriv か
NumRu::Derivative（対象はNArray）。また、これを GPhys に組み込ん
だものもつくる。

とりあえず、中央差分（dx_i = x_i+1 - x_i-1 タイプ）をつくる。
境界条件は取り替え可。デフォルトは dx_0 = x_1 - x_0, 
dx_N  = x_N - x_N-1。実装は境界条件にあわせて配列を拡張したのち、
x[2..-1] - x[0..-3] をとる。境界条件にあわせた拡張をメソッド化す
ることで、境界条件変更を容易にする。ちなみに今回の目的に照らすと、
極では上記のデフォルトのが妥当と考えられるので、とりあえずこの１
種類のみをつくる。（滞在時間が限られるので、基本方針は、「汎用化
のことを考えて設計しつつ、実装は必要最低限で突き進む。」）


> =begin
> = GPhys 用 EP-flux 計算モジュール作成作戦メモ
> 
> このメモは, 京大生存圏の堀之内さん部屋で行われた EP-flux モジュール作成に関する
> 堀之内さんと塚原の打ち合わせメモである.
> 
> == 履歴
> 2004/08/03(火) 塚原大輔 -- 新規作成.
> 
> 
> == 作成するモジュール
> 名称は仮である.
> 
> * ((<module NumRu::Spherical>)) / ((<module NumRu::GPhys::Spherical>))
>   * 球面座標での微分演算を行う関数群. NArray, GPhys の双方に対して作成.
> * ((<module NumRu::Cartesian>)) / ((<module NumRu::GPhys::Cartesian>))
>   * 直交直線座標での微分演算を行う関数群. NArray, GPhys の双方に対して作成.
> * ((<module NumRu::GPhys::EP_flux>))
>   * EP_flux を算出する関数. GPhys に対して作成.
> 
> === 微分演算 モジュール
> 
> * NumRu::Spherical 
>   * 球面座標での微分演算モジュール.
>     * xderiv, yderiv, rot, div, xdiv, ydiv, grad, xgrad, ygrad, lapla, ilapla, jacobian, .., ydiv_m0, ygrad_m0, yderiv_m0, etc.
>   * 等間隔グリッドデータを前提とする
>     * その他のグリッド(特にガウス緯度)についてはおいおい対処.
>     * ISPACK を導入しやすいよう意識して枠組みを作成
>   * スペクトル変換, 微分については DCL の SHTLIB を利用
>     * ((<URL|http://www.gfd-dennou.org/arch/ruby/products/ruby-dcl/ruby-dcl-doc/math2/node24.html>))
> * NumRu::Cartesian 
>   * デカルト座標での微分演算モジュール.
>     * スタッガードに対する関数はとりあえず作らない.
>   * とりあえず差分法
>   * cderiv(中央差分), fderiv(前方差分)..etc.
> 
> * NumRu::GPhys::Spherical
>   * NumRu::Spherical を利用
>   * データの定義領域のチェック(緯度は何度から何度までか, 緯度経度なのか緯度のみか, とか)
>     * メソッド化して, 下請けに応じて取り替えられるようにする
>       * ISPACK と SHTLIB でデータの定義領域が(ちょっと)異なる
>     * 軸の単位をチェック(角度かどうか)して, 判定
>       * 北から南にデータが入っている場合(NCEP はそう)は, 逆向きにしてあげる
>     * 将来的には, (可能な限り)領域を適合させられるようにする.
>   * 単位, 地球半径を考慮
> * NumRu::GPhys::Cartesian
>   * NumRu::Cartesian を利用
>   * 座標変数の自動更新
>     * 例 : 半整数グリッドと整数グリッド
>   * 単位を考慮
> 
> === EP_flux モジュール(NumRu::GPhys::EP_flux)
> * 基本は近似を行わない, いわゆるフルセットの EP_flux
>   * 近似を導入しようとも, 微分演算モジュールは作らねばならないため, コストは一緒.
>   * オプションで近似形を用いるかどうか切り替えられるようにする.
> * 入力変数
>   * 東西風, 南北風, 鉛直風(omega or w), 温度(temp or theta)
>   * 鉛直座標は圧力座標(P or logP)
> * 出力変数
>   * Fy, Fz, uv_bar, ...
> * divF 計算メソッドを作成
>   * Fy, Fz を引数にとる
> 
> == スケジュール(目標)
> 
> (1) NArray 版微分演算モジュールの作成     -- 8/4(水)
> (2) GPhys  版微分演算モジュールの作成     -- 8/5(木)
> (3) EP_flux 算出モジュールの作成          -- 8/6(金), 8/7(土)
> (4) チュートリアルの作成                  -- 8/8(日), 8/8(月)
> 
> === 留意点
> * コーディングとテストは同時進行
> * テストの方針
>   * 解析的解のデータと出力を比較する
> * テストがある程度通るようになったら, ドキュメントを作成(日本語 and RD)
> 
> =end