日行迹的计算

日行迹指的是在同一地点每天的同一时刻太阳在当地地平坐标系上位置变化的连线，通常可以使用固定位置的广角或者鱼眼相机进行记录。在地球上观察，太阳的日行迹呈“8”字型。由于对时间和位置精度的要求极高，据说目前世上只有几十张成功的日行迹照片。该现象的解释其实早有科普文章：

baijiahao.baidu.com/s?id=1643540164292027533 （此作者也有b站账号，但不知道为何没有将此文搬运至专栏）

但由于其文中只进行定性讨论，未有定量分析（一般科普文章也确实不这么干，大家不爱看），所以我特意进行了从零开始的计算模拟，以复现“8”字形的日行迹。

本文使用的基础数据

近日点：1.471e11 米

远日点：1.520e11 米

太阳质量：1.9891e30 千克

恒星年：31558150 秒

太阳日：86400 秒

恒星日：86164.1 秒

黄赤交角：23°26′11″=0.409022（弧度）

椭圆运动计算

以椭圆中心为坐标轴原点O，近日点方向为x轴正方向，近日点处地球运动方向为y轴正方向，x×y方向为z方向，建立三维直角坐标系。

a=1.49615e11（米）, b=1.4959494e11（米）, c=2.45e9（米）

离心率：0.0163754

太阳坐标(2.45e9,0,0)

地球所处位置与Ox方向的夹角为θ

参考椭圆运动计算方法得出时间t与夹角θ的函数关系：

(GM/a)^0.5 t = aθ - csinθ

由于此为超越方程，无解析解，根据实际的G、M、a、c值对其进行傅里叶展开后求得近似反函数，并针对一年周期进行优化，同时为方便后续计算，将t单位改为日，得到如下方程：

θ=t+0.000099sin0.00860106201t+0.951881sin0.01720212402t+0.007791sin0.03440424804t+0.000095sin0.05160637206t-0.00000536t

求出θ后，易得地球在t时刻的坐标

x=acosθ, y=bsinθ

在一年周期内该反函数与原函数的θ偏差小于0.000034°，时间t偏差不超过3秒。

时间选取

取近地点1月3日的正午12:00作为t=0时刻，此时θ=0，地球坐标为(1.49615e11,0,0)，地日向量为(-1.471e11,0,0)，地日方向向量(-1,0,0)。

向量P沿向量A旋转φ的求法

转轴向量A(xA, yA, zA)

待转动向量P(xP, yP, zP)

转动后向量P'(xP', yP', zP') = Pcosφ + A×Psinφ + [ A·P·(1-cosφ) ] A

xP' = xP cosφ + ( zP·yA - yP·zA ) sinφ + xA (xA·xP + yA·yP + zA·zP) (1-cosφ)

yP' = yP cosφ + ( xP·zA - zP·xA ) sinφ + yA (xA·xP + yA·yP + zA·zP) (1-cosφ)

zP' = zP cosφ + ( xP·yA - yP·xA ) sinφ + zA (xA·xP + yA·yP + zA·zP) (1-cosφ)

0时刻纬度为L的点的当地法向量L0

将时间拨至12月22日冬至，此时地日方向向量为s(-0.9763,0.2162,0)，由于此时太阳直射南半球，因此且地轴北方向在xOy平面上的投影与地日向量方向相反。再加上黄赤交角求出地轴的方向向量a(0.38907, -0.08254, 0.91750)

再将时间调回至0时刻，由于地球的公转，地轴在xOy平面的投影已经偏离了地日向量所在直线。将地轴向量×地日向量，求出太阳-地心-地轴所在平面F1的法向量a×s=F1，0时刻当地法向量L0可以通过将地轴方向向量a沿着F1旋转(90-纬度)°得到（取北纬为正，南纬为负）。

t时刻观测地法向量Lt

依据地球的恒星日86164.1、t时刻对应的日期和时间，求出t时刻地球沿着地轴转过的角度φ=86400/86164.1t。

将L0沿着地轴a转过t时刻的φ角得到Lt。

该法向量Lt同时也是观测地天顶（Zenith）所指向的方向。

t时刻太阳高度角E

通过当地法向量Lt和地日方向向量s的数量积可以求出t时刻太阳与天顶的夹角：

cos(90°-E)=Lt·s / |Lt|·|s|

继而可以求出高度角E

t时刻太阳方位角B

首先求出观测点-地心-地轴所在平面的法向量N，作为正北方向的参考

N=Lt×a

再求出观测点-地心-太阳所在平面的法向量B，作为太阳方位角的参考

B=Lt×s

依据向量N和向量B的夹角求出太阳的方位角B。

日行迹绘制

此时模型已经建立完成，输入纬度、日期和时间（以中午12:00为0:00，下午5点为5:00）可以求出当时的太阳方位角、高度角。

这里我选取北纬30°，下午1点至傍晚7点，每个时科各从1月3日起至次年1月3日取366个数据点，计算太阳的方位角和高度角。

以方位角和高度角作为横纵坐标绘制日行迹。

并挑选出特定节气（大寒、雨水、春分、谷雨、小满、夏至、大暑、处暑、秋分、霜降、小雪、冬至）单独做一张图，标记为黑点，叠加至日行迹上。

选取同一时刻（下午1时），绘制南纬80度到北纬80度的日行迹变化做成动图，可以发现整个日行迹始终保持小环（夏环）朝北、大环（冬环）朝南的方向。身处北半球时，日行迹在观察者的南方，看到的天球是北上南下，左东右西，因此小环朝上大环朝下。到了赤道地区，日行迹在西面，看到的天球是北右南左，下西上东，因此小环朝右大环朝左。到了南半球，日行迹在北面，看到的天球是北下南上，左西右冬，因此小环朝下大环朝上。大小环上下关系的颠倒来源于图像相对地面所见方向的颠倒。

选取同一地点不同时刻绘制日行迹变化动图，可以发现日落最晚的一天比夏至日略晚，日落最早的一天比冬至日略早。

重新选取北纬30°正午12点的日行迹，以仿鱼眼镜头的方式绘制日行迹，此时畸变最小，与人眼观察到的日行迹形状最接近。按图1类似的方法标记节气在日行迹的点上，我们可以发现交叉点在春分之后、谷雨之前，处暑之后、秋分之前。挨个输入日期尝试后发现，交叉点在4月12日和8月29日。冬至至上行交叉点为115日，上行交叉点到夏至为70日，夏至至下行交叉点为69日，下行交叉点到下一个冬至为111日，呈对称分布。

绘制地球绕日运行的椭圆轨迹，标记上各个时间点后如下图。日行迹的四个象限沿着冬至-夏至点连线呈对称状。在冬至、夏至日，自转轴和地日连线在同一平面上，自转与公转对太阳方位的影响相互抵消。受椭圆轨道不同区间公转速度不同，和地球自转的恒星年周期共同影响，在第一个四分之一环，地球自转稍微落后于公转，使得日行迹向东偏移；第二个四分之一环稍微超越，使得日行迹向西偏移；另两个四分之一环以此类推。

最后附上其中一张日行迹的计算过程