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の分母全体の括弧がありませんでしたので訂正しました．