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

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

# 実行環境
```zsh
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*}
$$




# ２体問題 実装

```c++: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`を使って座標位置を変更している。
<hr>

-  `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ファイルが巨大になってしまうので１周期分だけで書き出しをやめている。

<hr>

引数をデフォルトのfalseで実行します。なので、アニメーションgifは出力しない。gifをGUIアニメーション表示と同時に出力するとGUIアニメーションがスムーズに動かない。それはgifを出力しているときにアニメーションが止まるから。
```c++
root [0] .x n2body_anime.C()
root [1]
```
以下はGUIのアニメーションが終了したときのスクショです。

![n2body_screenshot.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/258007/2ce4cc0b-3367-c34e-f43e-47bf1c1e6a8e.png)

次は引数をtrueにしてアニメーションgifを出力し、GUIを表示させない。
```c++
root [1] n2body_anime(true)
Info in <TCanvas::Print>: gif file n2body.gif has been created
root [2]
```
![n2body.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/258007/05e08bf6-9e45-3896-cd8f-147e436d2ecd.gif)


# 関連

- [CERN ROOT (C++) グラフ作成 (３) for とsleepを使ってアニメーション (sinの進行波, sin cos描画)](https://qiita.com/h6591/items/15cdf5eefef725e26720)
- [CERN ROOT (C++) グラフ作成 (2) yをxの計算式として与え、x-yグラフを描画 (class TF1 )](https://qiita.com/h6591/items/a81b37c3c3b019db39c7)
- [CERN ROOT (C++) グラフ作成 (1) sin,cosのグラフの重ね表示](https://qiita.com/h6591/items/085668110a8a016622a5)

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