プログラミングの備忘録

プログラムをつくる過程を残すもの

processingの備忘録 -万有引力-

f:id:taq2777:20211022203726p:plain こんにちは。
今回はprocessingに「万有引力の法則」を導入して、惑星の運動のシミュレーション的なものをつくっていこうと思います。


万有引力の法則は、「全てのものは引き合っている」という物理学者のアイザック・ニュートンで有名な法則で、 $$ F=G \frac{Mm}{r^{2}} $$ と表されます。ここで、$F$は物体間にはたらく引力、$G$は万有引力定数、$M, m$は物体の重さ、$r$は物体間の距離を表します。


ではprocessingに導入していきます。
まずは、太陽と地球にあたるものを描いてみます。(実際は直径は100倍くらい違いますし、距離も遠いですが、適当に)

void setup(){
    size(500, 500);
}

void draw(){
    background(0);
    
    fill(255, 100, 0);
    ellipse(width/2, height/2, 40, 40);
    
    fill(0, 0, 255);
    ellipse(width*3/4, height/2, 20, 20);
}


これに動きをつけます。

float x, y, vx, vy;

void setup(){
    size(500, 500);
  
    x = width*3/4;
    y = height/2;
    vx = 0;
    vy = -3;
}

void draw(){
    background(0);
    
    fill(255, 100, 0);
    ellipse(width/2, height/2, 40, 40);
    
    x += vx;
    y += vy;    
    
    fill(0, 0, 255);
    ellipse(x, y, 20, 20);
}


このままでは円運動をせずに宇宙の彼方に飛んで行ってしまうので、万有引力を追加します。
LJポテンシャルを導入したときと同じやり方です。
力$F$がわかっているので運動方程式$ma=F$にあてはめて加速度$a$の式にして、これをx方向とy方向に分けて速度vxとvyに足していけばよいです。
やってみると加速度$a$は、 $$ a=-G \frac{M}{r^{2}} $$ と表せることがわかったので、プログラムに組み込みます。
万有引力はその名の通り引力なので、「-」符号をつけます。)

float xs, ys, ms;
float xe, ye, vxe, vye, me;
float G;

void setup(){
    size(500, 500);
  
    G = 10;
  
    xs = width/2;
    ys = height/2;
    ms = 100;
  
    xe = width*3/4;
    ye = height/2;
    vxe = 0;
    vye = -3;
    me = 1;
}

void draw(){
    background(0);
    
    fill(255, 100, 0);
    ellipse(xs, ys, 40, 40);
    
    float r = sqrt((xe-xs)*(xe-xs)+(ye-ys)*(ye-ys)); //距離の計算
    float a = -G*ms/(r*r); //加速度の計算
    float angle = atan2(ye-ys, xe-xs); //角度の計算
    
    vxe += a*cos(angle); //x方向に
    vye += a*sin(angle); //y方向に
    xe += vxe;
    ye += vye;    
    
    fill(0, 0, 255);
    ellipse(xe, ye, 20, 20);
}

atan2」はある点とある点の角度を求めることができる関数です。引数が(y方向の距離, x方向の距離)と他の関数とは逆になっていることに注意してください。

これで青丸に円運動をさせることができました。実際は赤丸も青丸の影響で少し運動していますが、ここでは無視しています。


ちなみに、赤丸にも万有引力の法則を適用すると以下のようなプログラムになります。
(msとmeを100、vxsを1、vysを2、vxeを-1、vyeを-2のように変えてみるとくるくる回ります。)

float xs, ys, vxs, vys, ms;
float xe, ye, vxe, vye, me;
float G;

void setup(){
    size(500, 500);
  
    G = 10;
  
    xs = width/2;
    ys = height/2;
    vxs = 0;
    vys = 0;
    ms = 100;
  
    xe = width*3/4;
    ye = height/2;
    vxe = 0;
    vye = -3;
    me = 1;
}

void draw(){
    background(0);
    
    float r = sqrt((xe-xs)*(xe-xs)+(ye-ys)*(ye-ys));
    float as = -G*me/(r*r);
    float angle = atan2(ye-ys, xe-xs);
    
    vxs += as*cos(angle);
    vys += as*sin(angle);    
    xs += vxs;
    ys += vys;    
    
    fill(255, 100, 0);
    ellipse(xs, ys, 40, 40);
    
    float ae = -G*ms/(r*r);
    
    vxe += ae*cos(angle);
    vye += ae*sin(angle);    
    xe += vxe;
    ye += vye;    
    
    fill(0, 0, 255);
    ellipse(xe, ye, 20, 20);
}


このまま何個も増やせますが、楽をするためにクラスに入れてみます。

float xs, ys, ms;
int n = 5;
Star[] stars;
int j = 0;

void setup(){
    size(500, 500);
    
    xs = width/2;
    ys = height/2;
    ms = 100;
    
    stars = new Star[n];
    for(int i = 0; i < n; i ++){
        stars[i] = new Star(i);
    }
}

void draw(){
    background(0);
    
    fill(255, 100, 0);
    ellipse(xs, ys, 40, 40);
    
    for(int i = 0; i < n; i ++){
        stars[i].display();
        stars[i].move();
    }
}

void mousePressed(){
    stars[j].data();
    j ++;
    if(j == n){
        j = 0;
    }
}

class Star{
    int id;
    float x, y, d, vx, vy, m;
    float G = 10;
    
    Star(int _id){
        _id = id;
    }
    
    void data(){
        x = mouseX;
        y = mouseY;
        d = 20;
        vx = random(-2, 2);
        vy = random(-2, 2);
        m = 1;
    }
    
    void display(){
        fill(0, 0, 255);
        ellipse(x, y, d, d);
    }
    
    void move(){
        float r = sqrt((x-xs)*(x-xs)+(y-ys)*(y-ys));
        float a = -G*ms/(r*r);
        float angle = atan2(y-ys, x-xs);
        
        vx += a*cos(angle);
        vy += a*sin(angle);    
        x += vx;
        y += vy;
    }
}

実際は惑星ですが、広義の星ということでstarにしました。
新たに「j」という変数をつくって、クリックするたびに星が追加されると同時にjも1ずつ足されるようにしています。そしてjがnと同じ値になったらjを0に戻しています。


これまでは星の速度はランダムに決めていましたが、引っ張ることで任意に速度を決められるようにしてみます。(モン○タース○ライク的に)

float xs, ys, ms;
int n = 5;
Star[] stars;
int j = 0;
boolean[] pause = new boolean[n];
float px, py;

void setup(){
    size(500, 500);
    noStroke();
    
    xs = width/2;
    ys = height/2;
    ms = 100;
    
    stars = new Star[n];
    for(int i = 0; i < n; i ++){
        stars[i] = new Star(i);
    }
}

void draw(){
    background(0);
    
    fill(255, 100, 0);
    ellipse(xs, ys, 40, 40);
    
    for(int i = 0; i < n; i ++){
        stars[i].display();
        if(pause[i] == false){
            stars[i].move();
        }
    }
    
    if(pause[j] == true){
        strokeWeight(1);
        stroke(255, 0, 0);
        line(px, py, px-1.5*(mouseX-px), py-1.5*(mouseY-py));
        noStroke();
    }
}

void mousePressed(){
    pause[j] = true;
  
    stars[j].x = mouseX;
    stars[j].y = mouseY;
    px = mouseX;
    py = mouseY;
    
    stars[j].data();
}

void mouseReleased(){
    pause[j] = false;
  
    stars[j].vx = -0.05*(mouseX-px); 
    stars[j].vy = -0.05*(mouseY-py);
    
    j ++;
    if(j == n){
        j = 0;
    }
}

class Star{
    int id;
    float x, y, d, vx, vy, m;
    float G = 10;
    
    Star(int _id){
        _id = id;
    }
    
    void data(){
        d = 20;
        m = 1;
    }
    
    void display(){
        fill(0, 0, 255);
        ellipse(x, y, d, d);
    }
    
    void move(){
        float r = sqrt((x-xs)*(x-xs)+(y-ys)*(y-ys));
        float a = -G*ms/(r*r);
        float angle = atan2(y-ys, x-xs);
        
        vx += a*cos(angle);
        vy += a*sin(angle);    
        x += vx;
        y += vy;
    }
}

mouseReleased()で、離したときに押したときの位置と離したときの位置の差を何倍かした値で星の速度を設定するようにしています。
「pause」という配列で押したときにポーズし(pause[j]=true)、離したときに動くようにする(pause[j]=false)ことと、star[i].move()の手前にif文で「i番目の星のpauseがfalseのとき」という処理を追加することで、今追加しようとしている星だけ止めて、他は動いたままにするという描画をできるようにしています。
j++以下の処理は都合上mouseReleased()に移しました。


さらに青丸同士での万有引力による相互作用も追加しようと思いましたが、うまくいきませんでした。

追記(2021/11/19)
一応、青丸同士での相互作用も追加できました。(「一応」の理由は後述)
コードは以下です。

float xs, ys, ms;
int n = 5;
Star[] stars;
int j = 0;
boolean[] pause = new boolean[n];
float px, py;
float G = 1;

void setup(){
    size(500, 500);
    noStroke();
    
    xs = width/2;
    ys = height/2;
    ms = 100;
    
    stars = new Star[n];
    for(int i = 0; i < n; i ++){
        stars[i] = new Star(i);
    }
}

void draw(){
    //background(0);
    fill(0, 50);
    rect(0, 0, width, height);
    
    fill(255, 100, 0);
    ellipse(xs, ys, 40, 40);
    
    for(int i = 0; i < n; i ++){
        stars[i].display();
        stars[i].gravitation();

        if(pause[i] == false){
            stars[i].move();
        }
    }
    
    if(pause[j] == true){
        strokeWeight(1);
        stroke(255, 0, 0);
        line(px, py, px-1.5*(mouseX-px), py-1.5*(mouseY-py));
        noStroke();
    }
}

void mousePressed(){
    pause[j] = true;
  
    stars[j].x = mouseX;
    stars[j].y = mouseY;
    px = mouseX;
    py = mouseY;
    
    stars[j].data();
}

void mouseReleased(){
    pause[j] = false;
  
    stars[j].vx = -0.05*(mouseX-px); 
    stars[j].vy = -0.05*(mouseY-py);
    
    j ++;
    if(j == n){
        j = 0;
    }
}

class Star{
    int id;
    float x, y, d, vx, vy, m, visible;
    
    Star(int _id){
        _id = id;
        visible = 0;
    }
    
    void data(){
        d = 20;
        m = 100;
        visible = 255;
    }
    
    void display(){
        if(visible > 0){
            fill(0, 0, 255, visible);
            ellipse(x, y, d, d);
        }
    }
    
    void move(){    
        x += vx;
        y += vy;
    }

    void gravitation(){
        if(visible > 0){
            float dx = x - xs;
            float dy = y - ys;
            float r = sqrt(dx*dx + dy*dy);
            float a = -G*ms/(r*r);
            float angle = atan2(dy, dx);
            
            vx += a*cos(angle);
            vy += a*sin(angle);
            
            for(int i = id+1; i < n; i ++){
                if(stars[i].visible > 0){
                    float dx2 = x - stars[i].x;
                    float dy2 = y - stars[i].y;
                    float r2 = sqrt(dx2*dx2 + dy2*dy2);
                    if(r2 == 0){
                        return;
                    }
                    float a2 = -G*stars[i].m/(r2*r2);
                    float angle2 = atan2(dy2, dx2);
                    
                    vx += a2*cos(angle2);
                    vy += a2*sin(angle2);
                    stars[i].vx -= a2*cos(angle2);
                    stars[i].vy -= a2*sin(angle2);
                }
            }
        }
    }
}

visibleという変数を追加して、見えているときは255、見えていないときは0としています。
新たにclass内にgravitation()という関数を追加して万有引力の計算を行うことにしました。ある青丸(x, y)と赤丸(xs, ys)、そしてある青丸(x, y)と他の青丸(stars[i].x, stars[i].y)(i=id+1としているため、番号でいうとある青丸stars[k]の次の青丸stars[k+1])との相互作用を計算する処理をそれぞれ書いています。(ひとつにまとめられそうだが断念)
そしてここが「一応」ポイントなのですが、変数を見ているとなぜかr2の値が0と正しい値を交互にとり、それに伴ってaの値も-Infinityと正しい値を交互にとることになってうまくいかないので、r2が0になったときは途中で止めて次の処理にいく(if(r2 == 0){return;}の部分)というように変更しました。(おそらく青丸をクリックで追加する処理が何か影響している)
というわけで、一応はうまくいったのですがまだまだ改善の余地はあります。
以上


ここに衝突など加えてもいいですが、惑星なら衝突=崩壊なので、ぶつかったら小さな丸に分ける処理なども追加してみるといいかもしれません。


ちなみに、惑星の運動に似たようなものに「ボーアモデル」というものもあります。
ボーアモデルは原子の構造を表したもので、太陽と地球のように、原子核の周りを電子が回っていると考えています。
(実際は電子は回っているわけではなく確率的に存在するとされていますが、ボーアモデルで考えるとうまく説明できることもあるのでよく使われます。)

ボーアモデルでの相互作用は万有引力ではなくクーロン力ですが、クーロン力は、 $$ F= \frac{1}{4 \pi \varepsilon_{0}} \frac{Qq}{r^{2}} $$ と表され、万有引力と同じ形なので同じプログラムで表すことができます。


というわけで、processingに万有引力クーロン力を導入して惑星の運動やボーアモデルを表せるようになりました。
いろいろいじって遊んでみてください。
では、また次回。