太陽運動 (3)：中心差

在【太陽運動 (2)：黃道與天球赤道，真近點角與平近點角】最後，我們做圖後發現了一個看似sine(正弦函數)的圖形，這也是【太陽運動 (1)：史書中暗藏的正弦函數】中我們意外發現史書中的數字所能畫出的圖形，以當時古代曆法的年代，當然不會知道這張圖代表的意義，但測量與天體運動的原理是不會改變的，我們來使用後代的數學來解釋為什麼看到了這樣子的圖形。

事實上，這張圖描述的即是近代天文物理學二體問題中有名的「中心差」(equation of the center)。

要推導接下來的式子其實步驟很多，以筆者的參考書籍【Methods of celestial mechanics】就花了$17$頁左右的篇幅，就科普文的來說不太適合，真要推導也只會變成抄書而已。筆者在此就簡單介紹幾個重要的步驟，以及會運用到的一些定理，若讀者有興趣可以直接前往閱讀。

我們的目的是要簡化「真近點角減平近點角」，但直接展開不容易，透過一些輔助會簡單得多。因此在介紹前，要先讓大家認識克卜勒方程式(Kepler's equation)，乍聽之下是大家熟知的克卜勒定律：橢圓軌道、角動量守恆、公轉周期平方與軌道半長軸立方成正比，但其實克卜勒方程式指的是另一條有名的方程式，也是一條與描述橢圓軌道有關的方程式：

為方便說明，地球軌道非真實比例

克卜勒在此引入了一個圓心O的輔助圓，假設在想像中的情況下行星以等速率圓周運動繞圓心O公轉，而當行星從「近日點」以等角速度行進至圖中的位置(圖中之「等角速度」)時，在相同時間內，若地球以按照克卜勒定律的方式運行，則會運行到圖中「地球」的位置。從地球向線段(O-近日點)做垂線並延伸至該輔助圓，可以得到點P。圖中的$\angle$地球-太陽-近日點稱為「真近點角」(true anomaly)，$\angle$等角速度-O-近日點稱為「平近點角」(mean anomaly)，$\angle$P-O-近日點則稱為「偏近點角」(eccentric anomaly)。

推導之後可得克卜勒方程式： \[ M = E - e \sin E \]

上式中的$M$是「平近點角」，$E$是「偏近點角」，$e$則是「離心率」即橢圓方程式：$\frac{x^2}{a^2} + \frac{y^2}{b^2} = 1 $ ，其中離心率即為 $ e = \sqrt{1-\frac{b^2}{a^2}} $。

這裡引入了一個「偏近點角」，透過偏近點角，將這些週期函數做傅立葉級數展開會比較容易些，同時我們也需要運用「真近點角$T$」與「偏近點角$E$」的關係：\[ \cos T = \frac{\cos E - e}{ 1 - e \cdot \cos E} \]

有了這幾個式子，我們可以對這些週期函數做傅立葉級數展開，即大家可能熟知的技巧： \[ \begin{align*} &f(x) = a_{0} + 2 \sum_{1}^{\infty} a_{p} \cos px + 2 \sum_{1}^{\infty} b_{p} \sin px \\ &a_{0} = \frac{1}{ 2\pi} \int_{0}^{2\pi} f(x) \,dx \\ &a_{p} = \frac{1}{ 2\pi} \int_{0}^{2\pi} f(x) \cos \,px \,dx \\ &b_{p} = \frac{1}{ 2\pi} \int_{0}^{2\pi} f(x) \sin \,px \,dx , p = +1, +2 ... \infty \\ \end{align*} \]

如果要繼續推導下去讀者可能會看得有點累，建議有興趣的讀者可以直接前往筆者的參考書籍【Methods of celestial mechanics】$60$頁至$76$頁查閱。 若展開到 $e^7$，偏近點角$E$可對平近點角$M$展開成： \[ \begin{align*} E= M &+ \left( e-\frac{1}{8}e^3 + \frac{1}{192}e^5 - \frac{1}{9216}e^7 \right) \sin M \\ &+ \left( \frac{1}{2}e^2 -\frac{1}{6}e^4 + \frac{1}{48}e^6 \right) \sin 2M + \left( \frac{3}{8}e^3 -\frac{27}{128}e^5 +\frac{243}{5120}e^7 \right) \sin 3M \\ &+ \left( \frac{1}{3}e^4 - \frac{4}{15}e^6 \right) \sin 4M + \left( \frac{125}{384}e^5 -\frac{3125}{9216}e^7 \sin 5M \right) \\ &+ \frac{27}{80} e^6 \sin 6M + \frac{16807}{46080}e^7 \sin 7M \\ \end{align*} \] 再根據偏近點角$E$與真近點角$T$的關係，最終真近點角$T$對平近點角$M$若展開到$e^7$則成： \[ \begin{align*} T = M &+ \left( 2e - \frac{1}{4}e^3 + \frac{5}{96}e^5 + \frac{107}{4608}e^7 \right) \sin M \\ &+ \left( \frac{5}{4}e^2 -\frac{11}{24}e^4 + \frac{17}{192}e^6 \right) \sin 2M + \left( \frac{13}{12}e^3 -\frac{43}{64}e^5 +\frac{95}{512}e^7 \right) \sin 3M \\ &+ \left( \frac{103}{96}e^4 - \frac{451}{480}e^6 \right) \sin 4M + \left( \frac{1097}{960}e^5 -\frac{5957}{4608}e^7 \sin 5M \right) \\ &+ \frac{1223}{960} e^6 \sin 6M + \frac{47273}{32256}e^7 \sin 7M \\ \end{align*} \]

因為真實的行星軌道的離心率都非常小接近0，因此上式中 $e^2$之後的數字都會非常小可以忽略不計，也就會變成： \[ T-M = 2e \sin M \] 即「真近點角減平近點角」會非常近似sine(正弦函數)，最大最小值皆為 $2e$。

為方便讀者想像，下圖是以假設地球公轉軌道為離心率0.833的軌道所做出的圖形：

因為離心率很接近1，這時候上面式子各e的次方項就不太能省略了，所以如果將「真近點角減平近點角」對時間做圖會畫出下面的圖形：

但地球真實軌道的離心率僅有0.0167086，很接近正圓，如下圖：

這時候$e^2$之後便能省略，也因此「真近點角減平近點角」對時間做圖就會畫出如上述所說的sine(正弦函數)的圖形了！

也可以看到最大最小值確實是$\pm 0.0167086 \times 2 =\pm 0.0334172 $。

這就是【太陽運動 (1)：史書中暗藏的正弦函數】中，所畫出的正弦函數的來源了！

花了不少時間，我們終於解惑了神祕的圖形是怎麼來的了，詳細的推導過程就請讀者們自行參閱參考書籍了，筆者列的參考書籍僅是其中一種推導過程，筆者看過各種不同的推導方式，當然最終結論都是相同的。

確認以及慢慢可以預測太陽運動的非等速率運動，對於曆法的制定也是非常重要的，除了顯而易見的節氣之外，朔望月也的長度也受此影響，不難想像，在近日點時的冬季，因為地球公轉稍快，月球要回到相同相位的時間需稍久些，因此對太陽運動有更多的認識後朔望月也終於可以慢慢改用「定朔」而非「平朔」了，這部分我們會在日後繼續說明。