next up previous
Next: 5 プログラム全体の構造 Up: 計算天文学 II 第10回 データ構造とアルゴリズム(2) Previous: 3 各関数のテスト

4 時間積分

時間積分には前にやったリープフロッグ公式を使うことにする。どんなふうに 書いてもいいが、今回は、分かりやすいと思われる


$\displaystyle x_{i+1}$ $\textstyle =$ $\displaystyle x_{i} + \Delta t v_{i} + \Delta t^2 a(x_i)/2$ (1)
$\displaystyle v_{i+1}$ $\textstyle =$ $\displaystyle v_{i} + \Delta t [a(x_i)+a(x_{i+1})]/2$ (2)

の形を使うことにしよう。以下の関数を particle クラスに入れる。

pt minus 1 pt

        void predict(real dt){
            real dt2 = dt*dt*0.5;
            pos += dt*vel + dt2*acc;
            vel +=  (dt*0.5)*acc;
        }
        void correct(real dt){
            vel +=  (dt*0.5)*acc;
        }

時間積分は以下のように、まず予測子を計算、それから加速度を計算し、最後に速度の修正子を計算する。 pt minus 1 pt

void integrate(bhnode * bn,
               int nnodes,
               particle * pp,
               int n,
               real eps2,
               real theta,
               real dt)
{
    for(int i = 0;i<n;i++) pp[i].predict(dt);
    calculate_gravity(bn,nnodes,pp,n,eps2,theta);
    for(int i = 0;i<n;i++) pp[i].correct(dt);
}

重力を計算する関数はこんな感じになる。

pt minus 1 pt

void calculate_gravity(bhnode* bn,
                       int nnodes,
                       particle * pp,
                       int n,
                       real eps2,
                       real theta)
{
    real rsize = calculate_size(pp,n);
    bn->assign_root(vector(0.0), rsize*2, pp, n);
    int heap_remainder = nnodes-1;
    bhnode * btmp = bn + 1;
    bn->create_tree_recursive(btmp,  heap_remainder);
    //    bn->dump();
    //    PRL(bn->sanity_check());
    bn->set_cm_quantities();
    clear_acc_and_phi(pp, n);
    
    particle * p = pp;
    for(int i = 0; i<n; i++){
        calculate_gravity_using_tree(p, bn, eps2, theta);
        //        PR(i); PR(p->pos);PR(p->phi); PRL(p->acc);
        p++;
    }
}

これは基本的には前のテスト用のメインプログラムから重力計算のところを切 り出してまとめただけである。



Jun Makino
平成15年1月27日