日行跡指的是在同一地點(diǎn)每天的同一時(shí)刻太陽(yáng)在當(dāng)?shù)氐仄阶鴺?biāo)系上位置變化的連線，通?？梢允褂霉潭ㄎ恢玫膹V角或者魚(yú)眼相機(jī)進(jìn)行記錄。在地球上觀察，太陽(yáng)的日行跡呈“8”字型。由于對(duì)時(shí)間和位置精度的要求極高，據(jù)說(shuō)目前世上只有幾十張成功的日行跡照片。該現(xiàn)象的解釋其實(shí)早有科普文章：

baijiahao.baidu.com/s?id=1643540164292027533 （此作者也有b站賬號(hào)，但不知道為何沒(méi)有將此文搬運(yùn)至專欄）

但由于其文中只進(jìn)行定性討論，未有定量分析（一般科普文章也確實(shí)不這么干，大家不愛(ài)看），所以我特意進(jìn)行了從零開(kāi)始的計(jì)算模擬，以復(fù)現(xiàn)“8”字形的日行跡。

本文使用的基礎(chǔ)數(shù)據(jù)

近日點(diǎn)：1.471e11 米

遠(yuǎn)日點(diǎn)：1.520e11 米

太陽(yáng)質(zhì)量：1.9891e30 千克

恒星年：31558150 秒

太陽(yáng)日：86400 秒

恒星日：86164.1 秒

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

橢圓運(yùn)動(dòng)計(jì)算

以橢圓中心為坐標(biāo)軸原點(diǎn)O，近日點(diǎn)方向?yàn)閤軸正方向，近日點(diǎn)處地球運(yùn)動(dòng)方向?yàn)閥軸正方向，x×y方向?yàn)閦方向，建立三維直角坐標(biāo)系。

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

離心率：0.0163754

太陽(yáng)坐標(biāo)(2.45e9,0,0)

地球所處位置與Ox方向的夾角為θ

參考橢圓運(yùn)動(dòng)計(jì)算方法得出時(shí)間t與夾角θ的函數(shù)關(guān)系：

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

由于此為超越方程，無(wú)解析解，根據(jù)實(shí)際的G、M、a、c值對(duì)其進(jìn)行傅里葉展開(kāi)后求得近似反函數(shù)，并針對(duì)一年周期進(jìn)行優(yōu)化，同時(shí)為方便后續(xù)計(jì)算，將t單位改為日，得到如下方程：

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

求出θ后，易得地球在t時(shí)刻的坐標(biāo)

x=acosθ, y=bsinθ

在一年周期內(nèi)該反函數(shù)與原函數(shù)的θ偏差小于0.000034°，時(shí)間t偏差不超過(guò)3秒。

時(shí)間選取

取近地點(diǎn)1月3日的正午12:00作為t=0時(shí)刻，此時(shí)θ=0，地球坐標(biāo)為(1.49615e11,0,0)，地日向量為(-1.471e11,0,0)，地日方向向量(-1,0,0)。

向量P沿向量A旋轉(zhuǎn)φ的求法

轉(zhuǎn)軸向量A(xA, yA, zA)

待轉(zhuǎn)動(dòng)向量P(xP, yP, zP)

轉(zhuǎn)動(dòng)后向量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時(shí)刻緯度為L(zhǎng)的點(diǎn)的當(dāng)?shù)胤ㄏ蛄縇0

將時(shí)間撥至12月22日冬至，此時(shí)地日方向向量為s(-0.9763,0.2162,0)，由于此時(shí)太陽(yáng)直射南半球，因此且地軸北方向在xOy平面上的投影與地日向量方向相反。再加上黃赤交角求出地軸的方向向量a(0.38907, -0.08254, 0.91750)

再將時(shí)間調(diào)回至0時(shí)刻，由于地球的公轉(zhuǎn)，地軸在xOy平面的投影已經(jīng)偏離了地日向量所在直線。將地軸向量×地日向量，求出太陽(yáng)-地心-地軸所在平面F1的法向量a×s=F1，0時(shí)刻當(dāng)?shù)胤ㄏ蛄?strong>L0可以通過(guò)將地軸方向向量a沿著F1旋轉(zhuǎn)(90-緯度)°得到（取北緯為正，南緯為負(fù)）。

t時(shí)刻觀測(cè)地法向量Lt

依據(jù)地球的恒星日86164.1、t時(shí)刻對(duì)應(yīng)的日期和時(shí)間，求出t時(shí)刻地球沿著地軸轉(zhuǎn)過(guò)的角度φ=86400/86164.1t。

將L0沿著地軸a轉(zhuǎn)過(guò)t時(shí)刻的φ角得到Lt。

該法向量Lt同時(shí)也是觀測(cè)地天頂（Zenith）所指向的方向。

t時(shí)刻太陽(yáng)高度角E

通過(guò)當(dāng)?shù)胤ㄏ蛄?strong>Lt和地日方向向量s的數(shù)量積可以求出t時(shí)刻太陽(yáng)與天頂?shù)膴A角：

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

繼而可以求出高度角E

t時(shí)刻太陽(yáng)方位角B

首先求出觀測(cè)點(diǎn)-地心-地軸所在平面的法向量N，作為正北方向的參考

N=Lt×a

再求出觀測(cè)點(diǎn)-地心-太陽(yáng)所在平面的法向量B，作為太陽(yáng)方位角的參考

B=Lt×s

依據(jù)向量N和向量B的夾角求出太陽(yáng)的方位角B。

日行跡繪制

此時(shí)模型已經(jīng)建立完成，輸入緯度、日期和時(shí)間（以中午12:00為0:00，下午5點(diǎn)為5:00）可以求出當(dāng)時(shí)的太陽(yáng)方位角、高度角。

這里我選取北緯30°，下午1點(diǎn)至傍晚7點(diǎn)，每個(gè)時(shí)科各從1月3日起至次年1月3日取366個(gè)數(shù)據(jù)點(diǎn)，計(jì)算太陽(yáng)的方位角和高度角。

以方位角和高度角作為橫縱坐標(biāo)繪制日行跡。

并挑選出特定節(jié)氣（大寒、雨水、春分、谷雨、小滿、夏至、大暑、處暑、秋分、霜降、小雪、冬至）單獨(dú)做一張圖，標(biāo)記為黑點(diǎn)，疊加至日行跡上。

選取同一時(shí)刻（下午1時(shí)），繪制南緯80度到北緯80度的日行跡變化做成動(dòng)圖，可以發(fā)現(xiàn)整個(gè)日行跡始終保持小環(huán)（夏環(huán)）朝北、大環(huán)（冬環(huán)）朝南的方向。身處北半球時(shí)，日行跡在觀察者的南方，看到的天球是北上南下，左東右西，因此小環(huán)朝上大環(huán)朝下。到了赤道地區(qū)，日行跡在西面，看到的天球是北右南左，下西上東，因此小環(huán)朝右大環(huán)朝左。到了南半球，日行跡在北面，看到的天球是北下南上，左西右冬，因此小環(huán)朝下大環(huán)朝上。大小環(huán)上下關(guān)系的顛倒來(lái)源于圖像相對(duì)地面所見(jiàn)方向的顛倒。

選取同一地點(diǎn)不同時(shí)刻繪制日行跡變化動(dòng)圖，可以發(fā)現(xiàn)日落最晚的一天比夏至日略晚，日落最早的一天比冬至日略早。

重新選取北緯30°正午12點(diǎn)的日行跡，以仿魚(yú)眼鏡頭的方式繪制日行跡，此時(shí)畸變最小，與人眼觀察到的日行跡形狀最接近。按圖1類似的方法標(biāo)記節(jié)氣在日行跡的點(diǎn)上，我們可以發(fā)現(xiàn)交叉點(diǎn)在春分之后、谷雨之前，處暑之后、秋分之前。挨個(gè)輸入日期嘗試后發(fā)現(xiàn)，交叉點(diǎn)在4月12日和8月29日。冬至至上行交叉點(diǎn)為115日，上行交叉點(diǎn)到夏至為70日，夏至至下行交叉點(diǎn)為69日，下行交叉點(diǎn)到下一個(gè)冬至為111日，呈對(duì)稱分布。

繪制地球繞日運(yùn)行的橢圓軌跡，標(biāo)記上各個(gè)時(shí)間點(diǎn)后如下圖。日行跡的四個(gè)象限沿著冬至-夏至點(diǎn)連線呈對(duì)稱狀。在冬至、夏至日，自轉(zhuǎn)軸和地日連線在同一平面上，自轉(zhuǎn)與公轉(zhuǎn)對(duì)太陽(yáng)方位的影響相互抵消。受橢圓軌道不同區(qū)間公轉(zhuǎn)速度不同，和地球自轉(zhuǎn)的恒星年周期共同影響，在第一個(gè)四分之一環(huán)，地球自轉(zhuǎn)稍微落后于公轉(zhuǎn)，使得日行跡向東偏移；第二個(gè)四分之一環(huán)稍微超越，使得日行跡向西偏移；另兩個(gè)四分之一環(huán)以此類推。

最后附上其中一張日行跡的計(jì)算過(guò)程