有關「現代曆法」——主要編者：songgz

06-14

　　幾個月前，增編寫個一個農曆計算javascript程序，那是我的第一個農曆程序。為了實現其中的演算法，曾花費了好幾個星期的時間研究天文計算相關的原理，當我算出結果之後，對程序的結果仍然沒有信心。不過，在「春光」老師介紹《天文演算法》一書之後，我認真閱讀並翻譯全書之後，我知道，我的演算法基本上沒有錯誤，但同時也認識到，天文學家們的演算法的確高明，比我當時的演算法好得多。如果你也對天文計算感興趣，那就很有必要學習他們的先進思想方法。　　當然，重要的是他們的思想方法以及相關的理論，而不一定是他們的計算結果。因為，《天文演算法》一書是早期出版的，有些數據比較舊，造成精度不一定很好。如果你對計算精度要球特別高，可能需要更換一些數據，但數據的處理的思想方法及技巧是一樣的。　　國內不少網友對天文計算比較感興趣，卻又不知如何下手。問題出在哪裡呢？為什麼苦苦研究幾個月甚至幾年也沒有進展？主要原因是國內有關的書籍太少，業餘條件下沒有機會學習到天文計算的理論。在我們國家，有不少機構對天文學有深入的研究，比如天文台、一些大學等，可他們不太願意出版天文演算法之類的書籍(這類書籍銷量少，出版要賠錢的)。既然如此，就讓我們自已想辦法解決問題吧。　　首先，我認為需要掌握一定的計算機程序設計技術，不要求很利害，但起碼也要有幾個月的程序設計經驗。如果你不會程序設計，那你只能用Excel或計算器之類的工具來處理計算問題，那簡直是在浪費你的生命。　　其次，應掌握《高等數學》中的一些知識。比如：極限、導數、微分、積分、極值問題、求根問題、最小二乘法、向量數學等。當然，我們更多使用高中的《立體幾何》、《解析幾何》、《函數》、《三角函數》等有關知識。還應了解《球面三角學》里的幾個公式。　　其三，《數值方法》這類書籍是必須讀的。　　其四，需對物理學有所了解，尤其是運動學相關知識。當然，如果你想用數值積分的方法解決天文計算問題，《理論力學》甚至是《天體力學》也是有必要了解的。　　對於多數具有大學學歷的人來說，基本具備以上知識，也就是說，只要你有興趣，就完全可以進行天文計算。從本質上講，日月運動、行星運動主要使用「牛頓力學」及數學方法(如微積分)，在牛頓那個年代，力學理論、數學理論、計算工具等都不可能和現在相比，在那個年代，就連「除法問題」主要是「教授」們才能掌握的！他們可以計算天體運動問題，我們為什麼就不可以呢？只要有信心，或多或少可以解決問題。　　相反，如果根本沒有學過《高中數學》《高中物理》《高等數學》，我建議你還是花點時間學習一下(弄不好要花費一兩年時間)，否則，即使算出了結果，也很難對你算出的結果形成理性的認識，甚至感性的認識也談不上。如果你不想學習這些「無用」的東西，能不能實現天文演算法，或許可以：通過閱讀別人的程序。　　在太陽系中，我們一致認為九大行星(八個主行星和1個冥王星)圍繞太陽系運動。幾百年前，哥白尼先生就已經這麼認為！現在的教科書中也都是這麼說的！因為太陽系中的星體圍繞太陽運動，所以我們常以太陽為中心來考慮問題。　　以太陽為中心考慮問題有不少好處：我們可以把行星圍繞太陽運動近似看做圓周運動(圓周運動是物理學中的一個重要名詞)。這樣，行星的運動顯得非常有規律。　　如果不以太陽為中心，而以地球為中心，我們會發現水星、金星的運動毫無規律，你可以充分發揮你的空間想像，你將發現，水星和金星在天球上的視運動不可能是圍繞地球做圓周運動，而是在太陽附近「擺動」。　　然而，我們不可能站在太陽上面進行天文觀測，通常是站在地球上進行天文觀測的，所以「地心說」也是非常科學的，也就是說，我們可以認為，太陽系的中心是地球。　　不管是「日心說」、「地心說」，只要你科學對待，不要隨意附加「神學的」或「迷信的」論調，都可以正確解天文現象。　　要記住，這個世界上沒有絕對的「中心」。太陽個頭大，質量大，所以行星圍繞它運轉。僅因它質量大，就能決定一切嗎？顯然不是。我們可以這麼想，某件事，是大人說了算還是小孩說了算，通常是大人說了算，但具體問題要具體分析，有時候小孩說了算。　　從純牛頓力學的角度看，不管是「日心」還是「地心」，通過計算，均能得到相同的數值計算結果。　　因此，根據實際需要，有時我們使用「地心」坐標，有時候使用「日心」坐標。　　為了計算方便以及更好的解釋天文現象，人們還引入了「太陽系質心」、「銀心」、地月系質心、月心、行星的「心」、純幾何的視中心，等等，一系列的「中心」。不管是什麼「中心」，我們必須正確理解這些「中心」的具體定義，必須明白這些「心」到底在哪裡，以便進行坐標變換。如果知道這些心的定義，我們未必原原本本的使用現成的天文坐標變換公式，也可適當根據程序設計的需要建立自已所需的坐標變換公式(比如在近似處理時)，使程序更加簡潔明了。

　　茫茫宇宙，無邊無際——宇宙大得讓我們難以想像。宇宙是球形的嗎？我們很難說清楚這個問題，但是，進入我們視野的天空，就象一個「球」，我們稱它為天球，我們位於天球的中心。這個球的半徑有多大？半徑很大很大，看成無窮大也無妨。太陽系內星體之間的距離根本不能與天球半徑相比，如果把天球半徑看為1，那麼太陽系內星體間的距離可以看做0。在我們的視野中，所有的恆星都在天球上。　　赤道的概念已被大家熟知。現在，我們擴展一下赤道的概念。赤道平面擴展到天球，與天球相交，交線為天球上的一個大圓，這個圓稱為天赤道。　　不管以太陽為中心，還是以地球為中心，我們看到的天赤道是相同的，所有恆星與天赤道的相對視位置關係是一樣的。為什麼呢？因為，從幾何學看來，太陽系的尺度太小，不能與天球半徑相比，在太陽系中的任何一個地方，都可以看做天球的中心，所以，恆星之間的視位置關係以及它們與天赤道之間的視覺位置關係不會因為中心的改變而發生改變。　　當然，恆星並不是真的在天球的表面上，個別恆星與我們之間的距離不是很遠，所以，當觀測點(即「中心」)改變時，它們的視位置也會有點變化，這就是我們常說的恆星周年視差，我們以後再討論這些問題吧。　　地球公轉的軌道面，同樣可以擴展到天球，得到的交線也是個大圓，稱為黃道。　　雖然赤道與天赤道是不同的兩個概念，但具體問題所涉及的是「天赤道」還是「赤道」往往是非常明確的(明顯的)，一般不會產生二意，所以通常「天赤道」也簡稱為「赤道」。　　天赤道與黃道必然在天球上產生兩個交點。一個稱為「春分點」(太陽經過升交點)，一個是「秋分點」(太陽經過降交點)。「春分點」是天球上重要的參考點。我們描述天體位置，通常是相對於這個參考點而言的。　　嚴格上說，黃道是地月質心公轉軌道在天球上的擴展。　　在天球上，並沒有一個恆星可以直接用來標定春分點的位置，它是看不見也摸不著的，我們如何確定它的位置？最重要的一種方法是：通過數以千計的恆星位置，反推出春風點在天球上的位置，我們常說的FK5天球坐標系統就與它有關。還有其它方法可以確定春分點，比如動力學方法。由於所用的方法不同，得到的春風點也會有一些微小的差別。難道有很多種春風點嗎？其實，「春風點」這個名詞最好不要亂用，如果說FK5中的那個「點」是真正的春分點，那麼動力學方法(如VSOP87)得到的那個點就不應該稱做春分點，其實在《天文演算法》一書中很少用「春風點」一詞，取代之「分點」一詞。當然，我們不一定計較那麼一點點微小的差別，統稱為「春分點」也不要緊。

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

　　當我們在天球上標定了赤道、黃道、分點之後，我們接著來看日、地在黃道坐標中的運動。以太陽系中的任意一點為天球的中心都不影響分點及黃道位置，當天球中心移到地心，地心直角坐標變為零，沒有經緯度，但可以得到太陽的經緯度，此時，如果跑到日心看地球，可以得到地球的經緯度，這兩個經緯度正好相反(日地連成的直線與天球相交的兩點正好完全相反)：　　地球經度 = 太陽經度+180度，地球緯度 = -太陽緯度。　　曆元、攝動、自球自轉軸、天北極、黃極、黃軸、黃經、黃緯、赤經、赤緯，這幾個概念不再講述，請在網路上找一些資料看看吧！　　春分點是黃道與赤道的交點！　　春分點是天球中的重要的參考點，遺憾的是，春分點會移動！　　在圖中，有兩條黃道、兩道赤道。

　　一條的曆元2000.0時刻黃道，一條是當日(Date)黃道。　　一條的曆元2000.0時刻赤道，一條是當日(Date)赤道。　　隨時間推移，黃道為什麼會發生變化？因為地球受至大行星的攝動。　　隨時間推移，赤道也會發生移動，因為受月球影響，地球自轉軸的方向發生移動，造成天北極圍繞黃軸緩慢移動。找一個直鐵棒，鐵棒放置的方向與黃軸相同，用電焊機把鐵棒焊在赤道面的中心。始終保持鐵棒與黃軸平行，並轉動鐵棒，這時你將看到，北極圍繞鐵棒做圓周運動，同時春分點也同步移動，當北極運行一周(約26000年)，春分點在黃道上也移動一周。　　春分點移動的後果：黃經、赤經的起算點是春分點，造成所有星體的黃經、赤經發生變化。假設，黃道不動(沒有發生變化)，赤道發生變化，引起春分點移動，如果春分點在黃道上向西移動了1度(在上圖中，往左下方移動)，那麼，所有星體的黃經增加了1度！　　春分點的移動量稱為歲差。歲差有很多種表達方法，最著名、最傳統的歲差是指黃道上的歲差。　　上圖中，有三個重要的交點：γo，γ1，γm，我們統稱它為分點。γo稱為「曆元2000.0標準分點」，其實就是曆元2000時刻的春風點，γm稱為「Date分點」或「當日分點」或「當日春分點」。　　圖中，兩個黃道的交點是N，定義N到γo與N到γo'的距離相等，那麼以γo'起算的黃經則與γo起算的黃經是相等的(除非天體過份靠近黃極)，對於行星，它在黃道附近，對於大部分恆星，也是遠離黃極的，所以γo'對於行星、月球、大部分恆星來說，是很重要的一個參考點。比如太陽，當它從γo出發，轉動數周后達到了γo'，所經歷的時間是整倍數的太陽繞地球運轉的周期。這種運行周期具有比較嚴格的力學意義。恆星在天空中的位置幾乎是不變的，以它作為參考系，具有「慣性」系的特徵，γo則是利用恆星確定的，也具有慣性特點。任意時刻，我們利用某一顆或幾顆恆星的位置，輕鬆的標定出γo的位置(當然，首次建立FK5系統時，需要大量恆星標定出γo，以提高精度，所以工程巨大)，當太陽運行到γo'，可以認為太陽相對於慣性參考點γo移動了一周或數個1周，即相對於慣性參考點移動了數個周期，這就是力學意義上的周期。γo與γo'相差很小，所以太陽到了γo與太陽到了γo'，太陽在星空的中相對位置幾乎相同，因此，太陽經過這兩點所經歷的時間是整倍數的恆星年(以恆星為參考，太陽在天球中運轉一周的時間)，由於這個原因，恆星年一般認為就是太陽(或地球)在慣性系中的運動周期，「在慣性系中的運動周期」常稱為真周期。

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

　　也許你會問，怎麼會有這麼多歲差參數，有必要嗎？回答，為了坐標變換方便而已。通常，只需只道3個歲數參數，就可利用坐標變換方法及多項式插值法導出其它的歲差參數。　　一組常用的歲差表達，取自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)　　上表是關於t的多項式係數，t等於J2000起算的儒略千年數，為儒略千年=365250天第一列對應t的0次方，第二列對應t的1次方，第三列對應t的2次方，其餘類推。如π=0.0*t^0 + 469.97560*t^1 - 3.35050*t^2 + ... 結果的單位是角秒。

　　以上歲差表達是非常精確的，誤差小於1毫角秒。　　不過t不能太大，t應小於1。　　如果t=10，那麼，會有什麼後果？比如計算θ，t的7次方項為0.00001*10000000=100角秒！實際上，表達式係數的舍入誤本身就是0.00001，所以誤差高達100角秒。　　如果由於某種需要，你非要計算t=10的情況，最好把t的5、6、7次方項忽略。這樣，如果t=5，誤差也可控制在1角秒之內。　　順便說一下，儒略千年數是一種時間計數方法。日常生活中，時間的計數法使用「年月日時分秒」，在天文計算過程，這種時間表達是很不方便的，原因是明顯的，需要6個數才能表達一個時間。　　天文學中，採用2000年1月1日12時0分0秒時作為計時的起點，並連續計時，單位是天。1天日86400秒。1秒是多少？1秒時間長度由原子鐘得到，現代天文學的時間基礎已不再使用星體運動位置來確定，改用原子鐘，原子鐘的精度要好得多。　　為什麼又稱為「儒略」千年數。這還要談到另一個問題，一年有多少天？你可能立刻回答到：365.2422日。然而，在現代天文學家眼裡，一年的長度是變化的，並不是固定為365.2422日，何況地球運動時還受到各種干擾，365.2422隻是一定時期內的平均值。因此，任意一年某時刻加上365.2422後，並不是下一年的同一時刻，這樣，把365.2422看做1個單位意義不大，天文學家更原意把365.25看做1個單位，它與公曆的置閏規則具有更好的「兼容」性，便於日曆轉換。365.25天是儒略年的長度。儒略1千年就是365250天。　　再次強調，1天=86400原子秒，與地球自轉沒有必然的關係。在歷史上，曾使用某一年的地球公轉與自轉來定義秒長，之後，地球自轉速度發生改變，1天不等於86400秒了。注意，此「天」非彼「天」，前面的天是時間單位，後面的天指平均自轉一周所用的時間長度。然而，人類的日常生活與地球自轉息息相關，日期與時間的計數必須儘可能與地球自轉同步，這就造成，我們有兩套基本的時間系統，一個是手錶時(或稱為協調世界時或稱為UTC時)，一個是力學時。力學時又時什麼東西呢？力學時具有鮮明的動力學色彩，假如我們使用2000年地球公轉一圈所用時間定義為1個單位，2萬年後，公轉周期可能發生了變化，如果不通過力學原理推算，我們將無法得到那時候一年的長度，但是，通過力學方法，我們可以得到2萬年後公轉的周期，比如是1.00001(亂寫一個)個單位。簡單的說，當我們定義了時間單位，我們就可以使用力學原理推算天體在某一位置所對應時間，這就夠成時間基礎(根據星體位置確定時間)。原子時就簡單多了，它不是通過天體位置來確定時間的，而是通過原子鐘內部的振蕩次數來來確定時間。原子鐘曾與力學時對準過，秒長是相同的，起點相差零點幾秒鐘。　　理論上，力學時、原子時是非常均勻的。手錶時與地球自轉(以恆星為參考點確定自轉速度)同步，而地球自轉時快時慢，所以不均勻。　　手錶時(去除地區的時差後的手錶時,即通過跳秒與UT1同步的時間)與原子時的差稱為ΔT，當然手錶時含有時區信息，如北京時間多了正好8個時。ΔT可由一組多項式表達式計算得到。（UT1同由自轉與恆星(春風點)確定的，是天文意義的時間系統，嚴格說，UT1與時角有關，即與角度有關；原子時則與「振蕩次數」有關；手錶時UTC的秒長與時角無關,它與原子時的秒長完全相同。不過，UTC通過年終的跳秒實現與UT1同步，這樣就造成UTC與原子時不同步。）　　手錶時的秒長是原子秒，為什麼手錶時不是原子時呢？因為國際上有人動了手腳，在每年的最後一天做了「正或負跳秒(正負閏秒)」處理，實現手錶時與地球自轉同步。

//原子時與UT1的差//部分參考了美國國家航空航天局網站里的演算法,www.nasa.gov//2050以後的外推參考了skymap10.5,在skymap的星曆表中分段找幾個數據,再用最小二乘法擬合double deltatT(double y){ //輸入年份(可帶小數點),如2000.78static double DTxs[15][12]={//{範圍下,範圍上,偏移,時間單位,多項各係數} {-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,力學時與世界時之差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;}我們在第(五)貼已經談到了赤道與黃道的變化將引起歲差。受到月球運動的影響，地球自轉軸發生圍繞黃軸的緩慢的長期的圓周運動，造成赤道發生變化，進而引起春風點變化。從上面的那個歲差示意圖中不難發現，春風點的變化不單單是赤道變化引起的，黃道的變化也會引起春風點的變化。不過，圖中的角π非常小，所以黃道變化引起和春風點移動是非常小的。從歲差示道意圖中的幾何關係：弧Nγ1 = 弧Nγm + χcos(ω)，注意，這是近似計算，但精度仍非常高兩邊同時減去П得： Nγ1 - П = (Nγm - П)+ χcos(ω)，即 ψ = p+χcos(ω) 其實ψ就是：黃道不動，赤道變化引起的黃經歲差，我們稱之是日月黃經歲差。其實χ就是：Date赤道上的行星歲差。如果把ψ = p+χcos(ω)改寫為p = ψ-χcos(ω) 那麼 -χcos(ω) 就是行星黃經歲差。同樣，從幾何關係，還可輕鬆得到行星黃緯歲差：χsin(ω)---------------- χ的歲差速度是：106"/千年，加速度為-238"/平方千年 χ歲差造成春風點東移：0".106*cos(23.5°)/年日月歲差(黃經)造成春風點西移：50.384.78750/年黃經總歲差速度約為二者之差：50".29/年---------------- [章動]：受到日月運動的影響，地球自轉軸發生振動，這種振動稱為章動。「章」字，從何講起？地球自轉的振周期主要與月球軌道升交點運動有關，升交點運動周期為18.6年，古人稱之為一章。在第(五)貼的圖中，已給出受到章動後的赤道位置。圖中的Δψ就是黃經章動值。章動的後果：春風點在黃道上移動了Δψ(稱為黃經章動)，黃赤交角變化了Δε(稱為交角章動) 考慮章動以後的赤道稱為真赤道，否則為平赤道。黃道的振動是非常微弱，所以一般不再區分「真」與「平」要得到真春點起算的黃經，必須加上Δψ，要繼續旋轉到真赤道坐標中，應使用真黃赤交角(即加上Δε)----------------- Δψ與Δε的計算前面已發了一個關於IAU1980與IAU2000B章動序列比較的貼子，裡面含有相關的所有計算，請閱讀程序。---------------- 精度控制完整的計算，計算量偏大，程序冗長。如果你只是為了計算農曆，根本不需要那麼高的精度。為什麼呢？因為，現代的紫金山中國農曆是基於北京手錶時的，這種時間與地球自轉有關，考慮到地球地轉的不規律性，無法對幾百年後的北京時間進行測算，誤差0.5分鐘是正常的，如果對5千年後的北京時間預測計算，誤差可達幾小時。0.5分鐘的誤差意味差什麼呢？意味著太陽位置誤差可達1".2，月亮位置誤差可達15"，因此沒有必要把Δψ算得很准，0".1已經足夠。以下javascript程序給出幾千年內精度為0".05的黃經章動及交角章動。<script language=javascript>/*********************對IAU1980章動計算的簡化處理許劍偉 2008年5月2日*********************/t=1; //J2000.0起算的儒略世紀數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); //黃經章動 dE3+=(nuB[i+5]+nuB[i+6]*t/10)*Math.cos(c); //交角章動}dL3/=10000; //黃經章動dE3/=10000; //交角章動document.write("[精簡的章動計算,精度0.05角秒,單位角秒]<br>");document.write("黃經章動：" + dL3 + " 交角章動：" + dE3 + "<br>");

</script>

--------------------- [問題]：程序中使用了「J2000.0起算的儒略世紀數」，你會計算它嗎？ [例題]：比如，2002年1月1日12:00:00 TD的儒略世紀數T是多少？TD表示力學時。 [解]： J2000.0對應2000年1月1日12:00:00 TD 由於2000年是閏年，含有366天，也就是說： 2001年1月1日12:00:00加上366天後變為2001年1月1日12:00:00 又因為2001年是平年，所以再加上365天就是2002年1月1日12:00:00 所以，所求時刻相對於J2000.0相差366+365=731天那麼，T = 731/36525 = 0.020013689

行星光行差計算起來比較麻煩，也沒有現成的簡單公式可以應用，我們總是根據具體的行星問題展開計算。用下文所述的方法得到的光行差是非常準確的，最後得到的行星視位置與skymap的計算結果分毫不差。當然，日月計算也會涉及行星光行差問題，但這隻以行星光行差中的最簡單的一個特例。

此主題相關圖片如下：

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

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

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

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

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

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

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

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

剛才說到「低精度」演算法，具體含義是什麼？我們的目標是準確計算出位置的差值，實際上是差值的精度達到3至4位有效數字即可。因為光行星光行差通常為幾十角秒，如果達到3位有效數字，誤差已小於0.05角秒。

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

如果你使用VSOP87方法進行星曆計算，建議使用這種方法進行光行差修正。因為，你已經擁有了星體位置計算的函數，調用該函數計算兩次位置即可輕鬆算出行星光行差。當然，要適當控制計算的項數，以免影響程序效率。

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

幾個注意點：

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

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

通過以下改正,可以大幅度提高VSOP87的精度1、以下對VSOP87日心球面地球序列擬合到DE406: 利用VSOP87計算出黃經L後,計算黃經的修正量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起算的儒略千年數最後 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" 黃經擬合時已考慮了最新的歲差速度訂正: -2.9965"/儒略千年黃緯與距離不用修正。3、使用J2000.0動力學分點,因此J2000黃道坐標中的黃經與DE405/406相差一個估定值0.04"。注意1，使用J2000.0分點(非常接近ICRF的赤經起標點)得到的星體坐標才能與最新的《中國天文年曆》相同。注意2,公元3000以後的擬合了「瑞士星曆表」4、許劍偉,地震後第12天《天文演算法》一書中給出了求各種天文現象發生時刻的方法，不過作者並沒有深入討論演算法的效率、簡潔性及具體原理。事實書中演算法的精度是固定的，不一定滿足現代天文計算的要求。比如試圖把農曆的精度計算到30秒,使用《天文演算法》一書中的演算法將不能實現。本貼討論高精度的、超高速的演算法。問題假設：設經度L(t) = Lo + v*t + f(t)，對於方程L(t)=W，求t的值。式中速度v值較大,f(t)對t求導數也得到速度量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)第二步,用第一步的方法重新執行一次,不過計算過程中可以考慮v使用v+dv代替。(3)第三步,重複第二步計算,但v一定要使用v+dv代替注意事項1:計算L(t)時不一定需要完全精度,第一、二步計算時只需計算幾個周期項即可。第三步則需多算幾個項。如果精度還不滿意，則重複第三步注意事項2:v的解析式可以使用VSOP87或ELP/MPP02的序列求導數得到,只需考慮速度量較大的幾個項即可。注意事項3:如果使用DE405/DE406星曆表或「瑞士星曆表」,可以考慮使用求變差的辦法求v+dv。「變差」法，指求相近兩個時間的坐標差值，除以時間差後得到速度。說明1：以上方法本質上就是《天文演算法》中的迭代逼近思想方法。說明2：上面講到的v的運動學本質是「平角速度」，用2*3.1415926535/v得到「平周期」，如地球運動的「平周期」是365.2422天說明3：對v+dv的取值還需與截斷誤差合併分析，這當中涉及「加速度」的分析，將在以後的貼子中討論。說明4：以上所述的是通用方法，對於實際問題，有必要進行更仔細的分析。

計算量：一次精確計算位置坐標的計算量為B，那麼該演算法的計算量小於1.3*B

演算法的極限精度取決於VSOP或ELP/MPP02理論本身《農曆天文演算法完全解——第一部分(分兩部分講述)》演算法本質：迭代技術一、日月黃經速度的表達 (1)地球黃經速度的表達(Date黃道坐標中的速度,含歲差速度) 以下表達式中t是千年數,誤差小於萬分之3,數據由VSOP87地球黃經序列求導數得到: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)月球黃經速度的表達(Date黃道坐標中的速度,不含歲差速度) 歲差速度遠小於月球黃經速度,所以下速度計算忽略了歲差速度。平速度是(單位：單位弧度/世紀,式中的t是世經數)：v=8400; //精確值是8399.68473007193 周期項速度(對ELP/MPP02黃經序列求導數得到)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的值可根據精度的需要從上面表達式中取一種。月球真速度為：v=v-dv; (3)月球平黃經與黃經的表達式(ELP/MPP02)w1 = 3.81034409083 +8399.68473007193*t -3.31895204255e-05*t2 +3.11024944911e-08*t3 -2.03282376489e-10*t4L=w1+周期項+泊松項，由於泊松項較小僅需考慮幾項就可以了，除非有高精度計算的需要。周期項與泊松項的具體演算法待敘。二、求月球經度為W時所對應的時間經度計算函數為 Lon(t,n); t為J2000.0起算的世紀數,t不超過+-30世紀;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); //為下一次計算準備速度,速度誤差5%左右以上兩行相當於解方程 W=3.81034 +8400*t,解為to = (W - 3.81034)/8400 to是月球黃經為W是的「平」時間 to與真時間的差值不超過18小時以上計算已經考慮了黃經的常數項(3.81034)和速度項(8400)。因為to的超不超過30，所以平黃經的高次項(加速度項)很小，不超過2度，黃經周期項不超過8度，所以黃經的的高次項及周期項總和不超過10度。因此如果把to再次回代到完整的序列中計算，上述方程最多只有10度(約0.17弧度)的誤差： W=3.81034 +8400*t+0.17，得to = (W -3.81034 -0.17)/8400，與原來的to比較不難發現，最大誤差就是月球運動0.17(即10度)所用的時間,約18小時。這18個小時的誤差主要是截斷引起的。(2)a.余項計算：dW = W - Lon(to,3); b.時間計算：to=to+dW/v; c.速度計算：v=v2(to); //為下一次計算準備速度,精度為0.3%，如果考慮由於時間不準確引起的額外誤差，精度仍可達到0.35%(取0.4%) 在步驟(1)中已經說到,誤差為18小時,我們稱之為余項(或稱之為截斷誤差),而本次計算已經截取了3個主要周期項，可修正這18小時的余項。由於速度的精度只有5%，所以殘留誤差為18*5%=0.9小時。不過，實際誤差不止這些，因為本次計算也存在截斷誤差。截斷部分誤差估計: 把被截斷的最大10項的振幅取和作為最大誤差估計,以後的項不用算(因為存在正負相消的作用),這10項取和為3000角秒=50角分，對應的最大時間誤差為1.7小時。與剛才討到的0.9小時綜合考慮，不妨把誤差估計為2小時。為什麼不把誤差估計為0.9+1.7=2.6小時？因為0.9與1.7的誤差估計值均指最大誤差，出現最大誤差的可能性已經很小了，二者同時出現最大值的可能性幾乎為0。當然，平均誤差僅為最大誤差的六分之一左右。第(1)步與第(2)步計算無需考慮光行時、章動的修正計算，因為誤差遠超過這些修正量。(3)a.余項計算：dW = W - Lon(to,21); b.時間計算：to = to + dW/v; c.無需再計算更精確的速度v 周期項取21位的理論精度是4',對應的時間誤差為8分鐘。速度不精確引起的誤差：2小時*0.004=0.008小時=0.48分鐘。總誤差不妨按8分鐘估計。(4)a.余項計算：dW = W - Lon(to,200); b.時間計算：to = to+dW/v; 計算200項的理論最大誤差約2角秒,對應的時間誤差為4秒鐘速度v不精確引起的誤差：8分*0.004=0.3秒總誤差不妨取4秒鐘。 to就是我們最終所需要計算的時間！(5)如果相要再提高精度,繼續計算： dW = W - Lon(to,全部); to = to+dW/v; 「全部」計算的精度為5毫角秒，相當於0.01秒速度v不精確引起的誤差：5秒*0.004=0.02秒總誤差不妨取0.03秒鐘。 ELP/MPP02在J2000前後幾十年內與DE406擬合得很好，差異為0.005角秒以內。(6)說明：由於都是在鄰近的時間內計算，可以認為速度不變(或者說加速度為0)。因為對ELP/MPP02求二次導數，我們發現加速度非常小。為什麼會這樣呢？其實，周期項的周期最小值為10天左右(個別為5天左右,但振幅非常小)，也就是說在10天範圍內速度大小才小幅波動一次，在幾小時或幾分鐘內，速度幾乎不變。速度大小波動幅度最大的周期項的周期是27.5天左右，它對速度的影響可達10%以上，即便是這種周期項，10小時對速度的影響最多只有：10%*10小時/27.5天 = 0.15%，2小時誤差對速度精度的影響為：0.03%(7)說明2：如果要計算視黃經的,應注意光行差、光行時、歲差、章動修正等。月球光行時約為1.2秒(約修正0.6角秒)。(8)計算量：以上針對月球的精確到5秒鐘(指最大可能的誤差為5秒,平均誤差不到1秒，標準差約為3秒)的總計算量約為250個周期項的計算量。三、地球經度為W時對應的時刻的演算法與月球的類似,但迭代次數要少得多。如果考慮了光行差，那麼該演算法可直接用於計算24節氣。四、日月距角為W時對應的時刻。方法同上。不同的是把上面的經度換為月球黃與太陽黃經的差值，速度換為月球黃經速度與太陽黃經速度的差。距角與月相密切相關，距角為0度對應新月，為90度是半滿上弦月，180度是滿月。由於章動效果同時影響太陽和月球坐標，所以不用計算章動。當然地球運動引起的恆星光行差是必須計算的，否則無法正確計算出月相，因為幾何意上的日月合朔不是光學(視覺)意義上的日月合朔。計算距角時，地球黃經的精度控制在1角秒已足夠，因為月球的精度也只控制在2角秒。五、這裡描述的演算法與本人在javascript農曆程序中提供的相應演算法有所不同：速度提高了4至5倍，精度提高了5倍。六、精度是根據程序員的意願決定的，但受到VSOP87或DE405/406理論的極限精度制約。七、如果讀者還時一知半解，無法寫出相應的程序，本考慮編寫具體的程序，以做示範。八、ELP/MPP02的黃經起算點與DE405/406星曆表有關應變換到J2000黃道坐標中，以免產生0.1秒的時間誤差。九、VSOP87是早期的理論，存在起點誤差以及微小的長期項誤差，同樣需做訂正，訂正後，幾百年內的誤差小於0.5秒。訂正方法在前面的貼子中已經講述。十、使用真正的DE405/DE406方法的計算時間，方法同上，只是在計算Lon()時使用DE405/406方法改進《農曆天文演算法完全解——第一部分》演算法本質：迭代技術一、日月黃經速度的表達 (1)地球黃經速度的表達(Date黃道坐標中的速度,含歲差速度) 以下表達式中t是千年數,誤差小於萬分之3,數據由VSOP87地球黃經序列求導數得到: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)月球黃經速度的表達(Date黃道坐標中的速度,不含歲差速度) 歲差速度遠小於月球黃經速度,所以下速度計算忽略了歲差速度。平速度是(單位：單位弧度/世紀,式中的t是世經數)：v=8399.68473007193 考慮歲差速度後,應把式中的8399.68473007193換為：8399.68473007193弧度/世紀+5028.79角秒/世紀=8399.70911033384弧度/世紀周期項速度(對ELP/MPP02黃經序列求導數得到)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的值可根據精度的需要從上面表達式中取一種。月球真速度為：v=v-dv; (3)月球平黃經與黃經的表達式(ELP/MPP02)w1 = 3.81034409083 +8399.68473007193*t -3.31895204255e-05*t2 +3.11024944911e-08*t3 -2.03282376489e-10*t4 注意:w1中不含歲差速度 (4)黃經表達式(相對於Date平分點): L=w1+周期項+泊松項，由於泊松項較小僅需考慮幾項就可以了，除非有高精度計算的需要。 (5)周期項與泊松項的具體演算法待敘。

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

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

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

本程序未經檢驗，可能有誤，有空的時候我會與《中國天文年曆》及「瑞士星曆表」比對一下。到時再做報告。

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; //項數是以周期項為基準,n是周期項的計算項數 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){ //月球經度計算,返回Date分點黃經,傳入世紀數 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; //月球平黃經 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){ //月球速度計算,傳入世經數 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){ //已知黃經求時間 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; //世紀數L=moonLon(testT,-1); //正算t=moon_t(L); //反算dt=(testT-t)*35625*86400;

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

</script>

為了能夠與《中國天文年曆》比對，我們須把以上程序中的月球坐標轉換為視坐標

這時我們應處理光行差及章動

月球光行差問題我們已討論過精密的光行差計算，但是，在農曆計算中根本無需精密計算，以下考慮簡化計算。從ELP的中序列截取5%精度的距離及速速表達式：日月距離：r=385001 +20905*sin( 0.7848 +8328.691425*t+ 1.523*t*t/10000 ); //最大誤差2% 黃經速度：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千米/世紀那麼有：r = r0+dr；v = v0+dv；τ= r/c； (式中τ光行時) 光行期間黃經的變化量： 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(φ)) 結論：月球的黃經光行差約為0.7角秒(最大誤差16%)，如果考慮上式中的周期，最大誤差5%

以下關於章動，取自IAU1980，最大誤差不超過0.1角秒function nutation(t){ //章動計算,t為世紀數,返回值為角秒單位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; //黃經章動this.dE=dE; //交角章動}

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

那麼：nu.dL為黃經章動，nu.dE為交角章動

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

rad = 180*3600/Math.PI;function rad2mrad(v){ //對超過0-2PI的角度轉為0-2PI v=v % (2*Math.PI); if(v<0) return v+2*Math.PI; return v;}function rad2str(d,tim){ //將弧度轉為字串//tim=0輸出格式示例: -23°59' 48.23"//tim=1輸出格式示例: 18h 29m 44.52svar 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; //世紀數L=moonLon(testT,-1); //正算t=moon_t(L); //反算dt=(testT-t)*35625*86400;

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

//以下是視位置計算function nutation(t){ //章動計算,t為世紀數,返回值為角秒單位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; //黃經章動this.dE=dE; //交角章動}

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

測試結果如下：

隨機插出幾個月球視黃經數據進行比較

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" 本程序我們不難發現，誤差與精密星曆相差無幾，黃經誤差在0.5"左右，對應的時應誤差為1秒左右。

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

如果用標準差表達誤差，誤差大約為3秒左右。

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

十、地球黃經為W時對應的時間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 方程右邊是地球黃經表達式，a是它的常數項，v是黃經速度，dW是黃經的高次項(t的1次方以上的)和周期項，稱dW為方程的余項。因dW較小，所以先忽略dW，解得「平」時間： to = (W - a)/v = (W-1.753469512)/6283.319653318 v = v(to) ▲to是地球黃經為W是的「平」時間，to與真時間的差值不超過48小時 ▲余項引起的誤差：或余項已知，to回代到原方程，得t = (W-a+dW)/v = to+dW/v 顯然dW/v就是to的誤差值。把to代入余項的表達式就可得到余項的值。此外，也可以Lon()函數覺察到余項：Lon(to,10) = a+v*to - dW = W + dW，即dW = W-Lon(to,10) ▲因為to不超過6儒略千年，所以平黃經的高次項引起的余項不超過1度，黃經周期項不超過2度，合計余項最多3度。地球每天運動1度左右，所以to最大誤差3天。 v(t)求導得最大可能的加速度為(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()函數的截斷誤差引起的。 v引起的最大可能誤差為 3天*0.17% = 7分鐘 Lon(to,10)截斷誤差為17角秒，地球每分鐘運動2.5角秒，所以to誤差7分鐘左右。合誤差可以估計為10分鐘(以上兩個誤差同時發生最大值的可能性非常小，所以誤差不必取7+7=14分鐘)。(3)執行to = to +(W-Lon(to,120)) /v; v引起的誤差為(這就是最後的迭代誤差): 10分鐘*0.17%=1秒，平均誤差應小於0.15秒 Lon(to,120)的截斷誤差是0.25角秒，對應時間誤差為6秒。

這是精簡的超高速的演算法

//以下是地球黃經數據,最大誤差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);

//地球黃緯數據,誤差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);

//地球向徑數據,誤差0.00001AUvar 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("地球黃經(弧度):"+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; //黃經章動this.dE=dE/100; //交角章動}function gxc_sun(t,L){ //太陽光行差(黃經),L為太陽黃經,t是千年數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){ //黃赤轉換(黃赤坐標旋轉) //黃道赤道坐標變換,赤到黃E取負 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是千年數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; //世紀數earthCoor(testT,z,-1); //地球黃經z[0]+=Math.PI,z[1]=-z[1]; //太陽黃經黃緯

//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); //轉到赤道坐標

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

</script>