一様密度の立方体惑星の重力ポテンシャル

ちょっと前に話題になった立方体惑星の話ですが、一辺 $2a$ で密度 $\rho$ の立方体型の惑星が作る重力ポテンシャルを計算してみました。座標系の原点は惑星の中心で、各辺は x,y,z 軸に平行としています。
$\Phi(x,y,z)=G\rho[A(x-a,y-a,z-a)-A(x-a,y-a,z+a)-A(x-a,y+a,z-a)+A(x-a,y+a,z+a)\cr\qquad\qquad-A(x+a,y-a,z-a)+A(x+a,y-a,z+a)+A(x+a,y+a,z-a)-A(x+a,y+a,z+a)]$
G は万有引力定数で、A は次のような関数です。
$A(X,Y,Z)=XY\log\left(Z+\sqrt{X^2+Y^2+Z^2}\right)+YZ\log\left(X+\sqrt{X^2+Y^2+Z^2}\right)+ZX\log\left(Y+\sqrt{X^2+Y^2+Z^2}\right)\cr\qquad\qquad+X^2\arctan\left(\frac{Y+Z+\sqrt{X^2+Y^2+Z^2}}{X}\right)+Y^2\arctan\left(\frac{Z+X+\sqrt{X^2+Y^2+Z^2}}{Y}\right)+Z^2\arctan\left(\frac{X+Y+\sqrt{X^2+Y^2+Z^2}}{Z}\right)$
惑星表面のポテンシャルを計算しようとすると、0 で割ったり 0 の対数を取ったりすることになりますが、そこは適当な極限値を使って下さい。