octaveでピンポン(?)

C言語とoctaveで、壁に跳ね返る球の軌跡を描いてみました。
なぜこんなことをしてみたかというと、統計力学の授業で、「ピンポン玉を壁にぶつけてはねかえり続けるとどうなるか？」と先生がおっしゃっていたからです。

アルゴリズム

というほどのこともないですが。
座標が箱のサイズ(今回は20*20)の範囲外になると、比例計算して箱の端っこにあるべき座標になります。そして、移動してきた方向を考えて次に進む向きを決めます。両端が半円の時もほぼ同じ仕組みです。

書いたコードと結果(四角形の場合)

▶コード(クリックで展開)

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

int main(int argc, const char * argv[]) {
    float x=0;
    float y=0;
    float x_new=0;
    float y_new=0;
    float v=3;
    float theta=40;
    
    float dif_x=1;
    float dif_y=1;
    
    int N=500;
    int countL=0;
    int countR=0;
    
    for (int i=0;i<=N;i++) {
        
        if (dif_x>0 && dif_y>0) {
            x_new=x+v*sin(theta*M_PI/180);
            y_new=y+v*cos(theta*M_PI/180);
        }else if (dif_x>0 && dif_y<=0){
            x_new=x+v*sin(theta*M_PI/180);
            y_new=y-v*cos(theta*M_PI/180);
        }else if (dif_x<=0 && dif_y>0){
            x_new=x-v*sin(theta*M_PI/180);
            y_new=y+v*cos(theta*M_PI/180);
        }else if (dif_x<=0 && dif_y<=0){
            x_new=x-v*sin(theta*M_PI/180);
            y_new=y-v*cos(theta*M_PI/180);
        }
        
        if(x<0){
            if (dif_y>=0) {
                y=y-((-x)*tan(theta*M_PI/180));
                y_new=y+v*cos(theta*M_PI/180);
            }else if (dif_y<0){
                y=y+((-x)*tan(theta*M_PI/180));
                y_new=y-v*cos(theta*M_PI/180);
            }
            x=0;
            theta=90-theta;
            x_new=x+v*sin(theta*M_PI/180);
        }else if (x>20){
            if (dif_y>=0) {
                y=y-((x-20)*tan(theta*M_PI/180));
                y_new=y+v*cos(theta*M_PI/180);
            }else if (dif_y<0){
                y=y+((x-20)*tan(theta*M_PI/180));
                y_new=y-v*cos(theta*M_PI/180);
            }
            x=20;
            theta=90-theta;
            x_new=x-v*sin(theta*M_PI/180);
        }
        if (y<0.0){
            if (dif_x>=0) {
                x=x-((-y)*tan(theta*M_PI/180));
                x_new=x+v*sin(theta*M_PI/180);
            }else if (dif_x<0){
                x=x+((-y)*tan(theta*M_PI/180));
                x_new=x-v*sin(theta*M_PI/180);
            }
            y=0;
            theta=90-theta;
            y_new=y+v*cos(theta*M_PI/180);
            
        }else if (y>20){
            if (dif_x>=0) {
                x=x-((y-20)*tan(theta*M_PI/180));
                x_new=x+v*sin(theta*M_PI/180);
            }else if (dif_x<0){
                x=x+((y-20)*tan(theta*M_PI/180));
                x_new=x-v*sin(theta*M_PI/180);
            }
            y=20;
            theta=90-theta;
            y_new=y-v*cos(theta*M_PI/180);
            
        }
        
        if (x<=10) {
            countL++;
        }else{
            countR++;
        }
        
        dif_x=x_new-x;
        dif_y=y_new-y;
        printf("%f,%f,%d,%d \n",x,y,countL,countR);
        x=x_new;
        y=y_new;
        }
}

書いたコードと結果(両端が丸い場合)

▶コード(クリックで展開)

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

double ran0() {
    return (double)rand()/((double)RAND_MAX+1);
}

int main(int argc, const char * argv[]) {
    float x=0;
    float y=0;
    float x_new=0;
    float y_new=0;
    float v=3;
    float theta=40;
    
    float dif_x=1;
    float dif_y=1;
    float dif;
    float dif_min;
    float dif_min_x;
    
    int N=500;
    int countL=0;
    int countR=0;
    
    srand( (unsigned)time( NULL ) );
    
    for (int i=0;i<=N;i++) {
        dif_min=20; //initialize
        dif_min_x=0; //initialize
        
        if (dif_x>0 && dif_y>0) {
            x_new=x+v*sin(theta*M_PI/180);
            y_new=y+v*cos(theta*M_PI/180);
        }else if (dif_x>0 && dif_y<=0){
            x_new=x+v*sin(theta*M_PI/180);
            y_new=y-v*cos(theta*M_PI/180);
        }else if (dif_x<=0 && dif_y>0){
            x_new=x-v*sin(theta*M_PI/180);
            y_new=y+v*cos(theta*M_PI/180);
        }else if (dif_x<=0 && dif_y<=0){
            x_new=x-v*sin(theta*M_PI/180);
            y_new=y-v*cos(theta*M_PI/180);
        }
        
        
        if(x<0){
            if ((y-10)*(y-10)>(100-x*x)) {
                for (int j=1;j<50;j++) {
                    x=x+0.1*j;
                    dif=fabsf((y-10)*(y-10)-(100-x*x));
                    if (dif_min>dif) {
                        dif_min=dif;
                        dif_min_x=x;
                    }
                }
                dif_min_x=dif_min_x+ran0();
                if (y>=10) {
                    y=10+sqrt(fabsf(100-dif_min_x*dif_min_x-dif_min));
                }else{
                    y=10-sqrt(fabsf(100-dif_min_x*dif_min_x-dif_min));
                }
                x=dif_min_x;
                theta=90-theta;
                
                if (dif_y>=0) {
                    y_new=y+v*cos(theta*M_PI/180);
                }else if (dif_y<0){
                    y_new=y-v*cos(theta*M_PI/180);
                }
                x_new=x+v*sin(theta*M_PI/180);
            }
        }else if (x>20){
                if ((y-10)*(y-10)>(100-(x-20)*(x-20))) {
                for (int j=0;j<=50;j++) {
                    x=x-0.1;
                    dif=fabsf((y-10)*(y-10)-(100-(x-20)*(x-20)));
                    if (dif_min>dif) {
                        dif_min=dif;
                        dif_min_x=x;
                    }
                }
                dif_min_x=dif_min_x+ran0();
                if (y>=10) {
                    y=10+sqrt(fabsf(100-(dif_min_x-20)*(dif_min_x-20)-dif_min));
                }else{
                    y=10-sqrt(fabsf(100-(dif_min_x-20)*(dif_min_x-20)-dif_min));
                }
                x=dif_min_x;
                theta=90-theta;
                
                if (dif_y>=0) {
                    y_new=y+v*cos(theta*M_PI/180);
                }else if (dif_y<0){
                    y_new=y-v*cos(theta*M_PI/180);
                }
                x_new=x-v*sin(theta*M_PI/180);
            }
            
        }
        
        if (y<0.0){
            if (dif_x>=0) {
                x=x-((-y)*tan(theta*M_PI/180));
                x_new=x+v*sin(theta*M_PI/180);
            }else if (dif_x<0){
                x=x+((-y)*tan(theta*M_PI/180));
                x_new=x-v*sin(theta*M_PI/180);
            }
            y=0;
            theta=90-theta;
            y_new=y+v*cos(theta*M_PI/180);
            
        }else if (y>20){
            if (dif_x>=0) {
                x=x-((y-20)*tan(theta*M_PI/180));
                x_new=x+v*sin(theta*M_PI/180);
            }else if (dif_x<0){
                x=x+((y-20)*tan(theta*M_PI/180));
                x_new=x-v*sin(theta*M_PI/180);
            }
            y=20;
            theta=90-theta;
            y_new=y-v*cos(theta*M_PI/180);
            
        }
        
        if (x<=10) {
            countL++;
        }else{
            countR++;
        }
        
        if (x_new==x) {
            if (dif_x>0 && dif_y>0) {
                x_new=x_new+v*sin(theta*M_PI/180);
                y_new=y_new+v*cos(theta*M_PI/180);
            }else if (dif_x>0 && dif_y<=0){
                x_new=x_new+v*sin(theta*M_PI/180);
                y_new=y_new-v*cos(theta*M_PI/180);
            }else if (dif_x<=0 && dif_y>0){
                x_new=x_new-v*sin(theta*M_PI/180);
                y_new=y_new+v*cos(theta*M_PI/180);
            }else if (dif_x<=0 && dif_y<=0){
                x_new=x_new-v*sin(theta*M_PI/180);
                y_new=y_new-v*cos(theta*M_PI/180);
            }
        }
        
        dif_x=x_new-x;
        dif_y=y_new-y;
        printf("%f,%f,%d,%d \n",x,y,countL,countR);
        x=x_new;
        y=y_new;
    }
}

なんか細かく見ると所々挙動がおかしいのですが、まあだいたいこんな感じになりました。
考察としては、ただ四角い箱だと同じ場所を跳ね返り続けるんですが、両端を丸くすると少しカオスっぽくなることがわかります。

それから、いつも思ってるんですけど全体的にコードが汚いですよね。と思っていてこの後の記事に繋がるんですが、if文の中に指示を入れまくるんじゃなくて関数としてまとめてからif文に入れた方が綺麗で良いらしい。というかそもそもswitch文を使うべきなのかもしれない。
まあこれについては次回の記事でまとめたいと思います。