==============
多体系
==============

.. toctree::
   :maxdepth: 1

　金星や火星あるいは月を加えた系を計算したくなる．
地球と太陽の二つしかないときは二体問題だが，月が増えれば三体問題．三体以上は多体と言い，解析的に軌道を得ることができないことが知られている．


数が増えると
=================

これまでのコードの書き方ではいくつかの問題が生じる．
一つは計算の主要部分，

.. code-block:: fortran

        fu_temp = fu(x,y)
        fv_temp = fv(x,y)
        x = x + fx(u,v)*delta_t + 0.5d0*fu_temp*delta_t**2
        y = y + fy(u,v)*delta_t + 0.5d0*fv_temp*delta_t**2
        u = u + 0.5d0*(fu_temp + fu(x,y))*delta_t
        v = v + 0.5d0*(fv_temp + fv(x,y))*delta_t

これまでの二体問題では，（太陽を中心に固定していたので）
地球の位置と速度の情報（x,y,u,v）だけを更新していればよかったが，
これにもう一つ天体が加わると，その分だけ変数が（2次元であれば）４つ増えることになる．
さらにもう一つ増やすと．．．と天体を増やすたびに変数がどんどん増えていく．
おまけに変数ごとに上記のような計算式を書いて，となると非常に面倒である．

次に，大切な相互作用の計算部分でも，

.. code-block:: fortran

        r = sqrt(x**2+y**2)
        xx = -y / r**3

これまでは，太陽との万有引力だけを計算していればよかったが，
新たに天体を加えたら，その天体との相互作用も計算しなくてはいけない．

そして，
運動方程式（再掲）も，


.. math:: 
    \frac{d^2 \boldsymbol{x}_i}{dt^2} = \sum_{j \neq i}Gm_j\frac{\boldsymbol{x}_j- \boldsymbol{x}_i}{|\boldsymbol{x}_j-\boldsymbol{x}_i|^3}
    :label: motion of eqation

太陽と地球だけだった場合は，無次元化により方程式の係数をすっきりさせることができたが，
それぞれ異なる質量をもっている天体が存在する場合はそうはいかない．
無次元化はするが，それぞれの質量に関する情報も必要である．

ということで，思い切って，どっと変更．

.. code-block:: fortran
    :linenos:
    :emphasize-lines: 77-114,230-237
    
    !+++++++++++++++++++++++++++++++++++++++++++++++
    module globals

    double precision,parameter::    pi = acos(-1.0d0)
    double precision,parameter::    G = 6.67384d-11
    
    integer,parameter::             ndim = 2
    integer,parameter::             nbody = 5

    end module globals
    !+++++++++++++++++++++++++++++++++++++++++++++++
    module var_trajectory
    
    use globals
    implicit none
    
    double precision                coord(ndim,nbody)
    double precision                veloc(ndim,nbody)

    end module var_trajectory
    !+++++++++++++++++++++++++++++++++++++++++++++++
    module param_planet

    use globals
    implicit none

    double precision                mass(nbody)

    end module param_planet
    !+++++++++++++++++++++++++++++++++++++++++++++++
    module cal_ene
    
    use globals
    use param_planet
    implicit none
    
    contains
        subroutine cal_kin(v,e)
            double precision,intent(in)::     v(ndim,nbody)
            double precision,intent(out)::    e
            integer     i,j
            double precision    v2

            e = 0.0d0
            do i=1,nbody
                v2 = 0.0d0
                do j=1,ndim
                    v2 = v2 + v(j,i)**2
                end do
                e = e + 0.5d0*mass(i)*v2
            end do
        end subroutine cal_kin
        !++++++++++++++++++++++++
        subroutine cal_pot(c,e)
            double precision,intent(in)::   c(ndim,nbody)
            double precision,intent(out)::    e
            integer     i,j,k
            double precision    rij(ndim)
            double precision    rij2,rij_abs

            e = 0.0d0
            do i=1,nbody-1
                do j=i+1,nbody
                    rij(:) = c(:,j)-c(:,i)
                    rij2 = 0.0d0
                    do k=1,ndim
                        rij2 = rij2 + rij(k)**2
                    end do
                    rij_abs = sqrt(rij2)
                    e = e + mass(i)*mass(j)/rij_abs
                end do
            end do
        end subroutine cal_pot

    end module cal_ene
    !+++++++++++++++++++++++++++++++++++++++++++++++
    module cal_force

    use globals
    use param_planet
    implicit none
    
    contains
        subroutine force(c,a)
            double precision,intent(in)::       c(ndim,nbody)
            double precision,intent(out)::      a(ndim,nbody)
            integer     i,j,k
            double precision    rij(ndim)
            double precision    rij2,rij_abs,inv_rij_abs
            double precision    ftmp
            double precision    fij(ndim),f(ndim,nbody)

            f = 0.0d0
            do i = 1,nbody-1
                do j = i+1,nbody
                    rij(:) = c(:,j)-c(:,i)
                    rij2 = 0.0d0
                    do k=1,ndim
                        rij2 = rij2 + rij(k)**2
                    end do
                    rij_abs = sqrt(rij2)
                    inv_rij_abs = 1.0d0/rij_abs
                    ftmp = mass(i)*mass(j)*(inv_rij_abs**3)
                    fij(:) = ftmp*rij(:)
                    f(:,i) = f(:,i) + fij(:)
                    f(:,j) = f(:,j) - fij(:)
                end do
            end do
            do i = 1,nbody
                a(:,i) = f(:,i) / mass(i)
            end do
        end subroutine force

    end module cal_force
    !+++++++++++++++++++++++++++++++++++++++++++++++
    program celestial
    
    use globals
    use var_trajectory
    use param_planet
    use cal_ene
    use cal_force

    implicit none

    integer                     i,dout
    integer,parameter::         n=10**6
    integer,parameter::         noutmax=10000
    integer,parameter::         nlive=10
    double precision            tmax
    double precision            t,dt,dt2
    double precision            cnor,vnor,mnor,tnor
    double precision            E,E0,Ev,Ep,dE,Eold,sum_dE
    double precision            start_time,finish_time
    double precision            accel(ndim,nbody)
    integer                     nliveout

    !----- CPU time
    call CPU_TIME(start_time)

    !----- mass
    mass(1) = 1.989d30          ! mass of sun[kg]
    mass(2) = 5.972d24          ! mass of earth[kg]
    mass(3) = 4.869d24          ! venus
    mass(4) = 6.419d23          ! mars
    mass(5) = 7.348d22          ! mass of moon[kg]

    !----- initial conditions
    coord(1,1) = 0.0d0
    coord(2,1) = 0.0d0
    coord(1,2) = 1.496d11       !average radius of revolution of earth[m]
    coord(2,2) = 0.0d0
    coord(1,3) = 1.082d11   
    coord(2,3) = 0.0d0
    coord(1,4) = 2.279d11
    coord(2,4) = 0.0d0
    coord(1,5) = 3.844d8        !average radius of revolution of moon[m]
    coord(1,5) = coord(1,5) + coord(1,2)
    coord(2,5) = 0.0d0

    veloc(1,1) = 0.0d0
    veloc(2,1) = 0.0d0
    veloc(1,2) = 0.0d0
    veloc(2,2) = 2.978d4        !average velocity of earth[m/s]
    veloc(1,3) = 0.0d0
    veloc(2,3) = 3.502d4
    veloc(1,4) = 0.0d0
    veloc(2,4) = 2.413d4
    veloc(1,5) = 0.0d0
    veloc(2,5) = 1.022d3        !average velocity of moon[m/s]
    veloc(2,5) = veloc(2,5) + veloc(2,2)


    !----- time
    t = 0.0d0
    tmax = 3.2d8              !365day*24hour*60min*60sec = 31536000sec = 3.2d7

    !----- normalizse
    cnor = sqrt((coord(1,2)-coord(1,1))**2 + (coord(2,2)-coord(2,1))**2)
    mnor = mass(1)
    tnor = sqrt(cnor**3/(G*mnor))
    vnor = cnor/tnor

    coord = coord/cnor
    veloc = veloc/vnor
    mass = mass/mnor
    t = t/tnor
    tmax = tmax/tnor

    !------ main
    dt = (tmax-t)/dble(n)
    dt2 = 0.5d0*dt

    open(unit=10,file='outputNvV.dat',status='replace')
    write(10,*)'t_NvV, E_NvV, dE_NvV, x1_NvV, y1_NvV, x2_NvV, y2_NvV, &
        &  x3_NvV, y3_NvV, x4_NvV, y4_NvV, x5_NvV, y5_NvV'

    nliveout = int(n/nlive)
    if(n .gt. noutmax) then
        dout = int(n/noutmax)
    else
        dout = 1
    end if
    
    do i = 1, n-1
        
        if(mod(i,nliveout) .eq. 0) then
            write(*,*)'current step:',i,'/total step:',n
        end if

        if(i .eq. 1) then
            call cal_kin(veloc,Ev)
            call cal_pot(coord,Ep)
            E0 = Ev-Ep
            Eold = E0
        end if

        call cal_kin(veloc,Ev)
        call cal_pot(coord,Ep)
        E = Ev-Ep

        dE = abs((E - E0)/E0)
        sum_dE = sum_dE + abs(E-Eold)
        
        if(mod(i-1,dout) .eq. 0) then
            write(10,*)t,',',E,',',dE,',',coord(1,1),',',coord(2,1),',',coord(1,2),',',coord(2,2),',',coord(1,3) &
            & ,',',coord(2,3),',',coord(1,4),',',coord(2,4),',',coord(1,5),',',coord(2,5)
        end if

        if(i .eq. 1) then
            call force(coord,accel)
        end if

        veloc = veloc + dt2 * accel
        coord = coord + dt * veloc
        call force(coord,accel)
        veloc = veloc + dt2 * accel

        t = t + dt
        Eold = E

    end do 

    close(10)

    write(*,*) 'Time step:',dt
    write(*,*) 'Average error of energy update:',sum_dE/dble(n)

    !-----Output CPU time
    call CPU_TIME(finish_time)
    write(*,*) 'CPU time:',finish_time - start_time

    end program celestial
    !+++++++++++++++++++++++++++++++++++++++++++++++

貼り付けるとかなり長いコードになってしまったが，中身はそれほど変わっていない．

大切な主要部分は書き方を変更．

.. code-block:: fortran

        if(i .eq. 1) then
            call force(coord,accel)
        end if

        veloc = veloc + dt2 * accel
        coord = coord + dt * veloc
        call force(coord,accel)
        veloc = veloc + dt2 * accel

以前の書き方よりもすっきりしているように見える．
ただし，座標（coord），速度（veloc）と加速度（accel）の変数はすべて配列で記述している．
モジュール「var_trajectory」の中で定義されているように，ndim=2,nbody=5，2×5の大きさをもつ配列である．
coord(1,1)は1番目の天体のx座標，veloc(2,5)は5番目の天体のy方向速度を格納する．

.. index::
    single: scheme; Velocity-Verlet

スキームはこれまでと同じく速度ベルレ法である．見た目が少々違うが，

.. math::
    \dot x &= u, \dot u &= \Phi (x) 

に対して，

次のようなアルゴリズムであり，中身は変わっていない．

.. math::
    &u^* = v(t) + \Phi (x(t))\frac{\Delta t}{2}\\
    &x(t+\Delta t) = x(t) + u^* \Delta t\\
    &u(t+\Delta t) = u^* + \Phi (x(t+\Delta t))\frac{\Delta t}{2}

まず，半歩進めた仮の速度v* を計算する．
最初のステップ（i=1）では，計算に必要なaccelの値がないので，この時のみif文内で計算している．
その後は，1ステップ前でaccelの計算が済んでいる．
この速度を使って位置の更新．
この位置を使って速度の更新である．

加速度の計算 :math:`\Phi (x)` は，サブルーチンを作って行うことにした．
このためのサブルーチン「force」はモジュール「cal_force」のなかに（contain）含まれている．

.. code-block:: fortran

    module cal_force

    use globals
    use param_planet
    implicit none
    
    contains
        subroutine force(c,a)
            double precision,intent(in)::       c(ndim,nbody)
            double precision,intent(out)::      a(ndim,nbody)
            integer     i,j,k
            double precision    rij(ndim)
            double precision    rij2,rij_abs,inv_rij_abs
            double precision    ftmp
            double precision    fij(ndim),f(ndim,nbody)

            f = 0.0d0
            do i = 1,nbody-1
                do j = i+1,nbody
                    rij(:) = c(:,j)-c(:,i)
                    rij2 = 0.0d0
                    do k=1,ndim
                        rij2 = rij2 + rij(k)**2
                    end do
                    rij_abs = sqrt(rij2)
                    inv_rij_abs = 1.0d0/rij_abs
                    ftmp = mass(i)*mass(j)*(inv_rij_abs**3)
                    fij(:) = ftmp*rij(:)
                    f(:,i) = f(:,i) + fij(:)
                    f(:,j) = f(:,j) - fij(:)
                end do
            end do
            do i = 1,nbody
                a(:,i) = f(:,i) / mass(i)
            end do
        end subroutine force

    end module cal_force

系内にあるすべての天体間の組合せに対して距離を計算し，万有引力を計算している．

.. code-block:: fortran 

    f(:,i) = f(:,i) + fij(:)
    f(:,j) = f(:,j) - fij(:)

天体jから天体iへの作用fijを計算したら，iに作用する力f(:,i)に加える，とともに，
天体iから天体jへの作用（-fij，反作用）としてjに作用する力f(:,j)にも加えている．
その他の部分はよく見れば何を計算しているか分かるだろう．

.. note::
    たとえばdoループを使って，二天体間の距離の二乗rij2を計算しているが，この時ループの直前で「rij2=0.0d0」と値をゼロにするのを忘れないように．これを忘れると，その他の天体間のための計算に用いたrij2の値を引き継いでしまう．

.. note::
    これまでの関数（function）による書き方と似ている．関数は，単一の値を返すのに対して，サブルーチンは引数に値を入れて返す．この場合はc=coordの値を受け取りa=accelに結果を入れて返している．

.. note::
    天体の質量，初期位置，初速度はWikipediaの情報を元に決めた．Wikipediaは理科年表の値などを載せているようである．初期位置はx軸上の平均公転半径に相当する場所，初期速度は平均軌道速度とした．

さて，以下が結果である．

.. figure:: ../img/337XvsY_n_vV.png
   :scale: 100 %
   :align: center
   :alt: X vs Y　of many-body

なんてことはない，それぞれの天体の円（にしか見えない）軌道が描かれる．

よく見ないと分からないが，例えば月と地球の軌道を見比べれば，月が地球の軌道の外側と内側を交互に移動しながら周回する様子も分かる．

あるいは，たとえば金星の速度を20%増しにして長時間の計算をしてみると，金星だけでなく火星の軌道も徐々にずれてくる．

.. figure:: ../img/339XvsY_n_vVm.png
   :scale: 100 %
   :align: center
   :alt: X vs Y　of many-body

.. figure:: ../img/338deltaE_n_VV.png
   :scale: 100 %
   :align: center
   :alt: delta_E of many body sim

一応，エネルギーの初期値からの変動も載せる．
時間刻みは :math:`\Delta t = 0.637` （赤）から :math:`\Delta t = 0.637^{-5}` （黒）まで
速度ベルレ法なので，初期値からの変動は小さく，また刻みを小さくすればするほど小さくなる．