講義（山本裕樹）

この文書は読取専用です。文書のソースを閲覧することは可能ですが、変更はできません。もし変更したい場合は管理者に連絡してください。
====== 万有引力の数値計算 ======

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

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

----
==== rbCanvas プログラム ====
プログラム rocket_orbit.rb は、地球上から発射した推進力のないロケット（[[https://www.spinlaunch.com|SpinLaunch]]）がどのように運動するかをシミュレーションする。

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

プログラムは [[https://rbcanvas.net/|rbCanvas]] で動作する。
rbCanvas editor に入力して[▶]をクリックしてみよう。

<file ruby rocket_orbit.rb[enable_line_numbers="true"]>
# coding: utf-8
# ロケットの運動を描画する
# 動作環境：rbCanvas
# 作者：山本裕樹

# 画像サイズ
WIDTH = 800
HEIGHT = 800

# 画像上の中心座標
IOX = WIDTH / 2
IOY = HEIGHT / 2

# 空間のスケール [pix/m]
SCALE = 5e-6

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

# 経過時間 [s]
t = 0

# $dt：時間の刻み [s]
$dt = 60

# ロケットの初期値
# (x0, y0)：ロケットの発射位置 [m]
x0 = RE + 50
y0 = 0

# (x, y)：ロケットの空間座標 [m]
x = x0
y = y0

# (u, v)：ロケットの速度 [m/s]
u = 0
v = 0

# ロケットの軌跡
tr = []
trmax = 1000

# 原点から (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) から画像座標 (ix,iy) へ変換
def x2ix(x) return (IOX + x * SCALE).round end
def y2iy(y) return (IOY - y * SCALE).round end

# 背景の設定
Window.width   = WIDTH
Window.height  = HEIGHT
Window.bgcolor = 'black'

# 発射角度のスライダーとラベル
label1_angle = UI.label(WIDTH + 20, 120, 100, 30, '発射角度')
label2_angle = UI.label(WIDTH + 370, 120, 180, 30, '')
slider_angle = UI.slider(WIDTH + 150, 120, 200, [-90, 90, 30, 1])

# 発射速度のスライダーとラベル
label1_speed = UI.label(WIDTH + 20, 170, 100, 30, '発射速度')
label2_speed = UI.label(WIDTH + 370, 170, 180, 30, '')
slider_speed = UI.slider(WIDTH + 150, 170, 200, [5, 15, 10, 0.1])

# 発射ボタンとラベル
$launch = false
button_launch = UI.button(WIDTH + 20, 220, 100, 30, '発射')
button_launch.on_action do
  $launch = true
end

flight = false

Window.loop do

  # 発射時
  if $launch
    # 発射角度 [rad]
    alpha = slider_angle.value.to_f * Math::PI / 180
    # 発射速度 [m/s]
    speed = slider_speed.value.to_f * 1e3
    
    u = speed * Math.cos(alpha)
    v = speed * Math.sin(alpha)
    x = x0
    y = y0
    t = 0
    tr = []
    $launch = false
    flight = true
  end

  # ロケットの空間座標から画像座標への変換
  ix = x2ix(x)
  iy = y2iy(y)

  # 飛行中
  if flight
    # 軌跡の記録
    if tr.length > trmax
      tr.shift
    end
    tr.push([ix, iy])
  end
  
  # ロケットの描画
  Window.draw_circle_fill(ix, iy, 3, 'green')

  # 軌跡の描画
  (1..(tr.length - 1)).each do |i|
    Window.draw_dot(tr[i - 1][0], tr[i - 1][1], 'cyan', weight: 1)
  end
    
  # 地球の描画
  Window.draw_circle_fill(x2ix(0), y2iy(0), RE * SCALE, 'blue')
  
  # ラベルの更新
  label2_angle.value = "#{slider_angle.value}度"
  label2_speed.value = "#{slider_speed.value}km/s"

  # 経過時間の表示
  Window.draw_font(0, 25, sprintf("経過時間:%d秒", t),
                   Font.default, color: 'yellow')
  
  if flight
    # 次のステップ
    u, v, x, y = next_step(u, v, x, y)
    t += $dt

    # 地球に衝突ならストップ
    if r(x, y) < RE
      flight =false
    end
  end
end
</file>

ロケットの発射速度と発射角度を色々変えてみよう。

^  発射速度 [m/s]  ^  発射角度 [度]  ^  結果                                                                                                  ^
|                8 |              90 | y方向に発射されたロケットは地上すれすれを周り続ける。この速度を__**第一宇宙速度**__という。            |
|             11.2 |                 | y方向に発射されたロケットは地球の引力を振り切って戻ってこない。この速度を__**第二宇宙速度**__という。  |

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

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

ロケットの速度を $(v_{x},v_{y})$ とすると、$(x,y)$ と以下の関係にある。
\begin{eqnarray*}
\frac{dx}{dt}&=&v_{x}\\
\frac{dy}{dt}&=&v_{y}
\end{eqnarray*}
ここで $\frac{dx}{dt}$ と $\frac{dy}{dt}$ は $x$ と $y$ をそれぞれ時間 $t$ で微分することを意味する。
 
ロケットは地球の引力を受けて運動するため、__**万有引力の法則**__から次の運動方程式が成り立つ。
\begin{eqnarray*}
m\frac{dv_{x}}{dt}&=&-\frac{GM_{\mathrm{E}}m}{r^{3}}x\\
m\frac{dv_{y}}{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{dv_{x}}{dt}&=&-\frac{GM_{\mathrm{E}}}{r^{3}}x\\
\frac{dv_{y}}{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}&=&v_{x}\\
\frac{dy}{dt}&=&v_{y}\\
\frac{dv_{x}}{dt}&=&f(x,y)\\
\frac{dv_{y}}{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*}
とした。

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