ドーナツ地球（３）

ドーナツ地球（２）からの続き 

12個の球で近似するくらいならば，いっそのことドーナツ型の内部に一様分布する無数の質点モデルをとったほうがよいかもしれない。$\bm{r}$の位置における単位質量に対する万有引力のポテンシャルは，$V(\bm{r}) = \sum_{i=1}^{N} \dfrac{G m_i }{\sqrt{(\bm{r}-\bm{r}_i )^2}}$である。$N$ 個の質点は，質量 $m_i = \frac{M}{N}$，位置ベクトル$\ \bm{r}_i\ $で与えられ，空間に一様分布している。

トーラス（ドーナツ）は，円環の半径が$R$で，断面が半径$a$の円形であるとする。$\ \bm{r}=(x,y,z)\ $とすると，$x=\cos\phi (R + r \cos\theta)$，$ y=\sin\phi (R + r \sin\theta)$，$z = r \sin\theta )$ と表わすことができる。ただし，$\phi$はトーラス中心に対し円環上で$x$軸から測った方位角，$\theta$はトーラス断面円の中心に対し$x-y$平面から測った仰角で，トーラス中心に遠い方から測ったものである。

トーラス断面円の中心からの距離が$r$であり，$(x-R \cos\phi)^2+(y-R \sin\phi)^2 +z^2 = r^2\ $である。そこで，トーラス面上の点は，$(x-R \cos\phi)^2+(y-R \sin\phi)^2 +z^2 = a^2\ $を，トーラス内部の点は不等式$\ (\sqrt{x^2+y^2}-R)^2 + z^2 < a^2\ $を満足する。

m = 1000000; R = 2.0; r = 1.0;

t = Table[{RandomReal[{-3, 3}], RandomReal[{-3, 3}], 

    RandomReal[{-3, 3}]}, m];

(* s = Select[t,(#[[1]]^2+#[[2]]^2+#[[3]]^2)<1 &]; *)

s = Select[t, ((Sqrt[#[[1]]^2 + #[[2]]^2] - R)^2 + #[[3]]^2) < 1 &]; 

smax = Length[s]

n[w_] := 1/Sqrt[w[[1]]^2 + w[[2]]^2 + w[[3]]^2 + 1/smax]

ListPointPlot3D[s, BoxRatios -> Automatic]

図１：トーラス内部に一様分布する質点

p = {4 - z, 0, 0};

d = Table[p - s[[i]], {i, 1, smax}];

w = 1/smax*Sum[n[d[[k]]], {k, 1, smax}];

v[z_] := w

g[1] = Plot[v[z], {z, 0, 8}]

p = {0, 0, z - 4};

d = Table[p - s[[i]], {i, 1, smax}];

w = 1/smax*Sum[n[d[[k]]], {k, 1, smax}];

v[z_] := w

g[2] = Plot[v[z], {z, 0, 8}]

図２：万有引力ポテンシャルの値（左x軸上，右z軸上）

不連続な質点モデルで計算しているため，たまたま観測点が質点と重なると大きなスパイクが現れる。この現象を緩和するため，分母にくる距離の項 n[d[[k]]] において，1/smaxという微小項を加えている。