« September 2010 | Main | December 2010 »

October 30, 2010

ヤマトが小惑星帯を突っ切る時、どれだけ小惑星に接近するか?

E225今までの「小惑星帯」シリーズでは、地球から木星へのホーマン軌道を航行する宇宙船だけを計算していた。

でも、それだけじゃなくて、色々な軌道を航行する宇宙船で計算したいよね。
既に書いたように、一連のプログラムは汎用的な宇宙船の軌道に対して計算できるように作ってある。それを説明しよう。

プログラム中、宇宙船を記述している部分は、$spaceCraft と言うオブジェクトだ。このオブジェクトを作るクラスは、SpaceCraft なんだが、このテンプレートは下記のようになっている。

class SpaceCraft
  def initialize()
   << 初期化 >>
  end
  def getStartEndDateTime()
   << シミュレーション開始時刻と終了時刻 >>
    return startDateTime,endDateTime
  end
  def getXV(t)
   << シミュレーション開始時刻から宇宙船の位置と速度を計算 >>
    return x,y,z,v
  end
end

$spaceCraft = SpaceCraft.new

このようにオブジェクト $spaceCraft は、たった3つの関数しか必要ない。
一つめの「initialize()」は初期化ルーチンで、もし、何も初期化する必要がなければ、空白でも構わない。
二つめの「getStartEndDateTime()」は、シミュレーションの開始時刻と終了時刻を返す関数だ。開始時刻、終了時刻ともに Ruby の DateTime 形式の変数である。
三つめの「getXV(t)」は、シミレーション時刻 t における宇宙船の位置と速度を返す関数だ。t はRuby の DateTime 形式の変数であり、位置は X,Y,Zの3値で単位は AU、速度は単位 m/s だ。(X軸は2000年の春分の日の太陽方向、Z軸は地球軌道面垂直の北方向だ)

と言っても、具体的じゃないと使い方は良く判らないよね。
と言うことで、宇宙戦艦ヤマトが小惑星帯を突っ切る時のシミュレーションを例にして説明しよう。

まず、ネットで調べると、ヤマトは2199年10月10日に地球を出発し、翌11日に木星に到着している。(細かく言うと、この間に一度月へ行き、その後火星にワープした後、木星に行っているのだが、複雑になるので地球から直接木星に行ったことにして計算している)
と言うことで、シミュレーションの開始時刻は2199年10月10日と終了時刻は2199年10月11日に決定。

次に出発日の地球の位置と、到着日の木星の位置を調べよう。
簡単なプログラムを作ってみた。planet.rbだ。
Keplerian Elements for Approximate Positions of the Major Planets のページにある惑星の軌道データである p_elem_t1.txtp_elem_t2.txtをダウンロードして同じディレクトリに入れておく。
使い方は簡単で、下記のようにコマンドすると、その日の各惑星の座標が表示される。下は LINUX での例だが、Windows でも動くはずだ。


$ ./planet.rb 2199/10/10
日時:2199-10-10T00:00:00+00:00
水星の座標:0.24981685038202 0.199517885753347 -0.00646163832547683
金星の座標:0.104758635915238 0.712548030132175 0.00410577758340108
地球の座標:0.969811669005803 0.242169528844897 -0.000159629098671379
火星の座標:0.208607999745768 1.5342134756758 0.0271095662169784
木星の座標:4.46690521161679 -2.25921162744416 -0.089720069022941
土星の座標:8.25727296133364 -5.33923031072755 -0.240218385841203
天王星の座標:2.01593671944826 19.0180999921417 0.0440323396842922
海王星の座標:27.8666085210368 10.5986082732312 -0.860653284622775
冥王星の座標:-27.3441244017523 22.564811905219 5.49766045953591
$ ./planet.rb 2199/10/11
日時:2199-10-11T00:00:00+00:00
水星の座標:0.22560894627618 0.221734367880399 -0.00243113800090403
金星の座標:0.0846423315647759 0.715105639376437 0.00530200612400009
地球の座標:0.965221462902767 0.258757577936856 -0.000167272010251388
火星の座標:0.195246255145545 1.53723608777542 0.0274973094274193
木星の座標:4.47021365430111 -2.25211110441074 -0.0898227559121188
土星の座標:8.25999112622408 -5.33456363399447 -0.240407888588168
天王星の座標:2.0120016479658 19.0183314422802 0.0440839994093446
海王星の座標:27.8654712607142 10.6015605792454 -0.860687867932851
冥王星の座標:-27.3453017930669 22.5621221647619 5.49828872449565

これで、出発日の地球の位置と到着日の木星の位置が判った。
後は、その間を結ぶ軌道を作るだけだ。でも、たった一日で地球〜木星を楕円軌道で結べない。あまりにも速くなるので、太陽の引力を振り切ってしまうのだ。
ここでは、最も単純に出発点と到着点を直線で結んで等速直線運動にしてしまう。
具体的には次のようなプログラムになる。

class SpaceCraft
  def initialize()
    @startDateTime = DateTime.new(2199,10,10,0,0,0)
    @x0 = 0.969811669005803
    @y0 = 0.242169528844897
    @z0 = -0.000159629098671379
    @endDateTime = DateTime.new(2199,10,11,0,0,0)
    @x1 = 4.47021365430111
    @y1 = -2.25211110441074
    @z1 = -0.0898227559121188
    @vx = (@x1 - @x0) / (@endDateTime - @startDateTime)
    @vy = (@y1 - @y0) / (@endDateTime - @startDateTime)
    @vz = (@z1 - @z0) / (@endDateTime - @startDateTime)
    @v = sqrt(@vx * @vx + @vy * @vy + @vz * @vz) * AU / (24.0 * 3600.0)
  end
  def getStartEndDateTime()
    return @startDateTime,@endDateTime
  end
  def getXV(t)
    x = @x0 + @vx * (t - @startDateTime)
    y = @y0 + @vy * (t - @startDateTime)
    z = @z0 + @vz * (t - @startDateTime)
    return x,y,z,@v
  end
end

上記を反映したプログラムはasteroidInterval3Y.rbだ。

計算してみると、最も接近した小惑星は 34.4 万キロ だった。これじゃわざわざ避けて通る必要もないねえ。

ところで、ネットで調べてみると、島大介が操縦して小惑星帯を抜けるシーンは火星〜木星間の小惑星帯じゃなくて、冥王星の外の第10惑星の成れの果てという設定だったんだ。30年も前のことだから、すっかり忘れていた。
まあ、そんな遠いところに小惑星があるか判らないし、あったとしても未だ地球から観測されていないので、火星〜木星間の小惑星帯のシミュレーションで許してもらおう。

訂正:
前に計算したときの「地球〜木星のホーマン軌道」に誤りがあった。出発日と到着日に肝心の地球と木星がその場所にない。
実は、今回の planet.rb を作る時に、以前の計算の時に間違いをしていたのに気が付いた。ユリウス日と修正ユリウス日を誤って使っていたのだ。だから、以前の計算は、「地球〜木星のホーマン軌道」ではなく、「地球と同じ太陽からの距離〜木星と同じ太陽からの距離へのホーマン軌道」と思って欲しい。
お詫びして訂正する。

| | Comments (1) | TrackBack (0)

October 23, 2010

クワッドコアの威力

E224前々回のコンテンツ「宇宙船で小惑星帯を通り抜けるとき、どのくらい小惑星に接近するの?」の計算では、Phenom X4 を使ったのだが、せっかくクワッドコアを持っているのに一つのCPUコアしか使わずに計算した。そのため、計算に約28時間(正確には1665.8分間)もかかった。

そこで、プログラムを修正して、クワッドコア対応で計算してみた。

プログラムはasteroidInterval3X.rbだ。前々回のプログラムと使い方はほとんど同じなのだが、引数が一つ増えている。
$ ./asteroidInterval3X.rb 0 4 >asteroidInterval3X.dat
とやると、4つのタスクを起動し、それぞれに53万個の小惑星を4分割して、13万2500個ずつ割り当てて計算させる。
それぞれ、別のコアが並列で計算した後、最後に4つの結果をシミュレーション時間毎に最短距離を選ぶようにマージする。6コアのCPUならば、最後の引数を6にすれば、全てのコアをフルに使える。さらに、ハイパースレッド対応のCPUの場合、例えば、6コアかつハイパースレッドならば、12を引数にすれば、全てのスレッドを使うことができる。

私の場合、ハイパースレッド非対応のクワッドコアなので、引数を4にして計算した結果が、冒頭のグラフだ。赤線でシングルコアで計算した結果を、緑線でクワッドコアで計算した結果を示しているが、ほとんど赤と緑が重なっているので、計算結果はあっていることがわかる。

プログラム以外は同一の条件だ。
・Phenom X4 9350e 2GHz
・8GBメモリ
・Ubuntu 10.04 32ビットPAEカーネル
・Ruby 1.8.7

で、計算時間は、
・シングルコア:約28時間(正確には1665.8分間)
・クワッドコア(引数:4):約3時間半(正確には207.3分間)
となり、ちょうど8分の1の時間で済んでる。

8分の1・・? クワッドコアなんだから、4分の1 じゃないか? (ちょっと、わざとらしい)

実は、私が作ったアルゴリズムでは、前回説明したように、小惑星の個数の1.5乗と計算時間が比例する。つまり、小惑星の個数を少なくすると計算時間が、それ以上に減る。だから、8分の1であっているのである。

じゃあ、さらに分割すれば、もっと短い時間で計算できるはずである。
実際、引数を4から16にすると、さらに半分の103.3分で計算が済む。この場合、各コアが4つずつのタスクが割り当てられているのだが、それでも計算時間が少なくなるのだ。

さらに引数を2倍の32にすると、2の平方根分短くなるはずで、実際77.0分で終わった。

調子に乗って、引数を64にすると、50分くらいで終わるはずだ・・・と、やってみたら、オペレーションシステム自体がアボートしてしまった。

まあ、全てを別タスクにすること自体が間違っていたので、本当のところは、タスクはコアもしくは、スレッドの数だけにして、その中で、小惑星の個数を分割して計算した方が良いに決まって居る。けど、まだプログラムを改修してない。

既に、まきのさんがコメントしているように、究極的には、全ての小惑星を個々に、つまり一個ずつ計算して、後からマージするのが、最も良いことになる。
ただし、これは、現状では、「小惑星の個数の1.5乗と比例する計算時間」が支配的だが、あまり細かく分割すると「ミュレーション時間毎に最短距離を選ぶようにマージする時間」の方が大きくなってしまう可能性もあるだろう。
本当のところ、何処に最適解があるかは、もう少し作業してみないと判らない。

けど、速度が欲しいなら、RubyじゃなくてCでコンパイルしろ・・とも言われる。

| | Comments (0) | TrackBack (0)

October 20, 2010

小惑星の数が増えると、どう宇宙船に再接近する距離が縮まるか、計算時間が増えるか

前回のコンテンツで、E223「小惑星の数の平方根の逆数に比例している。計算時間は、小惑星の数の1.5乗と比例」と書いたが、何故、そのような関係になるかを考察した。

ある瞬間における宇宙船と小惑星との距離は、三次元的ものである。ところが、宇宙船が小惑星帯を通り抜けるような旅程全体での最短距離は、概念的には、イラストに示したように、速度ベクトルに垂直な平面への垂線との交点との距離になる。つまり、小惑星との距離は二次元的な直線距離になる。

小惑星の個数と各小惑星が平面で占有できる面積は反比例の関係になり、従って最短距離は小惑星の個数の平方根と比例することになのだ。

計算時間の話だが、私のアルゴリズムでは宇宙船の出発から到着まで一定時間間隔で計算しているのではない。最短にある小惑星までの近い時は細かい間隔で、遠い時は荒い間隔で計算している。具体的には、どんな事があっても最短距離の小惑星が宇宙船にぶつかって通り過ぎることの無いように、小惑星の速度の最大値(そこの太陽からの距離での第三宇宙速度)と宇宙船の速度との単純和速度で最短距離の小惑星までの距離を割った時間を計算時間の間隔としている。
つまり、計算回数は小惑星の平方根と比例するのだ。
また、一回の計算時間は単純に小惑星の数に比例する。

計算回数は小惑星の平方根と比例し、一回の計算時間は単純に小惑星の数に比例するわけだから、全体の計算時間は小惑星の数の1.5乗に比例することになる。

ここまでに、実は私が断りもなしに大胆な仮定をおいていることを、聡明な読者は気づいているかも知れない。その一つが「小惑星の動きに比較して宇宙船の相対的な速度が大きいこと」、もう一つが「宇宙船の速度ベクトルが直線的であること」である。

実は、これら二つの仮定は、少なくとも前回のコンテンツで示した地球〜木星のホーマン軌道においては、成立しているようである。と言うのも、前回のコンテンツの最後に示したグラフのように、小惑星の個数の平方根と最短距離は比例しているし、小惑星の個数の1.5乗と計算時間が比例しているからである。どうやら、宇宙船の速度ベクトルは、各時点時点でのベクトルを直線で近似しても良いようだ。

しかし、小惑星帯の中を、小惑星と同じ速度で周回している宇宙船だと、こうならないだろう。私の予想では、最短距離は小惑星の個数の立方根と比例し、計算時間は小惑星の数の1.33乗に比例すると思うのだが、まだ試していない。

| | Comments (6) | TrackBack (0)

October 12, 2010

宇宙船で小惑星帯を通り抜けるとき、どのくらい小惑星に接近するの?

E2221先日のコンテンツ小惑星同士って、どのくらい離れているの?は、結構受けて、このブログへのコメントだけじゃなくて、 Twitter の方にも反響があった。
共通して一番聞かれたのが、「小惑星同士の距離が離れているのは判ったけど、小惑星帯の中を宇宙船で通り抜けるとき、どのくらい小惑星に接近するの?」って質問だ。

確かに、これには興味があるよね。
宇宙戦艦ヤマトにしろ、21エモンにしろ、小惑星帯の中を右へ左へ、小惑星を避けて通る宇宙船の操縦のシーンがある。あれが本当かどうか、知りたいよね。

さて、計算してみようと思ったら、ある瞬間での小惑星同士の距離を測るよりも、はるかにCPUパワーが要ることに気が付いた。

で、「CPUパワーが欲しい」って、呟いてみたんだが、呟いてみるもんだねえ。
「6コアのCPUに買い換えたので、4コアのCPUが余っているので、マザボ・メモリごと譲る」とのこと。
有難くいただくことにした。

で、うちで余っていたATX電源やHDD、キーボード、マウスなどと合わせて、Ubuntu をインストできるところまで、来た。残念ながら、ケースは無かったので、畳の上に裸で置いた状態で、計算開始だ。

体育の日を入れた三連休で、約28時間の大仕事だった。使った CPU は、Phenom X4 9350e 2GHz だ。Ruby で、その上、マルチコアに対応したプログラムではないので、3つのCPUコアが余った状態の、とっても無駄な時間を使った計算だったが、とにかく終わった。

その結果を最初のグラフに示す。
地球から木星までの約2年半のホーマン軌道で、その間に最も接近する小惑星は、なんと13万キロも離れている。13万キロまで小惑星が近づくのは一瞬だけで、ほとんどの期間は、200万キロ以上離れているのだ。

計算に使ったプログラムは、asteroidInterval3.rb だ。
例によって、JPL Solar System DynamicsEphemeridesにある情報で計算してみた。
JPL のAsteroidsSmall-Body Orbital Elements のページにある
ELEMENTS.NUMBR (圧縮版はELEMENTS.NUMBR.gz)と
ELEMENTS.UNNUM(圧縮版はELEMENTS.UNNUM.gz)が小惑星の軌道データを使う。

asteroidInterval3.rb と ELEMENTS.NUMBR と ELEMENTS.UNNUM を同じディレクトリに入れて、以下のように、コマンドすると、計算できる。
$ ./asteroidInterval3.rb 1000 >asteroidInterval.dat
なお、途中の「1000」は、小惑星の数を減らし、ランダムに 1000個選び出して計算するという意味だ。いきなり、軌道データの判っている全ての小惑星53万個で計算すると、大変な時間がかかるからねえ。「1000」の部分を「0」にすると、全ての小惑星の軌道データを用いて計算する。

E2222左に、選び出した小惑星の数と、最も小惑星に近づいた時の距離および計算時間との関係を示した。
このグラフは対数-対数の表示にしている。
グラフ上、赤い線は最も小惑星に近づいた時の距離だが、これは、ランダム性があるが、それでも細い青い線で示した小惑星の数の平方根の逆数に比例していることが判る。つまり、小惑星の個数が4倍になると、最も小惑星に近づいた時の距離が2分の1になるということ。
グラフ上、緑色の線は、計算時間だが、これは、紫色の細い線の小惑星の数の 1.5乗と比例している。

なぜ、こんなことを示したかと言うと、現在、小惑星は53万個見つかっているが、それだけじゃ無いってこと。あくまでも地球から観測できる小惑星が53万個だから、もっと小さくて見えにくい小惑星が、たくさんあるに違いない。
だから、本当は、どのくらい小惑星があったら、どのくらい接近するかの見当を付けておきたかったのだ。

で、接近する距離が、小惑星の数の平方根の逆数に比例するなら、もし、小惑星が見つかってない分も合わせて、現在の53万個の100倍の5300万個あったとすると、最短距離は、1.3万キロだって判るわけ。もし、1万倍の53億個なら、1300キロなわけだ。

小惑星って、53万個よりは多いが、100倍も無いような気がする。少なくとも1万倍って事はないだろう。
つまり、そんだけ小惑星が増えても、一番近付いたときでも1000キロも離れているってこと。

う~~ん、宇宙戦艦ヤマトの次から次に来る小惑星を避けてジグザグに通る操縦は必要ないんだなあ。

ところで、今回、公開したプログラムだけど、地球~木星間のホーマン軌道だけじゃなくて、汎用的にどんな軌道の宇宙船でも、小惑星まで最接近する距離を計算できるようにしてある。
その使い方については、いずれまた。

| | Comments (3) | TrackBack (0)

« September 2010 | Main | December 2010 »