SPHEREPACK

07/10/06 19:51

SPHEREPACKはLegendre変換をするFortranライブラリ.
Mac OS X上でg95を使ってコンパイルして使っている。

解像度変換のユーティリティルーチンtrssph()が配列サイズが大きいときにNaNを返すことがある。試した範囲では、288×145を1200×600にするのは大丈夫で、1440×720にするのはだめだった。

メモリ節約のため、trssph()は、Legendre陪関数を逐次計算するサブルーチンsphcom.fにあるlegin()でルジャンドル陪 関数が正しく計算されないようである。アルゴリズムとしては、robustな方法を使っているはずだし、ベクトルのLegendre変換をするルーチンは うまく動作する。また、SPHEREPACKを呼んでいるNCLでは、正しく動作する。

ソースを見たところ、ルジャンドル陪関数を計算するときに必要な係数abel, bbel, cbelが大きな整数のかけ算になっていたためであることが分かった。shagc.f, shsgc.f, shags.f, shsgs.f, shigc.f, shigs.fでは


- abel(imn)=sqrt(float((2*n+1)*(m+n-2)*(m+n-3))/
- 1 float(((2*n-3)*(m+n-1)*(m+n))))
- bbel(imn)=sqrt(float((2*n+1)*(n-m-1)*(n-m))/
- 1 float(((2*n-3)*(m+n-1)*(m+n))))
- cbel(imn)=sqrt(float((n-m+1)*(n-m+2))/
- 1 float(((n+m-1)*(n+m))))
+ abel(imn)=sqrt((2.*n+1.)*(m+n-2.)*(m+n-3.)/
+ 1 ((2.*n-3.)*(m+n-1.)*(m+n)))
+ bbel(imn)=sqrt((2.*n+1.)*(n-m-1.)*(n-m)/
+ 1 ((2.*n-3.)*(m+n-1.)*(m+n)))
+ cbel(imn)=sqrt((n-m+1.)*(n-m+2.)/
+ 1 ((n+m-1.)*(n+m)))



と修正 (diff -uの出力)したところ、問題なく動作した。またvhags.fとvhsgs.fは以下の通り修正した。



- abel = dsqrt(dble(float((2*n+1)*(m+n-2)*(m+n-3)))/
- 1 dble(float((2*n-3)*(m+n-1)*(m+n))))
- bbel = dsqrt(dble(float((2*n+1)*(n-m-1)*(n-m)))/
- 1 dble(float((2*n-3)*(m+n-1)*(m+n))))
- cbel = dsqrt(dble(float((n-m+1)*(n-m+2)))/
- 1 dble(float((m+n-1)*(m+n))))
+ abel = dsqrt((2.d0*n+1.d0)*(m+n-2.d0)*(m+n-3.d0)/
+ 1 ((2.d0*n-3.d0)*(m+n-1.d0)*(m+n)))
+ bbel = dsqrt((2.d0*n+1.d0)*(n-m-1.d0)*dble(n-m)/
+ 1 ((2.d0*n-3.d0)*(m+n-1.d0)*(m+n)))
+ cbel = dsqrt((n-m+1.d0)*(n-m+2.d0)/
+ 1 ((m+n-1.d0)*(m+n)))


- abel = dsqrt(dble(float((n+m)*(n-m+1))))/dcf
- bbel = dsqrt(dble(float((n-m)*(n+m+1))))/dcf
+ abel = dsqrt((n+m)*(n-m+1.d0))/dcf
+ bbel = dsqrt((n-m)*(n+m+1.d0))/dcf

2009/2/23 追記


パッチにバグがありました.vhags.fとvhsgs.fの修正の前半,abel, bbel, cbelの分母全体の括弧がありませんでしたので訂正しました.

コメント

このブログの人気の投稿

静的サイト生成ツールの導入

uptex

TeXのソースを分割してTeXShopでタイプセット