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(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++;
    }
}

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