使用 Skyfield 计算朔、望、弦
子曰：「有库即用矣。」（误）
发现 Skyfield 是个不错的库。第一个程序就拿来计算朔、望、弦的时刻吧。
| 月相 | 日月黄经差（模 360 度） |
| 朔 | \( 0^{\circ} \) |
| 上弦 | \( 90^{\circ} \) |
| 望 | \( 180^{\circ} \) |
| 下弦 | \( 270^{\circ} \) |
解决这个问题的关键，就是建立一个日月黄经差关于时间的函数 \( \Delta L(t) \)，只要这个函数有了，用牛顿法对 \( \Delta L(t) = 0^{\circ} \), \( \Delta L(t) = 90^{\circ} \) 等四个方程进行求根即可算出对应四个特殊月相的时刻。所以关键是搞出这个函数来。而 Skyfield 直接提供了计算给定时刻太阳和月球黄经的实现，所以简直得来全不费功夫。
首先，加载数据：
from skyfield.api import load, utc from datetime import datetime time_scale = load.timescale() planets = load('de421.bsp')
（PS: de421.bsp 这一个数据文件 (17M) 就包含了全部八大行星，拿来算月相有大炮打蚊子的嫌疑）
然后定义日月黄经差（角度，模 360）和时间（秒，UNIX timestamp）的函数关系：
def delta_l(t): T = time_scale.utc(datetime.utcfromtimestamp(t).replace(tzinfo=utc)) earth = planets['earth'] sun = planets['sun'] moon = planets['moon'] _, l1, __ = earth.at(T).observe(sun).ecliptic_latlon() _, l2, __ = earth.at(T).observe(moon).ecliptic_latlon() return (l2.degrees - l1.degrees) % 360
难以置信这就已经算完了（有库就是省事）。接下来的求根，就不是本文的重点了。通过改变迭代的初值，就可以求出各月的朔、望、上弦、下弦的时刻。至于精度，毕竟天文上变量太多情况太复杂，也不敢说这样算有多精确，但起码在分钟的数量级上是没有多少出入的。
顺带一提，有了这个库，用天文算法计算农历就不能成为问题。但是 JPL 的数据都太大，讲真还不如打表（事先算出节气和朔日的日期存下来）。
下面是使用上文代码计算出的 2017 年各月的月相数据（完整程序见 Gist）：
上弦 2017-01-06 03:47:39 望 2017-01-12 19:35:11 下弦 2017-01-20 06:14:12 朔 2017-01-28 08:07:01 上弦 2017-02-04 12:19:31 望 2017-02-11 08:34:11 下弦 2017-02-19 03:33:53 朔 2017-02-26 22:58:22 上弦 2017-03-05 19:33:00 望 2017-03-12 22:55:10 下弦 2017-03-20 23:58:56 朔 2017-03-28 10:57:13 上弦 2017-04-04 02:40:03 望 2017-04-11 14:09:33 下弦 2017-04-19 17:57:27 朔 2017-04-26 20:16:09 上弦 2017-05-03 10:47:30 望 2017-05-11 05:43:58 下弦 2017-05-19 08:33:29 朔 2017-05-26 03:44:27 上弦 2017-06-01 20:42:47 望 2017-06-09 21:11:04 下弦 2017-06-17 19:33:22 朔 2017-06-24 10:30:42 上弦 2017-07-01 08:51:48 望 2017-07-09 12:08:01 下弦 2017-07-17 03:26:17 朔 2017-07-23 17:45:34 上弦 2017-07-30 23:23:49 望 2017-08-08 02:12:01 下弦 2017-08-15 09:15:41 朔 2017-08-22 02:30:10 上弦 2017-08-29 16:13:43 望 2017-09-06 15:04:08 下弦 2017-09-13 14:25:38 朔 2017-09-20 13:29:52 上弦 2017-09-28 10:54:15 望 2017-10-06 02:41:23 下弦 2017-10-12 20:26:03 朔 2017-10-20 03:12:04 上弦 2017-10-28 06:22:48 望 2017-11-04 13:24:07 下弦 2017-11-11 04:37:03 朔 2017-11-18 19:42:08 上弦 2017-11-27 01:03:39 望 2017-12-03 23:48:09 下弦 2017-12-10 15:52:02 朔 2017-12-18 14:30:25 上弦 2017-12-26 17:20:49
