.. index:: law of universal gravitation

========================
万有引力
========================

.. toctree::
   :maxdepth: 1

　2つの物体の間には,物体の質量に比例し,2物体間の距離の2乗に反比例する引力が作用する.

.. math:: 
    F = -G\frac{Mm}{r^2}

太陽の周りを回る地球の軌道を計算してみよう．

運動方程式
===============

天体　:math:`i`　が複数の天体　:math:`j`　からの万有引力を受けて運動している場合，運動方程式は以下のように書ける．

.. 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

少々難しい書き方の気がするが，地球（:math:`i=e`)が太陽（:math:`j=s`）の周りをまわっているとすれば，地球の運動方程式は，

.. math:: 
    \frac{d^2 \boldsymbol{x}_e}{dt^2} = Gm_s\frac{\boldsymbol{x}_s- \boldsymbol{x}_e}{|\boldsymbol{x}_s-\boldsymbol{x}_e|^3}

太陽を中心に固定し，　:math:`\boldsymbol{x}_e- \boldsymbol{x}_s = \boldsymbol{r}`　と書けば，

.. math:: 
    \frac{d^2 \boldsymbol{r}}{dt^2} = - \frac{Gm_s}{|\boldsymbol{r}|^2} \frac{\boldsymbol{r}}{|\boldsymbol{r}|}

なじみの書き方になる．
xy座標系で計算することにして，

.. math::
    \frac{d^2 x}{dt^2} = - \frac{Gm_s}{|\boldsymbol{r}|^2} \frac{x}{|\boldsymbol{r}|} \\
    \frac{d^2 y}{dt^2} = - \frac{Gm_s}{|\boldsymbol{r}|^2} \frac{y}{|\boldsymbol{r}|}

である．

.. index:: normalize

無次元化（規格化？）
=============================

さて，運動方程式を数値計算すればよいのだが，万有引力定数 :math:`G = 6.67259 \times 10^{-11} \mathrm{m^3 s^{-2} kg^{-1}}` ，太陽の質量 :math:`m_s = 1.989 \times 10^{30} \mathrm{kg}` 
など非常に小さな数，大きな数が式中に現れている．これをそのまま計算するのは避けたい．
 
| 方程式を無次元化しよう．
| まずは距離を :math:`x^* = x/a`，時間を :math:`t^* = t/T` のように無次元化する． :math:`a` は，太陽と地球の距離（など）， :math:`T`　は後で決める．
 
すると，
  
.. math::
    \frac{dx}{dt}  = \frac{dx}{dt^*}\frac{dt^*}{dt}=\frac{1}{T}\frac{dx}{dt^*}=\frac{a}{T}\frac{dx^*}{dt^*}

同様に，

.. math::
    \frac{d^2x}{dt^2} =\frac{d}{dt}\left( \frac{a}{T}\frac{dx^*}{dt^*}\right)=\frac{d}{dt^*}\frac{dt^*}{dt}\left( \frac{a}{T}\frac{dx^*}{dt^*}\right)=\frac{a}{T^2}\frac{d^2 x^*}{d{t^*}^2}

これらを用いて運動方程式を書きなおすと，

.. math:: 
    \frac{a}{T^2}\frac{d^2x^*}{{dt^*}^2}=- \frac{Gm_s}{|a^2\boldsymbol{r}^*|^2} \frac{ax^*}{|a\boldsymbol{r}^*|} \\
    \frac{a^3}{Gm_sT^2}\frac{d^2x^*}{{dt^*}^2}=- \frac{1}{|\boldsymbol{r}^*|^2} \frac{x^*}{|\boldsymbol{r}^*|}

右辺からは単位をもった係数がなくなっている．距離がうまく無次元化できているならば大きさは1程度である．左辺の係数も大きさが1程度になるように選びたい．そこで，

.. math::
    \frac{a^3}{Gm_sT^2} = 1

になるように :math:`T` を決める．すなわち，

.. math::
    T = \sqrt{\frac{a^3}{Gm_s}}

とする．すると結局運動方程式は，

.. math::
    \frac{d^2x^*}{{dt^*}^2}=- \frac{1}{|\boldsymbol{r}^*|^2} \frac{x^*}{|\boldsymbol{r}^*|}

となる．元の式と比べると，大きな（あるいは小さな）値をもった係数がなく簡潔である．

計算しよう
========================

ボール投げの時と同様に，速度の変数を使う．（無次元変数を示す「*」は省略）

.. math::
    \dot x &= u, \dot u &= -\frac{x}{|\boldsymbol{r}^*|^3} \\
    \dot y &= v, \dot v &= -\frac{y}{|\boldsymbol{r}^*|^3}

ルンゲ－クッタ法を使って計算する．以下にコードを示す．

.. code-block:: fortran
    :linenos:
    
    !+++++++++++++++++++++++++++++++++++++++++++++++
    module globals

    double precision,parameter::    G = 6.67384d-11
    double precision                delta_t

    end module globals
    !+++++++++++++++++++++++++++++++++++++++++++++++
    module subfunc

    implicit none

    contains
        function fx(x,y,u) result(xx)

            double precision,intent(in)::    x,y,u
            double precision        xx
            
            xx = u
            
        end function fx
        function fu(x,y,u) result(xx)

            double precision,intent(in)::    x,y,u
            double precision        r,xx

            r = sqrt(x**2+y**2)
            xx = -x / r**3
            
        end function fu
        function fy(x,y,u) result(xx)

            double precision,intent(in)::    x,y,u
            double precision        xx
            
            xx = u
            
        end function fy
        function fv(x,y,u) result(xx)

            double precision,intent(in)::    x,y,u
            double precision        r,xx
            
            r = sqrt(x**2+y**2)
            xx = -y / r**3
            
        end function fv
    end module subfunc
    !+++++++++++++++++++++++++++++++++++++++++++++++
    program celestial
    
    use globals
    use subfunc

    implicit none

    integer                     i
    integer,parameter::         n=100
    double precision            tmax
    double precision            t(n),x(n),y(n),u(n),v(n)
    double precision            E(n)
    double precision            kx1,kx2,kx3,kx4,ku1,ku2,ku3,ku4
    double precision            ky1,ky2,ky3,ky4,kv1,kv2,kv3,kv4

    x(1) = 1.0d0
    y(1) = 0.0d0
    u(1) = 0.0d0
    v(1) = 1.0d0 
    t(1) = 0.0d0
    tmax = 10.0d0

    delta_t = (tmax-t(1))/dble(n)

    do i = 1, n-1
        kx1 = fx(x(i),y(i),u(i))*delta_t
        ky1 = fy(x(i),y(i),v(i))*delta_t
        ku1 = fu(x(i),y(i),u(i))*delta_t
        kv1 = fv(x(i),y(i),v(i))*delta_t
        kx2 = fx(x(i)+0.5d0*kx1,y(i)+0.5d0*ky1,u(i)+0.5d0*ku1)*delta_t
        ky2 = fy(x(i)+0.5d0*kx1,y(i)+0.5d0*ky1,v(i)+0.5d0*kv1)*delta_t
        ku2 = fu(x(i)+0.5d0*kx1,y(i)+0.5d0*ky1,u(i)+0.5d0*ku1)*delta_t
        kv2 = fv(x(i)+0.5d0*kx1,y(i)+0.5d0*ky1,v(i)+0.5d0*kv1)*delta_t
        kx3 = fx(x(i)+0.5d0*kx2,y(i)+0.5d0*ky2,u(i)+0.5d0*ku2)*delta_t
        ky3 = fy(x(i)+0.5d0*kx2,y(i)+0.5d0*ky2,v(i)+0.5d0*kv2)*delta_t
        ku3 = fu(x(i)+0.5d0*kx2,y(i)+0.5d0*ky2,u(i)+0.5d0*ku2)*delta_t
        kv3 = fv(x(i)+0.5d0*kx2,y(i)+0.5d0*ky2,v(i)+0.5d0*kv2)*delta_t
        kx4 = fx(x(i)+kx3,y(i)+ky3,u(i)+ku3)*delta_t
        ky4 = fy(x(i)+kx3,y(i)+ky3,v(i)+kv3)*delta_t
        ku4 = fu(x(i)+kx3,y(i)+ky3,u(i)+ku3)*delta_t
        kv4 = fv(x(i)+kx3,y(i)+ky3,v(i)+kv3)*delta_t

        x(i+1) = x(i) + (kx1+2.0d0*kx2+2.0d0*kx3+kx4)/6.0d0
        u(i+1) = u(i) + (ku1+2.0d0*ku2+2.0d0*ku3+ku4)/6.0d0
        y(i+1) = y(i) + (ky1+2.0d0*ky2+2.0d0*ky3+ky4)/6.0d0
        v(i+1) = v(i) + (kv1+2.0d0*kv2+2.0d0*kv3+kv4)/6.0d0

        t(i+1) = t(i) + delta_t
    end do 

    E = 0.5d0*(u**2 + v**2) - 1.0d0 / dsqrt(x**2+y**2)

    open(unit=10,file='outputRK.dat',status='replace')
    write(10,*)'t_rk, x_rk, y_rk, E_rk'
    do i = 1,n
        write(10,*)t(i),',',x(i),',',y(i),',',E(i)
    end do
    close(10)

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

.. warning::
    モジュールglobalはいらない．関数fuとfv内でのrの計算が重複している．関数の変数からtを（使わないので）削除した．またuだけでなくu,vとしたほうがよいのでは（使わないけど）．

座標(1.0,0.0)，速度(0.0,1.0)の初期条件で10秒(時間刻みは10/100=0.1)の計算を行っている．
ボール投げ問題と比べて関数fxとかfuの中身が少々変わっているだけで，その他はほとんど変更なし．

.. figure:: ../img/321XvsY_cl_RK.png
   :scale: 100 %
   :align: center
   :alt: X vs Y　of celestial body

| 軌道をプロットしてみると，円である．
| この問題は解析的に軌道を求めることができる．極座標( :math:`r, \theta` )表示で，初期条件が座標(1.0,0.0)，速度(0.0,v\ :sub:`0`\)の場合，

.. math::
    r = \frac{{v_0}^2}{1-(1-{v_0}^2)\cos\theta}

| である．ちょうど :math:`|v_0|=1.0` のとき円軌道， :math:`|v_0|=\sqrt{2}` のとき放物線軌道， :math:`|v_0|>\sqrt{2}` のとき双曲線軌道を描く．
| せっかくなので，試してみよう．コード68行目「v(1) = 1.0d0」の値を変えればよい．
| 以下は， :math:`|v_0|= 0.8,1.0,1.2,\sqrt{2},1.6` とした時の結果である．

.. figure:: ../img/322XvsY_cl_RK_V0.png
   :scale: 100 %
   :align: center
   :alt: v0 dependence of celestial body