關於作者提供的農曆編算方法檔案中的一些疑問
Closed this issue · 29 comments
將 t0 限制在5日和20日當然可以，您的計算不就已證實了嗎? 其實這方法接近用古代的平氣法算出的平節氣作為 t0。
至於收斂速度變化不大，就是牛頓-拉弗森法高明之處，此法收斂極速，是二階收斂算法。如果您詳細參閱此法的立法原則，就會知道如果欲求之方程是線性方程，用牛頓法求出的t1就已是正確答案。太陽黃經之變化可表示成線性函數加上修正項，這修正項古代稱為"日行盈縮"。日行盈縮是微小的修正項，所以求節氣的方程近似線性方程，即使 t0 並不很接近節氣時刻，用牛頓法算出的 t1 已經很準，算到 t2 基本上就已經收斂了，因此除非你的 t0 真的非常準確，能保證算出的 t1 就已收斂，否則如果仍是要算到 t2 才收斂，那就看不出收斂速度有大變化了。
我的其中一個算法是先算出某年1月1日之前的冬至時刻，然後以此作為下一節氣(小寒)的 t0，再以此節氣作為下一節氣(大寒)的 t0，按此方法便很容易編寫程式算出全年的節氣，用類此方法也輕易算出全年的月相，這就是我的 TDBtimes.txt 時刻的計算方法。此算法當然不是最有效快捷之法，但是容易編寫程式，亦不易出錯，而且需十幾秒就製成了1600年至3500年的月相節氣時刻表，如果簡化章動計算，或許可以快十倍。依此表便可編算現代農曆。
既然這些計算只須算一次，我認為不必太計較算法效率，更重要的是保證算出的數據正確。您算月相節氣也是為了編算農曆嗎? 還是另有用途? 您也是用 JPL 的 DE 曆表和 IAU 推薦的歲差和章動模型計算月相節氣嗎? 算出的結果是否與我的 TDBtimes.txt 數據相符?
閣下自稱晚輩，但不知是何方高士? 若不願在此公開，可否電郵告之?
拜讀過您給出的農曆編算方法檔案後晚輩先是實現了檔案中第六節地心視黃經的算法，我想到先通過紫金山天文臺給出的節氣數據進行反向驗證該步驟是否會有差錯，於是先將紫金山天文臺給出的2020年節氣發生的準確時間轉換為儒略日數然後通過實現的地心視黃經算法計出黃經角度，但是我得到的結果總是會與標準角度（例：立春的315度和立冬的225度等）相差0.1"，這與您檔案中給出的需要將角度的差異限定在0.04"內有很大的差異。當時我想到造成誤差的原因可能有二：其一是我實現的算法有誤，其二是我找尋得的驗證數據有誤。後者我無法驗證，但若是前者造成的差異可以縮減。
因為我讀取JPL插值數據的程式是早年根據Project Pluto C程式編寫的Java和Python程式，主要用於讀取405歷表，為了避免讀取程式造成的差異，我改用直接使用Project Pluto提供的C程式通過431歷表計算天體位置並重新編寫C程式計算。不過這次並沒有使用紫金山天文台的數據作為校驗，而是直接將您的TDBtimes檔案作為Test File。大致思路也比較簡單，首先讀取您給出的jd0然後使用我的程式向前計出冬至發生時刻根據該冬至時刻計出其後的全部24節氣，因為您檔案中的也是TDB時刻，所以我可以直接將計算結果進行比對而無需擔心時刻轉換造成的誤差。
但是我計出的結果發現在2000年前後500年內節氣時刻與您的結果誤差不超過1s，但是超過500年後誤差開始逐漸增大，在1600-3500年範圍內最大誤差超過8s。因為我使用的章動模型是IAU推薦的，想到您使用的是Vondr ́ak等人給出的章動算法，於是我又根據Vondr ́ak等人的文章修改了章動的算法，得到的結果依舊是在1600-3500年範圍內最大誤差多達幾秒。我想應該是我哪個計算步驟不對，但無法驗證。
想到月相的計算可以忽略掉∆ψ造成的誤差，於是調整程式開始計算月相，得到的結果在1600-3500年範圍內半數的時刻與您給出的結果在32個浮點位內完全一致，最大的時刻誤差不超過0.0002s。因此我確信應該是章動或歲差計算有差異。於是我重新審視了幾遍程式，但未發現有明顯的錯誤。以您給出的1600年的jd0(2305446.166666667)向前計算首個冬至發生時刻為例，根據該jd0得出的歲差矩陣為：
章動矩陣乘積為：
兩者與參考架偏差矩陣的乘積為：
因為缺乏驗證數據因此沒法確認歲差與章動的計算結果是否正確，但是我根據您的文檔校對了幾次歲差與章動的計算過程，沒有發現有差異，如若您有時間可否幫驗證下上述三個矩陣的結果是否與您的計算一致？
晚輩只是一個普通天文愛好者，非專業人士，以前使用VSOP87/2013、Swiss Ephemeris、ELP/MPP、Pluto99等歷表和算法計算太陽系內的天體運行軌跡，最開始只爲觀星用，使用的數據精度有過刪減，各天體的時刻精度在十到三十分鐘不等。後來也是慢慢發現可以用其計算歷表，因為十到三十分鐘的精度不太適合用於計算歷表，後來我改進了算法使用了全部的精確數據，雖然計算相對繁瑣但可以將時刻精度保持在五分鐘以內。再後來接觸到JPL歷表，最開始使用是的405版本的歷表，改進後的算法可將時刻精度在相當長的時間範圍內收斂在三分鐘以內；因為405的有效時間範圍不長，後準備使用431的歷表替換之，遂搜尋得您的這個Project。
除了月相和節氣的計算外，天體的其它天象數據我也有計算；精確計算這些數據的目的更多的需求是為了能得到一個精確歷表算法解算出的數據與我自行精簡得出的擬合數據進行一個比對，因為我所需的場景中無法承載像JPL歷表這樣如此龐大的檔案。如您所屬，日曆的編算僅需預先編排得到月相和節氣的 TDB時刻即可，但是其它的天象，如行星與天球上**（出入地平線、上下中天等）就需要依據測站的經緯度和海拔等等外部參數的變化不斷改變計算數據而無法通過預算的方式略去這一過程。
歪个楼请教一下，你这里计算过木星的数据吗？看古书上说太岁是木星。想把自己写的农历加入太岁。谢谢
找到問題所在了，黃經章動的第二項缺少了儒略世紀數的乘數，修正後1600年的jd0(2305446.166666667)向前計算首個冬至發生時刻的章動矩陣積為：
計得與您的TDBtimes中的節氣數據時刻1600-3500年內相差最大不超過0.2s，不知道這是否正常的誤差範圍？
這次檢視我又在計算月相時發現一個很有趣的現象，我所計得的月相時刻與您TDBtimes的時刻除了完全一致的數據外，其它數據都有一點點的誤差，不過最大誤差如上所述僅不超過0.0002s。不過有趣的是，該次計算我發現除了那些完全一致的時刻數據外，其餘有差異的所有數據的差值可分為六組一模一樣的數據，即僅有六個不同的差值數，具體為：
- 差值為 0.0000000000000000 的數據占了447項
- 差值為-0.0000402331352234 的數據占了648項
- 差值為-0.0000804662704468 的數據占了451項
- 差值為 0.0000402331352234 的數據占了164項
- 差值為-0.0001206994056702 的數據占了144項
- 差值為 0.0000804662704468 的數據占了 21項
- 差值為-0.0001609325408936 的數據占了 26項
看來我的歲差計算也存在問題
先回答Aquarian-Age的问题:
此处没有木星的资料，但我另有 github project 缯制星图或可供你参考。此 project 主要是绘制星图，木星资料须间接获取，最简单的方法是到 赤道星图网页，用 Date and Time 按键输入时间，然后从星图找出木星的符号(类似棕色的"4")，点击木星符号便弹出 popup box 显示木星的资料。popup box 的数据用 VSOP87 计算，据源码的说明，其精度在公元前1年至公元4000年约为1"。如果想要更准确或更长年限的数据，可用JPL horizons 网站，此网站用JPL DE431 历表计算行星位置，但岁差和章动的算法似乎还在用上世纪八十年代的公式计算。
gamccy:
閣下自言非專業人士，我在此重申:對於天文星曆表計算而言我也不是專業人士，您說以前使用過VSOP87/2013、Swiss Ephemeris、ELP/MPP、Pluto99等曆表，那您在這方面的經驗比我深，我的經驗止於VSOP87、ELP/MPP02及JPL的DE曆表，而且兩年前才開始玩的。說來慚愧，Pluto99曆表我是第一次聽到的。
您說月相時刻與我算出的偏差小於0.0002s，那就是一致了! 拙文說我定的收斂條件是1e-8日=0.00864s，偏差小於0.00864s 就完全一致了。而且我只用 double 來算，儒略日數整數部份已佔了七個數位，加上小數後8個位就已經15個數位了，9個位以後的數值不足信，那些只是 numerical garbage。
至於節氣時刻的差異可以探討一下，我相信是章動的計算有差異，有時間我會比較您的數據。您的章動怎樣算?是保留了一千多項嗎?還是有所簡化?
最後有一事不明，您把節氣時刻算準到0.2s有何用?須知太陽在0.2s平行0.008"，太陽視半徑約為半度或1800"，您已將太陽質心位置算準到視半徑的4e-6了，相信很多專業的天文研究亦無需如此精度，我聽說的地面天文研究項目嫌此精度不足者，就只有 pulsar timing array 一項。當然，派探測器探索太陽系天體亦需用高精度曆表，這也是JPL曆表的其中一個製表目的。
你的"章動矩陣積"應是算在 jd = 2305446.166666667 而不是之前的那個冬至吧，否則我們的差異會較大。從你給出的章動矩陣積可算出epsilon_A=23.491283425796度，Delta psi = 14.87451"。假設確是算在 jd = 2305446.166666667，則我們的 epsilon_A 在列出的十四個位都一致，我的Delta psi 是 14.87464"，我們在小數後第四個數位有偏差，但此偏差極微，相信之前那冬至時刻的偏差大致是0.003s，對嗎? 此偏差小於我預設的收斂條件，可以不理。
你說的 0.2s 偏差發生於何時? 若能找出來，我們再比較一下 epsilon_A 和 Delta psi 會有助於找出偏差根源。雖然追尋0.2s 偏差對於實際天文計算意義不大，但對於 debug 程式來說很有價值，你算出的不一定錯，可能是我算錯了。
我說:"你算出的不一定錯，可能是我算錯了" 竟然中了一半! 錯不在您，但我也沒有算錯，錯誤出現在拙文:公式(21)中F5的第三項7.74722" T^2 應改為 7.4722" T^2，(22)式F14最後一項應是0.00000538691 T^2，小數點後有五個0，我寫少了一個0。拙文會馬上改正，您看到這信息時或許已改正了。
我在編程式時是用 IERS 原文 Chapter 5 公式(5.43)及(5.44)，所以沒有算錯，撰寫拙文時不慎將數字寫錯了，我修改程式代入拙文錯誤的數字後得到的 Delta psi 與您算出的一致，可見錯不在您，對此向您致歉。請您修正F5和F14的計算，看看那 0.2s 偏差是否就消失了。
……
我在編程式時是用 IERS 原文 Chapter 5 公式(5.43)及(5.44)，所以沒有算錯，撰寫拙文時不慎將數字寫錯了，我修改程式代入拙文錯誤的數字後得到的 Delta psi 與您算出的一致，可見錯不在您，對此向您致歉。請您將數據修改，看看那 0.2s 偏差是否就消失了。
沒錯！是這個問題，修正該兩個輻角後我算出1600-3500年節氣時刻的差值與月相的差值幾乎完全一樣，與你TDBtimes檔案中數據的最大差值也是不超過0.0002s。
除了TDBtimes檔案外您還提供了-4000至8000年的EXTENDED檔案，我想既然如此應該也可以使用同樣的方式驗證該檔案的數據，於是直接解壓該數據文件執行我的程式，不過得到的結果是在-4000至8000年的時間跨度內我所計算的節氣時刻與您提供的數據差值最大已接近1500s，月相的最大差值已超過一個朔望周期，不過這並不奇怪，因為您編算檔案中說過IAU所薦的歲差模型僅在J2000前後千年內有效，超過這一區間精度會大減，只是我沒想到會有如此大的差額，因為首個朔日的計算依賴於首個冬至的時刻，如果冬至的時刻不夠準確，月相的時刻想必也會有很大偏差。
所以我也改用Vondr ́ak等人所推薦的歲差模型，計算後得到的數據差值有很大的收斂，與您EXTENDED檔案中的數據對比節氣的最大差值為3s，月相的差值不超過0.006s。
因為您的編算檔案中沒有提供Vondr ́ak等人推薦的歲差模型計算方式，我無法從該檔案中驗證我的計算結果是否與您一致。但是我參閱Vondr ́ak等人的文檔時發現其有提供詳細的計算程式以及測試用例。因為您有提及Vondr ́ak等人的歲差模型與IAU推薦的模型在計算上僅需替換三個角度的值即可，因此我僅需驗證這三個角度的計算結果即可。但是Vondr ́ak等人的計算程式中沒有直接的計算這三個角度的程式，不過其所提供的程式eg.1和 eg.2在計算過程中均有計算角度值的邏輯，因此我想如果實現這倆個程式計算出對應結果後再與其所提供的測試用例對比，如果結果一致則可以證明我計算角度的邏輯應該也是正確的。
編寫程式後我得到的結果與Vondr ́ak等人提供的測試用例結果在精度上達到15個浮點位的一致，因此我認為我計算角度邏輯應該是正確的。同時為了避免在拷貝表4/6的數據出錯，我也反復校對過兩次各位數字以及符號以及表4正餘弦角度值所對應的位置，確保表的數據無誤後我重新執行了一次計算，發現節氣和月相最大誤差依然與之前所計算的一致。以誤差值最大的7999年為例，該年所對應的jd0=4642633.166666667，以Vondr ́ak等人推薦的歲差模型計得歲差矩陣為：
該矩陣結果不知與您的結果是否一致呢？
我也想過傾角章動可能導致的誤差，但是我將傾角章動也計入後發現結果並沒有太大的差異。
至於節氣時刻的差異可以探討一下，我相信是章動的計算有差異，有時間我會比較您的數據。您的章動怎樣算?是保留了一千多項嗎?還是有所簡化?
章動的計算我也是使用的全部5.3a表的一千多項數據。
最後有一事不明，您把節氣時刻算準到0.2s有何用?須知太陽在0.2s平行0.008"，太陽視半徑約為半度或1800"，您已將太陽質心位置算準到視半徑的4e-6了，相信很多專業的天文研究亦無需如此精度，我聽說的地面天文研究項目嫌此精度不足者，就只有 pulsar timing array 一項。當然，派探測器探索太陽系天體亦需用高精度曆表，這也是JPL曆表的其中一個製表目的。
上面的Comment中我有提及，在我的應用場景中既需要時間跨度數千年的天體計算又不能佔用太大存儲空間，所以這需要在計算精度和資源佔用之間取一個平衡，當然我也嘗試過很多方式，比如最簡單的使用Jean Meeus在《Astronomical Algorithms》中提供的精簡過的VSOP87歷表、Swiss Ephemeris提供的精簡計算程式、以及我自行精簡的VSOP2013歷表等等。不管是何種精簡程式，都需要與一份權威的精確數據進行比對，以確保精簡後的精度能在一定範圍時間內。
最開始我使用 Swiss Ephemeris 所提供的檔案進行比對，但是其給出的檔案均為PDF文檔（其它的權威數據也大多是PDF文檔），無法很好地進行程式比對，我本想使用其提供的程式進行計算比對，但發現其計算依據依然是JPL的歷表，於是我就直接使用該歷表進行計算，說來慚愧，最開始使用JPL歷表計算時並沒有使用43x 的歷表，因為能找尋到的JPL歷表參考資料多為405歷表的，所使用的章動歲差模型也比較老舊，雖然得到的結果其實已經很精確。直至搜尋得您的Project，我開始根據您所提供的檔案以及章動歲差模型通過43x來計算數據。
除此之外，正如您編算檔案中所述，這些計算過程本身也很有趣，也能學到不少知識。
沒錯！是這個問題，修正該兩個輻角後我算出1600-3500年節氣時刻的差值與月相的差值幾乎完全一樣，與你TDBtimes檔案中數據的最大差值也是不超過0.0002s。
很好，這樣便說明我們的計算程式沒有 bug，如果我的計算公式沒錯，則整套算法應是沒有問題的。我的主要檢驗是 JPL horizons 網站及天文年曆的資料，得出的結果是一致的。不過 JPL horizons 網站還在用1980年代的歲差章動模型，得出的日月視黃經與我的算法相差約0.05"，天文年曆只給出準到1"的位置數據，所以沒有JPL horizons 網站有用，不過那裡有月相和二分二至時刻可比較。順便問一下，您提到紫金山天文臺的數據，未知在那裡可找到?我以前在網上看到了《2017年**天文年历》序言圖片，說歲差用IAU2006模型，即與我們用的一樣，章動似乎是用IAU2000B，是我們用的IAU2000A的簡化版，太陽系天體位置用 DE421/LE421 計算，未知2020年紫臺是否已改用新版本的DE曆表。
IAU所薦的歲差模型僅在J2000前後千年內有效，超過這一區間精度會大減，只是我沒想到會有如此大的差額
您可以看看Vondr ́ak 等人文章的 Fig. 5，就知道 IAU2006 模型過了其有效時段後或狂升或狂降，是完全不能用的。
以誤差值最大的7999年為例，該年所對應的jd0=4642633.166666667，以Vondrak等人推薦的歲差模型計得歲差矩陣為...
我的計算結果與您的稍有差異:
0.090165536831727 -0.918369572628242 -0.385314811600160
0.916615544695758 0.227820938922022 -0.328501998488269
0.389469022118506 -0.323565986870234 0.862333423306059
可否給出那三個歲差角，看看哪一個有偏差?
正如您編算檔案中所述，這些計算過程本身也很有趣，也能學到不少知識。
我有說過這樣的話嗎? 搜索拙文檔案方知原來在附錄那裡說的，自己早已忘記了。然則您連拙文的附錄也看了! 我撰該文本為方便自己日後參考，未想到竟然會有人仔細閱讀，實出乎意料之外。
可否給出那三個歲差角，看看哪一個有偏差?
依據式11/13及表4/6我計得jd0=4642633.166666667的三個歲差角為：
ψ = 1.4528501156194062
ω = 0.40301493466845972
χ = -0.030058050763585916
其中：
∑ψ = -0.054715068761488447
∑ω = -0.0054722701961134289
∑χ = -0.02998592307639926
ψ與ω有差異:
ψ = 1.452850720558418, ω = 0.40301525074386
其餘數值一致，這就是說公式(11)頭三項有偏差，這就好辦了，我怕是計算表4或表6項有偏差，那就要費多點力氣 debug 了。大家各自檢查一下計算，看看有沒有錯吧。
我檢查了我的程式，沒有發現錯誤，於是隨意玩玩，發現將(11)式的 T^3 項的正號改為負號後就得出您的數值。您的T^3項正負號是不是錯了?
我檢查了我的程式，沒有發現錯誤，於是隨意玩玩，發現將(11)式的 T^3 項的正號改為負號後就得出您的數值。您的T^3項正負號是不是錯了?
是這個問題！我太大意了！一直只覈對數據沒有注意符號的檢查。我直接拷貝的IAU的算式替換了各項但沒有同時修改符號。經更正後-4000至8000年內的差值與-1600至3500年內的一模一樣，也是最大不超過0.0002s，感謝！
順便問一下，您提到紫金山天文臺的數據，未知在那裡可找到?我以前在網上看到了《2017年**天文年历》序言圖片，說歲差用IAU2006模型，即與我們用的一樣，章動似乎是用IAU2000B，是我們用的IAU2000A的簡化版，太陽系天體位置用 DE421/LE421 計算，未知2020年紫臺是否已改用新版本的DE曆表。
紫金山天文台於2019年中印刷出版了一套新的《2020中國天文年歷》，我在最近的一些計算上有參考對比個中數據。不過我上面 Comment中提及的用於校驗節氣的2020節氣發生時刻數據來自于有趣天文館提供的2020年二十四节气精确时间。
紫金山天文台於2019年中印刷出版了一套新的《2020中國天文年歷》，我在最近的一些計算上有參考對比個中數據。
比較的結果怎樣?
我手上沒有《2020年**天文年历》，這裡的圖書館也沒有，我也不想花錢去買。不過我記起手上有《2015年**天文年历》，書中說用DE421曆表計算日月行星位置，歲差用IAU2006，章動用IAU2000B。於是拿來比較太陽視黃經，書中列出的黃經精度是0.1"，比較了幾天的黃經發現與書中的數據在0.1"精度範圍內一致。
不過我上面 Comment中提及的用於校驗節氣的2020節氣發生時刻數據來自于有趣天文館提供的2020年二十四節气精确時間。
感謝提供這網址。網頁列出了24節氣時刻準到0.01秒，還有準到1秒的72候的時刻，想是作者根據紫金山天文臺的資料計算的結果。
TDBtimes.txt的時刻經過TDB-UTC的修正後每一個2020年節氣時刻都比該網頁大2s，用DE421曆表重算也是這樣，太陽在2s平行0.08"，符合您發現的0.1"偏差。該網頁有2019年節氣時刻連結，驗算後發現我的2019年節氣時刻比他們大約1s，再檢驗2019年滿月時刻，發現我的時刻也比他們大約1s。這說明很可能我們的TDB-UTC數值不同，根據閏秒表，2019年和2020年的TT-UTC 都是69.184s，TDB-TT的差異小於0.002s應該可以不計，為慎重起見，我查了JPL horizons 網站提供的TDB-UTC 數據，2019-2020年的值在69.182311s 與 69.185680s 之間。也許不知何故他們的TDB-UTC在2019年取70s，在2020年取71s，這樣便可解釋偏差。我幾乎可以肯定這就是偏差的原因，稍後或會向那兩位作者請教。
您知道該文作者楊暘是男還是女嗎? 網上搜索了這名字，發現有男亦有女。
比較的結果怎樣?
結果還需本週我結束工作回到家中查閱書籍才能對比，結果我會在比較確認無誤後第一時間反饋到該 Comment。
您知道該文作者楊暘是男還是女嗎? 網上搜索了這名字，發現有男亦有女。
對方是先生還是女士我也不是太清楚，不過有趣天文觀的網站是由云南天文台抚仙湖观测站观测助手罗林先生搭建的，該網站也有其聯繫方式 MailTo:446469755@qq.com，大致可通過他聯絡到高良超、杨旸這兩位資料整理者。前些日時我也有通過E-Mail諮詢過其相關數據的來源及計算依據，不過到目前位置尚未有答複。
不過有趣天文觀的網站是由云南天文台撫仙湖觀測站觀測助手羅林先生搭建的，該網站也有其聯繫方式 MailTo:446469755@qq.com，大致可通過他聯絡到高良超、楊暘這兩位資料整理者。前些日時我也有通過E-Mail諮詢過其相關數據的來源及計算依據，不過到目前位置尚未有答複。
我是從楊暘女士/先生網頁找到聯絡電郵，我比您幸運，發電郵後13.5小時內得到了回覆:
我们这边用的是VSOP87星历，△T取值提前半年向紫金山天文台（以下简称"紫台"）问的，2020年为70.2s，还有就是舍入的关系，造成2秒不同原因，您那计算的2020/1/5/21:31:15.7是正确的，如果换算成紫台ΔT为2020/1/6/05:30:05.5。而我们采用的星历与DE431算的初始参数不同，会有时差个1秒钟左右。我们在ΔT取值方面为71.6秒，故而与您差2秒，但在力学时上时间一致。
如我所料，偏差主要是TDB-UTC取值不同，但是沒想到他們原來是用VSOP87來算"精確節氣時間"，然則所謂"紫金山天文臺數據"僅是△T而已，而且問了又不用。那70.2s 應該是2019年年中估計年末會有閏秒，結果沒有閏秒，TT-UTC在2020年仍是69.184s，但是不知何故他們不用 70.2s 而用莫明其妙的 71.6s，我已發電郵追問，且看何時會有回覆。
比較的結果怎樣?
花費了些許時日對比了《2020中國天文年曆》中的數據項，在該年度內的太陽視黃經誤差在[0.09, 0.16]"，月亮的視黃經誤差在[0.11, 0.17]"；我也嘗試計算對比了視赤經/緯的誤差發現也非常小。不過相較於太陽和月亮，其它系內大行星的視赤經/緯誤差稍微有所增大。
也許不知何故他們的TDB-UTC在2019年取70s，在2020年取71s，這樣便可解釋偏差。
如我所料，偏差主要是TDB-UTC取值不同，但是沒想到他們原來是用VSOP87來算"精確節氣時間"，然則所謂"紫金山天文臺數據"僅是△T而已，而且問了又不用。那70.2s 應該是2019年年中估計年末會有閏秒，結果沒有閏秒，TT-UTC在2020年仍是69.184s，但是不知何故他們不用 70.2s 而用莫明其妙的 71.6s
我手上沒有2019的天文年曆，遂不知紫金山天文台于2019年的天文年曆計算使用的ΔT為何值，不過我這裏有2020年的天文年曆，2020的天文年曆所用的ΔT數值卻為69.5s。而且該年在計算系內大行星時仍然使用的DE/LE421歷表至2014年來沒有變化。
先說ΔT。
還未收到楊暘老師的回覆，估計不會回覆了，不過我已猜到了那71.6s的來歷。
Espenak and Meeus Delta T 網頁有2005年-2050年的預推公式:
ΔT = 62.92 + 0.32217 * t + 0.005589 * t^2 其中 t = y-2000
代入t=20即得71.6s。那網頁 last update 是2007年10月3日，用這預推公式去算2020年的ΔT與實際數值有幾秒偏差不足為奇。有新資料不用，而去用十幾年前的預推公式，真使人費解。
還有一個嚴重問題:Espenak 和 Meeus 網頁的ΔT是用來計算日食，所以應是指 TT-UT1，因為天文觀測須用UT1而不是UTC，例如計算恆星時就要用UT1而不用UTC。但"精確節氣時刻"註明是用北京時間，北京時間是民用時UTC+8，雖說UTC與UT1相差不過0.9s，但"精確節氣時刻"列出了0.01s的精度，UTC與UT1的差異便不可忽略了。
根據IERS Bulletins的數據，2020年的UT1-UTC在 -0.25s 附近，故TT-UT1 = (TT-UTC) - (UT1-UTC) 在69.4s 附近。《**天文年历》的ΔT就是指TT-UT1，其2020年預推值 69.5s 與實測值頗為接近，但對於"精確節氣時刻"給出的北京時間 UTC+8 卻不相干。
由於《**天文年历》列出的日月行星位置都以TDB為準，比較數據時不須要考慮時間轉換。
不過相較於太陽和月亮，其它系內大行星的視赤經/緯誤差稍微有所增大。
這有點奇怪，我比較了《2015年**天文年历》八大行星的視赤經和視赤緯，我驗算的幾個例子偏差都很小(一般在0.01"內)。能否給出兩三個偏差比較大的例子讓我也算一算?
這有點奇怪，我比較了《2015年**天文年历》八大行星的視赤經和視赤緯，我驗算的幾個例子偏差都很小(一般在0.01"內)。能否給出兩三個偏差比較大的例子讓我也算一算?
以《2020年中國天文年歷》中天王星的視赤經/緯為例：
| 日期 | 視赤經 | 視赤緯 | 視赤經誤差 | 視赤緯誤差 |
|---|---|---|---|---|
| 4/9 | 2h 13m 49.325s | 12° 57' 47.63" | 0.24s | 0.31" |
| 7/1 | 2h 30m 25.125s | 14° 20' 28.91" | 0.30s | 0.19" |
| 11/3 | 2h 25m 27.292s | 13° 55' 02.68" | 0.27s | 0.30" |
| 11/19 | 2h 22m 57.053s | 13° 42' 08.96" | 0.15s | 0.22" |
以上則為天文年歷中數據與我所計得的結果誤差較大的幾項。
這有點奇怪，我比較了《2015年**天文年历》八大行星的視赤經和視赤緯，我驗算的幾個例子偏差都很小(一般在0.01"內)。能否給出兩三個偏差比較大的例子讓我也算一算?
以《2020年中國天文年歷》中天王星的視赤經/緯為例：
日期 視赤經 視赤緯 視赤經誤差 視赤緯誤差
4/9 2h 13m 49.325s 12° 57' 47.63" 0.24s 0.31"
7/1 2h 30m 25.125s 14° 20' 28.91" 0.30s 0.19"
11/3 2h 25m 27.292s 13° 55' 02.68" 0.27s 0.30"
11/19 2h 22m 57.053s 13° 42' 08.96" 0.15s 0.22"
以上則為天文年歷中數據與我所計得的結果誤差較大的幾項。
非天文年歷的數據項問題，而是我在合併代碼時誤將一個周期項的值修改成錯誤值，該操作同時也影響了之前與您核算的結果，經修正後誤差範圍收斂均在0.01s/"左右
日期 視赤經 視赤緯 視赤經誤差 視赤緯誤差
4/9 2h 13m 49.325s 12° 57' 47.63" 0.24s 0.31"
7/1 2h 30m 25.125s 14° 20' 28.91" 0.30s 0.19"
11/3 2h 25m 27.292s 13° 55' 02.68" 0.27s 0.30"
11/19 2h 22m 57.053s 13° 42' 08.96" 0.15s 0.22"
11/19的視赤緯恐怕有誤，列出的數值應是11/20的視赤緯。
以下是我用DE421及 IAU2006/2000A 的計算結果:
| 日期 | 視赤經 | 視赤緯 | 地心距離 |
|---|---|---|---|
| 4/9 | 2h 13m 49.3264s | 12°57'47.6355" | 20.76712127718 |
| 7/1 | 2h 30m 25.1252s | 14°20'28.898" | 20.28828319292 |
| 11/3 | 2h 25m 27.2948s | 13°55'02.66205" | 18.78876433087 |
| 11/19 | 2h 22m 57.0507s | 13°42'52.1846" | 18.84165958335 |
| 11/20 | 2h 22m 48.1406s | 13°42'08.93106" | 18.84753192604 |
可見與《2020年中國天文年曆》數據相差小於0.04"。我在檢查《2015年**天文年历》發現地心距離的計算值與天文年历的數據完全一致，這也是意料之中，因為距離直接用DE421曆表計算，視赤經和視赤緯要加上光行時、光行差、frame-bias matrix、歲差和章動修正，只要有一項計算與天文年曆計算有偏差都會影響計算結果。
《**天文年历》的歲差和章動計算和我們的計算除了章動稍有差異外是一樣的。《2015年**天文年历》並沒有說明行星光行差的計算方法，只是在"太陽表"說用(-20.49552"/地球向徑)來算太陽的光行差修正，這明顯是近似算法。據第三版 Exploratory Supplementary to the Astronomical Almanac 第267頁所述，此公式的誤差是0.001"，但公式沒有計及行星對地球軌道的重力攝動，實際誤差是多少我未有深究。書中也說到行星光行差有幾種近似算法，不知《**天文年历》採用什麼算法。
英美的天文年曆自2015年起改用DE430計算。以下是DE430及 IAU2006/2000A 的計算結果:
| 日期 | 視赤經 | 視赤緯 | 地心距離 |
|---|---|---|---|
| 4/9 | 2h 13m 49.322s | 12°57'47.6593" | 20.76712467528 |
| 7/1 | 2h 30m 25.1205s | 14°20'28.9229" | 20.28828630035 |
| 11/3 | 2h 25m 27.2898s | 13°55'02.68925" | 18.78876773606 |
| 11/19 | 2h 22m 57.0457s | 13°42'52.2117" | 18.84166307364 |
| 11/20 | 2h 22m 48.1356s | 13°42'08.95815" | 18.84753542145 |
可見兩曆表算出的數據有微小差異，我估計《2020年**天文年历》的地心距離完全符合DE421算出的結果。
此外，廣義相對論的光線重力偏折也會影響行星視位置。重力偏折就是1919年的日食觀測驗證了廣義相對論的推算，使愛因思坦一舉成名的那個相對論效應。下圖摘自第三版 Exploratory Supplementary to the Astronomical Almanac 第273頁圖7.5:
圖中可見當天王星與太陽角距離(elongation)小於90°時，重力偏折會達到0.01"以上，但《**天文年历》應該沒有計及這重力偏折效應。
可見要將行星位置列出準到0.01"的精度，必須說明計算用的曆表、時間標準、座標系統、包括了哪些修正項及計算方法，否則是沒有多大意義的。所以我對那些將月相節氣時刻列出準到0.01s而又不說明計算方法的做法很不以為然，太陽在0.01s平行0.0004"，月亮平行0.0055"。
非天文年歷的數據項問題，而是我在合併代碼時誤將一個周期項的值修改成錯誤值，該操作同時也影響了之前與您核算的結果，經修正後誤差範圍收斂均在0.01s/"左右
很好! 這說明我們的計算結果大致相同了。我也是經常將代碼挪動和更改，有時也會不慎犯錯，所以 unit testing 是很重要的。
Wow, I love your work. Your TDBtimes.txt file is fantastic. Thank you very much.
Wow, I love your work. Your TDBtimes.txt file is fantastic. Thank you very much.
Thanks. Glad to know that some people took a closer look at this project and appreciated the data.
I have always wanted to compute the 24 solar terms but it is a lot more complicated than I thought. I thought I could do this in Excel or Matlab. In the end, I gave up and just used the TDBtimes.txt file. I'm new to programming and coding, but I've published some feng shui app, although all are using Swift. Is there a way I could contact you directly for work collaboration?
Is there a way I could contact you directly for work collaboration?
You can simply email me, but I am not a fan of feng shui.
I have always wanted to compute the 24 solar terms but it is a lot more complicated than I thought. I thought I could do this in Excel or Matlab. In the end, I gave up and just used the TDBtimes.txt file.
TDBtimes.txt is only useful if you want to know the TDB times of 24 solar terms (defined according to the celestial coordinate system specified by the IERS 2010 convention) to about one second precision. Such a precision is of no astronomical significance, but is necessary for the calculation of modern Chinese calendar in accord with the GB/T 33661-2017 document. I will be amused if feng shui calculations require such an accuracy.
I bet you can find simple and fairly accurate formulae for the times of 24 solar terms. Without trying hard, I managed to derive a simple fitting formula for the solar term times accurate to within 20 minutes for years between 1900 and 2100. The root-mean-square error is 7 minutes. This is already more accurate than the best solar model used in the Qing dynasty. With a little more effort, I believe it's possible to derive a relatively simple fitting formula accurate to one minute. I'm just not interested in doing that since I've already got accurate times and the data are small, but you're welcome to try. Just take the data in TDBtimes.txt and use regressions or other statistical methods to fit the data to whatever precision you want.