4 時間積分

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

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

$$v_{i+1} = 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++;
    }
}

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