/**************************************************************
* ボールの衝突（完全弾性衝突）シミュレーション（２） *
* BALL01 ver 1.1 (1999/04/03) *
* Copyright(c) 1999 Yamada,K (1999.02.20) *
* *
* エネルギーも表示する *
* 個々のボールのエネルギー変化は、対応する色で表示し *
* 全エネルギーは黄色で表示する *
* スペースキーで描画を終了する *
***************************************************************
*/
#include <iostream.h>
#include <time.h>
#include "glibw32.h"// グラフィックス・クラスライブラリ GLIBW32
const int size_x = 640, size_y = 200, size_y2 = 200;
const double ZERO = 1.E-20; // ゼロと見なせる
const double PI = atan(1.0)*4.0; // 円周率
const int time_max = 1000;
// ボール・クラス
class Ball{
public:
Ball(); // デフォルト・コンストラクタ
~Ball() {--count;} // デストラクタ
void draw(GRAPH *); // ボールの描画
void non_coll(); // 衝突しないボールの次の位置を計算
void collision_w(); // 壁との衝突を調べ次の位置を計算
void collision(Ball &); // ボール同士の衝突を調べ次の位置を計算
int getcolor() const { return color;}
double getenergy() const { return m * v * v /2.0;}// エネルギーを計算
private:
static int count; // ボールの個数
int x_size; // ウィンドウサイズ（ｘ方向）
int y_size; // ウィンドウサイズ（ｙ方向）
int coll; // 次の１秒間に他のボールと衝突するかどうか
int color; // ボールの色
double x; // ボールのｘ座標
double y; // ボールのｙ座標
double v; // ボールの速度
double angle; // ボールの運動方向（度）
double r; // ボールの半径
double m; // ボールの質量
double x_next; // ボールの次のｘ座標
double y_next; // ボールの次のｙ座標
double dt; // 時間間隔（次のステップまでの時間）
int inside(double, double) const; // 領域の判定（0:外にある、1:内）
};
int Ball::count = 0; // 静的データメンバーの初期化
// デフォルト・コンストラクタ
Ball::Ball(): x_size(size_x), y_size(size_y), coll(0), dt(1.0)
{
++count;
if(count == 1){
time_t t;
srand((unsigned) time(&t)); //乱数の初期化
}
x =( rand() % (x_size/2) ) + x_size/4;// ｘ座標の初期値
y =( rand() % (y_size/2) ) + y_size/4;// ｙ座標の初期値
angle = rand() % 360; // 角度の初期値
v = rand() % 11 + 10; // 速度の初期値（１０〜２０）[m/s]
r = rand() % 6 + 5; // 半径（５〜１０）[m]
m = r * r * r * 4./3.* PI * 0.0001; // 質量 [Kg]
color = count % 15 + 1;
non_coll();
}
// ボールの描画
void Ball::draw(GRAPH *g)
{
// 更新
coll = 0;
dt = 1.0;
x = x_next;
y = y_next;
g->fcircle(x, y, r, color);// ボールを描く
}
// 衝突しない場合の次の位置を計算する
void Ball::non_coll()
{
x_next = x + v*cos(PI*angle/180.0)*dt;
y_next = y + v*sin(PI*angle/180.0)*dt;
}
// 壁との衝突
void Ball::collision_w()
{
double x_min = r, x_max = x_size - r;
double y_min = r, y_max = y_size - r;
// ２回衝突（壁の角で２回衝突する場合）
if( ( x_next <= x_min ) && ( y_next <= y_min ) ){
double tx = ( x - x_min) / (x - x_next);
double ty = ( y - y_min) / (y - y_next);
x_next = x_min - dt*(1.0-tx)*v*cos(PI*angle/180.0);
y_next = y_min - dt*(1.0-ty)*v*sin(PI*angle/180.0);
angle = 180.0 + angle;
}
if( ( x_next >= x_max ) && ( y_next <= y_min ) ){
double tx = ( x_max - x) / (x_next - x);
double ty = ( y - y_min) / (y - y_next);
x_next = x_max - dt*(1.0-tx)*v*cos(PI*angle/180.0);
y_next = y_min - dt*(1.0-ty)*v*sin(PI*angle/180.0);
angle = 180.0 + angle;
}
if( ( x_next <= x_min ) && ( y_next >= y_max ) ){
double tx = ( x - x_min) / (x - x_next);
double ty = ( y_max - y) / (y_next - y);
x_next = x_min - dt*(1.0-tx)*v*cos(PI*angle/180.0);
y_next = y_max - dt*(1.0-ty)*v*sin(PI*angle/180.0);
angle = 180.0 + angle;
}
if( ( x_next >= x_max ) && ( y_next >= y_max ) ){
double tx = ( x_max - x) / (x_next - x);
double ty = ( y_max - y) / (y_next - y);
x_next = x_max - dt*(1.0-tx)*v*cos(PI*angle/180.0);
y_next = y_max - dt*(1.0-ty)*v*sin(PI*angle/180.0);
angle = 180.0 + angle;
}
// １回衝突
if( ( x_next <= x_min ) && ( x >= x_min ) ){
double t = ( x - x_min) / (x - x_next);
angle = 180.0 - angle;
x_next = x_min + dt*(1.0-t)*v*cos(PI*angle/180.0);
}
if( ( x_next >= x_max ) && ( x <= x_max ) ){
double t = ( x_max -x ) / (x_next - x);
angle = 180.0 - angle;
x_next = x_max + dt*(1.0-t)*v*cos(PI*angle/180.0);
}
if( ( y_next <= y_min ) && ( y >= y_min ) ){
double t = ( y - y_min ) / (y - y_next);
angle = - angle;
y_next = y_min + dt*(1.0-t)*v*sin(PI*angle/180.0);
}
if( ( y_next >= y_max ) && ( y <= y_max ) ){
double t = ( y_max - y) / (y_next - y);
angle = - angle;
y_next = y_max + dt*(1.0-t)*v*sin(PI*angle/180.0);
}
}
// ボール同士の衝突
void Ball::collision(Ball &B)
{
if( this == &B) return; // 相手は自分自身
if( coll == 1 ) return; // 他のボールと衝突
// 重心と重心速度
double X1 = ( m * x + B.m * B.x ) / ( m + B.m);
double Y1 = ( m * y + B.m * B.y ) / ( m + B.m);
double X2 = ( m * x_next + B.m * B.x_next ) / ( m + B.m);
double Y2 = ( m * y_next + B.m * B.y_next ) / ( m + B.m);
double V_x = ( m * v * cos(PI*angle/180.0) +
B.m * B.v * cos(PI*B.angle/180.0) ) / ( m + B.m);
double V_y = ( m * v * sin(PI*angle/180.0) +
B.m * B.v * sin(PI*B.angle/180.0) ) / ( m + B.m);
if( (V_x * V_x + V_y * V_y) < ZERO ) return;
// 重心座標系
double v_p = (v*cos(PI*angle/180.0) - V_x)*(v*cos(PI*angle/180.0) - V_x)
+ (v*sin(PI*angle/180.0) - V_y)*(v*sin(PI*angle/180.0) - V_y);
v_p = sqrt(v_p); // 速さ
double x1 = x -X1, y1 = y - Y1, x2 = x_next - X2, y2 = y_next - Y2;
// 軌道直線の係数
double a = y2- y1, b = x2 - x1, c = x2*y1 - x1*y2;
double d = fabs(c)/sqrt(a*a + b*b); // 衝突径数
double R = B.m / ( m + B.m ) * (r + B.r) ;// 衝突距離
// 衝突条件（１）
if( d >= R) return;// 衝突しない
double theta, phi, psi;
if(fabs(b) < ZERO) theta = ( a*b > 0.0) ? 0.5*PI : -0.5*PI;
else theta = atan2(a,b);
phi = acos(d/R);
double x_c, y_c, x_d, y_d;// 衝突位置
if(fabs(b) < ZERO && x1 < 0 && y1 > 0 ) psi = PI - phi;
else if(fabs(b) < ZERO && x1 > 0 && y1 > 0 ) psi = phi;
else if(fabs(b) < ZERO && x1 < 0 && y1 < 0 ) psi = PI + phi;
else if(fabs(b) < ZERO && x1 > 0 && y1 < 0 ) psi = - phi;
else if(x1 < x2 && b*c > 0 ) psi = theta + phi + 0.5 * PI;
else if(x1 > x2 && b*c > 0 ) psi = theta - phi - 0.5 * PI;
else if(x1 < x2 && b*c < 0 ) psi = theta - phi - 0.5 * PI;
else if(x1 > x2 && b*c < 0 ) psi = theta + phi + 0.5 * PI;
else return;
psi = fmod(psi,360.0);
// 衝突位置
x_c = r*cos(psi);
y_c = r*sin(psi);
// 衝突時刻
double t = (fabs(b) < ZERO) ? ( y_c - y1 )/ a : ( x_c - x1 )/ b;
// 衝突条件（２）衝突は次のステップ（１秒後）の前に起こる
if(t<0.0 || t>1.0) return;
// 衝突位置（実験室系）
double xc = x + t*v*cos(PI*angle/180.0);
double yc = y + t*v*sin(PI*angle/180.0);
double xd = B.x + t*v*cos(PI*B.angle/180.0);
double yd = B.y + t*v*sin(PI*B.angle/180.0);
// 衝突条件（３） 衝突は領域内で起こる
if( inside(xc,yc) == 0 || B.inside(xd,yd) == 0 ) return;
// 衝突後の速度の計算
double v_x, v_y;
v_x = v_p * cos(2.0 * psi - theta + PI) + V_x;
v_y = v_p * sin(2.0 * psi - theta + PI) + V_y;
v = sqrt(v_x * v_x + v_y * v_y);
if(fabs(v_x) < ZERO) angle = (v_x * v_y > 0.0) ? 90.0 : -90.0;
else angle = (180.0 /PI ) * atan2(v_y , v_x);
// 現在位置を衝突位置に更新
x=xc;
y=yc;
dt = 1.0-t;
non_coll();
v_x = -( m / B.m ) * v_p * cos(2.0 * psi - theta + PI) + V_x;
v_y = -( m / B.m ) * v_p * sin(2.0 * psi - theta + PI) + V_y;
B.v = sqrt(v_x * v_x + v_y * v_y);
if(fabs(v_x) < ZERO) angle = (v_x * v_y > 0.0) ? 90.0 : -90.0;
else B.angle = (180.0 /PI ) * atan2(v_y , v_x);
// 現在位置を衝突位置に更新
B.x=xd;
B.y=yd;
B.dt = 1.0-t;
B.non_coll();
coll = 1;
B.coll = 1;
return;
}
// 領域の判定
int Ball::inside(double a, double b) const
{
double x_min = r, x_max = x_size - r;
double y_min = r, y_max = y_size - r;
if(a <= x_min || a >= x_max || b <= y_min || b >= y_max ) return 0;
else return 1;
}
main()
{
int no;
cout << "ボールの衝突運動のシミュレーションです。\n"
<< "スペースキーで描画を終了します。\n\n"
<< "ボールの個数を入力してください（２、３個がお薦め）：";
cin >> no;
Ball *ptr = new Ball[no]; // ボールの生成
double *e = new double[no];
double energy = 0.0;
for(int i = 0; i < no; ++i)
energy += (ptr+i)->getenergy();
cout << "開始時の全エネルギー： " << energy << " [N・m]" << endl;
// グラフィックスの初期化
ginit(size_x, size_y + size_y2);
GRAPH *g = new GRAPH[2];
// g[0]:ボールの描画、g[1]:エネルギーのグラフ
double axis_x = time_max / 10.0;
double axis_y = energy /10.0;
g[0].window(0, 0, size_x, size_y);
g[0].view(0, 0, size_x -1, size_y -1);
g[1].window(0, - 0.1 * energy , time_max, energy * 1.1);
g[1].view(0, size_y, size_x -1, size_y + size_y2 -1);
g[1].axis(axis_x, axis_y);
int count = 0; // 時刻のカウンタ
while( vkey() != VK_SPACE ){// スペースキーが押されるまで
g[0].cls(); // 画面消去
int i,j;
// まずは衝突しないとして全ボールの次の位置を計算
for(i = 0; i < no; ++i)
(ptr+i)->non_coll();
for(i = 0; i < no; ++i){
// ボール同士の衝突を調べて、次の位置を修正
for(j = i+1; j < no; ++j)
(ptr+i)->collision(*(ptr+j));
// 壁との衝突を調べて、次の位置を修正
(ptr+i)->collision_w();
}
// ボールの描画
for(i = 0; i < no; ++i)
(ptr + i) -> draw(&g[0]);
Sleep(20); // 実行を中断する時間をミリ秒単位で指定するWIN32 API関数
// エネルギーをプロットする
double energy_old = energy;
energy = 0.0;
for(i = 0; i < no; ++i){
double e_old = e[i];
g[1].moveto(count, e_old);
e[i] = (ptr+i)->getenergy();
int c = (ptr+i)->getcolor();
g[1].lineto(count, e[i], c);
energy += e[i];
}
g[1].moveto(count-1, energy_old);
g[1].lineto(count, energy, YELLOW);// 全エネルギー
++count;
count %= time_max;
if(count == 0){
g[1].cls();
g[1].axis(axis_x, axis_y);
}
}
cout << "終了時の全エネルギー： " << energy << " [N・m]" << endl;
delete [] g;
gend(); // グラフィックスの終了
delete [] e;
delete [] ptr; // ボールの削除
return 0;
}
| 戻る |
Copyright(c) Yamada,K