====== 万有引力の数値計算 ======

===== ロケットの軌道計算 =====

地球上からロケットを発射すると、ロケットは地球の引力の影響を受けて運動する。
宇宙空間は3次元空間だが、地球の引力による運動では2次元平面上で考えることができる。

----
==== Ruby プログラム ====
プログラム orbit_rocket.rb は、地球上から発射した推進力のないロケット（[[https://www.spinlaunch.com|SpinLaunch]]）がどのように運動するかを Curses で端末上に表示する。
初期値によってロケットの軌道が変化する。

**初期値**
  * $(x_{0},y_{0})$：ロケットの発射位置（単位 m（メートル））
  * $(u_{0},v_{0})$：ロケットの x 方向と y 方向の発射速度（単位 m/s（メートル毎秒））
  * $dt$：時間の刻み（単位 s（秒））

{{ :ruby:orbit_rocket.png?nolink |ロケットの初期値}}

<file ruby orbit_rocket.rb[enable_line_numbers="true"]>
#!/usr/koeki/bin/ruby
# -*- coding: utf-8 -*-
# Curses でロケットの運動を描画する
# 作者：山本裕樹

require "curses"

# 万有引力定数 G = 6.67430×10^-11 [m^3/kg/s^2]
# 地球の質量 ME = 5.9724×10^24 [kg]
# 地球の半径 RE = 6.3781x10^6 [m]
$G = 6.67430e-11
$ME = 5.9724e24
$RE = 6.3781e6

# 現実の座標から端末上の座標に変換するための x, y 方向のスケール
$xscale = 10.0
$yscale = 5.0

# 初期値
# ($x0, $y0)：ロケットの発射位置（単位 [m]）
# ($u0, $v0)：ロケットの発射速度（単位 [m/s]）
# $dt：時間の刻み(単位 [s])
$x0 = $RE + 10
$y0 = 0
$u0 = 6000
$v0 = 6000
$dt = 60

# 原点から (x, y) までの距離
def r(x, y)
  return Math.sqrt(x ** 2 + y ** 2)
end

# 運動方程式の右辺1
def f(x, y)
  return - $G * $ME * x / r(x, y) ** 3
end

# 運動方程式の右辺2
def g(x, y)
  return - $G * $ME * y / r(x, y) ** 3
end

# ku, kv, kx, ky の計算
def calc_k(u, v, x, y, du, dv, dx, dy)
  return f(x + dx, y + dy) * $dt,
         g(x + dx, y + dy) * $dt,
         (u + du) * $dt,
         (v + dv) * $dt
end

# ルンゲ=クッタ法で $dt 後の (x,y) と (u,v) を求める
def next_step(u, v, x, y)
  ku1, kv1, kx1, ky1 = calc_k(u, v, x, y, 0, 0, 0, 0)
  ku2, kv2, kx2, ky2 = calc_k(u, v, x, y,
                              ku1 * 0.5, kv1 * 0.5, kx1 * 0.5, ky1 * 0.5)
  ku3, kv3, kx3, ky3 = calc_k(u, v, x, y,
                              ku2 * 0.5, kv2 * 0.5, kx2 * 0.5, ky2 * 0.5)
  ku4, kv4, kx4, ky4 = calc_k(u, v, x, y, ku3, kv3, kx3, ky3)
  return u + (ku1 + 2 * ku2 + 2 * ku3 + ku4) / 6.0,
         v + (kv1 + 2 * kv2 + 2 * kv3 + kv4) / 6.0,
         x + (kx1 + 2 * kx2 + 2 * kx3 + kx4) / 6.0,
         y + (ky1 + 2 * ky2 + 2 * ky3 + ky4) / 6.0
end

# 現実の座標 (x,y) から端末上の座標 (cx,cy) へ変換
def x2cx(x) ($cox + x / $RE * $xscale).round end
def y2cy(y) ($coy - y / $RE * $yscale).round end

# 端末上の座標 (cx,cy) から現実の座標 (x,y) へ変換
def cx2x(cx) (cx - $cox) * $RE / $xscale end
def cy2y(cy) -(cy - $coy) * $RE / $yscale end

# (0, 0) を中心に地球を描く
def draw_earth
  # 地球の色の設定
  Curses.attrset(Curses.color_pair(2))

  # 端末上の全ての座標をスキャンして地球の半径内なら地球を出力
  0.upto(Curses.cols - 1) do |cx|
    0.upto(Curses.lines - 1) do |cy|
      if r(cx2x(cx), cy2y(cy)) < $RE
        Curses.setpos(cy, cx)
        Curses.addch(' ')
      end
    end
  end
end

# 現実の座標 (x,y) から端末上の座標へ変換して文字 ch を出力
def puts_char(x, y, ch)
  cx = x2cx(x)
  cy = y2cy(y)
  if (cx >= 0 && cx <= Curses.cols - 1) &&
     (cy >= 0 && cy <= Curses.lines - 1)
    Curses.setpos(cy, cx)
    Curses.addch(ch)
  end
end

# Cursesの初期化
Curses.init_screen

begin
  # 色の設定
  Curses.start_color
  Curses.init_pair(1, Curses::COLOR_WHITE, Curses::COLOR_BLACK)
  Curses.init_pair(2, Curses::COLOR_BLACK, Curses::COLOR_BLUE)
  Curses.init_pair(3, Curses::COLOR_CYAN, Curses::COLOR_BLACK)
  Curses.init_pair(4, Curses::COLOR_BLACK, Curses::COLOR_RED)

  # 端末の中央の座標
  $cox = ((Curses.cols - 1) * 0.5).round + 0.5
  $coy = ((Curses.lines - 1) * 0.5).round + 0.5

  # カーソルを非表示
  Curses.curs_set(0)

  # キー入力で待つ時間（ミリ秒）
  Curses.stdscr.timeout = 60
  
  # 初期値を代入
  x = $x0
  y = $y0
  u = $u0
  v = $v0

  # 時間
  i = 0

  # ロケットの軌道を描くループ
  while 1
    # 地球を出力
    draw_earth

    # デフォルトの色
    Curses.attrset(Curses.color_pair(1))
    
    # 時刻の出力
    Curses.setpos(Curses.lines - 1, 0)
    Curses.addstr(sprintf("t=%d[s]\n", i * $dt))
    
    if (r(x, y) > $RE)
      # ロケットを出力
      puts_char(x, y, 'R')
    else
      # 爆発を出力
      Curses.attrset(Curses.color_pair(4))
      puts_char(x, y, 'X')
      Curses.refresh
      # ループを抜ける
      break
    end

    # キー入力があるか
    if Curses.getch != nil
      # ループを抜ける
      break
    end

    # (x, y) にロケットの跡を出力
    Curses.attrset(Curses.color_pair(3))
    puts_char(x, y, '.')

    # 画面表示の更新
    Curses.refresh

    # $dt 後の (x,y) を求める
    u, v, x, y = next_step(u, v, x, y)
    i += 1
  end

  # キー入力があるまでずっと待つ
  Curses.stdscr.timeout = -1
  Curses.getch

ensure
  # Cursesの終了処理
  Curses.close_screen
end
</file>

ロケットの発射位置は地上から 10m の高さに設定している（23行目、24行目）。

ロケットの発射速度 $(u_{0},v_{0})$ を色々変えてみよう（25行目、26行目）。

^  No.  ^  $u_{0}$ [m/s]  ^  $v_{0}$ [m/s]  ^  結果                                                                                                  ^
|  1    |            6000 |            6000 | xy方向に発射されたロケットはしばらくは飛ぶが最後は地球に落下する。                                     |
|  2    |               0 |            7906 | y方向に発射されたロケットは地上すれすれを周り続ける。この速度を__**第一宇宙速度**__という。            |
|  3    |               0 |           10000 | y方向に発射されたロケットは地球の周りを回り続ける。                                                    |
|  4    |               0 |           11180 | y方向に発射されたロケットは地球の引力を振り切って戻ってこない。この速度を__**第二宇宙速度**__という。  |

----
==== 理論 ====

地球の中心を原点としたロケットの座標を $(x,y)$ とする。
ロケットは地球の引力を受けて運動するため、$(x,y)$ は時間 $t$ で変化する。

ロケットの速度を $(u,v)$ とすると、$(x,y)$ と以下の関係にある。
\begin{eqnarray*}
\frac{dx}{dt}&=&u\\
\frac{dy}{dt}&=&v
\end{eqnarray*}
ここで $\frac{dx}{dt}$ と $\frac{dy}{dt}$ は $x$ と $y$ をそれぞれ時間 $t$ で微分することを意味する。
 
ロケットは地球の引力を受けて運動するため、__**万有引力の法則**__から次の運動方程式が成り立つ。
\begin{eqnarray*}
m\frac{du}{dt}&=&-\frac{GM_{\mathrm{E}}m}{r^{3}}x\\
m\frac{dv}{dt}&=&-\frac{GM_{\mathrm{E}}m}{r^{3}}y
\end{eqnarray*}
ここで $G$ は万有引力定数、$M_{\mathrm{E}}$ は地球の質量、$m$ はロケットの質量である。
\[
G=6.67430\times 10^{-11}\mathrm{[m^{3}kg^{-1}s^{-2}]}
\]
\[
M_{\mathrm{E}}=5.9724\times 10^{24}\mathrm{[kg]}
\]
両辺を $m$ で割ると $m$ は消えるため、この計算にロケットの質量は不要である。
\begin{eqnarray*}
\frac{du}{dt}&=&-\frac{GM_{\mathrm{E}}}{r^{3}}x\\
\frac{dv}{dt}&=&-\frac{GM_{\mathrm{E}}}{r^{3}}y
\end{eqnarray*}
$r$ は地球の中心からロケットまでの距離で、三平方の定理（ピタゴラスの定理）から次の式で求められる。
\[
r=\sqrt{x^{2}+y^{2}}
\]

以上の式を解くとロケットの座標 $(x,y)$ が時間 $t$ とともにどのように変化するか、ロケットの軌道を求めることができる。

----
==== 数値計算 ====
以下の4つの方程式（__**微分方程式**__という）を解くとロケットの運動の軌道が求められる。
\begin{eqnarray*}
\frac{dx}{dt}&=&u\\
\frac{dy}{dt}&=&v\\
\frac{du}{dt}&=&f(x,y)\\
\frac{dv}{dt}&=&g(x,y)
\end{eqnarray*}
ここで
\begin{eqnarray*}
f(x,y)&=&-\frac{GM_{\mathrm{E}}}{r^{3}}x\\
g(x,y)&=&-\frac{GM_{\mathrm{E}}}{r^{3}}y
\end{eqnarray*}
とした。

微分方程式はオイラー法やルンゲ=クッタ法などを用いて数値的に計算できる。