昨天看了一个考古节目，说在距今4700年的山西省襄汾县陶寺城遗址发现一个观象台。在当地发掘中发现一处夯土台上有若干柱子的基础，考古人员把相关数据交给中国天文学者进行计算，认为这是一个古观象台。从一个固定观测点看，在不同的节气日，太阳升起时出现在不同的柱子缝隙之间。《史记-五帝本纪》记载尧帝【其仁如天，其知如神...乃命羲、和，敬顺昊天，数法日月星辰，敬授民时。】因此，考古人员大胆推断这个遗址可能正是帝尧的都城。（下图来自网络）

节目中提到，为了计算太阳升起时的角度，很费了一番周折。开始竟然是等到冬至去实测，结果天天气不作美。看到这，我想起了这篇《日夜的方程》中计算日出、日落时间的方程。稍加改变，就可以计算日出时太阳的角度了 -- 也就是阳光与地面正南方向的角度。如果地球自转轴没有倾斜而是与公转轨道平面垂直，那么太阳总是在正东方向升起，与正南方向夹角为90度。现在我们需要计算的是实际情况。为完整起见，我下面重复《日夜的方程》的部分内容。

要做这个计算，需要先确定坐标。那么就让我们以冬至这点的状态设置坐标，然后开始计算。下图中，O是地球的球心，N是北极，S为南极，地球绕 NS 轴自转，P为地球表面一点; A 为太阳所在，但应该想象太阳在右边无限远处，阳光是平行的，方向是从A到地心的连线；地球绕太阳的轨道是在水平面上，DOE 线垂直于轨道平面；地球自转轴与 DOE 线的夹角为 $\beta$，就像个倾斜的陀螺。冬至日的时候，NOAD 四点在同一个面上，垂直于地球轨道平面。此时，地球自转轴与阳光的夹角 <NOA 就是地球自转倾角 $\beta$ 加上90度。计算 OP 与阳光 的角度 $\theta$ 是一个相当直接的几何问题。坐标设置如下图中的 XYZ 与 X'Y'Z'。

上图中，设P点的纬度为 L，P在XY平面上的投影为 P', OP' 与 -Y 轴角度为 d 。考虑地球在自转的同时还在公转。我们希望得到的是一个含时间的方程，同时考虑自转与公转。在XYZ坐标系中，OP的方向为 $(\cos L \sin d, -\cos L\cos d, \sin L)$。而在X'Y'Z'坐标系中，AO'方向为 $(\sin\phi, -\cos\phi, 0)$。要计算两者的角度需要找出两个坐标之间变换关系。这是一个角度为 $\beta$的转动。我们有

$\hat{x} = \hat{x^\prime}\\ \hat{y}=\cos\beta \ \hat{y^\prime} +\sin\beta \ \hat{z^\prime}\\ \hat{z} = -\sin\beta \ \ \hat{y^\prime} +  \cos\beta \ \hat{z^\prime}$

代入计算OP 与 O'A 之间的角度，也就是阳光与垂直于地面的标杆的角度，我们得到

$\cos\theta = [\cos d \ \cos L \cos\beta + \sin L \sin\beta] \cos\phi + \sin\phi \sin d \cos L$

现在我们计算阳光与地面正南方向的夹角$\alpha$。XYZ 坐标中P处正南方向为 $(\sin L \sin d, -\sin L \cos d, -\cos L)$。进行坐标变换之后计算与阳光方向的余弦 ，我们得到

$\cos \alpha = \left[\cos\beta \sin L \cos d- \sin\beta \cos L\right] \cos\phi + \cos L \sin d \sin\phi$

太阳升起、落下的时候，阳光与垂直标杆垂直，因此 $\cos\theta =0$。所以，日出、日落时

$[\cos d \ \cos L \cos\beta + \sin L \sin\beta] \cos\phi + \sin\phi \sin d \cos L=0$

设地球自转角速度为 $\omega$, 公转角速度为 $\Omega$, t 为时间，则 $d = \omega \ t $, $\phi=\Omega \ t$。

因此太阳与直立标尺的角度 $\theta$ 的方程为

$\cos\theta = \left[ \sin L \sin\beta + \cos L \cos\beta \ \cos (\omega \ t)  \right]\ \cos (\Omega\ t) + \cos L \sin (\omega\ t) \sin(\Omega\ t)   $

至此，我们得到计算日出、日落时间的方程是:

$\left[ \sin L \sin\beta + \cos L \cos\beta \ \cos (\omega \ t)  \right]\ \cos (\Omega\ t) + \cos L \sin (\omega\ t) \sin(\Omega\ t)    =0$

或者：

$ \frac{ \sin L \sin\beta + \cos L \cos\beta \ \cos (\omega \ t)}{\cos L \sin (\omega\ t)}=-\tan\Omega t $

从上面的方程计算出日出时间，并代入阳光与正南方向夹角的公式

$\cos \alpha = \left[\cos\beta \sin L \cos (\omega t)- \sin\beta \cos L\right] \cos (\Omega t) + \cos L \sin (\omega t) \sin (\Omega t) $

即可算出太阳升起时阳光的角度。代入实际数据，倾角 \beta = 23.44 度，计算出冬至日北京（纬度取北纬40度）太阳升起阳光与正南方向的角度为121.27 度。夏至日太阳升起时角度 为58.7度。陶寺纬度为35.88度，算出冬至太阳升起角度为119.4度。如果考虑折射，这个角度将略有不同。另外，地球自转倾角以两万六千年为周期变化，因此4700年前倾角与现在也有不可忽略的差异。

附：Python 代码

http://www.zhenzhubay.com/home.php?mod=space&uid=2&do=blog&id=36911