惑星の運動のシミュレーション

惑星の運動のシミュレーション

先ほどの丸が移動してゆくプログラムをもとに、太陽系を動く惑星のような、重力によって速度が変わる系のシミュレーションをしてみましょう。

さきほどのプログラムを一部変更し、○の位置を表す変数(x, y)に加え、○の速度を表す変数(vx, vy)を定義します。

float x, y;
float vx, vy;

void setup(){
  size(400, 400);
  smooth();  //円の縁を滑らかにする
  stroke(255);
  background(0);
  x = 0;
  y = 0;
  vx = 1;
  vy = 1;
}

void draw(){
  ellipse(x, y, 5, 5);
  x = x + vx;
  y = y + vy;
}

このままだと、全く前と同じく、○が右下に向かって移動していきます。(等速直線運動)
これから、画面の中心(200,200)に太陽のような重力中心があり、○が重力によってその周りを回るような運動のシミュレーションを考えます。

vxやvyは○の速度を表していますが、その単位は m/s などではなく、
pixel/frame となっていることに注意してください。
pixelは画面上の画素を単位とした長さ、
frameはアニメーションの一コマ分(draw1回分)の時間です。

さて、○を宇宙空間を漂う物体とみなして、これに力が加わる時の運動について考えてみましょう。
物体の運動と力は、次の運動方程式で記述されます。

  F = ma

a は加速度で、単位時間当たりの速度の変化を表しています。
また、2物質間に働く重力は

  \displaystyle F = -G\frac{M m}{r^2}

以上2式から、

  \displaystyle a = -G\frac{M}{r^2}

であることがわかります。

r2 は べき乗を計算する命令 pow(a,b)を使い、
(pixel単位)で次のように求められます。
計算結果は rs という変数に入れておきましょう。
下の命令は、関数drawの中に入れます。

float rs;
rs = pow(x-200,2)+pow(y-200,2);

これで物体への加速度 a を求めることができます。
G と M を決めておかなくてはなりませんが、これは GM = 100 としておいてください。

float a;
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) として円運動をすれば正解です。

void draw(){
  ellipse(x, y, 5, 5);

  float rs, r, a, ax, ay;
  //ここでax, ayを求める

  vx = vx + ax;
  vy = vy + ay;

  x = x + vx;
  y = y + vy;
}

演習

  • 上のプログラムの //ここでax, ay を求める のところに正しい命令を入れて、○を円運動させてください。
  • 初期速度を少し変えて、惑星が楕円運動をすることなどを確かめてください。

サンプル