muttenz's blog

スイス星空だより

astrolibRでhelio_jdが誤った値を出す原因がわかった

まあ、こういう間違いの原因を私が探す義理は全くないのですが、どこか気持ち悪いのでastrolibRのスクリプトをよく見てみました。

原因はその中で使われているxyz()のスクリプトでした。

それも割に単純な間違いです。

スクリプトの最後の方でxyz座標をいろいろな三角関数を足したり引いたりして出すところです。

オリジナルのスクリプトでは

x = 0.99986 * cos(el)
-0.025127 * cos(g - el)
+0.008374 * cos(g + el)
+0.000105 * cos(g + g + el)
+6.3e-05 * t * cos(g - el)

-----------------

この後もずらずらと三角関数の項が並ぶのですが、

これをRで計算すると

> x = 0.99986 * cos(el)
> -0.025127 * cos(g - el)
[1] -0.007243908
> +0.008374 * cos(g + el)
[1] -0.007354704
> +0.000105 * cos(g + g + el)
[1] -0.0001032314
> +6.3e-05 * t * cos(g - el)
[1] 0.0002700628
> +3.5e-05 * cos(g + g - el)
[1] 2.896541e-05
このように最初の行の後はそれぞれの三角関数の値は計算されますがそこで完結し、最初の項にはそれらの値が足されたり引かれたりせずに全くそのままです。

ここは各三角関数の項の前のプラスやマイナスを前の行の最後に持ってこないと計算が継続されないです。このように

  x = 0.99986 * cos(el)-
0.025127 * cos(g - el)+
0.008374 * cos(g + el)+
0.000105 * cos(g + g + el)+
6.3e-05 * t * cos(g - el)+
3.5e-05 * cos(g + g - el)-

-----------------------------------

あるいはずらずらと長ーい1行にするか。

x = 0.99986 * cos(el)-0.025127 * cos(g - el)+0.008374 * cos(g + el)+0.000105 * cos(g + g + el)+6.3e-05 * t * cos(g - el)+3.5e-05 * cos(g + g - el)-   -----------------

このようにスクリプトをなおすとパイソンで書かれたスクリプトと全く同じように動きます。

これでも、このスクリプト、プロの人が書いたのでしょうか。。。。