新しく発言する EXIT インデックスへ
隕石の起動

  隕石の起動 けい 2004/11/12 02:56:35 
  けいさんがどの程度この問題に詳しくて、ま... 青木太一 2004/11/12 10:50:53 
   └返事遅くなって申し訳ありません。 けい 2004/11/15 02:29:00 

  隕石の起動 けい 2004/11/12 02:56:35  ツリーへ

隕石の起動 返事を書く
けい 2004/11/12 02:56:35
十進BASICで月と地球の間を通過する隕石の軌道を計算したいんですが、どのようなプログラムにすればよいのかわかりません。教えてください。

  けいさんがどの程度この問題に詳しくて、ま... 青木太一 2004/11/12 10:50:53  ツリーへ

Re: 隕石の起動 返事を書く
青木太一 2004/11/12 10:50:53
けいさんがどの程度この問題に詳しくて、また、このプログラムを何に使いたいかによって答えは変わってくると思います。
が、とりあえずてきとうに書いてみます。

軌道計算に限らず、運動方程式に従って物体がどう運動するかは、
たとえば以下のような方針で数値的に得られます。

ここからループ
 1.各物体にかかる力を計算
 2.力を質量で割って加速度を得る
 3.今までの速度に、その加速度(×時間刻み)を加算
 4.今までの位置に、その速度(×時間刻み)を加算
ここまでループ

ただし、この方法はあまり精度は高くありません。
以下のような用途なら十分でしょうが、
・なんとなく重力が働く物体間の挙動が見てみたい
・重力で運動している天体っぽい動画が見たい

以下のように高い精度が求められる用途には不十分だと思います。
・天体観測に使いたい
・地球に衝突しそうな小惑星の軌道計算したい
・太陽系が重力的に安定か確かめたい

より高い精度が必要な場合、たとえば3.と4.のところをもっと複雑に行うことで
精度を上げることができるようです(ルンゲ・クッタ法など)が、私はここに書き込めるほど
詳しく理解していません。(講義は受けたはずなんだけど...)

以下に、精度を重視しない方針で書いたコードを書きます。
実際の地球、月、隕石の位置、初速度、質量は面倒なので調べていません。ごめんなさい。
また、物体間が極端に近づくと精度が悪いです。(保存されるはずの力学的エネルギーが極端に変動しました)
けいさんの参考になったら幸いです。

!以下コード

randomize
set window 0,1,0,1
LET n=3
dim m(1 to n)!質量
dim x(1 to n),y(1 to n)!位置
dim vx(1 to n),vy(1 to n)!速度

!質量(てきとう)
LET m(1)=10
LET m(2)=3
LET m(3)=1.5

!初期位置、初期速度(てきとう)
for i= 1 to n
LET x(i)=rnd
LET y(i)=rnd
LET vx(i)=rnd-.5
LET vy(i)=rnd-.5
next i

LET dt=0.01 !時間刻み
LET G=6.67259e-11 !万有引力定数(他の値が適当なのであまり意味はない)
do

!sが絡むコメントアウトしてある部分は、
!力学的エネルギーが保存しているかどうか検査している。
! LET s=0 

!物体i,j間の万有引力(による速度変化(加速)を加算)
for i=1 to n
for j=i+1 to n
LET d=sqr((x(i)-x(j))^2 + (y(i)-y(j))^2) !i,j間の距離
! LET s=s-G*m(i)*m(j)/d
LET f=G*m(i)*m(j)/(d^2) !i,j間に働く引力 「1.各物体にかかる力を計算」

LET vx(i)=vx(i)+ f/m(i)*(x(j)-x(i))/d*dt!「2.力を質量で割って加速度を得る」
LET vx(j)=vx(j)+ f/m(j)*(x(i)-x(j))/d*dt!「3.今までの速度に、その加速度(×時間刻み)を加算」
LET vy(i)=vy(i)+ f/m(i)*(y(j)-y(i))/d*dt!を同時に行っている
LET vy(j)=vy(j)+ f/m(j)*(y(i)-y(j))/d*dt
next j
next i

!「4.今までの位置に、その速度(×時間刻み)を加算」
for i=1 to n
LET x(i)=x(i)+vx(i)*dt
LET y(i)=y(i)+vy(i)*dt

plot points:x(i),y(i)
next i

set draw mode explicit
set draw mode hidden
clear

!for i=1 to n
! LET s=s+1/2*m(i)*(vx(i)^2+vy(i)^2)
!next i
!print s

loop

END

   └返事遅くなって申し訳ありません。 けい 2004/11/15 02:29:00  ツリーへ

Re: けいさんがどの程度この問題に詳しくて、ま... 返事を書く
けい 2004/11/15 02:29:00
返事遅くなって申し訳ありません。
わざわざ有難うございます。
C言語は勉強中なんですが、Basicは始めたばかりなもので…。

//・なんとなく重力が働く物体間の挙動が見てみたい
これに当てはまります。とりあえず概要を把握したいんです。説明不足で申し訳ありません。
十分参考になりました。頑張ってみます。


インデックスへ EXIT
新規発言を反映させるにはブラウザの更新ボタンを押してください。