Lecture 5 萬有引力定律與克卜勒行星運動定律

建國高中特色選修課程 - 物理現象的程式設計與模擬
作者：賴奕帆
日期：2018/8/01

一、萬有引力定律

定律內容：牛頓於1687年提出 「任何兩個質點之間必有一吸引力 $F$ ，其量值和兩質量的乘積成正比，和它們之間的距離平方成反比，其方向在連心線上。」

\vec{F} = \frac{G M m}{r^{2}} \hat{r}

$\overrightarrow{F} = \dfrac{GMm}{r^{2}}\widehat{r}$

$r$ ：將兩球體視為質點時的兩球心距離，單位為公尺（ $\rm m$ ）。
$M$ 、 $m$ ：兩質點的質量單位為公斤（ $\rm kg$ ）。
$G$ ：萬有引力常數為 $\rm 6.67\times 10^{-11}\ N\cdot m^{2}/kg^{2}$

這堂課我們要將前三堂課所學的基本程式概念應用於行星運動。
讓我們先來寫一物體受到地球萬有引力作用，改變運動狀態與位置的Vpython程式。
請參考以下範例程式：

Example 1 : 衛星受地球萬有引力吸引落下










































"""
    建國中學 Vpython物理模擬
    作者: 物理科 賴奕帆老師
    特色課程 Lecture 05 萬有引力定律與克卜勒行星運動定律
    5_01_Gravity force.py
"""
from vpython import *  #引用視覺畫套件Vpython
"""
    1. 參數設定，設定變數及定義萬有引力公式
"""
G = 6.67*10**(-11) ; M = 6*10**24 ; m = 1000  
Re = 6.4*10**6 ; H = 20*Re ; t = 0 ; dt = 1
def Fg(x):                                 #定義公式
    return -G*M*m/(x**2)

"""
    2. 畫面設定
"""
scene = canvas(align = 'left',title ='4_01_Gravity force',  width=800, height=300, center=vec(0,0,0), background=vec(0.6,0.8,0.8)) #設定視窗
earth = sphere(pos=vec(0,0,0), radius=Re, texture=textures.earth) #放置物件地球
mater = sphere(pos=vec(H,0,0), radius=0.1*Re,color=color.red, make_trail=True) #放置物件衛星
materv = vec(0,0,0) #衛星速率=0
"""
    3. 執行迴圈
"""
while mater.pos.x-earth.pos.x >=Re:  #執行迴圈
    rate(5000)
    dist = ((mater.pos.x-earth.pos.x)**2+(mater.pos.y-earth.pos.y)**2+(mater.pos.z-earth.pos.z)**2)**0.5 #距離純量
    radiavector = (mater.pos-earth.pos)/dist #距離單位向量
    Fg_vector = Fg(dist)*radiavector # 萬有引力向量=萬有引力量值*單位向量
    
    materv += Fg_vector/m*dt   #Δv = F/m *dt
    mater.pos = mater.pos + materv*dt  # S = S0 + v *dt
  
    t = t+dt
"""
    4.當衛星碰觸地球表面，則衛星停止不動
"""
while True:  
    rate(1000)
    materv = vector(0,0,0)
    mater.pos = mater.pos

可見到一個紅色衛星緩緩受到地球引力吸引而掉落，但光看此動畫難以判斷是否正確，如何確認呢？讓我們加入圖表判觀察吧。

加入圖表觀察萬有引力







""" Fg_x圖(2. 畫面設定)"""
gd = graph(align='left',width=400,height=400,  #設定X-t繪圖視窗
              title='Fg', xtitle='R', ytitle='Fg(red)',
              foreground=color.black,background=color.white,
              xmax=20, xmin=-2, ymax=1.2, ymin=0)
f1 = gcurve(color=color.red)  #定義曲線
Fe = G*M*m/Re**2 #定義地球表面重力強度

Vpython的圖表可以即時觀察物件的各種狀態，包含受力、位置、速度、加速度、能量與時間等等，是一個方便有感覺的圖表工具。

二、克卜勒行星第一定律 (一) 圓形軌道衛星

討論「衛星質量 $m$ 以軌道半徑 $r$ 繞質量 $M$ 的行星作等速率圓周運動」。

$m$ 繞 $M$ 轉所需的向心力來自 $m$ 在該軌道上所受的萬有引力。
因萬有引力恆指向 $M$ 球心，所以其軌道平面通過球心。
$m$ 以萬有引力作為向心力作圓周運動：

F = \frac{G M m}{r^{2}} = m a_{c} = m \frac{v^{2}}{r} = m \frac{4 π^{2} r}{T^{2}}

$F =\dfrac{GMm}{r^{2}} = ma_{c} = m\dfrac{v^{2}}{r} = m\dfrac{4\pi ^{2}r}{T^{2}}$
移項整理後可得：

v = \sqrt{\frac{G M}{r}}

$v = \sqrt{\dfrac{GM}{r}}$ 和

T = 2 π \sqrt{\frac{r}{G M}}

$T = 2\pi \sqrt{\dfrac{r}{GM}}$

衛星 $m$ 的速率 : $v$ 與衛星的質量 $m$ 無關。
衛星 $m$ 的週期 : $T$ 與衛星的質量 $m$ 也無關。

在程式碼中加入衛星速度，觀察衛星運動



""" 加入圓周運動的衛星速度(1. 參數設定)"""
   H = 5*Re #更改衛星位置 
   V0 = _____________ #定義衛星初速度


   materv = vector(0,_______,0) #改變衛星初速度


   while True:  #原本是while mater.x-earth.x >=Re:

將公式中的速率置入程式內，執行出來果然是漂亮的正圓形，不過旁邊的圖表暫時沒用，可以先用Alt+3將程式碼刪除。

三、克卜勒行星第一定律 (二) 橢圓軌道衛星

又稱軌道定律，內容：『所有的行星各以太陽為焦點的橢圓形軌道上運行。』

橢圓的標準畫法：

把一條線的兩端 $F_{1}$ 和 $F_{2}$ 固定，用筆尖貼線拉緊，則筆尖所描繪出的軌跡即為一橢圓。兩固定點稱為焦點。從橢圓上任何一點 $P$ ，到兩焦點的距離之和，即「 $\overline{PF_{1}}+\overline{PF_{2}}$ 」，為一定值。
方程式： $\dfrac{x^{2}}{a^{2}}+\dfrac{y^{2}}{b^{2}}=1$ ， $a=$ 半長軸、 $b=$ 半短軸，其中 $a^{2}=b^{2}+c^{2}$
，橢圓面積 = $\pi ab$ ； (圓面積 = $\pi r^{2}$ )
偏心率 $e$ ： $e = \dfrac{c}{a}<1$ ，決定橢圓的扁平程度，偏心率愈大，軌道愈扁平。
近日點 $r_{min}= a+c = a (1+e)$ ，遠日點 $r_{max}= a-c = a (1-e)$ 。

如何把上節的正圓形軌道改成橢圓形軌道呢？

將Vpython中的正圓形軌道改成橢圓形軌道


   materv = vector(0,0.7*V0,0) #使衛星改為橢圓軌道運動

看起來很像橢圓形，但這軌跡真的遵守橢圓形方程式嗎？
讓我們找衛星運動軌跡的上下左右四個端點，以利得到橢圓形軌跡的半長軸 $a$ 與半短軸 $b$ 。

利用Lecture 2「三點記錄法」找橢圓四個端點


   """ 定義衛星前瞬間位置的名稱 (1. 參數設定)"""
   pre_mater_pos = vector(0,0,0)



   """ 在迴圈中記憶衛星前瞬間與前前瞬間的位置 (3. 執行迴圈，最前面)"""
   pre_pre_mater_pos = pre_mater_pos
   pre_mater_pos = vector(mater.pos.x , mater.pos.y, mater.pos.z)









   """ 找出衛星上下左右四個端點位置 (3.執行迴圈，最後面)"""
   if pre_mater_pos.x > pre_pre_mater_pos.x and pre_mater_pos.x > mater.pos.x :
      print (mater.pos) #印出右端點
   if pre_mater_pos.x < pre_pre_mater_pos.x and pre_mater_pos.x < mater.pos.x : 
      print (mater.pos) #印出左端點
   if ______________ and ______________  :
      print (mater.pos) #印出上端點
   if ______________ and ______________ :
      print (mater.pos) #印出下端點

由Print結果可知，在此條件下，橢圓運動 $x$ 方向最大值為 $3.2\times 10^{7}$ ， $x$ 方向最小值為 $-1.038\times 10^{7}$ ， $y$ 方向最大與最小值分別為 $\pm1.8228\times10^{7}$ 。
半長軸 $a=(3.2\times 10^{7}+1.038\times 10^{7})/2=2.119\times 10^{7}$ ，半短軸 $b=1.8228\times10^{7}$ 。

四、克卜勒行星第一定律 (三) 描繪出橢圓形軌道

上節最後我們找出了，要繪製橢圓形的半長軸 $a$ 與半短軸 $b$ ，讓我們利用vpython中的curve語法來繪製橢圓形吧。

利用curve繪製橢圓，看衛星是否在橢圓軌道上運行




"""利用curve畫橢圓 (2. 畫面設定)"""
oval = curve( color = color.black )
for N in range(0, 360, 1):
    oval.append( pos =(2.119*10**7*cos(N*pi/180)+0*(2.119**2-1.8228**2)**0.5*10**7,1.8228*10**7*sin(N*pi/180),0) )

在紅色衛星運動前，就先繪出了一個黑色橢圓軌道，紅色衛星果然在橢圓軌道上運行。
我們利用Vpython確認了克卜勒行星第一定律的正確性。

五、克卜勒行星第二定律 - 利用海龍公式計算面積

又稱等面積定律，內容：『太陽與行星的連線在相同的時間間隔內掃過相等的面積。』
定律詳細內容：如圖所示，當 $A\rightarrow B$ 與 $A{}'\rightarrow B{}'$ 時間間隔相同時，則行星與太陽連線在A至B的時間內掃過的面積等於在 $A{}'$ 至 $B{}'$ 的時間內掃過的面積。

面積速率公式：

\frac{△ A}{△ t} = \frac{1}{2} r v \sin ϕ = \frac{π a b}{T} = c o n s t .

$\dfrac{\bigtriangleup A}{\bigtriangleup t} = \dfrac{1}{2}rv\sin \phi = \dfrac{\pi ab}{T} = const.$

接下來要用Vpython來驗證『太陽與行星的連線在相同的時間間隔內掃過相等的面積。』是否正確，我們要如何找出每一時刻行星掃過的面積呢？
提示：我們有每一時刻行星在準確位置。

海龍公式：

將海龍公式，寫進前述橢圓軌道的程式之中











""" 計算每瞬間衛星與地球連線掃過的面積 (3.執行迴圈，最後面)"""
    a = mag(mater.pos-earth.pos)  #定義三角形的三邊長
    b = mag(earth.pos-pre_mater_pos)
    c = mag(pre_mater_pos-mater.pos)
    s = _______________        
    area = __________________________ #用海龍公式計算面積
    sum_area += area            #累加每dt時刻內的面積，計算總面積
    if t >= 2785 :
        print (t , sum_area)
        t = 0
        sum_area = 0

結果相當滿意。在模擬時間 $2785$ 秒內，衛星掃過的面積均為 $\rm 1.103\times 10^{14}m^{2}$ 。
但光數字好像感覺還不足，可以用動畫加以呈現嗎？

六、克卜勒行星第二定律 (一) 利用海龍公式計算面積

在if判別內，加入繪製直線


    if t >= 2785 :  #原有的程式碼
        cylinder(pos=____________, axis=________________, radius=0.01*Re, color=color.black)  #每2785秒將地球與衛星之間連線繪圖

我們利用Vpython確認了克卜勒行星第二定律的正確性。
那克卜勒行星第三定律呢 ?

Lecture 5 課堂作業

課堂作業 5-1

利用前述的圖表，討論衛星在橢圓運動中，是否有遵守力學能守恆。
註：在衛星空間尺度中，由於萬有引力不是定值，因此衛星與地球間的位能公式應為

U (r) = - \frac{G M m}{r}

$U(r)=-\frac{GMm}{r}$ 位能的零位面設定為無窮遠處，此公式會在高二物理第七章「功、動能與位能」中介紹。

課堂作業 5-1 參考解法步驟 :

在畫面設定處，更改圖表設定，準備畫出動能、位能與總力學能三條曲線。
在執行迴圈處，計算衛星任意瞬間的動能、位能與總力學能

課堂作業 5-2

再加入兩個衛星，讓三個衛星有不同的初速，觀察整體的運動狀況。

課堂作業 5-2 參考解法步驟 :

在畫面設定處，再加入兩個不同的衛星mater2與mater3，且讓每個衛星的初速度均不相同
在執行迴圈處，分別計算兩個不同的衛星mater2與mater3受太陽的引力，並使開始其運動

課堂作業 5-3

利用三點紀錄法找出衛星運動的週期。

課堂作業 5-3 參考解法步驟 :

在參數設定處，除了原有的t=0外，再加上一個T=0，以計算運動週期
在執行迴圈處，加入T=T+dt，使T一樣會隨著迴圈累加時間。
在前述利用pre找四端點位置處，將其中三個位置用ALT+3暫時刪除後，並修改print，使可印出運轉週期T，並在印出後使T歸零(T=0)

進階作業 5-1

請在衛星運動中，加入速度與加速度(分解為 $a_{T}$ 與 $a_{N}$ )向量的箭頭，我們可以利用功能定理，討論為何近地點(近日點)速率較快，而遠地點(遠日點)速率較慢。

進階作業 5-1 參考解法步驟 :

在畫面設定處，加入兩個向量箭頭，標示衛星速度向量(紅色)與加速度向量(黃色)，並在執行迴圈處，讓兩個向量箭頭與衛星一起運動，讓其長度代表衛星的速度與加速度。(注意為了讓動畫向量顯著，可自行調整量值的大小)
再繪製兩個箭頭，並利用內積與外積的數學概念將加速度分解為 $a_{T}$ 與 $a_{N}$ 。



"""讓切線加速度箭頭與法線加速度箭頭，隨著衛星一起運動(3. 執行迴圈)"""
    aT_arrow.axis = dot(Fg/m,materv)*(norm(materv)/mag(materv))*5*10**6
    aN_arrow.axis = (mag(cross(Fg/m,materv))/mag(materv))*cross(norm(materv),vector(0,0,-1))*5*10**6

進階作業 5-2

請嘗試驗證克卜勒行星第三定律。

克卜勒行星第三定律又稱週期定律。
定律內容：『各行星若繞同一恆星時，其公轉週期的平方與行星到太陽平均距離的三次方成正比。』
物理公式： $\dfrac{a^{3}}{T^{2}}=const.$ (半長軸 $\ a=$ 行星到太陽平均距離 $\ =\dfrac{r_{近}+r_{遠}}{2}$ )
此為克卜勒以大量觀察數據歸納出的，後來牛頓以行星受到萬有引力當作向心力來驗證其結果。

進階作業 5-2 參考解法步驟 :

繪製太陽、火星與地球，火星與地球以圓形軌道繞行太陽運行。
利用課堂作業5-3計時器概念，分別找出火星與地球繞行太陽軌道的週期。
將火星與地球的軌道半徑與週期代入克卜勒行星第三定律的公式中，驗證在正圓形軌道下，克卜勒行星第三定律是否正確。
略為改變火星與地球的初速度，使其改為橢圓形軌道。
利用本節內容，找出橢圓形軌道的半長軸 $a$ 與運行軌道週期 $T$ 。

Lecture 5 作業完整參考解法

同學們，在看參考解之前，請先想清楚自己是否真的有花時間想過如何寫課堂作業和進階作業了呢?
學習無他，但求自己覺得好玩而已。
若真的想把一件事情學好，務必花時間才能掌握。
希望這邊的作業參考解法，真的都只是讓同學們拿來參考，而不是為了繳交作業的複製貼上。
想過了嗎? 那請點下面的網址吧。

# Lecture 5 作業完整參考解法