octaveでピンポン その2

こんにちは。
以前、ピンポンのような動きをシミュレートするコードを書きましたが、今回はそれに速度を加えたものを作ってみました。さらに、玉を二つにしました。

アルゴリズム

今回は2つの球を扱い、さらに速度を考えるので速度交換の計算が必要になります。
完全弾性衝突だと、衝突した時に衝突面と垂直な成分のみが交換します。(図のような感じです)

また、角度の扱いに苦労したので、0から2πの範囲で計算するようにしました。

コード

▶クリックで展開

#include <stdio.h>
#include <math.h>

double angle(double theta, double x, double y){
    if (x*y >= 0) {
        if (x < 0) {
            theta += M_PI;
        }
    } else {
        if (x < 0) {
            theta = M_PI - theta;
        } else {
            theta = 2*M_PI - theta;
        }
    }
    return theta;
}

int main() {
    double x1,y1,x2,y2;
    double vx1,vy1,vx2,vy2;
    double t = 0;
    double dt = 0.001;
    double R = 3.0;
    double r = 0.14;
    double distanceFromOrigin1;
    double distanceFromOrigin2;
    double distanceOfStones;
    double theta,theta1,theta2,v1,v2,v1_ver,v2_ver,v1_pal,v2_pal,v1_verx,v1_very,v1_palx,v1_paly,v2_verx,v2_very,v2_palx,v2_paly;
    double dif_x,dif_y;
    int count = 0;
    
    scanf("%lf",&x1);
    scanf("%lf",&y1);
    scanf("%lf",&vx1);
    scanf("%lf",&vy1);
    scanf("%lf",&x2);
    scanf("%lf",&y2);
    scanf("%lf",&vx2);
    scanf("%lf",&vy2);
    
    for (int i=0;i<=2000;i++) {
        
        distanceFromOrigin1 = x1*x1+y1*y1;
        distanceFromOrigin2 = x2*x2+y2*y2;
        distanceOfStones = pow(x1-x2, 2)+pow(y1-y2, 2);
        
        if (pow(R-r, 2)-0.05 <= distanceFromOrigin1 && distanceFromOrigin1 <= pow(R-r, 2)+0.05) {
            theta = fabs(atan(y1/x1));
            theta1 = fabs(atan(vy1/vx1));
            v1 = sqrt(vx1*vx1 + vy1*vy1);
            v1_ver = v1 * cos(fabs(angle(theta1, vx1, vy1) - angle(theta, x1, y1)));
            v1_pal = v1 * sin(fabs(angle(theta1, vx1, vy1) - angle(theta, x1, y1)));
            
            v1_verx = v1_ver * cos(angle(theta, x1, y1));
            v1_very = v1_ver * sin(angle(theta, x1, y1));
            if (angle(theta1, vx1, vy1) >= angle(theta, x1, y1)) {
                v1_palx = v1_pal * cos(angle(theta, x1, y1) + M_PI/2);
                v1_paly = v1_pal * sin(angle(theta, x1, y1) + M_PI/2);
            } else {
                v1_palx = v1_pal * cos(angle(theta, x1, y1) - M_PI/2);
                v1_paly = v1_pal * sin(angle(theta, x1, y1) - M_PI/2);
            }
            
            vx1 = -v1_verx + v1_palx;
            vy1 = -v1_very + v1_paly;
            printf("1,%lf, %lf,%lf,%lf,%lf\n",t,x1,y1,x2,y2);
        }
        
        if (pow(R-r, 2)-0.05 <= distanceFromOrigin2 && distanceFromOrigin2 <= pow(R-r, 2)+0.05) {
            theta = fabs(atan(y2/x2));
            theta2 = fabs(atan(vy2/vx2));
            v2 = sqrt(vx2*vx2 + vy2*vy2);
            v2_ver = v2 * cos(fabs(angle(theta2, vx2, vy2) - angle(theta, x2, y2)));
            v2_pal = v2 * sin(fabs(angle(theta2, vx2, vy2) - angle(theta, x2, y2)));
            v2_verx = v2_ver * cos(angle(theta, x2, y2));
            v2_very = v2_ver * sin(angle(theta, x2, y2));
            if (angle(theta2, vx2, vy2) >= angle(theta, x2, y2)) {
                v2_palx = v2_pal * cos(angle(theta, x2, y2) + M_PI/2);
                v2_paly = v2_pal * sin(angle(theta, x2, y2) + M_PI/2);
            } else {
                v2_palx = v2_pal * cos(angle(theta, x2, y2) - M_PI/2);
                v2_paly = v2_pal * sin(angle(theta, x2, y2) - M_PI/2);
            }
            vx2 = -v2_verx + v2_palx;
            vy2 = -v2_very + v2_paly;
            printf("1,%lf, %lf,%lf,%lf,%lf\n",t,x1,y1,x2,y2);
        }
        
        if (pow(2*r, 2)-0.003 <= distanceOfStones && distanceOfStones <= pow(2*r, 2)+0.003) {
            theta = fabs(atan((y2-y1) / (x2-x1)));
            theta1 = fabs(atan(vy1/vx1));
            theta2 = fabs(atan(vy2/vx2));
            if (y1 >= y2) {
                dif_x = x1 - x2;
                dif_y = y1 - y2;
            } else {
                dif_x = x2 - x1;
                dif_y = y2 - y1;
            }
            
            v1 = sqrt(vx1*vx1 + vy1*vy1);
            v2 = sqrt(vx2*vx2 + vy2*vy2);
            v1_ver = v1 * cos(fabs(angle(theta1, vx1, vy1) - angle(theta, dif_x, dif_y)));
            v1_pal = v1 * sin(fabs(angle(theta1, vx1, vy1) - angle(theta, dif_x, dif_y)));
            v2_ver = v2 * cos(fabs(angle(theta2, vx2, vy2) - angle(theta, dif_x, dif_y)));
            v2_pal = v2 * sin(fabs(angle(theta2, vx2, vy2) - angle(theta, dif_x, dif_y)));
            
            v1_verx = v1_ver * cos(angle(theta, dif_x, dif_y));
            v1_very = v1_ver * sin(angle(theta, dif_x, dif_y));
            v2_verx = v2_ver * cos(angle(theta, dif_x, dif_y));
            v2_very = v2_ver * sin(angle(theta, dif_x, dif_y));
            
            if (angle(theta1, vx1, vy1) >= angle(theta, dif_x, dif_y)) {
                v1_palx = v1_pal * cos(angle(theta, dif_x, dif_y) + M_PI/2);
                v1_paly = v1_pal * sin(angle(theta, dif_x, dif_y) + M_PI/2);
            } else {
                v1_palx = v1_pal * cos(angle(theta, dif_x, dif_y) - M_PI/2);
                v1_paly = v1_pal * sin(angle(theta, dif_x, dif_y) - M_PI/2);
            }
            
            if (angle(theta2, vx2, vy2) >= angle(theta, dif_x, dif_y)) {
                v2_palx = v2_pal * cos(angle(theta, dif_x, dif_y) + M_PI/2);
                v2_paly = v2_pal * sin(angle(theta, dif_x, dif_y) + M_PI/2);
            } else {
                v2_palx = v2_pal * cos(angle(theta, dif_x, dif_y) - M_PI/2);
                v2_paly = v2_pal * sin(angle(theta, dif_x, dif_y) - M_PI/2);
            }
            
            vx1 = v2_verx + v1_palx;
            vy1 = v2_very + v1_paly;
            vx2 = v1_verx + v2_palx;
            vy2 = v1_very + v2_paly;
            
            count++;
            printf("1,%lf, %lf,%lf,%lf,%lf\n",t,x1,y1,x2,y2);
        } else if (i % 100 == 0) {
            printf("0,%lf, %lf,%lf,%lf,%lf\n",t,x1,y1,x2,y2);
        }
        
        t = t + dt;
        x1 = x1 + vx1 * dt;
        y1 = y1 + vy1 * dt;
        x2 = x2 + vx2 * dt;
        y2 = y2 + vy2 * dt;
        
        if (count == 2) {
            break;
        }
    }
}

できたもの

pic.twitter.com/m9nQjzhGrB

— たけのこ (@kabostrawberry) December 2, 2018

これを使えばビリヤードもできそうですね。