夏日炎炎。在网上偶尔看到一篇英文提到中国古代天文学，里面有个词“Winter Extreme"。顾名思义这是中文的”冬至”，英文叫 Winter Solstice。我之前一直想当然地认为“冬至”是冬天来到的意思。难道这个“至”确实是至极之意？查阅之后发现确实如此。《后汉书-律历志》注解中写道：“冬至之为极有三意焉：昼漏极短，去极极远，晷景极长。极者，至而还之辞也。” 冬至是白天时间最短的一天。还对极点给出了数学直觉极好的定义：到了然后返回之点 。Solstice 一词来自拉丁文，意思是太阳在天球上的纬度看起来不变之意。用数学术语来说，相当于一个曲线顶部的平坡处（导数为零）。殊途同归。冬至一词这个至字的含义我这才弄清楚，还是看了英文翻译才幡然醒悟。感觉需要继续思考一下相关问题了。

且让我们根据基本天文知识--地球自转的同时绕太阳公转---来写个数学公式，计算每天太阳升起与落下的时间，阳光投射角度随时间的变化等等。要做这个计算，需要先确定坐标。那么就让我们以冬至这点的状态设置坐标，然后开始计算。下图中，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 。阳光方向 为AO 。在XYZ坐标系里，OP的方向为 $(\cos L \sin d, -\cos L\cos d, \sin L)$，而OA方向为$(0,-\cos\beta, \sin\beta)$。

我们有

$\cos\theta =  \cos\beta \ \cos L \ \cos d + \sin \beta \ sin L $

太阳升起与落下时，OP 与阳光角度垂直，$cos\theta=0$. 所以, 日出、日落的方程是：

$ \cos\beta \ \cos L \ \cos d + \sin \beta \ sin L=0$

也就是说

$\cos d = -\tan \beta \tan L$

从上面方程求解出 d 也就知道了这天太阳升起的时间。以北京为例， 纬度L约为 40度，代入地球轨道倾角 23.43688 度，解出 d = 111.33 度。地球每小时转15度，因此这天的太阳升起的时间是 111.33/15 = 7.42 （小时），也就是7点25 分。日落时 角度为 270-21.33，日落时间因此是 (270-21.33)/15 =16.6，也就是下午4点多。

另外，上面的方程并不见得总是有解。对于高纬度地区，右边的绝对值完全可以大于1。这对应于两极地区极昼或者极夜的情况。

冬至点之后，地球继续绕太阳转动，到了O' 的位置。跳出太阳系看世界，地球自转轴NS的方向不变，但是阳光方向变成了 AO' 方向，NO'AD 四点不再在一个平面上，地球自转轴与阳光的夹角 <NO'A 也变了。怎么办？

如果想继续使用上面的公式，我们需要首先计算出夹角 <NO'A，然后用这个角度减去90度, 得出的 $\beta^\prime$，然后代入上面的$\sin d$ 公式如法炮制。这相当于把地球翻转成上图冬至时的情形（但 DOE不再垂直于轨道平面）。从上图中可以看出

$\sin\beta^\prime = \cos\phi\ \sin\beta$

从这个公式可见，从冬至点起，地球在绕太阳轨道运行90度之后，地轴与阳光方向垂直 (公式右边因为 cos 项为零，因此 $\beta^\prime =0$) 。日出方程成为 $\sin d=0$。此日昼夜长度相等，古代中国称为春分。西方称为 Equinox （意思是昼夜相等）。

综合起来，我们的日出方程为

$\cos d = -\tan[\sin^{-1}(\cos\phi\sin\beta)]\tan L = -\frac{\sin\beta\cos\phi\tan L}{\sqrt{1-\sin^2\beta\cos^2\phi}}$

上面的结果看起来不错，但有一个问题。计算是把地球移到轨道的某一处，然后假设地球在原处转动，计算从夜中开始到凌晨的转角。但是正确的计算应该考虑地球在自转的同时还在公转。我们希望得到的是一个含时间的方程，同时考虑自转与公转。如上所述，在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$

设地球自转角速度为 $\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 $

其中 L 为 纬度，$\beta$ 为行星自转倾角（地球约23.5 度），$\omega$ 为行星自转角速度，$\Omega$ 为行星公转角速度，时间是从冬至算起。当然上述计算中假定轨道是圆形，如果是实际椭圆轨道，那方程将复杂些。另外，空气折射也没有考虑。

要解这个方程一时也看不出有什么可以简化，似乎只能用蛮力。这正是计算机的长处。对比两种方法，以北纬40度为例计算从冬至半夜 (d=0) 到天明的时间，后一种晚约一分19秒，这是一个不可忽略的差别。这也说明第一种计算存在相当的缺陷。

文末PYTHON代码计算了北纬40度365天的日出、日落时间。

代码链接：

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