有關(guān)“現(xiàn)代歷法”——主要編者：songgz

白發(fā)殘翁 2012-12-05

展開全文

有關(guān)“現(xiàn)代歷法”——主要編者：songgz

　　幾個月前，增編寫個一個農(nóng)歷計算javascript程序，那是我的第一個農(nóng)歷程序。為了實現(xiàn)其中的算法，曾花費了好幾個星期的時間研究天文計算相關(guān)的原理，當(dāng)我算出結(jié)果之后，對程序的結(jié)果仍然沒有信心。不過，在“春光”老師介紹《天文算法》一書之后，我認(rèn)真閱讀并翻譯全書之后，我知道，我的算法基本上沒有錯誤，但同時也認(rèn)識到，天文學(xué)家們的算法的確高明，比我當(dāng)時的算法好得多。如果你也對天文計算感興趣，那就很有必要學(xué)習(xí)他們的先進思想方法。

　　當(dāng)然，重要的是他們的思想方法以及相關(guān)的理論，而不一定是他們的計算結(jié)果。因為，《天文算法》一書是早期出版的，有些數(shù)據(jù)比較舊，造成精度不一定很好。如果你對計算精度要球特別高，可能需要更換一些數(shù)據(jù)，但數(shù)據(jù)的處理的思想方法及技巧是一樣的。

　　國內(nèi)不少網(wǎng)友對天文計算比較感興趣，卻又不知如何下手。問題出在哪里呢？為什么苦苦研究幾個月甚至幾年也沒有進展？主要原因是國內(nèi)有關(guān)的書籍太少，業(yè)余條件下沒有機會學(xué)習(xí)到天文計算的理論。在我們國家，有不少機構(gòu)對天文學(xué)有深入的研究，比如天文臺、一些大學(xué)等，可他們不太愿意出版天文算法之類的書籍(這類書籍銷量少，出版要賠錢的)。既然如此，就讓我們自已想辦法解決問題吧。

　　首先，我認(rèn)為需要掌握一定的計算機程序設(shè)計技術(shù)，不要求很利害，但起碼也要有幾個月的程序設(shè)計經(jīng)驗。如果你不會程序設(shè)計，那你只能用Excel或計算器之類的工具來處理計算問題，那簡直是在浪費你的生命。

　　其次，應(yīng)掌握《高等數(shù)學(xué)》中的一些知識。比如：極限、導(dǎo)數(shù)、微分、積分、極值問題、求根問題、最小二乘法、向量數(shù)學(xué)等。當(dāng)然，我們更多使用高中的《立體幾何》、《解析幾何》、《函數(shù)》、《三角函數(shù)》等有關(guān)知識。還應(yīng)了解《球面三角學(xué)》里的幾個公式。

　　其三，《數(shù)值方法》這類書籍是必須讀的。

　　其四，需對物理學(xué)有所了解，尤其是運動學(xué)相關(guān)知識。當(dāng)然，如果你想用數(shù)值積分的方法解決天文計算問題，《理論力學(xué)》甚至是《天體力學(xué)》也是有必要了解的。

　　對于多數(shù)具有大學(xué)學(xué)歷的人來說，基本具備以上知識，也就是說，只要你有興趣，就完全可以進行天文計算。從本質(zhì)上講，日月運動、行星運動主要使用“牛頓力學(xué)”及數(shù)學(xué)方法(如微積分)，在牛頓那個年代，力學(xué)理論、數(shù)學(xué)理論、計算工具等都不可能和現(xiàn)在相比，在那個年代，就連“除法問題”主要是“教授”們才能掌握的！他們可以計算天體運動問題，我們?yōu)槭裁淳筒豢梢阅?？只要有信心，或多或少可以解決問題。

　　相反，如果根本沒有學(xué)過《高中數(shù)學(xué)》《高中物理》《高等數(shù)學(xué)》，我建議你還是花點時間學(xué)習(xí)一下(弄不好要花費一兩年時間)，否則，即使算出了結(jié)果，也很難對你算出的結(jié)果形成理性的認(rèn)識，甚至感性的認(rèn)識也談不上。如果你不想學(xué)習(xí)這些“無用”的東西，能不能實現(xiàn)天文算法，或許可以：通過閱讀別人的程序。

　　在太陽系中，我們一致認(rèn)為九大行星(八個主行星和1個冥王星)圍繞太陽系運動。幾百年前，哥白尼先生就已經(jīng)這么認(rèn)為！現(xiàn)在的教科書中也都是這么說的！因為太陽系中的星體圍繞太陽運動，所以我們常以太陽為中心來考慮問題。

　　以太陽為中心考慮問題有不少好處：我們可以把行星圍繞太陽運動近似看做圓周運動(圓周運動是物理學(xué)中的一個重要名詞)。這樣，行星的運動顯得非常有規(guī)律。

　　如果不以太陽為中心，而以地球為中心，我們會發(fā)現(xiàn)水星、金星的運動毫無規(guī)律，你可以充分發(fā)揮你的空間想象，你將發(fā)現(xiàn)，水星和金星在天球上的視運動不可能是圍繞地球做圓周運動，而是在太陽附近“擺動”。

　　然而，我們不可能站在太陽上面進行天文觀測，通常是站在地球上進行天文觀測的，所以“地心說”也是非常科學(xué)的，也就是說，我們可以認(rèn)為，太陽系的中心是地球。

　　不管是“日心說”、“地心說”，只要你科學(xué)對待，不要隨意附加“神學(xué)的”或“迷信的”論調(diào)，都可以正確解天文現(xiàn)象。

　　要記住，這個世界上沒有絕對的“中心”。太陽個頭大，質(zhì)量大，所以行星圍繞它運轉(zhuǎn)。僅因它質(zhì)量大，就能決定一切嗎？顯然不是。我們可以這么想，某件事，是大人說了算還是小孩說了算，通常是大人說了算，但具體問題要具體分析，有時候小孩說了算。

　　從純牛頓力學(xué)的角度看，不管是“日心”還是“地心”，通過計算，均能得到相同的數(shù)值計算結(jié)果。

　　因此，根據(jù)實際需要，有時我們使用“地心”坐標(biāo)，有時候使用“日心”坐標(biāo)。

　　為了計算方便以及更好的解釋天文現(xiàn)象，人們還引入了“太陽系質(zhì)心”、“銀心”、地月系質(zhì)心、月心、行星的“心”、純幾何的視中心，等等，一系列的“中心”。不管是什么“中心”，我們必須正確理解這些“中心”的具體定義，必須明白這些“心”到底在哪里，以便進行坐標(biāo)變換。如果知道這些心的定義，我們未必原原本本的使用現(xiàn)成的天文坐標(biāo)變換公式，也可適當(dāng)根據(jù)程序設(shè)計的需要建立自已所需的坐標(biāo)變換公式(比如在近似處理時)，使程序更加簡潔明了。

　　茫茫宇宙，無邊無際——宇宙大得讓我們難以想象。宇宙是球形的嗎？我們很難說清楚這個問題，但是，進入我們視野的天空，就象一個“球”，我們稱它為天球，我們位于天球的中心。這個球的半徑有多大？半徑很大很大，看成無窮大也無妨。太陽系內(nèi)星體之間的距離根本不能與天球半徑相比，如果把天球半徑看為1，那么太陽系內(nèi)星體間的距離可以看做0。在我們的視野中，所有的恒星都在天球上。

　　赤道的概念已被大家熟知?，F(xiàn)在，我們擴展一下赤道的概念。赤道平面擴展到天球，與天球相交，交線為天球上的一個大圓，這個圓稱為天赤道。

　　不管以太陽為中心，還是以地球為中心，我們看到的天赤道是相同的，所有恒星與天赤道的相對視位置關(guān)系是一樣的。為什么呢？因為，從幾何學(xué)看來，太陽系的尺度太小，不能與天球半徑相比，在太陽系中的任何一個地方，都可以看做天球的中心，所以，恒星之間的視位置關(guān)系以及它們與天赤道之間的視覺位置關(guān)系不會因為中心的改變而發(fā)生改變。

　　當(dāng)然，恒星并不是真的在天球的表面上，個別恒星與我們之間的距離不是很遠，所以，當(dāng)觀測點(即“中心”)改變時，它們的視位置也會有點變化，這就是我們常說的恒星周年視差，我們以后再討論這些問題吧。

　　地球公轉(zhuǎn)的軌道面，同樣可以擴展到天球，得到的交線也是個大圓，稱為黃道。

　　雖然赤道與天赤道是不同的兩個概念，但具體問題所涉及的是“天赤道”還是“赤道”往往是非常明確的(明顯的)，一般不會產(chǎn)生二意，所以通?！疤斐嗟馈币埠喎Q為“赤道”。

　　天赤道與黃道必然在天球上產(chǎn)生兩個交點。一個稱為“春分點”(太陽經(jīng)過升交點)，一個是“秋分點”(太陽經(jīng)過降交點)?！按悍贮c”是天球上重要的參考點。我們描述天體位置，通常是相對于這個參考點而言的。

　　嚴(yán)格上說，黃道是地月質(zhì)心公轉(zhuǎn)軌道在天球上的擴展。

　　在天球上，并沒有一個恒星可以直接用來標(biāo)定春分點的位置，它是看不見也摸不著的，我們?nèi)绾未_定它的位置？最重要的一種方法是：通過數(shù)以千計的恒星位置，反推出春風(fēng)點在天球上的位置，我們常說的FK5天球坐標(biāo)系統(tǒng)就與它有關(guān)。還有其它方法可以確定春分點，比如動力學(xué)方法。由于所用的方法不同，得到的春風(fēng)點也會有一些微小的差別。難道有很多種春風(fēng)點嗎？其實，“春風(fēng)點”這個名詞最好不要亂用，如果說FK5中的那個“點”是真正的春分點，那么動力學(xué)方法(如VSOP87)得到的那個點就不應(yīng)該稱做春分點，其實在《天文算法》一書中很少用“春風(fēng)點”一詞，取代之“分點”一詞。當(dāng)然，我們不一定計較那么一點點微小的差別，統(tǒng)稱為“春分點”也不要緊。

　　地平線——天地的交線，它又是一個圓。圓中心位于觀測點。考慮到天球很大，我們可以把圓(以及圓中心)平移到地心或日心，這樣，黃道、赤道、地平線的中心就相同了。也許你會問，平移后，會不會影響日月在地平坐標(biāo)中的位置？當(dāng)然會有一點影響，不過我們可以通過視差修正的方法來解決問題。

　　當(dāng)我們在天球上標(biāo)定了赤道、黃道、分點之后，我們接著來看日、地在黃道坐標(biāo)中的運動。以太陽系中的任意一點為天球的中心都不影響分點及黃道位置，當(dāng)天球中心移到地心，地心直角坐標(biāo)變?yōu)榱?，沒有經(jīng)緯度，但可以得到太陽的經(jīng)緯度，此時，如果跑到日心看地球，可以得到地球的經(jīng)緯度，這兩個經(jīng)緯度正好相反(日地連成的直線與天球相交的兩點正好完全相反)：

　　地球經(jīng)度 = 太陽經(jīng)度+180度，地球緯度 = -太陽緯度。

　　歷元、攝動、自球自轉(zhuǎn)軸、天北極、黃極、黃軸、黃經(jīng)、黃緯、赤經(jīng)、赤緯，這幾個概念不再講述，請在網(wǎng)絡(luò)上找一些資料看看吧！

　　春分點是黃道與赤道的交點！

　　春分點是天球中的重要的參考點，遺憾的是，春分點會移動！

　　在圖中，有兩條黃道、兩道赤道。

　　一條的歷元2000.0時刻黃道，一條是當(dāng)日(Date)黃道。

　　一條的歷元2000.0時刻赤道，一條是當(dāng)日(Date)赤道。

　　隨時間推移，黃道為什么會發(fā)生變化？因為地球受至大行星的攝動。

　　隨時間推移，赤道也會發(fā)生移動，因為受月球影響，地球自轉(zhuǎn)軸的方向發(fā)生移動，造成天北極圍繞黃軸緩慢移動。找一個直鐵棒，鐵棒放置的方向與黃軸相同，用電焊機把鐵棒焊在赤道面的中心。始終保持鐵棒與黃軸平行，并轉(zhuǎn)動鐵棒，這時你將看到，北極圍繞鐵棒做圓周運動，同時春分點也同步移動，當(dāng)北極運行一周(約26000年)，春分點在黃道上也移動一周。

　　春分點移動的后果：黃經(jīng)、赤經(jīng)的起算點是春分點，造成所有星體的黃經(jīng)、赤經(jīng)發(fā)生變化。假設(shè)，黃道不動(沒有發(fā)生變化)，赤道發(fā)生變化，引起春分點移動，如果春分點在黃道上向西移動了1度(在上圖中，往左下方移動)，那么，所有星體的黃經(jīng)增加了1度！

　　春分點的移動量稱為歲差。歲差有很多種表達方法，最著名、最傳統(tǒng)的歲差是指黃道上的歲差。

　　上圖中，有三個重要的交點：γo，γ1，γm，我們統(tǒng)稱它為分點。γo稱為“歷元2000.0標(biāo)準(zhǔn)分點”，其實就是歷元2000時刻的春風(fēng)點，γm稱為“Date分點”或“當(dāng)日分點”或“當(dāng)日春分點”。

　　圖中，兩個黃道的交點是N，定義N到γo與N到γo'的距離相等，那么以γo'起算的黃經(jīng)則與γo起算的黃經(jīng)是相等的(除非天體過份靠近黃極)，對于行星，它在黃道附近，對于大部分恒星，也是遠離黃極的，所以γo'對于行星、月球、大部分恒星來說，是很重要的一個參考點。比如太陽，當(dāng)它從γo出發(fā)，轉(zhuǎn)動數(shù)周后達到了γo'，所經(jīng)歷的時間是整倍數(shù)的太陽繞地球運轉(zhuǎn)的周期。這種運行周期具有比較嚴(yán)格的力學(xué)意義。恒星在天空中的位置幾乎是不變的，以它作為參考系，具有“慣性”系的特征，γo則是利用恒星確定的，也具有慣性特點。任意時刻，我們利用某一顆或幾顆恒星的位置，輕松的標(biāo)定出γo的位置(當(dāng)然，首次建立FK5系統(tǒng)時，需要大量恒星標(biāo)定出γo，以提高精度，所以工程巨大)，當(dāng)太陽運行到γo'，可以認(rèn)為太陽相對于慣性參考點γo移動了一周或數(shù)個1周，即相對于慣性參考點移動了數(shù)個周期，這就是力學(xué)意義上的周期。γo與γo'相差很小，所以太陽到了γo與太陽到了γo'，太陽在星空的中相對位置幾乎相同，因此，太陽經(jīng)過這兩點所經(jīng)歷的時間是整倍數(shù)的恒星年(以恒星為參考，太陽在天球中運轉(zhuǎn)一周的時間)，由于這個原因，恒星年一般認(rèn)為就是太陽(或地球)在慣性系中的運動周期，“在慣性系中的運動周期”常稱為真周期。

　　p，γo'到γm的角距離就是Date黃道上的歲差
　　ψ，γo到γ1的角距離是J2000黃道上的歲差
　　χ，γ1到γm是Date赤道上的歲差
　　εo，是歷元2000.0的黃赤交角
　　ε，是當(dāng)日的黃赤交角
　　ω， Date平赤道與J2000黃道的夾角
　　π，是兩個黃道之前的夾角
　　П，是N到γo的角距離，它等于N到γo'的角距離
　　
　　另一組常用的歲差計算參數(shù)
　　90°-ζ，γo到Q的角距離，是一個常用的歲差計算參數(shù)
　　90°+z， γm到Q的角距離，是一個常用的歲差計算參數(shù)
　　θ是兩個赤道面的夾角

　　也許你會問，怎么會有這么多歲差參數(shù)，有必要嗎？回答，為了坐標(biāo)變換方便而已。通常，只需只道3個歲數(shù)參數(shù)，就可利用坐標(biāo)變換方法及多項式插值法導(dǎo)出其它的歲差參數(shù)。

　　一組常用的歲差表達，取自IAU2000
　　P:(     0.0,        41.99604, 19.39715, -0.22350, -0.01035, 0.00019, 0.0,     0.0)
　　Q:(    0.0,      -468.09550,   5.10421, 0.52228, -0.00569, -0.00014, 0.00001, 0.0)
　　π:(     0.0,       469.97560, -3.35050, -0.12370, 0.00030, 0.0,      0.0,     0.0)
　　II:(629543.988,   -8679.218,    15.342,    0.005,   -0.037,   -0.001,    0.0,     0.0)
　　p:(     0.0,     50287.92262, 111.24406, 0.07699, -0.23479, -0.00178, 0.00018, 0.00001)
　　θ:(     0.0,     20041.90936, -42.66980,-41.82364, -0.07291, -0.01127, 0.00036, 0.00009)
　　ζ:(     2.72767, 23060.80472, 30.23262, 18.01752, -0.05708, -0.03040, -0.00013, 0.0)
　　z :(    -2.72767, 23060.76070, 109.56768, 18.26676, -0.28276, -0.02486, -0.00005, 0.0)
　　ε:( 84381.40880, -468.36051, -0.01667, 1.99911, -0.00523, -0.00248, -0.00003, 0.0)
　　ω:( 84381.40880,    -0.26501,   5.12769, -7.72723, -0.00492, 0.03329, -0.00031,-0.00006)
　　ψ:(     0.0,     50384.78750,-107.19530, -1.14366, 1.32832, -0.00940, -0.00350, 0.00017)
　　χ:(     0.0,       105.57686,-238.13769, -1.21258, 1.70238, -0.00770, -0.00399, 0.00016)
　　上表是關(guān)于t的多項式系數(shù)，t等于J2000起算的儒略千年數(shù)，為儒略千年=365250天
第一列對應(yīng)t的0次方，第二列對應(yīng)t的1次方，第三列對應(yīng)t的2次方，其余類推。
如π=0.0*t^0 + 469.97560*t^1 - 3.35050*t^2 + ...
結(jié)果的單位是角秒。

　　以上歲差表達是非常精確的，誤差小于1毫角秒。

　　不過t不能太大，t應(yīng)小于1。

　　如果t=10，那么，會有什么后果？比如計算θ，t的7次方項為0.00001*10000000=100角秒！實際上，表達式系數(shù)的舍入誤本身就是0.00001，所以誤差高達100角秒。

　　如果由于某種需要，你非要計算t=10的情況，最好把t的5、6、7次方項忽略。這樣，如果t=5，誤差也可控制在1角秒之內(nèi)。

　　順便說一下，儒略千年數(shù)是一種時間計數(shù)方法。日常生活中，時間的計數(shù)法使用“年月日時分秒”，在天文計算過程，這種時間表達是很不方便的，原因是明顯的，需要6個數(shù)才能表達一個時間。

　　天文學(xué)中，采用2000年1月1日12時0分0秒時作為計時的起點，并連續(xù)計時，單位是天。1天日86400秒。1秒是多少？1秒時間長度由原子鐘得到，現(xiàn)代天文學(xué)的時間基礎(chǔ)已不再使用星體運動位置來確定，改用原子鐘，原子鐘的精度要好得多。

　　為什么又稱為“儒略”千年數(shù)。這還要談到另一個問題，一年有多少天？你可能立刻回答到：365.2422日。然而，在現(xiàn)代天文學(xué)家眼里，一年的長度是變化的，并不是固定為365.2422日，何況地球運動時還受到各種干擾，365.2422只是一定時期內(nèi)的平均值。因此，任意一年某時刻加上365.2422后，并不是下一年的同一時刻，這樣，把365.2422看做1個單位意義不大，天文學(xué)家更原意把365.25看做1個單位，它與公歷的置閏規(guī)則具有更好的“兼容”性，便于日歷轉(zhuǎn)換。365.25天是儒略年的長度。儒略1千年就是365250天。

　　再次強調(diào)，1天=86400原子秒，與地球自轉(zhuǎn)沒有必然的關(guān)系。在歷史上，曾使用某一年的地球公轉(zhuǎn)與自轉(zhuǎn)來定義秒長，之后，地球自轉(zhuǎn)速度發(fā)生改變，1天不等于86400秒了。注意，此“天”非彼“天”，前面的天是時間單位，后面的天指平均自轉(zhuǎn)一周所用的時間長度。然而，人類的日常生活與地球自轉(zhuǎn)息息相關(guān)，日期與時間的計數(shù)必須盡可能與地球自轉(zhuǎn)同步，這就造成，我們有兩套基本的時間系統(tǒng)，一個是手表時(或稱為協(xié)調(diào)世界時或稱為UTC時)，一個是力學(xué)時。力學(xué)時又時什么東西呢？力學(xué)時具有鮮明的動力學(xué)色彩，假如我們使用2000年地球公轉(zhuǎn)一圈所用時間定義為1個單位，2萬年后，公轉(zhuǎn)周期可能發(fā)生了變化，如果不通過力學(xué)原理推算，我們將無法得到那時候一年的長度，但是，通過力學(xué)方法，我們可以得到2萬年后公轉(zhuǎn)的周期，比如是1.00001(亂寫一個)個單位。簡單的說，當(dāng)我們定義了時間單位，我們就可以使用力學(xué)原理推算天體在某一位置所對應(yīng)時間，這就夠成時間基礎(chǔ)(根據(jù)星體位置確定時間)。原子時就簡單多了，它不是通過天體位置來確定時間的，而是通過原子鐘內(nèi)部的振蕩次數(shù)來來確定時間。原子鐘曾與力學(xué)時對準(zhǔn)過，秒長是相同的，起點相差零點幾秒鐘。

　　理論上，力學(xué)時、原子時是非常均勻的。手表時與地球自轉(zhuǎn)(以恒星為參考點確定自轉(zhuǎn)速度)同步，而地球自轉(zhuǎn)時快時慢，所以不均勻。

　　手表時(去除地區(qū)的時差后的手表時,即通過跳秒與UT1同步的時間)與原子時的差稱為ΔT，當(dāng)然手表時含有時區(qū)信息，如北京時間多了正好8個時。ΔT可由一組多項式表達式計算得到。（UT1同由自轉(zhuǎn)與恒星(春風(fēng)點)確定的，是天文意義的時間系統(tǒng)，嚴(yán)格說，UT1與時角有關(guān)，即與角度有關(guān)；原子時則與“振蕩次數(shù)”有關(guān)；手表時UTC的秒長與時角無關(guān),它與原子時的秒長完全相同。不過，UTC通過年終的跳秒實現(xiàn)與UT1同步，這樣就造成UTC與原子時不同步。）

　　手表時的秒長是原子秒，為什么手表時不是原子時呢？因為國際上有人動了手腳，在每年的最后一天做了“正或負(fù)跳秒(正負(fù)閏秒)”處理，實現(xiàn)手表時與地球自轉(zhuǎn)同步。

//原子時與UT1的差
//部分參考了美國國家航空航天局網(wǎng)站里的算法,www.
//2050以后的外推參考了skymap10.5,在skymap的星歷表中分段找?guī)讉€數(shù)據(jù),再用最小二乘法擬合
double deltatT(double y){ //輸入年份(可帶小數(shù)點),如2000.78
static double DTxs[15][12]={//{范圍下,范圍上,偏移,時間單位,多項各系數(shù)}
   {-9999,-500,1820,100,-20,      0,        32},
   {-500, 500, 0,   100, 10583.6,-1014.41, 33.78311, -5.952053, -0.1798452,   0.022174192, 0.0090316521},
   {500, 1600,1000,100, 1574.2, -556.01,   71.23472, 0.319781, -0.8503463, -0.005050998, 0.0083572073},
   {1600, 1700,1600,1,   120,    -0.9808, -0.015320, 1.0/7129},
   {1700, 1800,1700,1,   8.83,    0.1603, -0.0059285, 0.00013336,-1.0/1174000},
   {1800, 1860,1800,1,   13.72, -0.332447, 0.0068612, 0.0041116, -0.00037436, 1.21272e-5, -1.699e-7, 8.75e-10},
   {1860, 1900,1860,1,   7.62,    0.5737, -0.2517540, 0.01680668,-0.0004473624,1.0/233174},
   {1900, 1920,1900,1, -2.79,    1.494119,-0.0598939, 0.0061966, -0.000197},
   {1920, 1941,1920,1,   21.20,   0.84493, -0.0761000, 0.0020936},
   {1941, 1961,1950,1,   29.07,   0.407,   -1.0/233,   1.0/2547},
   {1961, 1986,1975,1,   45.45,   1.067,   -1.0/260, -1.0/718},
   {1986, 2005,2000,1,   63.86,   0.3345, -0.060374, 0.0017275, 0.000651814, 0.00002373599},
   {2005, 2150,2000,100, 64.723,-11.050,    232.182, -86.08},
   {2150, 9999,2000,100, 59.033, 103.4335, 28.980314, 0.000276}
};
//計算TD-UT,力學(xué)時與世界時之差
double *p,dt=0,u;
int i,j;
for(i=0;i<15;i++){
   p=DTxs[i];
   if(y<p[0]||y>p[1]) continue;
   y=(y-p[2])/p[3];
   for(j=4,u=1;j<12;j++) dt+=u*p[j],u*=y;
   return dt;
}
return 0;
}

我們在第(五)貼已經(jīng)談到了赤道與黃道的變化將引起歲差。受到月球運動的影響，地球自轉(zhuǎn)軸發(fā)生圍繞黃軸的緩慢的長期的圓周運動，造成赤道發(fā)生變化，進而引起春風(fēng)點變化。從上面的那個歲差示意圖中不難發(fā)現(xiàn)，春風(fēng)點的變化不單單是赤道變化引起的，黃道的變化也會引起春風(fēng)點的變化。不過，圖中的角π非常小，所以黃道變化引起和春風(fēng)點移動是非常小的。
從歲差示道意圖中的幾何關(guān)系：
    弧Nγ1 = 弧Nγm + χcos(ω)，注意，這是近似計算，但精度仍非常高
兩邊同時減去П得：
    Nγ1 - П = (Nγm - П)+ χcos(ω)，即 ψ = p+χcos(ω)
其實ψ就是：黃道不動，赤道變化引起的黃經(jīng)歲差，我們稱之是日月黃經(jīng)歲差。
其實χ就是：Date赤道上的行星歲差。
如果把ψ = p+χcos(ω)改寫為p = ψ-χcos(ω)
那么 -χcos(ω) 就是行星黃經(jīng)歲差。
同樣，從幾何關(guān)系，還可輕松得到行星黃緯歲差：χsin(ω)
----------------
χ的歲差速度是：106"/千年，加速度為-238"/平方千年
χ歲差造成春風(fēng)點東移：0".106*cos(23.5°)/年
日月歲差(黃經(jīng))造成春風(fēng)點西移：50.384.78750/年
黃經(jīng)總歲差速度約為二者之差：50".29/年
----------------
[章動]：
受到日月運動的影響，地球自轉(zhuǎn)軸發(fā)生振動，這種振動稱為章動?！罢隆弊郑瑥暮沃v起？地球自轉(zhuǎn)的振周期主要與月球軌道升交點運動有關(guān)，升交點運動周期為18.6年，古人稱之為一章。
在第(五)貼的圖中，已給出受到章動后的赤道位置。圖中的Δψ就是黃經(jīng)章動值。
章動的后果：春風(fēng)點在黃道上移動了Δψ(稱為黃經(jīng)章動)，黃赤交角變化了Δε(稱為交角章動)
考慮章動以后的赤道稱為真赤道，否則為平赤道。
黃道的振動是非常微弱，所以一般不再區(qū)分“真”與“平”
要得到真春點起算的黃經(jīng)，必須加上Δψ，要繼續(xù)旋轉(zhuǎn)到真赤道坐標(biāo)中，應(yīng)使用真黃赤交角(即加上Δε)
-----------------
Δψ與Δε的計算
前面已發(fā)了一個關(guān)于IAU1980與IAU2000B章動序列比較的貼子，里面含有相關(guān)的所有計算，請閱讀程序。
----------------
精度控制
完整的計算，計算量偏大，程序冗長。如果你只是為了計算農(nóng)歷，根本不需要那么高的精度。
為什么呢？
因為，現(xiàn)代的紫金山中國農(nóng)歷是基于北京手表時的，這種時間與地球自轉(zhuǎn)有關(guān)，考慮到地球地轉(zhuǎn)的不規(guī)律性，無法對幾百年后的北京時間進行測算，誤差0.5分鐘是正常的，如果對5千年后的北京時間預(yù)測計算，誤差可達幾小時。0.5分鐘的誤差意味差什么呢？意味著太陽位置誤差可達1".2，月亮位置誤差可達15"，因此沒有必要把Δψ算得很準(zhǔn)，0".1已經(jīng)足夠。
以下javascript程序給出幾千年內(nèi)精度為0".05的黃經(jīng)章動及交角章動。
<script language=javascript>
/*********************
對IAU1980章動計算的簡化處理
許劍偉 2008年5月2日
*********************/
t=1; //J2000.0起算的儒略世紀(jì)數(shù)
var nuB=new Array(//章動表
2.1824,   -33.75705, 36e-6,-171996,-1742,92025, 89,
3.5069, 1256.66393, 11e-6, -13187, -16, 5736,-31,
1.3375, 16799.41822,-51e-6, -2274,   -2, 977, -5,
4.3649,   -67.51409, 72e-6,   2062,    2, -895, 5,
0.0431, -628.30196, 3e-6, -1426,   34,   54, -1,
2.3556, 8328.69143,155e-6,    712,    1,   -7, 0,
3.4638, 1884.96589, 8e-6,   -517,   12, 224, -6,
5.4383, 16833.17527,-87e-6,   -386,   -4, 200, 0,
3.6931, 25128.10965,103e-6,   -301,    0, 129, -1,
3.5501,   628.36198, 13e-6,    217,   -5, -95, 3
);
var dL3=0,dE3=0;
for(i=0;i<nuB.length;i+=7){
c=nuB[i] +nuB[i+1]*t +nuB[i+2]*t*t;
dL3+=(nuB[i+3]+nuB[i+4]*t/10)*Math.sin(c); //黃經(jīng)章動
dE3+=(nuB[i+5]+nuB[i+6]*t/10)*Math.cos(c); //交角章動
}
dL3/=10000; //黃經(jīng)章動
dE3/=10000; //交角章動
document.write("[精簡的章動計算,精度0.05角秒,單位角秒]<br>");
document.write("黃經(jīng)章動：" + dL3 + " 交角章動：" + dE3 + "<br>");

</script>

---------------------
[問題]：程序中使用了“J2000.0起算的儒略世紀(jì)數(shù)”，你會計算它嗎？
[例題]：比如，2002年1月1日12:00:00 TD的儒略世紀(jì)數(shù)T是多少？TD表示力學(xué)時。
[解]：
   J2000.0對應(yīng)2000年1月1日12:00:00 TD
   由于2000年是閏年，含有366天，也就是說：
     2001年1月1日12:00:00加上366天后變?yōu)?001年1月1日12:00:00
   又因為2001年是平年，所以再加上365天就是2002年1月1日12:00:00
   所以，所求時刻相對于J2000.0相差366+365=731天
   那么，T = 731/36525 = 0.020013689

行星光行差計算起來比較麻煩，也沒有現(xiàn)成的簡單公式可以應(yīng)用，我們總是根據(jù)具體的行星問題展開計算。用下文所述的方法得到的光行差是非常準(zhǔn)確的，最后得到的行星視位置與skymap的計算結(jié)果分毫不差。當(dāng)然，日月計算也會涉及行星光行差問題，但這只以行星光行差中的最簡單的一個特例。

此主題相關(guān)圖片如下：

t時刻，地球在A點，天體在B點。此時，我們在地球上看到的星光并不是并非t時刻天體發(fā)出來的光。由于光速是有限值，星光傳到地球需要一定的時間τ，所以t時刻我們看到的星光是天體在t-τ時刻發(fā)出來的。在圖中，O和S分別表示t-τ時刻地球和天體的直位置。我們稱τ為光行時。

可以把τ看作一個時間單位，圖中矢量的既可表示速度大小和方向，也可表示位移的大小和方向。

把t時刻看到的光線（圖中的SA矢量），可分角為SO矢和SC矢。顯然SC矢與地球的速度(或位移)OA是相同的，這樣，我們站在地心看光波，感覺不到SC分量，光的速度只剩下SO分量（這好比兩個人并行走路，感覺對方是相對自已是靜止的）。SO矢量與CA矢量是相同的，考慮到光線最終落到A點，所以星體的視覺方向為AC。因此，視方向與真方向的差為圖中的角CAB。以上推理是基于伽俐略變換的，我們總結(jié)如下：

t時刻星體的視方向是由t-τ時刻星體發(fā)出的并傳到觀測者的光線決定的，在地心坐標(biāo)中，視方向由該光線的SO分量決定。換句話說，在地心坐標(biāo)中，t時刻星體的視方向是t-τ是刻星體的真位置方向。

我們定義：t時刻，行星光行差等于視方向與真方向的角度差。這個角度差就是圖中的角CAB。通常，角度差（行星光行差）用經(jīng)緯度修正量來表度。反過來說，真方向加上行星光行差就得到視位置。

我們已經(jīng)知道，行星光行差與t-τ=時刻的真方向有關(guān)，由于光行時τ是未知的，這給計算帶來一點小麻煩。不過，多數(shù)情況下，B點的地心坐標(biāo)(經(jīng)度U、緯度V、距離Δ)已經(jīng)精確算出，那么我們可以用τ=Δ/c來估算光行時，式中c為光速，Δ的單位中千米，τ的單位是秒。

接下來我們有兩種計算方法求解修正量：

（1）使用低精度算法重新計算t-τ時刻與t時刻天體的地心坐標(biāo)。然后算出這兩個坐標(biāo)的經(jīng)緯度差值即可，這個差值就是行星光行差。應(yīng)注意，用t-τ的坐標(biāo)值減t時刻的坐標(biāo)值。還應(yīng)注意，應(yīng)使用完全相同的算法算出這兩個坐標(biāo)，并且計算過程中請保留所有的數(shù)字(16位)，不要做任何舍入，以保證能夠分辨出差值。

剛才說到“低精度”算法，具體含義是什么？我們的目標(biāo)是準(zhǔn)確計算出位置的差值，實際上是差值的精度達到3至4位有效數(shù)字即可。因為光行星光行差通常為幾十角秒，如果達到3位有效數(shù)字，誤差已小于0.05角秒。

當(dāng)使用VSOP87算法時，周期項的表達為P=A*(t^n)*COS(B+C*t)，進行差值計算時，P對t的導(dǎo)數(shù)P'越大，對差值的影響就越大，每個項的影響量=P'*τ/(365250*86400)弧度。除了幾個A的值比較大的中短周期項的P’比較大，其它的都非常小，實際應(yīng)用中每個表取4項已足夠，這相當(dāng)于把行星運動看作橢圓運動并加上1兩個主要攝動項。

如果你使用VSOP87方法進行星歷計算，建議使用這種方法進行光行差修正。因為，你已經(jīng)擁有了星體位置計算的函數(shù)，調(diào)用該函數(shù)計算兩次位置即可輕松算出行星光行差。當(dāng)然，要適當(dāng)控制計算的項數(shù)，以免影響程序效率。

（2）先計算出因地球運動引起的偏差，在天文學(xué)中，這個偏差稱為“光行差”，對應(yīng)圖中的角ASO。地球運動造成光線傳播期間地球發(fā)生位移(圖中的OA)，它對S點的張角(即角ASO)就是恒星周年光行差，除了使用上述的“速度或位移的分解法”，我們也可以用速度合成的方法去理解它，恒星周年光行差有現(xiàn)成的公式可能使用。同樣，由于行星運動，也形成了一個張角SAB，恒星周年光行差減去這個張角就得到行星光行差。行星光行差加上真位置就得到視位置。

幾個注意點：

（1）雖然星體B的運動改變了光的方向(伽俐略原理)，但對于觀測者，只能接收到到來自SA路徑的光線，所以觀測者并沒有覺察到運動方向的改變。

（2）計算光行時τ的時候，我們使用τ=Δ/c計算，式中Δ=|AB|。嚴(yán)格的說，τ=|CA|/c，然而CA是未知的。不過，由于行星的運動速度遠小于光速c，所以|CA|與|AB|相差很小。因此，用τ=Δ/c計算產(chǎn)生的誤差最多為v/c，式中v是行星的地心速度，v/c只有萬分之幾，引起的誤差僅是毫角秒數(shù)量級的。另外，如果使用伽利變換，在地心坐標(biāo)中，行星發(fā)出的光線相對地球的速度最大可達|c±v|，而我們計算τ時使用c，這也將造成一定的誤差。其實，如果考慮相對論效應(yīng)，伽利略修正的誤差與Δ、c取值不嚴(yán)格所造成的誤差是同一數(shù)量級的，所以沒有必要把τ計算得很嚴(yán)格。

通過以下改正,可以大幅度提高VSOP87的精度
1、以下對VSOP87日心球面地球序列擬合到DE406:
   利用VSOP87計算出黃經(jīng)L后,計算黃經(jīng)的修正量dL如下:
      dL = 0.2 -0.22*t -0.054*t2 + 0.1*cos(t)*t - 0.11*cos(1.5*t)*t -0.13*cos(3*t)
   式中t是J2000.0起算的儒略千年數(shù)
   最后 L = L+dL
   修正前的最大誤差:
     J2000+-3000年1.00"
          +-2000年0.67"
          +-1000年0.39"
          +- 500年0.28"
          +- 250年0.17"
   修正后的最大誤差:
     J2000+-3000年0.6"
          +-2000年0.25"
          +-1000年0.12"
          +- 500年0.02"
          +- 250年0.012"
   距離與黃緯不再修正,其最大誤差為:
     距離：+-250年3E-8AU, +-1000年1E-7AU, +-3000年1E-7AU
     黃緯：+-250年0.01" , +-1000年0.024", +-3000年0.08"
2、以下對VSOP87的date球面地球序列擬合到DE406:
   dL = 0.2 +2.95*t -0.15*t2 - 0.02*t3
        - 0.01*sin(5*t) -0.13*cos(3*t)- 0.1*cos(1*t)*t
   最后L = L + dL
   修正后的最大誤差:
     J2000+-3000年0.8"
          +-2000年0.4"
          +-1000年0.1"
          +- 500年0.02"
          +- 250年0.014"
黃經(jīng)擬合時已考慮了最新的歲差速度訂正: -2.9965"/儒略千年
黃緯與距離不用修正。
3、使用J2000.0動力學(xué)分點,因此J2000黃道坐標(biāo)中的黃經(jīng)與DE405/406相差一個估定值0.04"。
   注意1，使用J2000.0分點(非常接近ICRF的赤經(jīng)起標(biāo)點)得到的星體坐標(biāo)才能與最新的《中國天文年歷》相同。
   注意2,公元3000以后的擬合了“瑞士星歷表”
4、許劍偉,地震后第12天

《天文算法》一書中給出了求各種天文現(xiàn)象發(fā)生時刻的方法，不過作者并沒有深入討論算法的效率、簡潔性及具體原理。事實書中算法的精度是固定的，不一定滿足現(xiàn)代天文計算的要求。
比如試圖把農(nóng)歷的精度計算到30秒,使用《天文算法》一書中的算法將不能實現(xiàn)。本貼討論高精度的、超高速的算法。
問題假設(shè)：設(shè)經(jīng)度L(t) = Lo + v*t + f(t)，對于方程L(t)=W，求t的值。
式中速度v值較大,f(t)對t求導(dǎo)數(shù)也得到速度量dv,如果它相對于v很小,那么:
解法(針對VSOP87或ELP2000-82或MPP02)：
(1)第一步,初步確定t的值：
t1 = (W-Lo)/v
這樣，L的偏差為 dL = L(t1) - W
那么t1的偏差估計為dt = dL/v
對t1修正 t2 = t1 + dt
(2)第二步,用第一步的方法重新執(zhí)行一次,不過計算過程中可以考慮v使用v+dv代替。
(3)第三步,重復(fù)第二步計算,但v一定要使用v+dv代替
注意事項1:計算L(t)時不一定需要完全精度,第一、二步計算時只需計算幾個周期項即可。第三步則需多算幾個項。如果精度還不滿意，則重復(fù)第三步
注意事項2:v的解析式可以使用VSOP87或ELP/MPP02的序列求導(dǎo)數(shù)得到,只需考慮速度量較大的幾個項即可。
注意事項3:如果使用DE405/DE406星歷表或“瑞士星歷表”,可以考慮使用求變差的辦法求v+dv。“變差”法，指求相近兩個時間的坐標(biāo)差值，除以時間差后得到速度。
說明1：以上方法本質(zhì)上就是《天文算法》中的迭代逼近思想方法。
說明2：上面講到的v的運動學(xué)本質(zhì)是“平角速度”，用2*3.1415926535/v得到“平周期”，如地球運動的“平周期”是365.2422天
說明3：對v+dv的取值還需與截斷誤差合并分析，這當(dāng)中涉及“加速度”的分析，將在以后的貼子中討論。
說明4：以上所述的是通用方法，對于實際問題，有必要進行更仔細(xì)的分析。

計算量：一次精確計算位置坐標(biāo)的計算量為B，那么該算法的計算量小于1.3*B

算法的極限精度取決于VSOP或ELP/MPP02理論本身
《農(nóng)歷天文算法完全解——第一部分(分兩部分講述)》
算法本質(zhì)：迭代技術(shù)
一、日月黃經(jīng)速度的表達
(1)地球黃經(jīng)速度的表達(Date黃道坐標(biāo)中的速度,含歲差速度)
以下表達式中t是千年數(shù),誤差小于萬分之3,數(shù)據(jù)由VSOP87地球黃經(jīng)序列求導(dǎo)數(shù)得到:
v = 6283.32 //誤差小于萬分3
+210 *Math.sin(1.527+f)
+4.4 *Math.sin(1.48+f*2)
+12.9*Math.sin(5.82+f)*t
+0.55*Math.sin(4.21+f)*t*t;
式中f=6283.07585*t;
(2)月球黃經(jīng)速度的表達(Date黃道坐標(biāo)中的速度,不含歲差速度)
歲差速度遠小于月球黃經(jīng)速度,所以下速度計算忽略了歲差速度。
平速度是(單位：單位弧度/世紀(jì),式中的t是世經(jīng)數(shù))：
v=8400; //精確值是8399.68473007193
周期項速度(對ELP/MPP02黃經(jīng)序列求導(dǎo)數(shù)得到)
dv=914*sin( 0.785 +8328.69142*t); //誤差5%
dv=914*sin( 0.785 +8328.69142*t) //誤差不超過0.8%
+179*sin( 2.54 +15542.7543*t )
+160*sin( 0.19 +7214.0629*t )
+ 62*sin( 3.14 +16657.383 *t )
+ 34*sin( 4.8 +16866.932 *t )
+ 22*sin( 4.9 +23871.45 *t )
+ 12*sin( 2.6 +14914.45 *t );
dv=914*sin( 0.7848 +8328.691425*t+ 1.523*t2 ) //誤差小于0.3%
+179*sin( 2.543 +15542.7543*t )
+160*sin( 0.1874 + 7214.0629*t )
+62 *sin( 3.14   +16657.3828*t )
+34 *sin( 4.827 +16866.9323*t )
+22 *sin( 4.9    +23871.4457*t )
+12 *sin( 2.59   +14914.4523*t )
+7 *sin(0.23    + 6585.7609*t )
+5 *sin(0.9     +25195.624 *t )
+5 *sin(2.32    - 7700.3895*t )
+5 *sin(3.88    + 8956.9934*t )
+5 *sin(0.49    + 7771.3771*t )
+4 *sin(5.5     +24986.074 *t )
+4 *sin(4.3     +22756.817 *t )
+2 *sin(1       +32200.137 *t )
+2 *sin(1.9     +14428.126 *t )
+2 *sin(0.4     +31085.509 *t )
+2 *sin(1.528   + 628.302 *t );
dv的值可根據(jù)精度的需要從上面表達式中取一種。
月球真速度為：
v=v-dv;
(3)月球平黃經(jīng)與黃經(jīng)的表達式(ELP/MPP02)
w1 = 3.81034409083
    +8399.68473007193*t
    -3.31895204255e-05*t2
    +3.11024944911e-08*t3
    -2.03282376489e-10*t4
L=w1+周期項+泊松項，由于泊松項較小僅需考慮幾項就可以了，除非有高精度計算的需要。
周期項與泊松項的具體算法待敘。
二、求月球經(jīng)度為W時所對應(yīng)的時間
經(jīng)度計算函數(shù)為 Lon(t,n); t為J2000.0起算的世紀(jì)數(shù),t不超過+-30世紀(jì);n為序列截取數(shù)量，計算步驟如下:
數(shù)據(jù)準(zhǔn)備,表達式準(zhǔn)備:
to=0;
v0 = 8400; //平速度
v1(t) = 8400 + 914*sin(0.785+8328.69142*t); //低精度表達,5%精度
v2(t) = 中精度速度表達(見上文),0.3%精度
v = v0
(1)a.余項計算：dW = W - L(to,0)
   b.時間計算：to = to+ dW/v
   c.速度計算：v=v1(to); //為下一次計算準(zhǔn)備速度,速度誤差5%左右
   以上兩行相當(dāng)于解方程 W=3.81034 +8400*t,解為to = (W - 3.81034)/8400
   to是月球黃經(jīng)為W是的“平”時間
   to與真時間的差值不超過18小時
   以上計算已經(jīng)考慮了黃經(jīng)的常數(shù)項(3.81034)和速度項(8400)。因為to的超不超過30，所以平黃經(jīng)的高次項(加速度項)很小，不超過2度，黃經(jīng)周期項不超過8度，所以黃經(jīng)的的高次項及周期項總和不超過10度。因此如果把to再次回代到完整的序列中計算，上述方程最多只有10度(約0.17弧度)的誤差：
W=3.81034 +8400*t+0.17，得to = (W -3.81034 -0.17)/8400，與原來的to比較不難發(fā)現(xiàn)，最大誤差就是月球運動0.17(即10度)所用的時間,約18小時。這18個小時的誤差主要是截斷引起的。
(2)a.余項計算：dW = W - Lon(to,3);
   b.時間計算：to=to+dW/v;
   c.速度計算：v=v2(to); //為下一次計算準(zhǔn)備速度,精度為0.3%，如果考慮由于時間不準(zhǔn)確引起的額外誤差，精度仍可達到0.35%(取0.4%)
   在步驟(1)中已經(jīng)說到,誤差為18小時,我們稱之為余項(或稱之為截斷誤差),而本次計算已經(jīng)截取了3個主要周期項，可修正這18小時的余項。由于速度的精度只有5%，所以殘留誤差為18*5%=0.9小時。不過，實際誤差不止這些，因為本次計算也存在截斷誤差。
   截斷部分誤差估計: 把被截斷的最大10項的振幅取和作為最大誤差估計,以后的項不用算(因為存在正負(fù)相消的作用),這10項取和為3000角秒=50角分，對應(yīng)的最大時間誤差為1.7小時。與剛才討到的0.9小時綜合考慮，不妨把誤差估計為2小時。
   為什么不把誤差估計為0.9+1.7=2.6小時？因為0.9與1.7的誤差估計值均指最大誤差，出現(xiàn)最大誤差的可能性已經(jīng)很小了，二者同時出現(xiàn)最大值的可能性幾乎為0。
   當(dāng)然，平均誤差僅為最大誤差的六分之一左右。
   第(1)步與第(2)步計算無需考慮光行時、章動的修正計算，因為誤差遠超過這些修正量。
(3)a.余項計算：dW = W - Lon(to,21);
   b.時間計算：to = to + dW/v;
   c.無需再計算更精確的速度v
   周期項取21位的理論精度是4',對應(yīng)的時間誤差為8分鐘。
   速度不精確引起的誤差：2小時*0.004=0.008小時=0.48分鐘。
   總誤差不妨按8分鐘估計。
(4)a.余項計算：dW = W - Lon(to,200);
   b.時間計算：to = to+dW/v;
計算200項的理論最大誤差約2角秒,對應(yīng)的時間誤差為4秒鐘
速度v不精確引起的誤差：8分*0.004=0.3秒
總誤差不妨取4秒鐘。
to就是我們最終所需要計算的時間！
(5)如果相要再提高精度,繼續(xù)計算：
   dW = W - Lon(to,全部);
   to = to+dW/v;
“全部”計算的精度為5毫角秒，相當(dāng)于0.01秒
速度v不精確引起的誤差：5秒*0.004=0.02秒
總誤差不妨取0.03秒鐘。
ELP/MPP02在J2000前后幾十年內(nèi)與DE406擬合得很好，差異為0.005角秒以內(nèi)。
(6)說明：由于都是在鄰近的時間內(nèi)計算，可以認(rèn)為速度不變(或者說加速度為0)。因為對ELP/MPP02求二次導(dǎo)數(shù)，我們發(fā)現(xiàn)加速度非常小。為什么會這樣呢？其實，周期項的周期最小值為10天左右(個別為5天左右,但振幅非常小)，也就是說在10天范圍內(nèi)速度大小才小幅波動一次，在幾小時或幾分鐘內(nèi)，速度幾乎不變。速度大小波動幅度最大的周期項的周期是27.5天左右，它對速度的影響可達10%以上，即便是這種周期項，10小時對速度的影響最多只有：10%*10小時/27.5天 = 0.15%，2小時誤差對速度精度的影響為：0.03%
(7)說明2：如果要計算視黃經(jīng)的,應(yīng)注意光行差、光行時、歲差、章動修正等。月球光行時約為1.2秒(約修正0.6角秒)。
(8)計算量：以上針對月球的精確到5秒鐘(指最大可能的誤差為5秒,平均誤差不到1秒，標(biāo)準(zhǔn)差約為3秒)的總計算量約為250個周期項的計算量。
三、地球經(jīng)度為W時對應(yīng)的時刻的算法與月球的類似,但迭代次數(shù)要少得多。如果考慮了光行差，那么該算法可直接用于計算24節(jié)氣。
四、日月距角為W時對應(yīng)的時刻。方法同上。不同的是把上面的經(jīng)度換為月球黃與太陽黃經(jīng)的差值，速度換為月球黃經(jīng)速度與太陽黃經(jīng)速度的差。
距角與月相密切相關(guān)，距角為0度對應(yīng)新月，為90度是半滿上弦月，180度是滿月。由于章動效果同時影響太陽和月球坐標(biāo)，所以不用計算章動。當(dāng)然地球運動引起的恒星光行差是必須計算的，否則無法正確計算出月相，因為幾何意上的日月合朔不是光學(xué)(視覺)意義上的日月合朔。計算距角時，地球黃經(jīng)的精度控制在1角秒已足夠，因為月球的精度也只控制在2角秒。
五、這里描述的算法與本人在javascript農(nóng)歷程序中提供的相應(yīng)算法有所不同：速度提高了4至5倍，精度提高了5倍。
六、精度是根據(jù)程序員的意愿決定的，但受到VSOP87或DE405/406理論的極限精度制約。
七、如果讀者還時一知半解，無法寫出相應(yīng)的程序，本考慮編寫具體的程序，以做示范。
八、ELP/MPP02的黃經(jīng)起算點與DE405/406星歷表有關(guān)應(yīng)變換到J2000黃道坐標(biāo)中，以免產(chǎn)生0.1秒的時間誤差。
九、VSOP87是早期的理論，存在起點誤差以及微小的長期項誤差，同樣需做訂正，訂正后，幾百年內(nèi)的誤差小于0.5秒。訂正方法在前面的貼子中已經(jīng)講述。
十、使用真正的DE405/DE406方法的計算時間，方法同上，只是在計算Lon()時使用DE405/406方法
改進《農(nóng)歷天文算法完全解——第一部分》
算法本質(zhì)：迭代技術(shù)
一、日月黃經(jīng)速度的表達
(1)地球黃經(jīng)速度的表達(Date黃道坐標(biāo)中的速度,含歲差速度)
以下表達式中t是千年數(shù),誤差小于萬分之3,數(shù)據(jù)由VSOP87地球黃經(jīng)序列求導(dǎo)數(shù)得到:
v = 6283.32
    +210 *sin(1.5274 + 6283.07585*t) //誤差小于萬分3
    +4.4 *sin(1.4845 + 12566.15170*t)
    +12.9*sin( 2.6782 + 6283.07585*t)*t
    +0.55*sin( 1.0721 + 6283.07585*t)*t2;
(2)月球黃經(jīng)速度的表達(Date黃道坐標(biāo)中的速度,不含歲差速度)
歲差速度遠小于月球黃經(jīng)速度,所以下速度計算忽略了歲差速度。
平速度是(單位：單位弧度/世紀(jì),式中的t是世經(jīng)數(shù))：v=8399.68473007193
考慮歲差速度后,應(yīng)把式中的8399.68473007193換為：8399.68473007193弧度/世紀(jì)+5028.79角秒/世紀(jì)=8399.70911033384弧度/世紀(jì)
周期項速度(對ELP/MPP02黃經(jīng)序列求導(dǎo)數(shù)得到)
dv=914*sin( 0.785 +8328.69142*t); //誤差5%
dv=914*sin( 0.785 +8328.69142*t) //誤差不超過0.8%
+179*sin( 2.54 +15542.7543*t )
+160*sin( 0.19 +7214.0629*t )
+ 62*sin( 3.14 +16657.383 *t )
+ 34*sin( 4.8 +16866.932 *t )
+ 22*sin( 4.9 +23871.45 *t )
+ 12*sin( 2.6 +14914.45 *t );
dv=914*sin( 0.7848 +8328.691425*t+ 1.523*t2 /10000) //誤差小于0.3%
+179*sin( 2.543 +15542.7543*t )
+160*sin( 0.1874 + 7214.0629*t )
+62 *sin( 3.14   +16657.3828*t )
+34 *sin( 4.827 +16866.9323*t )
+22 *sin( 4.9    +23871.4457*t )
+12 *sin( 2.59   +14914.4523*t )
+7 *sin(0.23    + 6585.7609*t )
+5 *sin(0.9     +25195.624 *t )
+5 *sin(2.32    - 7700.3895*t )
+5 *sin(3.88    + 8956.9934*t )
+5 *sin(0.49    + 7771.3771*t )
+4 *sin(5.5     +24986.074 *t )
+4 *sin(4.3     +22756.817 *t )
+2 *sin(1       +32200.137 *t )
+2 *sin(1.9     +14428.126 *t )
+2 *sin(0.4     +31085.509 *t )
+2 *sin(1.528   + 628.302 *t );
dv的值可根據(jù)精度的需要從上面表達式中取一種。
月球真速度為：v=v-dv;
(3)月球平黃經(jīng)與黃經(jīng)的表達式(ELP/MPP02)
w1 = 3.81034409083
    +8399.68473007193*t
    -3.31895204255e-05*t2
    +3.11024944911e-08*t3
    -2.03282376489e-10*t4
注意:w1中不含歲差速度
(4)黃經(jīng)表達式(相對于Date平分點): L=w1+周期項+泊松項，由于泊松項較小僅需考慮幾項就可以了，除非有高精度計算的需要。
(5)周期項與泊松項的具體算法待敘。

二、求月球經(jīng)度為W時所對應(yīng)的時間
經(jīng)度計算函數(shù)為 Lon(t,n); t為J2000.0起算的世紀(jì)數(shù),t不超過+-30世紀(jì);n為序列截取數(shù)量，計算步驟如下:
數(shù)據(jù)準(zhǔn)備,表達式準(zhǔn)備:
to=0;
v0 = 8399.70911033384; //平速度(含歲差速度)
v1(t) = 速度表達,0.3%精度(見上文)
(1)a.余項計算：dW = W - L(to,0)
   b.時間計算：to = to+ dW/v0
   c.速度計算：v=v1(to); //為下一次計算準(zhǔn)備速度,速度誤差0.3%左右
   以上兩行相當(dāng)于解方程 W=3.81034 +8399.70911033384*t,解為to = (W - 3.81034)/8399.70911033384
   ▲to是月球黃經(jīng)為W是的“平”時間，to與真時間的差值不超過18小時
   ▲以上計算已經(jīng)考慮了黃經(jīng)的常數(shù)項(3.81034)和速度項(8399.70911033384)。因為to不超過30世紀(jì)，所以平黃經(jīng)的高次項(加速度項)很小，不超過2度。而黃經(jīng)周期項不超過8度，所以黃經(jīng)的的高次項及周期項總和不超過10度。因此如果把to再次回代到完整的序列中計算，上述方程最多只有10度(約0.17弧度)的誤差：
       W=3.81034 +8400*t+0.17，得to = (W -3.81034 -0.17)/8399.70911033384
   與原來的to比較不難發(fā)現(xiàn)，最大誤差就是月球運動0.17(即10度)所用的時間,約18小時。這18個小時的誤差就是截斷高次項及周期項引起的。
   ▲其實最后得到的v的最大誤差大于0.3%，下文還會討論到,由于to誤差18小時,也將造成v產(chǎn)生0.3%的誤差。v的總誤差估計為0.6%
   ▲請務(wù)必注意，“平速度”的精度要求比“速度”的精度要求要高得多。原因如下：
   t=(W-3.81034)/(v0+dv)=to+dt，式中dv是平速度的誤差；利用二項式定理展開后得：
   t=[(W-3.81034)/vo] * (1-*dv/Vo) = to*(1-dv/vo)
   當(dāng)W很大時(即to很大時)：dt = t - to = to*dv/vo，式中dt就是dv引起的時間誤差
   設(shè) to=20世紀(jì)，dv=1弧度/世紀(jì)，我們得到：dt=0.0024世紀(jì)=2087小時！而v的相對誤差只有dv/v=1/8400=0.012%
   所以“平速度”的精度要求非常高。
   ▲不含歲差的平速度v=8399.68473007193，當(dāng)然這指的是角速度，那么月亮平周期就是：
     T=2*3.1415926535897933/v=0.00074802632587923106世紀(jì)=27.32166155273891天(1世經(jīng)=36525天)，這就是恒星月
    還應(yīng)注意，月球黃經(jīng)還有“加速度”項
(2)a.余項計算：dW = W - Lon(to,21);
   b.時間計算：to=to+dW/v;
   在步驟(1)中已經(jīng)說到,誤差為18小時,我們稱之為余項(或稱之為截斷誤差),而本次計算已經(jīng)截取了21個主要周期項，可修正這18小時的余項。由于速度的精度只有0.6%，所以殘留誤差為18*0.6%小時=7分鐘。
   不過，實際誤差不止這些，因為本次計算也存在截斷誤差。
   ▲截斷部分誤差估計: 把被截斷的最大10項的振幅取和作為最大誤差估計,以后的項不用算(因為存在正負(fù)相消的作用),這10項取和為130角秒，對應(yīng)的最大時間誤差為5分鐘。與剛才討到的7分鐘綜合考慮，不妨把誤差估計為12*0.7=10分鐘。
   為什么不把誤差估計為7+5=12分鐘？因為7與5的均指最大誤差，出現(xiàn)最大誤差的可能性已經(jīng)很小了，二者同時出現(xiàn)最大值的可能性幾乎為0。
   當(dāng)然，平均誤差僅為最大誤差的六分之一左右。
   ▲第(1)步與第(2)步計算無需考慮月球光行時、地球章動的修正計算，因為誤差遠超過這些修正量。
(3)a.余項計算：dW = W - Lon(to,159);
   b.時間計算：to = to + dW/v;
   周期項取159項的理論精度是3",對應(yīng)的時間誤差為6秒。
   速度不精確引起的誤差：10分鐘*0.6%=3.6秒，平均誤差按1/6計算為0.6秒，不過這個估值比實測值要大許多。
   總誤差不妨按 6+0.6=7秒鐘估計。
(4)如果相要再提高精度,繼續(xù)計算：
   dW = W - Lon(to,全部);
   to = to+dW/v;
“全部”計算的精度為5毫角秒，相當(dāng)于0.01秒
速度v不精確引起的誤差：5秒*0.006=0.03秒
總誤差不妨取0.03秒鐘。
ELP/MPP02在J2000前后幾十年內(nèi)與DE406擬合得很好，差異為0.005角秒以內(nèi)。
(5)說明：由于都是在鄰近的時間內(nèi)計算，可以認(rèn)為速度不變(或者說加速度為0)。因為對ELP/MPP02求二次導(dǎo)數(shù)，我們發(fā)現(xiàn)加速度非常小。為什么會這樣呢？其實，周期項的周期最小值為10天左右(個別為5天左右,但振幅非常小)，也就是說在10天范圍內(nèi)速度大小才小幅波動一次，在幾小時或幾分鐘內(nèi)，速度幾乎不變。速度大小波動幅度最大的周期項的周期是27.5天左右，它對速度的影響可達10%以上，即便是這種周期項，18小時對速度的影響最多只有：10%*18小時/27.5天 = 0.3%，2小時誤差對速度精度的影響為：0.03%
(6)說明2：如果要計算視黃經(jīng)的,應(yīng)注意光行差、光行時、歲差、章動修正等。月球光行時約為1.2秒(約修正0.6角秒)。
(7)計算量：以上針對月球的精確到7秒鐘(指最大可能的誤差為7秒,平均誤差約1秒，標(biāo)準(zhǔn)差約為3秒)的總計算量約為200個周期項的計算量。
三、地球經(jīng)度為W時對應(yīng)的時刻的算法與月球的類似,但迭代次數(shù)要少得多。如果考慮了光行差，那么該算法可直接用于計算24節(jié)氣。
四、日月距角為W時對應(yīng)的時刻。方法同上。不同的是把上面的經(jīng)度換為月球黃與太陽黃經(jīng)的差值，速度換為月球黃經(jīng)速度與太陽黃經(jīng)速度的差。
距角與月相密切相關(guān)，距角為0度對應(yīng)新月，為90度是半滿上弦月，180度是滿月。由于章動效果同時影響太陽和月球坐標(biāo)，所以不用計算章動。當(dāng)然地球運動引起的恒星光行差是必須計算的，否則無法正確計算出月相，因為幾何意上的日月合朔不是光學(xué)(視覺)意義上的日月合朔。計算距角時，地球黃經(jīng)的精度控制在1角秒已足夠，因為月球的精度也只控制在2角秒。
五、這里描述的算法與本人在javascript農(nóng)歷程序中提供的相應(yīng)算法有所不同：速度提高了4至5倍，精度提高了5倍。
六、精度是根據(jù)程序員的意愿決定的，但受到VSOP87或DE405/406理論的極限精度制約。
七、如果讀者還時一知半解，無法寫出相應(yīng)的程序，本考慮編寫具體的程序，以做示范。
八、ELP/MPP02的黃經(jīng)起算點與DE405/406星歷表有關(guān)應(yīng)變換到J2000黃道坐標(biāo)中，以免產(chǎn)生0.1秒的時間誤差。
九、VSOP87是早期的理論，存在起點誤差以及微小的長期項誤差，同樣需做訂正，訂正后，幾百年內(nèi)的誤差小于0.5秒。訂正方法在前面的貼子中已經(jīng)講述。

這是一個超高速的算法，基于ELP/MPP02，最大理論誤差7秒，理論平均誤差1秒，程序大小8.5k。

請修改程序中的testT，以便計算其它時刻的情形。

本程序未經(jīng)檢驗，可能有誤，有空的時候我會與《中國天文年歷》及“瑞士星歷表”比對一下。到時再做報告。

xjw01于十中，2008年5月27日晚編寫

M_L1=new Array(
1.677,4.669,628.30196,-0.03,0,0,
0.516,3.372,6585.7609,-2.16,-19,0.1,
0.414,5.728,14914.4523,-0.64,6,0,
0.371,3.969,7700.3895,1.55,25,0,
0.276,0.74,8956.9934,1.5,25,0,
0.246,4.23,-2.3012,1.5,25,0,
0.071,0.14,7842.365,-2.2,-19,0,
0.061,2.50,16171.056,-0.7,6,0,
0.045,0.44,8399.679,-0.4,0,0,
0.040,5.77,14286.150,-0.6,6,0,
0.037,4.63,1256.604,-0.1,0,0,
0.037,3.42,5957.459,-2.1,-19,0,
0.036,1.80,23243.144,0.9,31,0,
0.024,0,16029.081,3,50,0,
0.022,1.0,-1742.931,-4,-44,0,
0.019,3.1,17285.685,3,50,0,
0.017,1.3,0.329,2,25,0,
0.014,0.3,8326.390,3,50,0,
0.013,4.0,7072.088,2,25,0,
0.013,4.4,8330.993,0,0,0,
0.013,0.1,8470.667,-2,-19,0,
0.011,1.2,22128.515,-3,0,0,
0.011,2.5,15542.754,-1,0,0,
0.008,0.2,7214.06,-2,0,0,
0.007,4.9,24499.75,1,0,0,
0.007,5.1,13799.82,-4,0,0,
0.006,0.9,-486.33,-4,0,0,
0.006,0.7,9585.30,1,0,0,
0.006,4.1,8328.34,2,0,0,
0.006,0.6,8329.04,2,0,0,
0.005,2.5,-1952.48,1,0,0,
0.005,2.9,-0.71,0,0,0,
0.005,3.6,30457.21,-1,0,0);

M_L2=new Array(
0.0049,4.67,628.3020,0,0,0,
0.0023,2.67,-2.301,1.5,25,0,
0.0015,3.37,6585.761,-2.2,-19,0,
0.0012,5.73,14914.452,-0.6,6,0,
0.0011,3.97,7700.389,2,25,0,
0.0008,0.7,8956.993,1,25,0);

function Mnn(F,t,t2,t3,t4,n){ //計算M_L0或M_L1或M_L2
n = Math.floor(n * F.length / M_L0.length+0.5)*6; //項數(shù)是以周期項為基準(zhǔn),n是周期項的計算項數(shù)
if(n>F.length) n=F.length;
var i,v=0; t2/=1e4,t3/=1e8,t4/=1e8;
for(i=0;i<n;i+=6)
    v+=F[i]*Math.cos(F[i+1] +t*F[i+2] +t2*F[i+3] +t3*F[i+4] +t4*F[i+5]);
return v;
}
function moonLon(t,n){ //月球經(jīng)度計算,返回Date分點黃經(jīng),傳入世紀(jì)數(shù)
if(n==-1) n = Math.floor(M_L0.length/6+0.1);
var t2=t*t, t3=t2*t, t4=t3*t, t5=t4*t;
var v = Mnn( M_L0, t,t2,t3,t4, n );
      v += Mnn( M_L1, t,t2,t3,t4, n )*t;
      v += Mnn( M_L2, t,t2,t3,t4, n )*t2;
var w = 3.81034409 + 8399.684730072*t -3.319e-05*t2 + 3.11e-08*t3 - 2.033e-10*t4; //月球平黃經(jīng)
var p =5028.792262*t + 1.1124406*t2 + 0.00007699*t3 - 0.000023479*t4 -0.0000000178*t5; //歲差
return w + (v+p)/rad;
}

function moonV(t,jing){ //月球速度計算,傳入世經(jīng)數(shù)
var v=8399.70911033384;
if(jing==0) return v;
v-=914*Math.sin( 0.7848 +8328.691425*t+ 1.523*t*t/10000 );
if(jing==1) return v;                    //誤差小于5%
v-=179*Math.sin( 2.543 +15542.7543*t ) //誤差小于0.3%
    +160*Math.sin( 0.1874 + 7214.0629*t )
    +62 *Math.sin( 3.14   +16657.3828*t )
    +34 *Math.sin( 4.827 +16866.9323*t )
    +22 *Math.sin( 4.9    +23871.4457*t )
    +12 *Math.sin( 2.59   +14914.4523*t )
    +7 *Math.sin(0.23    + 6585.7609*t )
    +5 *Math.sin(0.9     +25195.624 *t )
    +5 *Math.sin(2.32    - 7700.3895*t )
    +5 *Math.sin(3.88    + 8956.9934*t )
    +5 *Math.sin(0.49    + 7771.3771*t )
    +4 *Math.sin(5.5     +24986.074 *t )
    +4 *Math.sin(4.3     +22756.817 *t )
    +2 *Math.sin(1       +32200.137 *t )
    +2 *Math.sin(1.9     +14428.126 *t )
    +2 *Math.sin(0.4     +31085.509 *t )
    +2 *Math.sin(1.528   + 628.302 *t );
   return v;
}

function moon_t(W){ //已知黃經(jīng)求時間
var t,dW,v= 8399.70911033384;
t = ( W - 3.81034 )/v; v=moonV(t,2); //v的精度0.6%，詳見原文
t += ( W - moonLon(t,20) )/v;
t += ( W - moonLon(t,159))/v;
return t;
}

testT=30; //世紀(jì)數(shù)
L=moonLon(testT,-1); //正算
t=moon_t(L); //反算
dt=(testT-t)*35625*86400;

document.write("高速迭代法求指定Date平分點黃經(jīng)的發(fā)生時刻。測試如下：<br>");
document.write("正算時間T=" + testT + "(世紀(jì)) 時的黃經(jīng)是:" + L + "(弧度)<br>");
document.write("反算黃經(jīng)L=" + L + "(弧度) 時的時間是:" + t + "(世紀(jì))<br>");
document.write("迭代誤差(不是實際誤差):" +dt +"秒<br>");

</script>

為了能夠與《中國天文年歷》比對，我們須把以上程序中的月球坐標(biāo)轉(zhuǎn)換為視坐標(biāo)

這時我們應(yīng)處理光行差及章動

月球光行差問題
我們已討論過精密的光行差計算，但是，在農(nóng)歷計算中根本無需精密計算，以下考慮簡化計算。
從ELP的中序列截取5%精度的距離及速速表達式：
    日月距離：r=385001 +20905*sin( 0.7848 +8328.691425*t+ 1.523*t*t/10000 ); //最大誤差2%
    黃經(jīng)速度：v=8400   - 914*sin( 0.7848 +8328.691425*t+ 1.523*t*t/10000 ); //最大誤差5%
為了求證方便，令：
    r0=385001，v0=8400
    φ= 0.7848 +8328.691425*t+ 1.523*t*t/10000
    dr= 20905*sin(φ)
    dv=   914*sin(φ)
    光速：c=300000千米/秒=9.467E14千米/世紀(jì)
那么有：r = r0+dr；v = v0+dv；τ= r/c； (式中τ光行時)
光行期間黃經(jīng)的變化量：
    dL=v*τ=v*r/c，代入得：
    dL = (v0+dv)*(r0+dr)/c = (v0*r0/c)*(1+dv/v+dr/r) (略去二階小量)
    dL = (3.416E-6弧度)*(1+(20905/385001-914/8400)*sin(φ))
       = (0.7角秒)*(1-0.055*sin(φ))
結(jié)論：月球的黃經(jīng)光行差約為0.7角秒(最大誤差16%)，如果考慮上式中的周期，最大誤差5%

以下關(guān)于章動，取自IAU1980，最大誤差不超過0.1角秒
function nutation(t){ //章動計算,t為世紀(jì)數(shù),返回值為角秒單位
var c, dL=0,dE=0, t2=t*t, y0=-17.20-0.01742*t;
c= 2.1824    -33.75705*t +36e-6*t2; dL = y0 *Math.sin(c); dE = 9.20*Math.cos(c);
c= 3.5069 +1256.66393*t +11e-6*t2; dL+= -1.32*Math.sin(c); dE+= 0.57*Math.cos(c);
c= 1.3375 +16799.4182*t -51e-6*t2; dL+= -0.23*Math.sin(c); dE+= 0.10*Math.cos(c);
c= 4.3649    -67.5141*t +72e-6*t2; dL+= 0.21*Math.sin(c); dE+=-0.09*Math.cos(c);
c= 0.04   -628.302*t; dL+= -0.14*Math.sin(c);
c= 2.36 +8328.691*t; dL+= 0.07*Math.sin(c);
c= 3.46 +1884.966*t; dL+= -0.05*Math.sin(c); dE+= 0.02*Math.cos(c);
c= 5.44 +16833.175*t; dL+= -0.04*Math.sin(c); dE+= 0.02*Math.cos(c);
c= 3.69 +25128.110*t; dL+= -0.03*Math.sin(c);
c= 3.55   +628.362*t; dL+= 0.02*Math.sin(c);
this.dL=dL; //黃經(jīng)章動
this.dE=dE; //交角章動
}

調(diào)用方法：nu = new nutation(testT);

那么：nu.dL為黃經(jīng)章動，nu.dE為交角章動

我們修改第八貼的測試程序的末尾部分，以便能夠與《天文年歷》比對

rad = 180*3600/Math.PI;
function rad2mrad(v){   //對超過0-2PI的角度轉(zhuǎn)為0-2PI
v=v % (2*Math.PI);
if(v<0) return v+2*Math.PI;
return v;
}
function rad2str(d,tim){ //將弧度轉(zhuǎn)為字串
//tim=0輸出格式示例: -23°59' 48.23"
//tim=1輸出格式示例: 18h 29m 44.52s
var s="+";
var w1="°",w2="’",w3="”";
if(d<0) d=-d,s='-';
if(tim){ d*=12/Math.PI; w1="h ",w2="m ",w3="s "; }
else     d*=180/Math.PI;
var a=Math.floor(d); d=(d-a)*60;
var b=Math.floor(d); d=(d-b)*60;
var c=Math.floor(d); d=(d-c)*100;
d=Math.floor(d+0.5);
if(d>=100) d-=100, c++;
if(c>=60) c-=60, b++;
if(b>=60) b-=60, a++;
a="   "+a, b="0"+b, c="0"+c, d="0"+d;
s+=a.substr(a.length-3,3)+w1;
s+=b.substr(b.length-2,2)+w2;
s+=c.substr(c.length-2,2)+".";
s+=d.substr(d.length-2,2)+w3;
return s;
}

testT=0.08+(1-1.5)/36525; //世紀(jì)數(shù)
L=moonLon(testT,-1); //正算
t=moon_t(L); //反算
dt=(testT-t)*35625*86400;

document.write("高速迭代法求指定Date平分點黃經(jīng)的發(fā)生時刻。測試如下：<br>");
document.write("正算時間T=" + testT + "(世紀(jì)) 時的黃經(jīng)是:" + L + "(度)<br>");
document.write("反算黃經(jīng)L=" + L + "(度) 時的時間是:" + t + "(世紀(jì))<br>");
document.write("迭代誤差(不是實際誤差):" +dt +"秒<br>");

//以下是視位置計算
function nutation(t){ //章動計算,t為世紀(jì)數(shù),返回值為角秒單位
var c, dL=0,dE=0, t2=t*t, y0=-17.20-0.01742*t;
c= 2.1824    -33.75705*t +36e-6*t2; dL = y0 *Math.sin(c); dE = 9.20*Math.cos(c);
c= 3.5069 +1256.66393*t +11e-6*t2; dL+= -1.32*Math.sin(c); dE+= 0.57*Math.cos(c);
c= 1.3375 +16799.4182*t -51e-6*t2; dL+= -0.23*Math.sin(c); dE+= 0.10*Math.cos(c);
c= 4.3649    -67.5141*t +72e-6*t2; dL+= 0.21*Math.sin(c); dE+=-0.09*Math.cos(c);
c= 0.04   -628.302*t; dL+= -0.14*Math.sin(c);
c= 2.36 +8328.691*t; dL+= 0.07*Math.sin(c);
c= 3.46 +1884.966*t; dL+= -0.05*Math.sin(c); dE+= 0.02*Math.cos(c);
c= 5.44 +16833.175*t; dL+= -0.04*Math.sin(c); dE+= 0.02*Math.cos(c);
c= 3.69 +25128.110*t; dL+= -0.03*Math.sin(c);
c= 3.55   +628.362*t; dL+= 0.02*Math.sin(c);
this.dL=dL; //黃經(jīng)章動
this.dE=dE; //交角章動
}

nu=new nutation(testT);
L2=rad2mrad(L+nu.dL/rad-0.7/rad); //補章動與光行差
Ls=rad2str(L2,0);
document.write("視黃經(jīng)：" + Ls+"<br>");

測試結(jié)果如下：

隨機插出幾個月球視黃經(jīng)數(shù)據(jù)進行比較

2008年01月01日0h TD
+197°19' 24.43" 中國天文年歷
+197°19’24.91" 本程序

2008年01月06日0h TD
+256°54' 36.32" 中國天文年歷
+256°54' 36.11" 本程序

2008年01月18日0h TD
+ 56°04' 29.83" 中國天文年歷
+ 56°04' 29.68" 本程序

2100年01月01日0h TD
+157°24' 01.183" 瑞士星歷表
+157°24' 01.96" 本程序

2100年01月18日0h TD
+ 22°14' 39.400" 瑞士星歷表
+ 22°14' 40.47" 本程序

2200年01月02日0h TD
+108°26' 45.916" 瑞士星歷表
+108°26’46.12" 本程序

我們不難發(fā)現(xiàn)，誤差與精密星歷相差無幾，黃經(jīng)誤差在0.5"左右，對應(yīng)的時應(yīng)誤差為1秒左右。

不過，請務(wù)必注意，這只是平均誤差，誤差分析表明，最大可能誤差是3"左右，對應(yīng)的時間誤差是6秒左右。

如果用標(biāo)準(zhǔn)差表達誤差，誤差大約為3秒左右。

這已足已證明它是個精確的月歷程序，程序不到10k

十、地球黃經(jīng)為W時對應(yīng)的時間to的求解
v(t) = 6283.32 //誤差小于萬分3
+210 *Math.sin(1.527+f)
+4.4 *Math.sin(1.48+f*2)
+12.9*Math.sin(5.82+f)*t
+0.55*Math.sin(4.21+f)*t*t;
式中f=6283.07585*t;
(1)解方程 W = Lon(t) = a+v*t-dW
方程右邊是地球黃經(jīng)表達式，a是它的常數(shù)項，v是黃經(jīng)速度，dW是黃經(jīng)的高次項(t的1次方以上的)和周期項，稱dW為方程的余項。
因dW較小，所以先忽略dW，解得“平”時間：
to = (W - a)/v = (W-1.753469512)/6283.319653318
v = v(to)
▲to是地球黃經(jīng)為W是的“平”時間，to與真時間的差值不超過48小時
▲余項引起的誤差：或余項已知，to回代到原方程，得t = (W-a+dW)/v = to+dW/v
顯然dW/v就是to的誤差值。
把to代入余項的表達式就可得到余項的值。
此外，也可以Lon()函數(shù)覺察到余項：Lon(to,10) = a+v*to - dW = W + dW，即dW = W-Lon(to,10)
▲因為to不超過6儒略千年，所以平黃經(jīng)的高次項引起的余項不超過1度，黃經(jīng)周期項不超過2度，合計余項最多3度。地球每天運動1度左右，所以to最大誤差3天。
v(t)求導(dǎo)得最大可能的加速度為(3.6弧度/千年)/天，如果to誤差3天，那么v(t)算得的v值誤差為3.6*3=10弧度/千年。
這樣，v的相對誤差為 10/6283=0.17%
(2)在步驟(1)留下的主要余項：dW = W-Lon(to,10);
   那么修正to為：
   to = to+ (W-Lon(to,10)) / v;
   ▲to的最大誤差由v及Lon()函數(shù)的截斷誤差引起的。
   v引起的最大可能誤差為 3天*0.17% = 7分鐘
   Lon(to,10)截斷誤差為17角秒，地球每分鐘運動2.5角秒，所以to誤差7分鐘左右。
   合誤差可以估計為10分鐘(以上兩個誤差同時發(fā)生最大值的可能性非常小，所以誤差不必取7+7=14分鐘)。
(3)執(zhí)行to = to +(W-Lon(to,120)) /v;
   v引起的誤差為(這就是最后的迭代誤差): 10分鐘*0.17%=1秒，平均誤差應(yīng)小于0.15秒
   Lon(to,120)的截斷誤差是0.25角秒，對應(yīng)時間誤差為6秒。

這是精簡的超高速的算法

//以下是地球黃經(jīng)數(shù)據(jù),最大誤差0.25"
var E_L0=new Array(
33416565,4.6692568,6283.07584999,
348943,4.626102,12566.1517,34971,2.74412,5753.38488,34176,2.82887,3.52312,31359,3.62767,77713.77147,
26762,4.41808,7860.41939,23427,6.13516,3930.2097,13243,0.74246,11506.76977,12732,2.0371,529.69097,
11992,1.10963,1577.34354,9903,5.2327,5884.9268,9019,2.0451,26.2983,8572,3.5085,398.149,
7798,1.1788,5223.6939,7531,2.5334,5507.5532,5053,4.5829,18849.2275,4924,4.2051,775.5226,
3567,2.9195,0.0673,3171,5.849,11790.6291,2841,1.8987,796.298,2710,0.3149,10977.0788,
2428,0.3448,5486.7778,2062,4.8065,2544.3144,2054,1.8695,5573.1428,2023,2.4577,6069.7768,
1555,0.8331,213.2991,1322,3.4112,2942.4634,1262,1.083,20.7754,1151,0.6454,0.9803,
1029,0.636,4694.003,1019,0.9757,15720.8388,1017,4.2668,7.1135,992,6.21,2146.165,
976,0.681,155.42,858,5.983,161000.686,851,1.299,6275.962,847,3.671,71430.696,
796,1.808,17260.155,788,3.037,12036.461,747,1.755,5088.629,739,3.503,3154.687,
735,4.679,801.821,696,0.833,9437.763,624,3.978,8827.39,611,1.818,7084.897,
570,2.784,6286.599,561,4.387,14143.495,556,3.47,6279.553,520,0.189,12139.554,
516,1.333,1748.016,511,0.283,5856.478,490,0.487,1194.447,410,5.368,8429.241,
409,2.399,19651.048,392,6.168,10447.388,368,6.041,10213.286,366,2.57,1059.382,
360,1.709,2352.866,356,1.776,6812.767,333,0.593,17789.846,304,0.443,83996.847,
300,2.74,1349.867,254,3.165,4690.48,247,0.215,3.59,237,0.485,8031.092,
236,2.065,3340.612,228,5.222,4705.732,219,5.556,553.569,214,1.426,16730.464,
211,4.148,951.718,203,0.371,283.859,199,5.222,12168.003,199,5.775,6309.374,
191,3.822,23581.258,189,5.386,149854.4,179,2.215,13367.973,175,4.561,135.065,
162,5.988,11769.854,151,4.196,6256.778,144,4.193,242.729,143,3.724,38.028,
140,4.401,6681.225,136,1.889,7632.943,125,1.131,5.523,121,2.622,955.6,
120,1.004,632.784,113,0.177,4164.312,108,0.327,103.093,105,0.939,11926.254,
105,5.359,1592.596,103,6.2,6438.496,100,6.029,5746.271,98,1,11371.7,
98,5.24,27511.47,94,2.62,5760.5,92,0.48,522.58,92,4.57,4292.33,
90,5.34,6386.17,86,4.17,7058.6,84,3.3,7234.79,84,4.54,25132.3,
81,6.11,4732.03,81,6.27,426.6,80,5.82,28.45,79,1,5643.18,
78,2.96,23013.54,77,3.12,7238.68,76,3.97,11499.66,73,4.39,316.39,
73,0.61,11513.88,72,4,74.78,71,0.32,263.08,68,5.91,90955.55,
66,3.66,17298.18,65,5.79,18073.7,63,4.72,6836.65,62,1.46,233141.31,
61,1.07,19804.83,60,3.32,6283.01,60,2.88,6283.14,55,2.45,12352.85);

var E_L1=new Array(
2060589,2.6782346,6283.07585,43034,2.63513,12566.1517,4253,1.5905,3.5231,1193,5.7956,26.2983,
1090,2.9662,1577.3435,935,2.592,18849.228,721,1.138,529.691,678,1.875,398.149,
673,4.409,5507.553,590,2.888,5223.694,560,2.175,155.42,454,0.398,796.298,
364,0.466,775.523,290,2.647,7.114,208,5.341,0.98,191,1.846,5486.778,
185,4.969,213.299,173,2.991,6275.962,162,0.032,2544.314,158,1.43,2146.165,
146,1.205,10977.079,125,2.834,1748.016,119,3.258,5088.629,118,5.274,1194.447,
115,2.075,4694.003,106,0.766,553.569,100,1.303,6286.599,97,4.24,1349.87,
95,2.7,242.73,86,5.64,951.72,76,5.3,2352.87,64,2.65,9437.76,
61,4.67,4690.48,58,1.77,1059.38,53,0.91,3154.69,52,5.66,71430.7,
52,1.85,801.82,50,1.42,6438.5,43,0.24,6812.77,43,0.77,10447.39,
41,5.24,7084.9,37,2,8031.09,36,2.43,14143.5,35,4.8,6279.55,
34,0.89,12036.46,34,3.86,1592.6,33,3.4,7632.94,32,0.62,8429.24,
32,3.19,4705.73,30,6.07,4292.33,30,1.43,5746.27,29,2.32,20.36,
27,0.93,5760.5,27,4.8,7234.79,25,6.22,6836.65,23,5,17789.85,
23,5.67,11499.66,21,5.2,11513.88,21,3.96,10213.29,21,2.27,522.58,
21,2.22,5856.48,21,2.55,25132.3,20,0.91,6256.78,19,0.53,3340.61,
19,4.74,83996.85,18,1.47,4164.31,18,3.02,5.52,18,3.03,5753.38,
16,4.64,3.29,16,6.12,5216.58,16,3.08,6681.22,15,4.2,13367.97,
14,1.19,3894.18,14,3.09,135.07,14,4.25,426.6,13,5.77,6040.35,
13,3.09,5643.18,13,2.09,6290.19,13,3.08,11926.25,12,3.45,536.8);

var E_L2=new Array(
87198,1.0721,6283.07585,3091,0.8673,12566.1517,273,0.053,3.523,163,5.188,26.298,
158,3.685,155.42,95,0.76,18849.23,89,2.06,77713.77,70,0.83,775.52,
51,4.66,1577.34,41,1.03,7.11,38,3.44,5573.14,35,5.14,796.3,
32,6.05,5507.55,30,1.19,242.73,29,6.12,529.69,27,0.31,398.15,
25,2.28,553.57,24,4.38,5223.69,21,3.75,0.98,17,0.9,951.72,
15,5.76,1349.87,14,4.36,1748.02,13,3.72,1194.45,13,2.95,6438.5,
12,2.97,2146.17,11,1.27,161000.69,10,0.6,3154.69,10,5.99,6286.6,
9,4.8,5088.63,9,5.23,7084.9,8,3.31,213.3,8,3.42,5486.78,
7,6.19,4690.48,7,3.43,4694,6,1.6,2544.31,6,1.98,801.82,
6,2.48,10977.08,5,1.44,6836.65,5,2.34,1592.6,5,1.31,4292.33,
5,3.81,149854.4,4,0.04,7234.79,4,4.94,7632.94,4,1.57,71430.7,
4,3.17,6309.37,3,0.99,6040.35,3,0.67,1059.38,3,3.18,2352.87,
3,3.55,8031.09,3,1.92,10447.39,3,2.52,6127.66,3,4.42,9437.76,
3,2.71,3894.18,3,0.67,25132.3,3,5.27,6812.77,3,0.55,6279.55);

var E_L3=new Array(
2892,5.8438,6283.0758,168,5.488,12566.152,30,5.2,155.42,
13,4.72,3.52,7,5.3,18849.23,6,5.97,242.73,4,3.79,553.57,1,4.3,6286.6,1,0.91,6127.66);

var E_L4=new Array(77,4.13,6283.08,8,3.84,12566.15,4,0.42,155.42);
var E_L5=new Array(2,2.77,6283.08,1,2.01,155.42);

//地球黃緯數(shù)據(jù),誤差0.2"
var E_B0=new Array(
2796,3.1987,84334.6616,1016,5.4225,5507.5532,804,3.88,5223.694,438,3.704,2352.866,
319,4,1577.344,227,3.985,1047.747);

var E_B1=new Array(90,3.9,5507.55,62,1.73,5223.69);
var E_B2=new Array(17,1.63,84334.66);

//地球向徑數(shù)據(jù),誤差0.00001AU
var E_R0=new Array(
1000139888,0,0,16706996,3.09846351,6283.07584999,
139560,3.055246,12566.1517,30837,5.19847,77713.77147,
16285,1.17388,5753.38488,15756,2.84685,7860.41939,
9248,5.4529,11506.7698,5424,4.5641,3930.2097,
4721,3.661,5884.9268,3460,0.9637,5507.5532,
3288,5.8998,5223.6939,3068,0.2987,5573.1428,
2432,4.2735,11790.6291,2118,5.8471,1577.3435,
1858,5.0219,10977.0788,1748,3.0119,18849.2275);

var E_R1=new Array(1030186,1.1074897,6283.07585,17212,1.06442,12566.1517,7022,3.1416,0);
var E_R2=new Array(43594,5.78455,6283.07585,1236,5.5793,12566.1517,123,3.142,0,88,3.63,77713.77);
var E_R3=new Array(1446,4.2732,6283.0758,67,3.92,12566.15);
var E_R4=new Array(39,2.56,6283.08,3,2.27,12566.15);
var E_R5=new Array(1,1.22,6283.08);

//迭代算法測試
t=1;
L=earthLon(t,-1);
t2=earth_t(L);
document.write("迭代算法測試<br>");
document.write("輸入時間(千年):"+t+"<br>");
document.write("地球黃經(jīng)(弧度):"+L+"<br>");
document.write("返算時間(千年):"+t2+"<br>");
document.write("迭代誤差(秒):"+(t2-t)*365250*86400+"<br><br>");

//=======================
// 以下是視位置計算
//=======================
var nutB=new Array(
2.1824, -33.75705, 36e-6,-1720,920,
3.5069, 1256.66393, 11e-6,-132, 57,
1.3375,16799.4182, -51e-6, -23, 10,
4.3649, -67.5141, 72e-6, 21, -9,
0.04,   -628.302,   0,     -14, 0,
2.36,   8328.691,   0,       7, 0,
3.46,   1884.966,   0,      -5, 2,
5.44, 16833.175,   0,      -4, 2,
3.69, 25128.110,   0,      -3, 0,
3.55,    628.362,   0,       2, 0);
function nutation(t){ //章動計算
var i,c,a, t2=t*t,dL=0,dE=0;
for(i=0;i<nutB.length;i+=5){
c = nutB[i]+nutB[i+1]*t+nutB[i+2]*t2;
if(i==0) a=-1.742*t; else a=0;
dL+=(nutB[i+3]+a)*Math.sin(c);
dE+= nutB[i+4]   *Math.cos(c);
}
this.dL=dL/100; //黃經(jīng)章動
this.dE=dE/100; //交角章動
}
function gxc_sun(t,L){ //太陽光行差(黃經(jīng)),L為太陽黃經(jīng),t是千年數(shù)
var v = L-(282.93735+17.1946*t+ 0.046*t*t)/180*Math.PI; //真近點角
var e = 0.016708634-0.00042037*t-0.00001267*t*t;
return -20.49552 * (1+e*Math.cos(v));
}
function HCconv(JW,E){ //黃赤轉(zhuǎn)換(黃赤坐標(biāo)旋轉(zhuǎn))
//黃道赤道坐標(biāo)變換,赤到黃E取負(fù)
var HJ=rad2mrad(JW[0]), HW=JW[1];
var sinE =Math.sin(E),cosE =Math.cos(E);
var sinW=cosE*Math.sin(HW)+sinE*Math.cos(HW)*Math.sin(HJ);
var J=Math.atan2( Math.sin(HJ)*cosE-Math.tan(HW)*sinE, Math.cos(HJ) );
JW[0]=rad2mrad(J);
JW[1]=Math.asin(sinW);
}

function hcjj (t){ //返回黃赤交角,t是千年數(shù)
var t2=t*t, t3=t2*t,t4=t3*t;
return (84381.4088 -468.36051*t -0.01667*t2 -1.99911*t3-0.00523*t4)/rad;
}

z=new Array();
da=1;
testT=0.008+(da-1.5)/365250;    //世紀(jì)數(shù)
earthCoor(testT,z,-1);         //地球黃經(jīng)
z[0]+=Math.PI,z[1]=-z[1];      //太陽黃經(jīng)黃緯

//L=z[0]-50287.9*(testT-0.008)/rad;
//L=rad2mrad(L);
//document.write(rad2str(L,0)+"<br>");

gx= gxc_sun(testT,z[0]);       //太陽光行差
nu= new nutation(testT*10);    //章動計算
z[0] += (nu.dL+gx)/rad;        //補章動與光行差
E = hcjj(testT)+nu.dE/rad;     //得到真黃赤交角
HCconv(z,E); //轉(zhuǎn)到赤道坐標(biāo)

document.write("坐標(biāo)精度測試<br>");
document.write("testT：" + testT+"<br>");
document.write("視赤經(jīng)：" + rad2str(z[0],1)+"<br>");
document.write("視赤緯：" + rad2str(z[1],0)+"<br>");
document.write("距離：" + z[2]+"<br>");
document.write("請與《2008中國天文年歷》太陽視赤經(jīng)赤緯比對。該年的第 "+da+" 天0h TD。注意：幾乎分秒不差！<br>");
document.write("本程序中VSOP87地球數(shù)據(jù)序列以做了長期項擬合到DE405/406<br>");

</script>