九色国产,午夜在线视频,新黄色网址,九九色综合,天天做夜夜做久久做狠狠,天天躁夜夜躁狠狠躁2021a,久久不卡一区二区三区

農(nóng)歷24節(jié)氣算法

 [摘要] 古老而又現(xiàn)代的中國農(nóng)歷，是一種天文學(xué)性質(zhì)的陰陽歷。本文采用VSOP87星歷算法并結(jié)合運(yùn)動學(xué)方法以及牛頓求根法得到24節(jié)氣的精確時間，揭開農(nóng)歷計算的神秘面紗。

[關(guān)鍵字] 農(nóng)歷算法、星歷、節(jié)氣

[正文] 計算中國農(nóng)歷，首先要計算出二十四節(jié)氣時刻。在計算機(jī)問世之前，二十四節(jié)氣的許算是非常復(fù)雜的。隨著計算機(jī)及互聯(lián)網(wǎng)的普及，美國航空航天局、法國巴黎天文臺各自在網(wǎng)絡(luò)上發(fā)布了精密星歷表的計算方法，這使得民間計算農(nóng)歷成為可能。本文以法國巴黎天文臺的VSOP87算法為基礎(chǔ)，給出中國農(nóng)歷的二十四節(jié)氣算法。

在農(nóng)歷中，太陽黃經(jīng)為0度時，對應(yīng)春風(fēng)節(jié)氣。相鄰節(jié)氣對應(yīng)的太陽黃經(jīng)相差15度。一周年內(nèi)，太陽黃經(jīng)從0度變化到360度，共有24個節(jié)氣。

一、時間標(biāo)尺——儒略日數(shù)計算

計算星歷之前首先要解決時間尺問題。公歷規(guī)定平年365日，閏年366日。1582年10月4日以前，公歷規(guī)定每4年設(shè)置一個閏年，平均年長度365.25天，這期間的公歷稱為儒略歷。在1582年10月15日之后實行格里高利歷，規(guī)定每400年97閏，平均年長度為365.2425天。

由于儒略歷存在嚴(yán)重的“多閏”問題，到了1582年，公歷跑快了10天左右，當(dāng)時就人為調(diào)整了10天，并從此實行格里歷。因此務(wù)必注意1582年10月4日（儒略歷）的下一日為1582年10月15日（格里歷）。就是說1582年10月份少了10天。

在儒略歷中，能被4整除的年份為閏年，這一年有366天，其它年份為平年（365天）。 如900年和1236年為閏年，而750年和1429年為平年。

格里高利歷法也采用這一規(guī)則，但下列年份除外：不能被100整除的年份為平年，如1700年，1800年，1900年和2100年。其余能被400整除的年份則為閏年，如1600年，2000年和2400年。

儒略日數(shù)（簡稱儒略日）：

儒略日數(shù)是指從公元 -4712 年開始連續(xù)計算日數(shù)得出的天數(shù)及不滿一日的小數(shù)，通常記為 JD (**)。傳統(tǒng)上儒略日的計數(shù)是從格林尼治平午，即世界時12點(diǎn)開始的。若以力學(xué)時（或歷書時）為標(biāo)尺，這種計數(shù)通常表達(dá)為“儒略歷書日”，即JDE (**)，其中E只是一種表征，即按每天86400個標(biāo)準(zhǔn)秒長嚴(yán)格地計日。例如：

1977年4月26.4日 UT = JD 2443259.9

1977年4月26.4日 TD = JDE 2443259.9

儒略日的計算：

設(shè)Y為給定年份，M為月份，D為該月日期（可以帶小數(shù)）。

若M > 2，Y和M不變，若 M =1或2，以Y–1代Y，以M+12代M，換句話說，如果日期在1月或2月，則被看作是在前一年的13月或14月。

對格里高利歷有 ：A = INT（Y/100） B = 2 - A + INT(A/4)

對儒略歷，取 B = 0

儒略日即為：

JD = INT(365.25(Y+4716))+INT(30.6001(M+1))+D+B-1524.5

使用數(shù)值30.6取代30.6001才是正確的，但我們?nèi)允褂?/font>30.6001，以確保總能取得恰當(dāng)?shù)恼麛?shù)。事實上可用30.601甚至30.61來取代30.6001。例如，5乘30.6精確等于153，然而大多數(shù)計算機(jī)不能精確表示出30.6，這導(dǎo)致得出一個152.999 9998的結(jié)果，它的整數(shù)部分為152，如此算出的JD就不正確了。

由儒略日推算歷日：

將JD加上0.5，令 Z 為其整數(shù)部分，F 為尾數(shù)（小數(shù)）部分。

若 Z < 2299161，取A = Z

若 Z 大于等于2299 161，計算

α=INT((Z-1867216.25)/36524.25)

A=Z+1+α-INT(α/4)

然后計算

B = A+1524

C = INT((B-122.1)/365.25)

D = INT(365.25C)

E = INT((B-D)/30.6001)

該月日期（帶小數(shù)部分）則為：

d = B - D - INT(30.6001E) + F

月份m為：

IF E < 14 THEN m = E – 1

IF E=14 or E=15 THEN m = E – 13

年份為y：

IF m>2 THEN y = C – 4716

IF m =1 or m=2 THEN y = C – 4715

這個公式里求E時用的數(shù)30.6001不能代之以30.6，哪怕計算機(jī)沒有先前所說的問題。否則，你得到的結(jié)果會是2月0日而不是1月31日，或者4月0日而不是3月31日。

值得記住的一個常數(shù)是：2000年1月1日12:00:00的儒略日數(shù)是J2000 = 2451545

二、力學(xué)時與世界時的差值（deltat T）計算

一般的，可以把手表時（UTC）近似看作世界時（UT），二者的主要差別在于時區(qū)。如北京手表時8點(diǎn)對應(yīng)世界時0點(diǎn)。世界時與地球自轉(zhuǎn)嚴(yán)格同步，但有趣的是，我們的手表時實際上稱為協(xié)調(diào)世界時，它的秒長是原子鐘的秒長，由于地球自轉(zhuǎn)速度不均勻，時快時慢，這就注定手表時與地球自轉(zhuǎn)不完全同步?，F(xiàn)在，地球自轉(zhuǎn)速度正在變慢，我們不得不在某些年份的年末把手表撥慢1秒，使得手表時更好的與地球自轉(zhuǎn)同步，并美言為“跳秒”。力學(xué)時是根據(jù)太陽系的動力學(xué)原理導(dǎo)出的，是一種均勻的時間系統(tǒng)，其秒長與原子鐘的秒長相同。因此，協(xié)調(diào)世界時（UTC）與世界時（記為UT）其本同步，但力學(xué)時（記作TD）與世界時不太同步，二者的差值記作deltat T或記作△T。利用直接的天文觀測可以得知每年的△T，利用古代的日月食觀測資料可以反推古代的△T。所有年份的△T計算出來后，可以擬合出以下多項式表達(dá)，使得△T的計算更快捷，計算結(jié)果的單位是秒。

我們利用下表可以嚴(yán)格計算△T（即△T =TD - UT）

年份     a        b      c      d

-4000,108371.7,-13036.80,392.000, 0.0000

 -500, 17201.0,  -627.82, 16.170,-0.3413

 -150, 12200.6,  -346.41,  5.403,-0.1593

  150,  9113.8,  -328.13, -1.647, 0.0377

  500,  5707.5,  -391.41,  0.915, 0.3145

  900,  2203.4,  -283.45, 13.034,-0.1778

 1300,   490.1,   -57.35,  2.085,-0.0072

 1600,   120.0,    -9.81, -1.532, 0.1403

 1700,    10.2,    -0.91,  0.510,-0.0370

 1800,    13.4,    -0.72,  0.202,-0.0193

 1830,     7.8,    -1.81,  0.416,-0.0247

 1860,     8.3,    -0.13, -0.406, 0.0292

 1880,    -5.4,     0.32, -0.183, 0.0173

 1900,    -2.3,     2.06,  0.169,-0.0135

 1920,    21.2,     1.69, -0.304, 0.0167

 1940,    24.2,     1.22, -0.064, 0.0031

 1960,    33.2,     0.51,  0.231,-0.0109

 1980,    51.0,     1.29, -0.026, 0.0032

 2000,    63.87,    0.1,   0,     0,

 2005

表中每一行適用一定的年代范圍，如第1行適用于公元-4000年到-500年，第2行適用于公元-500到-1500年，其它類推。每行的起始年份記作Y1，終止年份記作Y2，如果年份y在Y1到Y2之間，那么該年的deltat T表達(dá)為：

△T = a + b*t1 + c*t2 + d*t3，單位是秒

其中t1 = （y-Y1）/（Y2-Y1）*10， t2 = t1*t1， t3 = t1*t1*t1

對于2005年以后的deltat T是未知的，要做外推計算：

2005至2014年建議使用1995到2005年期間△T的平均增速計算，即：

△T = F(y) = 64.7 + (y-2005) * b， 其中速度 b = 0.4

2114年以后可以使用二次曲線外推

△T = f(y) = -20+ a * [(y-1820)/100]^2 ，其中加速度a = 31

2114年到2014年之間的外推，可以在上面兩個外推算式的基礎(chǔ)上做一次的曲線連接，使之連續(xù)即可。比如可以這么計算：

△T  = f(y) + (y-2114) * [f(2014) – F(2014)] /100

以下數(shù)值可供程序驗證參考

2008年△T = 66.0秒

1950年△T = 29秒

500年 △T = 5710秒

三、太陽視黃經(jīng)（真分點(diǎn)視坐標(biāo)）

算法基于VSOP87半解析法。

力學(xué)時t為J2000起算的儒略世紀(jì)數(shù)，t2 = t*t，t3 = t2*t，t4 = t3*t

A、低精度算法

L0(t) = 48950621.66 + 6283319653.318*t 弧度

B、中精度算法

  L1(t) =  [ 48950621.66 + 6283319653.318*t + 53*t*t

 + 334116*cos( 4.67+628.307585*t)

 + 2061*cos( 2.678+628.3076*t)*t ] / 10000000 弧度

C、高精度算法

L2(t) = [ 48950621.66 + 6283319653.318*t

+ 52.9674*t2 + 0.00432*t3 - 0.001124*t4

 +334166 * cos( 4.669257+  628.307585*t)

   +3489 * cos( 4.6261  + 1256.61517*t )

   + 350 * cos( 2.744   +  575.3385*t)

   + 342 * cos( 2.829   +    0.3523*t)

   + 314 * cos( 3.628   + 7771.3771*t)

   + 268 * cos( 4.418   +  786.0419*t)

   + 234 * cos( 6.135   +  393.021*t )

   + 132 * cos( 0.742   + 1150.677*t )

   + 127 * cos( 2.037   +   52.9691*t)

   + 120 * cos( 1.11    +  157.7344*t)

   +  99 * cos( 5.23    +  588.493*t )

   +  90 * cos( 2.05    +    2.63*t  )

   +  86 * cos( 3.51    +   39.815*t )

   +  78 * cos( 1.18    +  522.369*t )

   +  75 * cos( 2.53    +  550.755*t )

   +  51 * cos( 4.58    + 1884.923*t )

   +  49 * cos( 4.21    +   77.552*t )

   +  36 * cos( 2.92    +    0.07*t  )

   +  32 * cos( 5.85    + 1179.063*t )

   +  28 * cos( 1.9     +   79.63*t  )

   +  27 * cos( 0.31    + 1097.71*t  )

 +2060.6 * cos( 2.67823 +  628.307585*t ) * t

   +43.0 * cos( 2.635   + 1256.6152*t   ) * t

   +8.72 * cos( 1.072   +  628.3076*t   ) * t2

   -994 – 834 * sin(2.1824-33.75705*t)

- 64 * sin(3.5069+1256.66393*t) ] / 10000000 弧度

  最后兩行分別為光行差和章動

四、太陽黃經(jīng)速度

平速度：  v0 = 628.3319653318

即時速度：v1 = 628.332 +21 * sin(1.527+628.307585*t)

速度的單位是“弧度/儒略世紀(jì)”即“弧度/36525天”

注意，平速度比即時速度的精度要高得多，務(wù)必保留足夠的有效數(shù)字，否則將帶來嚴(yán)重的計算誤差。

五、節(jié)氣時刻計算

以上天體黃經(jīng)時間的函數(shù)，即L = f(t)，所謂的求節(jié)氣時刻就是已知L求t，顯然這是在求解一個關(guān)于t的方程。偉大的英國天文學(xué)家物理學(xué)家牛頓給出了一種非常有效的迭代算法：牛頓求根法。用這種方法，求t所花費(fèi)的時間僅是求f(t)花費(fèi)時間的1.2——1.3倍。設(shè)某個節(jié)氣對應(yīng)的黃經(jīng)為W，那么算法如下。

牛頓迭代算法設(shè)計：

第1步迭代：t = 0

第2步迭代：t = t + ( W – L0(t) ) / v0

第3步迭代：t = t + ( W – L1(t) ) / v1(t)

第4步迭代：t = t + ( W – L2(t) ) / v1(t)

誤差：算法誤差2分鐘以內(nèi)，實際找到的誤差一般在30秒以內(nèi)，平均15秒

注意：W指的是太陽黃經(jīng)。1999年春分對應(yīng)W=0，以后每W每增加15度對應(yīng)下一個節(jié)氣。迭代的的結(jié)果是力學(xué)時，單位是儒略世紀(jì)數(shù)。最后結(jié)果還應(yīng)轉(zhuǎn)換為北京時間，即：JD = J2000 + t*36525 - △T/86400 + 8/24

最后使用“儒略日數(shù)轉(zhuǎn)公歷”所述方法得到節(jié)氣的日期和時間。

六、計算結(jié)果比較

為了進(jìn)行誤差比較，下文列出2007年的24節(jié)氣，并與《壽星天文歷》比對?！秹坌翘煳臍v》是筆者制作的一款精度優(yōu)于1秒的農(nóng)歷工具，已發(fā)布于互聯(lián)網(wǎng)上，其算法與本文類似。

[參考文獻(xiàn)]

1、Pierre Bretagnon與Gerard Francou 《VSOP87行星運(yùn)動理論》

2、Jean Meeus《Astronomical.Algorithms》

3、NASA下屬JPL實驗室的《DE405/406星歷表》