LoginSignup
1
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@h6591

CERN ROOT (C++) グラフ作成 (4) 2体問題シミュレーション　アニメーション

Last updated at Posted at 2024-11-05

概要

２体問題の物理シミュレーションをTGraphを使ってアニメーション化します。

実行環境

sw_vers
ProductName:		macOS
ProductVersion:		14.6.1
BuildVersion:		23G93

root --version
ROOT Version: 6.32.06
Built for macosxarm64 on Sep 21 2024, 18:21:53
From tags/6-32-06@6-32-06

２体問題 理論

２体(body1, body2)の万有引力は一般に次の式で表される

$$
\vec{F}_{12} = -G\frac{m_1m_2}{|\vec{r}_1-\vec{r}_2|^3}(\vec{r}_1-\vec{r}_2)
$$

$\vec{F}_{12}$ : body2がbody1に及ぼす力
$m1,m2$ : body1,body2の質量
$\vec{r}_1,\vec{r}_2$ : body1,body2の位置ベクトル
よって、body1の運動方程式と加速度$\vec{a}_1$は
$$
\begin{align*}
m_1\vec{a}_1&=-G\frac{m_1m_2}{|\vec{r}_1-\vec{r}_2|^3}(\vec{r}_1-\vec{r}_2) \\
\vec{a}_1&=-G\frac{m_2}{|\vec{r}_1-\vec{r}_2|^3}(\vec{r}_1-\vec{r}_2)
\end{align*}
$$
body2が受ける力はbody1とは反対向きなのでbody2の加速度$\vec{a}_2$はbody1と同様に
$$
\vec{a}_2=G\frac{m_1}{|\vec{r}_1-\vec{r}_2|^3}(\vec{r}_1-\vec{r}_2)
$$

Leapfrog法を使ってシミュレーションが出来るようにする。

$$
\begin{align*}
\vec{v}_{n+1/2}&=\vec{v}_n + \vec{a}_n \Delta t/2 \\
\vec{r} _{n+1}&=\vec{r} _n + \vec{v} _{n+1/2}\Delta t\\
\vec{v} _{n+1} &= \vec{v} _{n+1/2} + \vec{a} _{n+1} \Delta t/2
\end{align*}
$$

２体問題 実装

n2body_anime.C
void set_xyrange(TGraph* gr){
    //TCanvasを直接引数にしたほうが簡単
    //TCanvasを直接アクセスできないときの例として記述した
    TCanvas* c = (TCanvas*)gROOT->FindObject("canvas");
    if (c == nullptr) return;  //エラー処理が必要
    auto height = c->GetWindowHeight();
    auto width = c->GetWindowWidth();
    constexpr double x_range = 2.0;
    // x軸とy軸の目盛りに対するピクセル数を同じにする
    auto y_range = x_range * height / width; 
    gr->SetMaximum(y_range);   //y軸の最大値
    gr->SetMinimum(-y_range);  //y軸の最小値
    gr->GetXaxis()->SetLimits(-x_range, x_range);  //x軸の範囲    
}
void n2body_anime(bool gif=false) {
    gROOT->SetBatch(gif); //true: GUIを表示させない

    //アニメーションgifを書き出すタイミング
    constexpr int gif_max_interval = 10;
    //アニメーションgifのフレーム数
    constexpr int gif_max_print = 25;
    //Sleep時間s
    constexpr double timer = 1000/120.0; //ミリ秒 1秒間に約120回描画
    //実行時間
    constexpr int exec_time = 10;  //秒
    //実行ステップ回数
    constexpr int repeat_count = (int)1000 * exec_time / timer; //ミリ秒
    //time step
    constexpr double dt =  0.01;
    constexpr double dt_half =  dt / 2.0;
    //万有引力定数
    constexpr double G = 6.67430e-11;    
    //body1の質量
    constexpr double body1_m = 1.0e+11;
    //body2の質量
    constexpr double body2_m = 2.0e+11;

    constexpr double body1_init_vel_y = 1.1;
    //body1とbody2の運動量の和が0になるように決める
    constexpr double body2_init_vel_y = -body1_init_vel_y * body1_m / body2_m;

    //body1の位置
    ROOT::Math::XYVector body1_pos(-1.5, 0.0);
    //body1の速度
    ROOT::Math::XYVector body1_vel(0.0, body1_init_vel_y);
    //body2の位置
    ROOT::Math::XYVector body2_pos(1.0, 0.0);
    //body2の速度
    ROOT::Math::XYVector body2_vel(0.0, body2_init_vel_y);
    
    auto canvas = new TCanvas("canvas");
    auto body1 = new TGraph();
    auto body2 = new TGraph();


    body1->AddPoint(body1_pos.x(),body1_pos.y());
    body1->SetMarkerStyle(20);
    body1->SetMarkerColor(kRed);
    body1->SetMarkerSize(2);

    body2->AddPoint(body2_pos.x(),body2_pos.y());
    body2->SetMarkerStyle(20);
    body2->SetMarkerColor(kGreen);
    body2->SetMarkerSize(4);
    
    set_xyrange(body1);

    body1->Draw("ap");
    body2->Draw("p");
    canvas->ModifiedUpdate();

    auto r_vec = body1_pos - body2_pos;  //rベクトル
    auto r3  = pow(r_vec.mag2(),1.5);    //r^3
    auto force1_vec = -G * body2_m / r3 * r_vec; //body1に作用する力
    auto force2_vec =  G * body1_m / r3 * r_vec; //body2に作用する力
   
    int gif_interval = gif_max_interval;
    int gif_count = 0;

    for(int i = 0; i < repeat_count && !gSystem->ProcessEvents(); i++) {
        if (!gif)
            gSystem->Sleep(timer);
        auto body1_vel_half =  body1_vel + dt_half * force1_vec;
        auto body2_vel_half =  body2_vel + dt_half * force2_vec;
        body1_pos +=  dt * body1_vel_half;
        body2_pos +=  dt * body2_vel_half;

        r_vec = body1_pos - body2_pos;    
        r3  = pow(r_vec.mag2(),1.5);      
        force1_vec = -G * body2_m / r3 * r_vec; 
        force2_vec =  G * body1_m / r3 * r_vec;
        
        body1_vel = body1_vel_half + force1_vec * dt_half; 
        body2_vel = body2_vel_half + force2_vec * dt_half; 

        body1->SetPoint(0,body1_pos.x(), body1_pos.y());
        body2->SetPoint(0,body2_pos.x(), body2_pos.y());
        
        set_xyrange(body1);
        canvas->ModifiedUpdate();
        if (gif) {
            if (gif_interval == gif_max_interval) {
                if (gif_count <= gif_max_print) {
                    if (gif_count == gif_max_print) {
                        canvas->Print("n2body.gif++30++"); //永久ループ
                        break;
                    } else {
                        if (gif_count == 0)
                            canvas->Print("n2body.gif");   
                        else
                            canvas->Print("n2body.gif+30");  //30*10ms
                    }
                    gif_count++;
                } 
                gif_interval = 0;
            }
            gif_interval++; 
        }
    }
}

永久ループにもできますが、10秒間だけ動くようにプログラムしました。途中でやめるときはCTRL+Cで中止してください。

body1とbody2の初期速度はbody1とbody2の運動量(質量x速度)の和がゼロになるように決めて、重心位置が移動しないようにする。
上記例では、速度のx成分はどちらもゼロでy成分の運動量の和がゼロになるように計算して設定している。
body1,body2の位置の移動は->SetPointを使って座標位置を変更している。

  • gROOT->SetBatch(gif); //true: GUIを表示させない
    • アニメーションgifを出力させるときはGUIを表示させないように、逆にgifを出力しないときはGUIを表示させるようにする
  • canvas->Print("n2body.gif++30++"); //永久ループ
    • アニメーションgifが永久ループするようにする。最後のフレームの書き出しのときに指定する
  • canvas->Print("n2body.gif");
    • gifファイルの最初のフレームの書き出し
    • n2body.gif+30はappendなので、n2body.gifがもともと存在すると意図したアニメーションgifにならない。n2body.gifが存在しないときは無くてもOK。
  • canvas->Print("n2body.gif+30"); //30*10ms
    • 30*10msだけ停止してフーレムを表示する

アニメーションgifを毎回書き出すたとgifファイルが巨大になってしまうので１周期分だけで書き出しをやめている。

引数をデフォルトのfalseで実行します。なので、アニメーションgifは出力しない。gifをGUIアニメーション表示と同時に出力するとGUIアニメーションがスムーズに動かない。それはgifを出力しているときにアニメーションが止まるから。

root [0] .x n2body_anime.C()
root [1]

以下はGUIのアニメーションが終了したときのスクショです。

n2body_screenshot.png

次は引数をtrueにしてアニメーションgifを出力し、GUIを表示させない。

root [1] n2body_anime(true)
Info in <TCanvas::Print>: gif file n2body.gif has been created
root [2]

n2body.gif

関連

終わりに

そんなに難しくなく出来たと思います。やはり、アニメーションgifで見るよりGUIのアニメーションで見たほうが全然よいです。実際はゲームエンジンで作ったほうがもっとスムーズで綺麗かもしれませんが。

1
2
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
1
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?