next up previous
Next: 5 プログラム全体の構造 Up: 計算天文学 II 第11回 データ構造とアルゴリズム(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(double dt){
            double dt2 = dt*dt*0.5;
            pos += dt*vel + dt2*acc;
            vel +=  (dt*0.5)*acc;
        }
        void correct(double dt){
            vel +=  (dt*0.5)*acc;
        }

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

void integrate(bhnode * bn,
               int nnodes,
               particle * pp,
               int n,
               double eps2,
               double theta,
               double 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,
                       double eps2,
                       double theta)
{
    double rsize = calculate_size(pp,n);
    bn->assign_root(myvector(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
平成17年1月23日