惑星の運動のシミュレーション
先ほどの丸が移動してゆくプログラムをもとに、太陽系を動く惑星のような、重力によって速度が変わる系のシミュレーションをしてみましょう。
さきほどのプログラムを一部変更し、○の位置を表す変数(x, y)に加え、○の速度を表す変数(vx, vy)を定義します。
x = 0.0 y = 0.0 vx = 1.0 vy = 1.0 def setup(): size(400, 400) smooth() #円の縁を滑らかにする stroke(255) background(0) def draw(): global x, y, vx, vy ellipse(x, y, 5, 5) x = x + vx y = y + vy
このままだと、全く前と同じく、○が右下に向かって移動していきます。(等速直線運動)
これから、画面の中心(200,200)に太陽のような重力中心があり、○が重力によってその周りを回るような運動のシミュレーションを考えます。
vxやvyは○の速度を表していますが、その単位は m/s などではなく、
pixel/frame となっていることに注意してください。
pixelは画面上の画素を単位とした長さ、
frameはアニメーションの一コマ分(draw1回分)の時間です。
さて、○を宇宙空間を漂う物体とみなして、これに力が加わる時の運動について考えてみましょう。
物体の運動と力は、次の運動方程式で記述されます。
a は加速度で、単位時間当たりの速度の変化を表しています。
また、2物質間に働く重力は
以上2式から、
であることがわかります。
r2 は べき乗を計算する命令 pow(a,b)を使い、
(pixel単位)で次のように求められます。
計算結果は rs という変数に入れておきましょう。
下の命令は、関数drawの中に入れます。
rs = pow(x-200,2)+pow(y-200,2);
これで物体への加速度 a を求めることができます。
G と M を決めておかなくてはなりませんが、これは GM = 100 としておいてください。
a = -100/rs;
a はベクトルなので、x 方向と y 方向に分けておかなくてはなりません。
例えば a の x成分(ax) は a に r と x の比(x/r)をかけることで求めることができます。
(ここでの x はx軸方向の○と中心の間の距離)
r は sqrt(rs),
x は (x-200)
なので、ax は
ax = a * x / r;
ですね。符号が中心に向かう方向になるように 注意しなくてはなりません。
あとは次のように、draw中に 加速度を加える処理をつければ惑星の運動がシミュレートできます。
x, y の初期値を (200, 100)、
vx, vyの初期値を (1, 0) として円運動をすれば正解です。
def draw(): ellipse(x, y, 5, 5) #ここでax, ayを求める vx = vx + ax vy = vy + ay x = x + vx y = y + vy
演習
- 上のプログラムの #ここでax, ay を求める のところに正しい命令を入れて、○を円運動させてください。
- 初期速度を少し変えて、惑星が楕円運動をすることなどを確かめてください。