前言

這篇文章的內容最早是在我高一入學軍訓(2019年9月)午休預習高中地理時的一些奇思妙想，在一段時間的思考後的確收獲瞭不少的知識：

高一上學期獨立發現瞭日出日落時間計算公式(本片文章未提及，但讀者可以很容易地推導)，太陽赤經與黃經的換算公式，開普勒方程(幾何法)，並學習瞭球面三角基礎；後來高一下學期對微積分、線性代數等數學知識有瞭一定的基礎後，推導出瞭時間關於黃經的積分表示；高二上學期瞭解到瞭Bessel函數給出的開普勒方程級數解；高二下學期的一次偶然有興趣利用它來計算學校對面山的高差....

高中這麼幾年來反復思考這個課題，致使這篇文章也改版瞭很多次，文章質量自然也隨之不斷上漲。這篇文章首發於疫情期間(2020年4月)，之後在2021年1月有過大的改版。這次(22年元旦)也是一次改版，融入瞭很多新的內容，主要包括一下內容：

• 天體測量學的入門概念
• 開普勒三定律的推導
• 開普勒方程
• 測量學校對面山的高差

1.一點地理知識鋪墊

剛進入高中的時候，在地理課上學習瞭太陽直射點緯度的概念，老師的ppt上出現瞭類似於這樣的一張圖片：

太陽直射點的回歸運動

我就比較好奇這函數曲線的表達式會是什麼樣的(￣▽￣)"，畢竟當時數學課也正在學習函數嘛，看到函數圖像老想著去求表達式(●ˇ∀ˇ●)

不妨設太陽直射點維度記為 delta (以下均稱為太陽赤緯)，以春分日起的天數記為 tau ，那麼該圖對應的函數即為 delta=delta(tau) 。

既然要計算 delta(tau) ，那麼就得理解 delta 為何要隨著時間 tau 變化，對吧？(￣︶￣*))

在地理課上我們學習到：由於地球繞太陽公轉，這會導致太陽直射點在南北回歸線上回歸運動，其運動周期為一個回歸年^[1]。

為瞭比較好理解為什麼太陽直射點會在南北回歸線之間運動，我們不妨引入幾個簡單的概念：

可以參照下圖1來理解這幾個概念

9279ed0bd5ede772390fc03f0c55fb79圖1 圖源自百度

相信你已經理解它們的意思瞭吧||ヽ(*￣▽￣*)ノ，既然理解瞭那就重新回到我們的思路。我們現在的目的是理解：為何一個回歸年中，太陽直射點會在南北回歸線之內回歸運動？

既然黃道是地球公轉軌道在天球上的投影，那麼根據相對運動，太陽在天球上的投影一定也位於黃道上(可以參考上圖)，對吧？(｡･∀･)ﾉ。那麼地球在一個回歸年中繞太陽公轉，對應的就是太陽在天球上的投影沿著黃道運動一周。

所以不難思考：太陽在天球上的投影(即太陽直射點)在一個回歸年中是沿著黃道自西向東運動一周的。

其中與黃道相切的兩個平行於天赤道的小圓就稱為南北回歸線，所以在一個回歸年中太陽直射點會在南北回歸線之內回歸運動。是不是很好理解呢！(*≧︶≦)

既然太陽直射點(一下簡稱為點 S )會運動，那麼現在我們就有必要引入坐標來描述它的運動瞭。球面上最常見的坐標莫屬經緯度瞭：

地理教科書又告訴我們，黃道和天赤道的二面角稱為黃赤交角，其數值 varepsilonapprox 23°26’ ，且南北回歸線為於平行於天赤道並於黃道相切的小圓，那麼可見南北回歸線的赤緯自然就是 pm varepsilon 。

而點 S 的運動是沿著黃道自西向東的，所以我們隻需要知道 S 與春分點 gamma 關於天球球心 O 的夾角 lambda:=angle gamma OS (如下圖2)即可完全掌握 S 的運動。

ed3624437b3134385b8bb2be308bc50f圖2 Geogebra繪制(後同)

顯見 lambda 是關於 tau 的一個單調函數 lambda=lambda(tau) 。且根據球面上一點幾何關系，可以確定 lambda 與 delta 的換算關系。所以為求得 delta(tau) ，我們得先求得 lambda 與 delta 的換算方程，以及 lambda(tau) 。

2.太陽赤緯與太陽黃經的換算方程

參照圖2，我們可以以點 O 為原點，天赤道所在平面為 x-y 平面，射線 Ogamma 為 x 軸正方向建立右手系 {x,y,z} 。(如下圖3)

圖3

因為我們研究的問題是天球上的投影，與天球本身半徑並無關系，所以為簡便運算，不妨定義天球的半徑為單位 1 。

考慮坐標變換：使 y,z 軸繞 x 軸逆時針旋轉 varepsilon ，得到新的坐標系 {x,y',z'} ，使得 z' 軸垂直於黃道面。那麼自然有：

begin{pmatrix} coslambda\ sinlambda\ 0\ end{pmatrix}= begin{pmatrix} 1&0& 0 \ 0&cosvarepsilon&sinvarepsilon \ 0&-sinvarepsilon&cosvarepsilon\ end{pmatrix} begin{pmatrix} cosdeltacosalpha \ cosdeltasinalpha \ sindelta\ end{pmatrix}\

於是我們得到：

begin{align} sinlambda&=cosvarepsiloncosdeltasinalpha+sinvarepsilonsindeltatag{1}\ 0&=-sinvarepsiloncosdeltasinalpha+cosvarepsilonsindeltatag{2}\ end{align}\

對公式(1)同乘一個 sinvarepsilon後有： sinvarepsilonsinlambda=sinvarepsiloncosvarepsiloncosdeltasinalpha+sin^2varepsilonsindelta\

再將公式(2)代入此式得： sinvarepsilonsinlambda=cosvarepsilonleft(cosvarepsilonsindeltaright)+sin^2varepsilonsindelta=sindelta\

所以我們得到瞭 lambda,delta 的換算方程^{[3]： color{blue}{sindelta=sinvarepsilonsinlambdatag{3}\}}

接下來就還剩下求解 lambda(tau) 的問題瞭。

//馬上要進入物理領域瞭喲，同學們系好安全帶// (/▽＼)

註意到 lambda 這個角度與地球繞太陽公轉的徑角是相等的，如圖4所示

圖4 其中γ代表春分日地球經過地點，P代表地球所在位置，S代表太陽

可見，想要求得 lambda(tau) ，不牽扯到萬有引力或者開普勒定律是不可能的瞭。

開普勒第二定律無非是一個角動量守恒，作圖如下：

圖源自百度

取空間內一點 O 作為參考點，定義太陽(質量為 m_{1} )、地球(質量為 m_{2} )的位矢分別為 boldsymbol{x}_{1},boldsymbol{x}_{2} ，且 boldsymbol{r}:=boldsymbol{x}_{2}-boldsymbol{x}_{1}\

我們可以發現，若以太陽為參考系，那麼有： begin{align} &boldsymbol{ddot{r}}=boldsymbol{ddot{x}}_{2}-boldsymbol{ddot{x}}_{1}=left(-frac{Gm_{1}}{r^3}boldsymbol{r}right)-left(frac{Gm_{2}}{r^3}boldsymbol{r}right)=-frac{Gleft(m_{1}+m_{2}right)}{r^3}boldsymbol{r}\[2ex] &boldsymbol{rtimesboldsymbol{ddot{r}}}=boldsymbol{rtimes}left(-frac{G(m_{1}+m_{2})}{r^3}boldsymbol{r}right)=0\[2ex] Rightarrow&frac{text d}{text dt}left(boldsymbol{rtimesboldsymbol{dot{r}}}right)=underbrace{boldsymbol{dot{r}timesdot{r}}}_{0}+boldsymbol{rtimesddot{r}}=0\[2ex] end{align}\

我們可以定義如下守恒量 boldsymbol{h}:=boldsymbol{rtimesdot{r}}\

這是萬有引力為有心力^{[4]的必然結果。在此守恒量的基礎上，很容易導出開普勒定律，不想瞭解開普勒定律如何推導的也不要緊，隻需知道下文公式(*)為一個守恒量即可跳過一下選讀內容[5]：}

3.開普勒三定律簡單推導[選讀]

取極坐標 {r,theta} ，極點為 S ，極軸為近日點所在橢圓長軸射線。可以定義兩個正交歸一基底^{[6]： boldsymbol{hat{r}}:=frac{partial}{partial r}~,~boldsymbol{hat{theta}}:=frac{1}{r}frac{partial}{partialtheta}\}

再取一平面直角坐標系 {x,y} ，原點位於 S ， x 正半軸與極軸相同，同樣可以定義正交歸一基底： boldsymbol{hat{i}}=frac{partial}{partial x}~,~boldsymbol{hat{j}}=frac{partial}{partial y}\

那麼根據基底變換式^[7]：

frac{partial}{partial x^{mu}}=frac{partial {x'}^{nu}}{partial x^{mu}}frac{partial}{partial {x'}^{nu}}~,~~~~small{mu,nu=1,2}\

我們可以將極坐標基底展開為直角坐標基底： begin{align} boldsymbol{hat{r}}:&=frac{partial}{partial r}=frac{partial x}{partial r}frac{partial}{partial x}+frac{partial y}{partial r}frac{partial}{partial y}=frac{partial}{partial r}left(rcosthetaright)boldsymbol{hat{i}}+frac{partial}{partial r}left(rsinthetaright)boldsymbol{hat{j}}\[2ex] &=costhetaboldsymbol{hat{i}}+sinthetaboldsymbol{hat{j}}\[2ex] boldsymbol{hat{theta}}:&=frac{1}{r}frac{partial}{partial theta}=frac{1}{r}left(frac{partial x}{partial theta}frac{partial}{partial x}+frac{partial y}{partial theta}frac{partial}{partial y}right)=frac{1}{r}left(frac{partial}{partial theta}left(rcosthetaright)boldsymbol{hat{i}}+frac{partial}{partial theta}left(rsinthetaright)boldsymbol{hat{j}}right)\[2ex] &=-sinthetaboldsymbol{hat{i}}+costhetaboldsymbol{hat{j}}\ end{align}\ 那麼求導有：

begin{align} &boldsymbol{dot{hat{r}}}=-dot{theta}sinthetaboldsymbol{hat{i}}+dot{theta}costhetaboldsymbol{hat{j}}=dot{theta}boldsymbol{hat{theta}}\[2ex] &boldsymbol{dot{hat{theta}}}=-dot{theta}costhetaboldsymbol{hat{i}}-dot{theta}sinthetaboldsymbol{hat{j}}=-dot{theta}boldsymbol{hat{r}}\[2ex] end{align}\

計算

begin{align} boldsymbol{h}&=boldsymbol{rtimesdot{r}}=boldsymbol{rtimes}frac{text d}{text dt}left(rboldsymbol{hat{r}}right)=boldsymbol{rtimes}left(dot{r}boldsymbol{hat{r}}+rdot{theta}boldsymbol{hat{theta}}right)\[2ex] &=r^2dot{theta}left(boldsymbol{hat{r}timeshat{theta}}right)=r^2dot{theta}left(left(costhetaboldsymbol{hat{i}}+sinthetaboldsymbol{hat{j}}right)boldsymbol{times}left(-sinthetaboldsymbol{hat{i}}+costhetaboldsymbol{hat{j}}right)right)\[2ex] &=r^2dot{theta}left(boldsymbol{hat{i}timeshat{j}}right)tag{4}\[2ex] end{align}\

顯然在已確定的地球軌道(即已知長半軸長 a 以及離心率 e )上，地球質量 m 對lambda(tau) 的是沒有任何影響的。由此，根據上式(4)我們不妨定義一個新的守恒量： h:=r^2dot{theta}tag{*}\

(1)開普勒第二定律：frac{text d}{text dt}left(frac{1}{2}r^2dot{theta}right)=frac{1}{2} frac{text dh}{text dt}=0\

(2)開普勒第一定律^{[8]：對引力加速度(定義 mu:=G(m_{1}+m_{2}) 與 boldsymbol{h} 做叉積有begin{align} &left(ddot{boldsymbol{r}}+frac{mu}{r^3}boldsymbol{r}right)boldsymbol{times h}=0\[2ex] Rightarrow~&left(frac{text ddot{boldsymbol{r}}}{text dt}boldsymbol{times h}+underbrace{dot{boldsymbol{r}}boldsymbol{times}frac{text dboldsymbol{h}}{text dt}}_{0}right)+frac{mu}{r^3}boldsymbol{rtimes}left(boldsymbol{rtimes}dot{ boldsymbol{r}}right)=0\[2ex] Rightarrow~&frac{text d}{text dt}left(dot{boldsymbol{r}}boldsymbol{times h}right)+frac{mu}{r^3}left(boldsymbol{r}left(boldsymbol{rcdot}dot{boldsymbol{r}}right)-dot{boldsymbol{r}}left(boldsymbol{rcdot r}right)right)=0\[2ex] Rightarrow~&frac{text d}{text dt}left(dot{boldsymbol{r}}boldsymbol{times h}right)+frac{mu}{r^3}left(rdot{r}boldsymbol{r}-r^2dot{boldsymbol{r}}right)=0\[2ex] Rightarrow~&frac{text d}{text dt}left(dot{boldsymbol{r}}boldsymbol{times h}right)+muleft(frac{dot{r}boldsymbol{r}-dot{boldsymbol{r}r}}{r^2}right)=0\[2ex] Rightarrow~&frac{text d}{text dt}left(dot{boldsymbol{r}}boldsymbol{times h}-mufrac{boldsymbol{r}}{r}right)=0\[2ex] end{align}\}

所以我們可以定義一個守恒量，稱為龍格-楞次矢量(Runge-Lenz vector)： color{red}{boldsymbol{B}:=dot{boldsymbol{r}}boldsymbol{times h}-mufrac{boldsymbol{r}}{r}tag{5}\}

再將該矢量與 boldsymbol{r} 點乘有：begin{align} &bm{r}cdot(bm{dot{r}}timesbm{h})-mufrac{bm{r}cdotbm{r}}{r}=bm{r}cdotbm{B}\[2ex] Rightarrow~&(bm{r}timesbm{dot{r}})cdotbm{h}-mu r=bm{r}cdotbm{B}\[2ex] Rightarrow~&h^2=rBcos{theta}+mu r \[2ex] Rightarrow~&r=frac{h^2/mu}{1+left(B/muright)cos{theta}}\[2ex] end{align}\。(3)開普勒第三定律：對比橢圓極坐標方程以及軌跡方程，自然有下式成立： color{blue}{frac{b^2}{a}=frac{h^2}{mu}tag{6}\}

既然地球的運動軌跡為橢圓，那麼結合第二定律可以知道有如下方程成立： frac{pi ab}{T}=frac{1}{2}htag{7}\

其中 T 為一個恒星年， pi ab 代表橢圓面積。那麼可以驗證： frac{a^3}{T^2}=a^3left(frac{h}{2pi ab}right)^2=frac{a}{4pi^2b^2}left(frac{mu b^2}{a}right)=frac{mu}{4pi^2}\

其中第二個等號用到公式(6)。

註意到 theta 與 lambda 實際上隻相差一個常數 theta_{gamma} : theta=theta_{gamma}+lambda\

這個 theta_{gamma} 就是春分日地球所在位置對應的極角(圖4中點 gamma 的極角)，目前它大約為 theta_{gamma}approx76°40′28.09″ 。

4.開普勒方程：太陽黃經與時間的換算

為計算 lambda(tau) 隻需要計算 theta(tau) 即可。註意到守恒量 h 含有 dot{theta} ,可以考慮分離變量積分(不要忘記把時間單位制換為天數，弧度制改為角度制哦~)： begin{align} h&=r^2dot{theta}\[2ex] Rightarrow~tau&=int_{theta_{gamma}}^{theta_{gamma}+lambda}frac{r^2}{h}text dtheta=int_{theta_{gamma}}^{theta_{gamma}+lambda}frac{b^4/a^2h}{left(1+ecosthetaright)^2}text dtheta\[2ex] &=int_{theta_{gamma}}^{theta_{gamma}+lambda}frac{left(b^4/a^2right)left(T/360 abright)}{left(1+ecosthetaright)^2}text dtheta= color{red}{frac{left(1-e^2right)^{frac{3}{2}}T}{360}int_{theta_{gamma}}^{theta_{gamma}+lambda}frac{text dtheta}{left(1+ecosthetaright)^2}}tag{8}\[2ex] end{align}\

其中用到橢圓方程以及公式(7)。

可見，公式(3)和公式(8)即為 delta(tau) 曲線的參數方程。由參數方程繪制 delta(tau) 曲線如下：

圖5 太陽直射點緯度關於時間的函數圖像(節選) 圖源自Desmos

雖說是把參數方程得到瞭，但想要得到 delta(tau) 的表達式還是需要比較強的數學能力，下面說明一下：

首先解決積分： begin{align} &intfrac{text dtheta}{(1+ecostheta)^2}=frac{1}{e^2}intfrac{e^2cos^2theta}{(1+ecostheta)^2}text dtheta+frac{1}{e^2}intfrac{e^2sin^2theta}{(1+ecostheta)^2}text dtheta\[2ex] &~~~~~~~~~~~~~~~~~~~~~~~~~~~~=frac{1}{e^2}intfrac{e^2cos^2theta}{(1+ecostheta)^2}text dtheta+frac{1}{e^2}int esinthetatext dleft(frac{1}{1+ecostheta}right)\[2ex] &~~~~~~~~~~~~~~~~~~~~~~~~~~~~=frac{1}{e^2}intfrac{e^2cos^2theta}{(1+ecostheta)^2}text dtheta-frac{1}{e^2}intfrac{ecostheta}{1+ecostheta}text dtheta+frac{sintheta}{eleft(1+ecosthetaright)}\[2ex] &~~~~~~~~~~~~~~~~~~~~~~~~~~~~=-frac{1}{e^2}intfrac{ecostheta}{(1+ecostheta)^2}text dtheta+frac{sintheta}{eleft(1+ecosthetaright)}\[2ex] &~~~~~~~~~~~~~~~~~~~~~~~~~~~~=-frac{1}{e^2}intfrac{text{d}theta}{1+ecostheta}+frac{1}{e^2}intfrac{text dtheta}{(1+ecostheta)^2}+frac{sintheta}{eleft(1+ecosthetaright)}\[2ex] &~~~~~~~~~~~~~~~~~~~~~~~~~~~~=frac{1}{1-e^2}intfrac{text dtheta}{1+ecostheta}-frac{esintheta}{(1-e^2)(1+ecostheta)}\[2ex] &intfrac{text dtheta}{1+ecostheta}=intfrac{2text dtheta}{(e+1)(1+ecostheta)-(e-1)(1-ecostheta)}\[2ex] &~~~~~~~~~~~~~~~~~~~~~~~~=frac{2}{1+e}intfrac{text dleft(tancfrac{theta}{2}right)}{1+cfrac{1-e}{1+e}cfrac{1-ecostheta}{1+ecostheta}}{xlongequal{0leq thetaleqpi}}frac{2}{sqrt{1-e^2}}intfrac{text dleft(sqrt{cfrac{1-e}{1+e}}tancfrac{theta}{2}right)}{1+left(sqrt{cfrac{1-e}{1+e}}tancfrac{theta}{2}right)^2}\[2ex] &~~~~~~~~~~~~~~~~~~~~~~~~=frac{2}{sqrt{1-e^2}}tan^{-1}left(sqrt{frac{1-e}{1+e}}tancfrac{theta}{2}right)+c\[2ex] end{align}\

那麼根據公式(8)，若令 tau_{gamma}:=frac{left(1-e^2right)^{frac{3}{2}}T}{360}int_{0}^{theta_{gamma}}frac{text dtheta}{(1+ecostheta)^2}approx75.8231\

則有： color{blue}{frac{2pi left(tau+tau_{gamma}right)}{T}=(1-e^2)^{frac{3}{2}}int_{0}^{theta}frac{text domega}{(1+ecosomega)^2}=2tan^{-1}left(sqrt{frac{1-e}{1+e}}tanfrac{theta}{2}right)-frac{esqrt{1-e^2}sintheta}{1+ecostheta}}\

此時令

M=frac{2pileft(tau+tau_{gamma}right)}{T}~,~E=2tan^{-1}left(sqrt{frac{1-e}{1+e}}tanfrac{theta}{2}right)\

自然有

begin{align} sin E&=frac{2sincfrac{E}{2}coscfrac{E}{2}}{cos^2cfrac{E}{2}+sin^2cfrac{E}{2}}=frac{2tancfrac{E}{2}}{1+tan^2cfrac{E}{2}}=frac{2(1+e)tancfrac{E}{2}}{(1+e)+(1-e)cfrac{1+e}{1-e}tan^2cfrac{E}{2}}\[2ex] &=frac{2sqrt{1-e^2}left(sqrt{cfrac{1+e}{1-e}}tancfrac{E}{2}right)}{1+left(sqrt{cfrac{1+e}{1-e}}tancfrac{E}{2}right)^2+eleft(1-left(sqrt{cfrac{1+e}{1-e}}tancfrac{E}{2}right)^2right)} =frac{sqrt{1-e^2}left({2tancfrac{theta}{2}}right)}{left(1+tan^2cfrac{theta}{2}right)+eleft(1-tan^2cfrac{theta}{2}right)}\[2ex] &=frac{sqrt{1-e^2}sintheta}{1+ecostheta}\[2ex] end{align}\

所以得到瞭橢圓形式的開普勒方程(Kepler's equation)^[9]:

bbox[tan]{boxed{E-esin E=Mtag{9}}}\

想要求得 lambdaleft(tauright) 隻需求解該方程即可(我的能力到這瞭...( ＿ ＿)ノ｜)，得到 lambda(tau) 後代入式(3)即可得到 delta(tau) 。至於開普勒方程的解析解有很多種形式的，下面這篇文章總結瞭很多，大傢感興趣可以閱讀一下：

醬紫君：x = cos x 的解析形式

其中我個人比較喜歡貝塞爾級數解結合多重對數函數得到的的取整積分形式^[10]：

boxed{E=int_0^pi left(leftlfloor{phi-esinphi+Mover2pi}rightrfloor-leftlfloor{phi-esinphi-Mover2pi}rightrfloorright)text dphitag{10}}\

醬紫君的文中提及到的這篇文獻(An exact analytical solution of Kepler's Equation)更是將 E 的定義延拓至復數集，給出瞭更一般的解法，這篇知乎回答^{[11]也給該論文給出瞭詳細講解：}

農夫三拳有點疼：方程 cosx＝x 是否有解析解？

5.地平坐標系的參數的關系與太陽方位測量

太陽高度 h 在一天中達到最大值時稱太陽上中天(正午)，達到最小值時稱太陽下中天(子午)。

對於非地球兩極的觀測點，我們稱過兩次中天太陽所處位置的天球上的大圓稱為子午線。子午線與地平線的交點為南北點，靠近北方的點稱為北點，靠近南方的點稱為南點。

過天體所在位置和天頂的半大圓與地平線的交點到北點的大圓弧度稱為該天體此時的方位角 A 。

由坐標 (A,h) 刻畫天體位置的坐標系就稱為地平坐標系。如下圖：

圖6 地平坐標系與太陽周日視運動示意圖，以夏季北半球中緯度為例

其中地理緯度為 phi ，角度 t 稱為時角。空間直角坐標系 {x,y,z} 原點位於天球圓心，沿正北方位 x 軸正方向，正東方為 y 軸正方向，天底方向為 z 軸正方向建立。 boldsymbol{r} 為太陽位矢， boldsymbol{n} 為下中天位矢，boldsymbol{d} 為太陽視運動軌跡圓的圓心位矢。

那麼在此直角坐標系下各向量可表示如下： begin{align} &boldsymbol{r}=begin{pmatrix}cos hcos A&cos hsin A &-sin hend{pmatrix}^{mathrm{T}}\[2ex] &boldsymbol{n}=begin{pmatrix}sin(phi+delta)&0&cos(phi+delta)end{pmatrix}^{mathrm{T}}\[2ex] &boldsymbol{d}=begin{pmatrix}sindelta cosphi &0&-sindeltasinphiend{pmatrix}^{mathrm{T}}\[2ex] &boldsymbol{tilde{r}}:=boldsymbol{r}-boldsymbol{d}=begin{pmatrix}cos hcos A-sin delta cosphi&cos hsin A&-sin h+sindeltasinphi end{pmatrix}^{mathrm{T}}\[2ex] &boldsymbol{tilde{n}}:=boldsymbol{n}-boldsymbol{d}=begin{pmatrix}sinphicosdelta&0&cosphicosdeltaend{pmatrix}^{mathrm{T}}\[2ex] end{align}\

那麼可以計算：

begin{align} boldsymbol{tilde{r}timestilde{n}}&= begin{vmatrix} boldsymbol{hat{i}}&boldsymbol{hat{j}}&boldsymbol{hat{k}}\[2ex] cos hcos A-sin delta cosphi&cos hsin A&-sin h+sindeltasinphi\[2ex] sinphicosdelta&0&cosphicosdelta\[2ex] end{vmatrix}\[2ex]&= begin{pmatrix} cos hsin Acosphicosdelta\[2ex] -sin hsinphi -cos h cos Acosphicosdelta+sindeltacosdelta\[2ex] -cos h sin Asinphicosdelta\[2ex] end{pmatrix}=sin tcos^{2}deltabegin{pmatrix} cosphi\[2ex] 0\[2ex] -sinphi end{pmatrix}\[2ex] end{align}\ 所以我們得到瞭方程組： begin{cases} cos hsin A=cosdeltasin t\[2ex] sindelta=sin hsinphi +cos hcos phi cos A\[2ex] end{cases}\tag{11}\

消去方位角A 得： sin^2 h-2sin delta sinphisin h+cos^2deltacos^2phileft(tan^2deltatan^2phi-cos^2 tright)=0\

解二次方程得

sin h=sindeltasinphi pmcos tcos deltacosphi\

由於 t=0,t=pi 時分別代表太陽下中天和上中天，即 sin h 取的最小值和最大值，則 pm 取負。即： color{red}{sin h=sin deltasinphi -cos tcosdelta cosphitag{12}}\ 式(12)即可使我們計算太陽高度，而對於太陽方位角，隻需對方程組(11)中的第一個方程變一下形即可：

color{blue}{sin A=frac{cos delta sin t}{cos h}tag{13}}\

6.測量學校對面山的高差

先介紹一下我學校周圍的環境吧。

圖7 學校附近區域鳥瞰圖(黃線一端是待測山脈，另一端為學校教學樓所在地)，此圖向上大致朝西北，圖片源自Google Earth(後同)

重慶第二外國語學校^{[12]坐落於重慶市南岸區的南山山頂的一片背斜谷地[13]中央。地勢中間低，東西兩側高，即學校東西兩側是兩端近乎平行的山脈，在俯瞰圖中對應 "wedge"型的山嶺。這次需要測量的主角就是"wedge"型的山嶺更靠近兩江的一側山嶺。如下圖：}

3743d3d9deb00286385c55083549b314學校對面小山嶺的正視圖

實際景象呢可以參考一下圖片：

圖片來源於我的地理老師(下同) 由於拍攝時期和我觀測的時期不同，所以Google Earth上標記的日落點略有不同

觀測時間：2021年3月30日和5月28日傍晚；

• 觀測地點：先後有兩次觀測，一次在學校籃球場附近；一次在教室窗邊(5樓)；
• 觀測所需數據來源：

1. 觀測地大致經緯度：通過手表的GPS，獲取學校所在地大致經緯度；
2. 落山時刻太陽的方位與高度：先後有兩次在不同地點觀測太陽日落，並記錄下手表上顯示的時間。然後通過前文所撰寫的內容來運算出兩次觀測太陽的方位與高度；
3. 教室觀測點到教學樓末端的長度：使用學生直尺測量教學樓一塊地磚(正方形)的邊長，並計數從教室門口到走廊盡頭的地磚個數，從而得到觀測點教學樓教室觀測點到教學樓末端的長度；
4. 教室觀測點到籃球場觀測點的高差：由於教學樓下是平地，即籃球場與教學樓一樓等高，通過樓層數加上人的正常身高可計算出高差；
5. 教學樓朝向的方位角：通過"立筆見影"粗略測得(詳情見後)；

• 數據處理方法：運用前文所撰寫的知識以及一些簡單的解三角形知識。

來看看兩個觀測地內部環境，如圖：

圖8 教學樓(O2)與籃球場(O1)周邊鳥瞰圖 點G為教學樓走廊盡頭

為瞭簡便運算，得假設幾個條件：

1. O_1Gbot O_{2}G
2. 被測山脈兩次落山點的海拔相等

通過這兩個假設條件，就可以把這個測量問題完全轉化為數學計算題目瞭。作圖如下：

圖9 兩觀測地簡圖

根據學校經緯度 (29°32'38'' mathrm{N},106°36'23''mathrm{E}) ；3月30日日落時間 18^{mathrm{h}}40^{mathrm{m}} ,5月28日落時間 19^{mathrm{h}}29^{mathrm{m}} ；按照前文內容計算出一下數據： boxed{ begin{array}{c|c|c|c|c|c|c|c|c|c} text{Date} & text{Time } &tauleft(text {days}right) & E(text{rad}) & theta(text {rad}) &lambda(°) &delta (°) &t(°) &h(°) &A(°)\ hline 3.30 &18^{text{h}}40^{text{m}} & 10 &1.49433 &1.51129 &10.76744&4.26078&266.60639&5.04803&272.04397\ hline {5.28} &18^{text{h}}32^{text{m}} & 69 &2.50315 & 2.51321&68.17323&21.66461& 264.60639&14.95350&286.72895\ &19^{text{h}}40^{text{m}} & &&&&&278.85639&3.299570&293.10239\ end{array}}\

現在來計算圖9中的角 beta:=angle GO_{2}W 。如圖10：

圖10 5月28日18時32分太陽光打在窗臺上，通過"立筆見影"測量教學樓暢向方向角 俯視圖

用手扶著一支筆使其立在窗外的窗臺上，控制好筆的位置使其影子尖端恰好位於窗戶邊緣 C ,記錄下此時筆的位置 A 。使用學生直尺測出點 A 到窗戶的距離 AB=18text{cm} ，以及筆的長度 l=14 text{ cm} 。

由表可知，此時太陽高度 h=14.95350° ，那麼可知 AC=ltan h ，進而有 angle ACB=sin^{-1}left(frac{AB}{l tan h}right)approx 20.08359\

由表可知，此時太陽方位角 A=286.72895° ,即可計算得 color{blue}{beta=3.118411°} ;

由表可知，當日日落方位角 A=293.10239° ,結合 beta 即可計算得圖9中 color{red}{alpha:=angle S_{2}O_{2}G=26.220801°} ;

如圖9，過點 S_{2} 作 S_{1}O_{1} 平行線交 angle GO_{1} 延長線為點 {O_{1}}' 。由表可知3月30日日落方位角為 A=272.04397° ，做差得圖9中角 color{blue}{mu:=angle O_{2}S_{2}{O_{1}}'=21.05842°} ；

設圖9中 nu:=angle {O_{1}}'O_{2}G , L=O_{1}G 。過點 {O_{1}}' 作 O_{2}S_{2} 的垂線，垂足為點 D ，記 L'=O_{2}D 。

由於 O_{1}Gbot O_{2}G , {O_{1}}'Dbot O_{2}D ,故點 G,D 在以 O_{2}{O_{1}}' 為直徑的圓上，如圖11：

圖11

其中點 O 為圓心，設 varphi:=angle O_{1}{O_{1}}'G ，圓的半徑為 R 。由圓周角定理顯然得到 angle DOG=2alpha 。

根據幾何關系有： L=2Rsin frac{varphi}{2}~,~L'=2Rsin left(frac{varphi}{2}+alpharight)\

所以得到

L'=Lfrac{sinleft(cfrac{varphi}{2}+alpharight)}{sin cfrac{varphi}{2}}=Lfrac{cos nu}{cosleft(nu+alpharight)}tag{14}

如圖9，設 d_{1}={O_{1}}'S_{2},d_{2}=O_{2}S_{2} 。那麼對於 Delta S_{2}O_{2}{O_{1}}' ,有如下關系： d_{1}=L'frac{tannu}{sinmu}~,~d_{2}=L'left(1+frac{tannu}{tanmu}right)tag{15}\

設兩次落日時的太陽高度分別為 h_{1},h_{2} ，教學樓到籃球場的高差 h_{0} ，那麼在豎直方向上有：

d_{1}tan h_{1}-d_{2}tan h_{2}=h_{0}tag{16}\

聯立式(14)(15)(16)整理得關於 nu 的方程： frac{tannu}{sinmu}tan h_{1}-left(1+frac{tannu}{tanmu}right)tan h_{2}=frac{h_{0}}{L}frac{cosleft(v+alpharight)}{cosnu}tag{17}\

由於教室位於五樓，考慮人身高大致為 2 text {m} ，一層樓大約高 3 text {m} ，所以有 h_{0}=3times 4+2 text{m}=14 text{m} ;

用學生直尺測出教室一塊地磚(正方形)的邊長為 s=61text{cm} ， 計數 O_{2}G 線段方向上一共有 108 塊地磚，所以有 L=108times 0.61 text{m}=65.88text{ m} ；

根據表格，有 h_{1}=5.04803°~,~h_{2}=3.299570° 。將以上數據代入方程(17)解得 color{red}{nu=52.57647°} ;

再將 nu 代回式(14)(15)得 color{blue}{d_{2}=905.48437text{m}} ;故山頂到籃球場的高差為 color{red}{H=d_{2}tan h_{2}+h_{0}=68.47978text{m}}\

大功告成！！！

用谷歌地球試瞭一下，它給出的 d_{2} 大約為七百多米，比我算出來的數據小瞭不少，看來測量誤差還是蠻大的，想瞭想，大概誤差主要是 角beta ，高差 h_{0} ，走廊長 L 帶來的吧。怎麼說，至少知道那座山的高差在 45sim70 text{m} 之間[]~(￣▽￣)~*。

好瞭，這篇改版瞭不少次的文章，也是我來知乎寫的第一篇文章，就到此吧！感謝閱讀！

參考

1. ^註意是回歸年！不是恒星年！恒星年是指地球繞太陽公轉一周的周期；回歸年是指太陽直射點連續兩次由南半球經過赤道進入北半球的間隔時間。回歸年和恒星年是不相等的，因為地球自轉軸在由東向西進動，導致回歸年略小於恒星年
2. ^參考“天球“的百度百科 http://baike.baidu.com/item/%E5%A4%A9%E7%90%83/494605
3. ^其實這一步用幾何法更加簡單，直接設出黃道法向量以及太陽直射點徑矢做點乘即可；或者學過一點球面幾何的同學直接運用球面三角餘弦定理均可。但鑒於幾何法要重新作圖，作一張圖又比較麻煩，所以我轉而用旋轉矩陣來代替幾何法
4. ^有心力是指力方向的延長線始終經過一個固定點的力
5. ^參考 http://zhuanlan.icpchaxun.com/p/20684078
6. ^矢量空間的基底是正交歸一的，是指其中任二基矢u,v滿足g(u,v)=0且g(e_{mu},e_{nu})=pm 1，g是定義在該矢量空間上的一個度規
7. ^這裡的矢量是指滿足線性性以及萊佈尼茨律的一個從光滑函數到實數的映射。這裡取的特殊基底變換公式可以根據多元函數微積分快速推導
8. ^其實可以不用龍格-楞次矢量來推導，還有比耐方程法、復數法等等方法都可以推證
9. ^最早是由開普勒提出，並用幾何方法證明的，該文章的舊版有寫過幾何證明，這次改版就不寫瞭；在開普勒之後先後有拉格朗日、貝塞爾等著名數學傢參與過該方程的求解研究
10. ^這個解的化簡由知乎答主"夏草"在這篇回答中給出的 http://www.icpchaxun.com/question/65449609/answer/1703213535
11. ^我也沒看得懂/我是菜(。>︿<)_θ/
12. ^沒錯，就是不久前B站知名科普工作者畢導分享視頻動態"為君翻作「琵琶行」"中的學校 http://www.bilibili.com/video/BV1uT4y1R7ca?share_source=copy_web
13. ^背斜構造是指巖層發生折曲時，其形狀向上凸起的巖層。根據板塊構造學說，處在運動之中的板塊在發生相對運動時，板塊交界處會相互碰撞形成褶皺，就可能發育為背斜。在板塊內部，強大的擠壓作用亦會令到巖層發生彈性變形而形成背斜。背斜構造由於因巖石拉張易被侵蝕而形成谷地，即背斜谷