こんにちは。引き続き、processingの話題です。
前回、「クラス」の話をしました。今回は「衝突」の話をします。
「衝突」は、「ものとものとのぶつかり」のことですが、これをプログラム上で表してみたいと思います。
そのためには「運動量保存則」と「エネルギー保存則」が必要になります。
まずは簡単のために一直線上に丸1と丸2があるときを考えます。
一方の丸を質量$m_1$、速度$v_1$、もう一方の丸を質量$m_2$、速度$v_2$として、衝突後にそれぞれの速度が$v_1'$、$v_2'$になるとすると、
運動量保存則は、
$$
m_1 v_1 + m_2 v_2 = m_1 v_1' + m_2 v_2'
$$
エネルギー保存則は、
$$
\frac{1}{2} m_1 {v_1}^{2} + \frac{1}{2} m_2 {v_2}^{2} = \frac{1}{2} m_1 {v_1}'^{2} + \frac{1}{2} m_2 {v_2}'^{2}
$$
この二つの式から、衝突後のそれぞれの丸の速度を求めると、
$$
\begin{align}
v_1' = \frac{m_1 v_1 + 2 m_2 v_2 - m_2 v_1}{m_1 + m_2}\\
v_2' = \frac{2 m_1 v_1 + m_2 v_2 - m_1 v_2}{m_1 + m_2}
\end{align}
$$
となります。
これを導入する前に、二つの丸を描くプログラムを書いておきます。
float x1, x2; float v1, v2; float m1, m2; void setup(){ size(500, 500); x1 = 50; x2 = width/2; v1 = 3; v2 = 0; m1 = 1; m2 = 1; } void draw(){ background(0); ellipse(x1, 400, 50, 50); ellipse(x2, 400, 50, 50); x1 += v1; x2 += v2; }
どちらの丸も質量は同じとし、丸2は静止しているとしておきました。
これを実行しても、丸1が丸2を通り抜けてしまい、衝突しません。
ここで、先ほど求めた$v_1'$、$v_2'$の式を組み込んで、衝突できるようにします。
（プログラム中では、$v_1'$、$v_2'$はそれぞれv1a、v2aと表しました。）
float x1, x2; float v1, v2; float v1a, v2a; float m1, m2; float distance; void setup(){ size(500, 500); x1 = 50; x2 = width/2; v1 = 3; v2 = 0; m1 = 1; m2 = 1; } void draw(){ background(0); distance = abs(x1-x2); if(distance < 50){ //丸同士が接しているとき、衝突の処理をする v1a = (m1*v1+2*m2*v2-m2*v1)/(m1+m2); v2a = (2*m1*v1+m2*v2-m1*v2)/(m1+m2); v1 = v1a; v2 = v2a; } x1 += v1; x2 += v2; ellipse(x1, 400, 50, 50); ellipse(x2, 400, 50, 50); }
さらに、壁を追加します。
float x1, x2; float v1, v2; float v1a, v2a; float m1, m2; float distance; void setup(){ size(500, 500); x1 = 50; x2 = width/2; v1 = 3; v2 = 0; m1 = 1; m2 = 1; } void draw(){ background(0); distance = abs(x1-x2); if(distance < 50){ v1a = (m1*v1+2*m2*v2-m2*v1)/(m1+m2); v2a = (2*m1*v1+m2*v2-m1*v2)/(m1+m2); v1 = v1a; v2 = v2a; } x1 += v1; x2 += v2; ellipse(x1, 400, 50, 50); ellipse(x2, 400, 50, 50); if(x1 < 25){ x1 = 25; v1 *= -1; } if(x2 > width-25){ x2 = width-25; v2 *= -1; } }
壁を追加したことで、丸の質量が同じのためニュートンのゆりかごのような動きをするようになりました。
ではここに、以前に説明したような手順でクラスと配列を導入してみます。
（細かい手順は省き、最終的な結果を示します。）
int n; Circle[] circles; void setup(){ size(500, 500); n = 2; circles = new Circle[n]; for(int i = 0; i < n; i ++){ float xi = 50+100*i; float vi = random(-3, 3); circles[i] = new Circle(xi, 400, 1, vi); } } void draw(){ background(0); for(int i = 0; i < n; i ++){ circles[i].display(); circles[i].move(); circles[i].wall(); circles[i].collide(); } } class Circle{ float x, y, m, v; Circle(float _x, float _y, float _m, float _v){ x = _x; y = _y; m = _m; v = _v; } void display(){ ellipse(x, y, 50, 50); } void collide(){ for(int i = 0; i < n; i ++){ float distance = abs(circles[i].x-x); if(distance < 50){ float v1a, v2a; v1a = (circles[i].m*circles[i].v+2*m*v-m*circles[i].v)/(circles[i].m+m); v2a = (2*circles[i].m*circles[i].v+m*v-circles[i].m*v)/(circles[i].m+m); circles[i].v = v1a; v = v2a; } } } void move(){ x += v; } void wall(){ if(x < 25){ x = 25; v *= -1; } if(x > width-25){ x = width-25; v *= -1; } } }
これでクラスと配列を導入できました。
クラスの中の「collide」の部分は、自分の位置(x, 400)と他の粒子の位置(circles[i].x, 400)から丸同士の距離を計算し、その距離が50より小さければ衝突の処理を行う、という形になっています。
これまでは一直線上にある場合を考えましたが、同一平面上にある場合に変えてみます。
画像で表すとこんな感じです。
二次元になるため、y軸方向に関する情報も加え、さらに1回の衝突で複数回処理が行われることを防ぐため、衝突の処理を行う条件に「距離が縮んでいる」という条件を追加します。
int n = 2; Circle[] circles = new Circle[n]; void setup(){ size(500, 500); for(int i = 0; i < n; i ++){ circles[i] = new Circle(i); } } void draw(){ background(0); for(int i = 0; i < n; i ++){ circles[i].display(); circles[i].move(); circles[i].wall(); circles[i].collide(); } } class Circle{ int id; float x, y, m, d, vx, vy, elast; float[] odistance = new float[n]; Circle(int _id){ id = _id; d = 50; x = id%10*(width/10)+d/2; y = id/10*(height/10)+d/2; m = 1; vx = random(-3, 3); vy = random(-3, 3); for(int i = 0; i < n; i ++){ odistance[i] = width+height; } } void display(){ ellipse(x, y, d, d); } void move(){ x += vx; y += vy; } void collide() { for (int i = id+1; i < n; i ++){ float dx = circles[i].x - x; float dy = circles[i].y - y; float distance = sqrt(dx*dx + dy*dy); if (distance < d && distance < odistance[i]){ //距離が縮んでいるという条件を追加 float angle = atan2(dy, dx); //2つの丸を結んだ線のなす角を計算 float v1h = vx*cos(angle)+vy*sin(angle); //丸の速度を衝突方向に水平な成分と垂直な成分に分解 float v1v = vx*sin(angle)-vy*cos(angle); float v2h = circles[i].vx*cos(angle)+circles[i].vy*sin(angle); float v2v = circles[i].vx*sin(angle)-circles[i].vy*cos(angle); float v1ha = (m*v1h+2*circles[i].m*v2h-circles[i].m*v1h)/(m+circles[i].m); //水平方向は衝突を計算 float v2ha = (2*m*v1h+circles[i].m*v2h-m*v2h)/(m+circles[i].m); vx = v1ha*cos(angle)+v1v*sin(angle); //x軸方向、y軸方向に書き直す vy = v1ha*sin(angle)-v1v*cos(angle); circles[i].vx = v2ha*cos(angle)+v2v*sin(angle); circles[i].vy = v2ha*sin(angle)-v2v*cos(angle); } odistance[i] = distance; } } void wall(){ if(x < d/2){ x = d/2; vx *= -1; } if(x > width-d/2){ x = width-d/2; vx *= -1; } if(y < d/2){ y = d/2; vy *= -1; } if(y > height-d/2){ y = height-d/2; vy *= -1; } } }
これで二次元に拡張することができました。
さて、なぜ二次元に拡張したかというと、最終的には単原子分子の運動のシミュレータのようなものをつくりたいからです。
このままでも、分子間に相互作用が無いと考えた場合の運動になりますが、完成形は分子間力も加えたものにしたいと考えています。
というわけで、processingで物体同士の衝突の処理をできるようにしてみました。
区切りがよいので今回はここまでにします。また次回。
