• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    《皇極歷》外行星算法及精度分析

    2022-08-15 07:52:40辛佳岱
    中國科技史雜志 2022年2期
    關(guān)鍵詞:黃經(jīng)外行星土星

    辛佳岱 唐 泉

    (西北大學(xué)科學(xué)史高等研究院,西安 710127)

    自西漢《太初歷》起,行星理論就是歷法的重要內(nèi)容。隋代之前的歷法家根據(jù)行星和太陽的平運動來推算行星的見、伏等特征天象發(fā)生的時刻及其位置。公元560年前后,張子信發(fā)現(xiàn)太陽與行星運動不均勻現(xiàn)象之后,歷法家不斷探索行星不均勻運動修正算法,以達到提高行星位置計算精度的目的。劉焯的《皇極歷》(600)是現(xiàn)存中國第一部設(shè)計算法修正太陽與行星不均勻運動的歷法,并且影響了隋代及初唐歷法中步五星術(shù)算法的編制,如張胄玄《大業(yè)歷》(607)、傅仁均《戊寅元歷》(619)和李淳風(fēng)《麟德歷》(665)等三部歷法都沿襲了《皇極歷》中五星入氣加減差算法,以修正行星不均勻運動。僧一行在《大衍歷》(724)中,將《皇極歷》等歷法中的五星入氣加減差發(fā)展為五星爻象歷,基本奠定了中國古代行星不均勻運動修正算法的框架,這個算法在楚衍和宋行古等人編制的《崇天歷》(1023)中逐漸趨于成熟,此后一直到元代《授時歷》和明代《大統(tǒng)歷》,幾乎未發(fā)生實質(zhì)性變化。

    現(xiàn)有研究涉及行星運動的理論模型[1,2]、隋代及初唐歷法中的“五星入氣加減差”算法及其精度[3—6]、行星理論中一些重要算法和天文數(shù)表[7—10],以及中印行星理論比較研究[11,12]等方面內(nèi)容。關(guān)于行星計算精度的研究主要集中于宋元歷法[13—17],對明代《大統(tǒng)歷》和《回回歷法》也有相關(guān)考察[18],而對宋代之前諸歷的行星精度較少涉及[19,20]。這對我們?nèi)嬖u價中國古代行星計算精度造成了很大困難。

    在中國古代行星理論發(fā)展史上,隋代劉焯所制《皇極歷》占據(jù)著重要地位,故本文以《皇極歷》行星算法為研究對象,在對算法術(shù)文全面解讀的基礎(chǔ)上,利用Python語言模擬其算法并討論行星視位置計算精度,試圖回答如下問題:張子信發(fā)現(xiàn)太陽視運動與行星公轉(zhuǎn)不均勻現(xiàn)象后,《皇極歷》是如何將這些天文發(fā)現(xiàn)融入行星運動理論中的?《皇極歷》推算行星視位置的誤差如何?

    1 《皇極歷》五星推步思路

    《皇極歷》記載了五星的天文常數(shù)、入氣加減差、動態(tài)描述和五個推步算法([21],頁1963—1970),其行星理論的核心目標是推算行星自始見至始伏期間每日的黃道宿度。隋代歷法中通常不考慮行星伏行段的運行情況,這是因為伏行段的行星位置不易觀測,無法準確給出其運動狀態(tài),當時歷法中有“伏不書度”或“伏不注度”等記載([21],頁1847、1919)。此外,隋代歷法家對內(nèi)行星運動規(guī)律掌握還不夠好,如《皇極歷》中水星的“入氣加減差”僅給出一些經(jīng)驗觀測的描述性言語概述,沒有量化數(shù)據(jù),無法計算任意時刻內(nèi)行星的位置[3]。故本文的討論僅限于外行星,不涉及內(nèi)行星。

    外行星包含五個天文常數(shù)。以土星為例,依次是土數(shù)、伏半平、復(fù)日及余、殘日及余和見去日。土數(shù)(H)意為土星的會合周期,單位為“分”;以“氣日法”(A=46644)除之得復(fù)日及余(h),即將其換算為以日為單位的會合周期;伏半平(F)意為土星伏行時間的一半,單位為“分”;殘日及余是土星會合周期與回歸年日數(shù)之差。見去日則是晨始見時刻土星與太陽的角距(1)《皇極歷》記載木星見去日為14度,火星為16度,土星為16.5度。。

    五個推步算法即“推星平見術(shù)”“求常見日”“求定見日”“求星見所在度”和“求次日”。其中,求常見日涉及修正行星不均勻運動的五星“入氣加減差”;求次日需要借助行星自晨始見至夕始伏的動態(tài)描述(五星動態(tài)表)。這五個算法是為了計算行星特征天象發(fā)生時刻和任意時刻的行星位置而設(shè)計。具體而言,前三個算法推算行星晨始見發(fā)生的時刻,后兩個算法推算在一個會合周期內(nèi)行星自始見至始伏的每日位置。下面對各個算法的名稱、目的及其算法思路進行釋讀。

    (1)推星平見術(shù)。按照太陽和行星的平運動推算行星平見時刻與其前一個冬至?xí)r刻之間的時距。

    據(jù)術(shù)文,假設(shè)Nn為上元時刻距某一年的積年數(shù),t為回歸年的日數(shù),h為行星會合周期的日數(shù),f為行星伏行一半的日數(shù),則以會合周期推算行星平見時刻與其前一個冬至?xí)r刻的時距γ0,由下式得到:

    tNn-f≡h-γ0(modh)

    (1)

    《皇極歷》定義上元O點起于五星與太陽同度的合時刻,而該算法計算平見時刻距離其前一個冬至?xí)r刻的時距γ0。因此,從上元至所求年冬至?xí)r刻的積日數(shù)及余tNn中,減去行星伏行時間一半的日數(shù)f,余為A0P,即將起算點從上元行星的合時刻O點推后到上元之后的行星始見時刻A0點。再從A0P中減去整數(shù)倍的行星會合周期日數(shù)h,不足部分為AnP,以行星會合周期AnAn+1減之,所得為PAn+1,即術(shù)文所求γ0。如圖1所示。

    圖1 《皇極歷》行星平見時刻與所求年天正冬至?xí)r刻的時距

    (2)求常見日??紤]行星中心差,對平見時刻進行修正,得到修正后的行星常見時刻距離其前一個冬至?xí)r刻的時距。

    如圖2,根據(jù)公式(1)推算出行星平見時刻與其前一個冬至?xí)r刻的時距γ0,該算法求經(jīng)行星中心差修正后的行星常見時刻距離其前一個冬至?xí)r刻的時距γ1,由下式得到:

    (2)

    圖2 《皇極歷》行星常見時刻與所求年天正冬至?xí)r刻的時距

    平見時刻與常見時刻之間的修正值Δm,其數(shù)值及符號取決于行星平見時刻所入節(jié)氣。換言之,由于修正行星中心差的結(jié)果取決于“平見”所值節(jié)氣,故稱為“入氣加減差”算法。

    在隋代及初唐的一些歷法中,“入氣加減差”是為了修正行星中心差對五星始見時刻的影響而設(shè)計的算法。雖然五星“入氣加減差”算法帶有一定的經(jīng)驗色彩,但對于外行星運動不均勻性的改正是有效的、科學(xué)的[3]。《皇極歷》木星和土星的“入氣加減差”算法前人已解釋清楚[3—6],火星的“入氣加減差”術(shù)文記載如下:

    平見,在雨水前,以十九乘去大寒日;清明前,又十八乘去雨水日,增雨水所乘者;夏至后,以十六乘去處暑日;小滿后,又十五日;寒露前,以十八乘去白露日;小雪前,又十七乘去寒露所乘者(2)嚴敦杰承擔(dān)??薄端鍟ぢ蓺v志》工作時,在《皇極歷》術(shù)文后的“??庇洝盵二三]指出:“又十七乘去寒露所乘者 當作[又十七乘去寒露日,增寒露所乘者]?!眳⒁妳⒖嘉墨I[21],第1973頁。;大雪后,二十九乘去大寒日,為減,小雪至大雪減二十五日。([21],頁1964)

    經(jīng)分析,這段文字應(yīng)該存有脫訛,陳美東也指出這一點(3)陳美東的校勘是:“……清明前,又十八乘去雨水日,增雨水所乘者;[清明至夏至加二十七日;]夏至后,以十六乘去處暑日;(小滿)[處暑]后,又(十五日)[二十八乘去白露日。減處暑所乘者];寒露前,……”參見參考文獻[3]。,然而筆者按照??备膭幼钌僭瓌t,對術(shù)文??比缦拢?/p>

    平見,在雨水前,以十九乘去大寒日;清明前,又十八乘去雨水日,增雨水所乘者;[清明至小滿加二十七日;]夏至后,以十六乘去處暑日;小滿后,又十五[乘去夏至]日;寒露前,以十八乘去白露日;小雪前,又十七乘去寒露日,增寒露所乘者;大雪后,二十九乘去大寒日,為減,小雪至大雪減二十五日。

    上述文字表明,若依據(jù)火星會合周期推算的始見時刻“平見”所入節(jié)氣不同,則火星運動不均勻性影響的改正值就會改變?;鹦请S著“平見”所入節(jié)氣而引起時間提前或是滯后的改正曲線,如圖3所示。

    圖3 《皇極歷》火星“平見日”至“常見日”的改正值曲線

    圖3中橫坐標為黃經(jīng)值,春分對應(yīng)黃經(jīng)0°,每經(jīng)過一氣黃經(jīng)值的改變量為15°,據(jù)此可換算到其余對應(yīng)的節(jié)氣??v坐標為時間改正值,以“日”為單位,正值表示經(jīng)修正行星不均勻運動得到的“常見日”較“平見日”滯后(如圖2中的B″點),負值則表示提前(如圖2中的B′點)。

    (3)求定見日??紤]太陽中心差,對常見時刻進行修正,得到修正后的行星定見時刻距離其前一個冬至?xí)r刻的時距。

    與“求常見日”類似,此處要對γ1進行太陽中心差修正,求修正后的行星定見時刻距離其前一個冬至?xí)r刻的時距γ2,由下式得到:

    (3)

    常見時刻與定見時刻之間的修正值為太陽中心差Cs,該算法在《皇極歷》步日躔術(shù)中記載([21],頁1937—1939)。其中,“先后數(shù)”即太陽中心差,以“先后數(shù)”除以“轉(zhuǎn)法”(52)得到真太陽與平太陽的行度之差。日躔表中羅列二十四氣各氣初始時刻的太陽先后數(shù),劉焯設(shè)計了等間距二次內(nèi)插算法以求每日先后數(shù)([22],頁152—155)。

    (4)求星見所在度。由于外行星運行速度較太陽慢,故以定見時刻太陽的位置減去“見去日”,即可得到定見時刻行星的位置。

    術(shù)文中記載定見時刻太陽位置的方法是:計算行星定見時刻所在晨前夜半的太陽位置,再加上定見日余對應(yīng)的太陽行度即可([21],頁1969—1970)。為了便于利用Python語言模擬其算法,此處以一種等價方法替換。根據(jù)《皇極歷》步月離中的“推日度術(shù)”([21],頁1950—1951),計算所求年冬至?xí)r刻的太陽位置,再由公式(3)得到行星定見時刻距離所求年冬至?xí)r刻的時距γ2以步日躔中的太陽中心差算法推算出在γ2內(nèi)太陽運行的度數(shù),即可得到定見時刻的太陽位置。

    以土星定見時刻的位置為例,若令λ0為太陽在所求年冬至?xí)r刻的黃道宿度,d為在γ2內(nèi)太陽的實際運行度數(shù),并且根據(jù)土星“見去日”常數(shù)得知定見時刻土星落后太陽16.5度。故土星定見時刻的黃道宿度λ土,由下式得到:

    λ土=λ0+d-16.5

    (4)

    (5)求次日。計算行星在任意一日的位置。由公式(4)推算出行星定見時刻的位置,在此基礎(chǔ)上累加行星各日行度,即可得到行星定見之后至始伏期間每一日的位置。

    《皇極歷》記載了行星在一個會合周期中自始見至始伏的逐日行度。根據(jù)土星的動態(tài)描述([21],頁1966),羅列出土星自晨始見至夕始伏的運動狀態(tài),如表1所示。

    表1 土星動態(tài)表

    由表1可見,土星運動狀態(tài)劃分為五個段目:前順行、前留、逆行、后留和后順行(4)前順行段目以土星晨始見為起點,后順行段目以土星夕始伏為結(jié)束?!痘蕵O歷》術(shù)文中并未出現(xiàn)“前順行、前留”等這些名稱,是我們?yōu)榱朔奖忝枋?,依?jù)唐代歷法五星動態(tài)表中的名稱,結(jié)合《皇極歷》土星動態(tài)描述給出相應(yīng)段目的名稱。,并未給出伏行段目行星的動態(tài)描述,僅以“伏半平”和“見去日”兩常數(shù)籠統(tǒng)刻畫行星伏行的時間和度數(shù)。木星的運行狀態(tài)與土星類似,相對容易理解。

    《皇極歷》中火星的運動劃分為七個段目:前疾、前遲、前留、逆行、后留、后遲和后疾?;鹦窃谇凹才c后疾段目的動態(tài)描述與以往歷法大有不同,表現(xiàn)為這兩段目的火星運動狀態(tài)并非直接給出,而是要先判斷火星前疾和后疾初日距離其前一個冬至的日數(shù),以此推算火星在該段目上的運行時間(日率)和度數(shù)(度率)等([21],頁1964—1965)。這種處理方式與印度6世紀中葉《五大歷書匯編》(Pacasiddhāntikā)中的《寶莉莎歷數(shù)書》(Paulia)將火星動態(tài)表與行星的黃經(jīng)相聯(lián)系的設(shè)計方案非常類似[23]。

    火星在前疾與后疾兩段目動態(tài)描述首先根據(jù)火星晨始見時刻或后疾初始時刻所值節(jié)氣,對應(yīng)算出火星在前疾段或后疾段的日率和度率,之后在某些節(jié)氣上,還要對日率和度率進行二次修正。其中,對某些節(jié)氣的日率、度率二次修正時,有“半度之行”的說法,如術(shù)文言“白露至寒露,初日行半度,四十日行二十度”([21],頁1964)。以往學(xué)者的理解是火星晨見初日運動速度均為半度(5)張培瑜等學(xué)者認為白露至寒露時,火星晨見初日運動速度均為半度,見參考文獻[24]。劉洪濤認為火星在白露至寒露間初見的平均速度是日行半度,見后勻速行40日移20度,見參考文獻[25]。王應(yīng)偉認為術(shù)文疑有脫字或脫句,見參考文獻[26]。。唐泉依據(jù)初唐《戊寅元歷》《麟德歷》中對火星前疾段目的修正法則時的論述,認為此處“半度之行”,并非火星晨始見在白露至寒露這個時段內(nèi),其初日視行速為0.5度/日,而是若在這個時段內(nèi)要對原有日率和度率進行二次修正,即日率減去40日,度率減去20度[19]。

    中國傳統(tǒng)歷法家計算行星位置的基本思路在《皇極歷》中已雛形初現(xiàn)。概言之,其五星推步的思路是:首先,依據(jù)行星會合周期推算平見日;其次,以行星中心差修正平見日,得到常見日,再以太陽中心差對常見日加以修正,得到定見日;再次,推算定見時刻的行星位置;最后,在此基礎(chǔ)上再累加五星動態(tài)表中的每日行度,得到行星的每日位置。

    2 《皇極歷》行星算法的計算機模擬與算例

    為了加深對《皇極歷》五星推步算法的理解,此處以公元600年冬至之后的土星在一個會合周期內(nèi)自晨始見至夕始伏每日夜半的極黃經(jīng)推算過程為例,給出計算實例。

    《皇極歷》中與土星位置相關(guān)的基本天文常數(shù)包含:上元距開皇二十年(600)積年數(shù)為N600=1008837年(6)此處的積年數(shù)包含所求年,歷法中通常稱為“算上”。。歲數(shù)T=17036466.5,土數(shù)H=17635594,伏半平F=864995,它們均以氣日法A=46644為分母,由此得到:

    (5)

    參照上一節(jié)《皇極歷》計算行星位置的基本思路,為了得到土星每日夜半的歷取黃經(jīng)值,根據(jù)五星推步算法過程將其分為五個步驟執(zhí)行。具體如下:

    第1步:求公元600年冬至?xí)r刻距離后一個土星平見時刻的時距γ0。

    依據(jù)“推星平見術(shù)”,將上元積年N600、回歸年t、土星會合周期h等常數(shù)代入公式(1),即

    由此得到γ0約是338.2460日。

    第2步:求公元600年冬至?xí)r刻距離后一個土星常見時刻的時距γ1。

    依據(jù)“求常見日”,可知平見與常見時刻之間的修正值取決于行星平見時刻所入節(jié)氣?!痘蕵O歷》平氣日數(shù)約為15.2185日,由第1步所得γ0推算土星平見時刻入于冬至后的第22氣,即入于小雪后約第3.4385日。按照土星入氣加減差修正并代入公式(2),即

    由此得到γ1約是337.6508日。

    第3步:求公元600年冬至?xí)r刻距離后一個土星定見時刻的時距γ2,以及定見日期。

    依據(jù)“求定見日”,可知常見與定見時刻之間的修正值取決于行星常見時刻所入節(jié)氣。土星常見時刻入冬至后第22氣約第6.7417日,根據(jù)步日躔術(shù),計算出由于太陽運動不均勻性而修正土星常見時刻的數(shù)值約是0.7765日。因此,結(jié)合公式(3),即

    γ2=337.6508+0.7765(日)。

    由此得到γ2約是338.4274日。

    推算土星定見日期,需要得知公元600年歷取冬至?xí)r刻的日期。依據(jù)“推氣術(shù)”算法([21],頁1936—1937),以如下公式求得歷取冬至?xí)r刻:

    tN600≡r(mod 60)

    (6)

    代入常數(shù),得到冬至大小余r:

    《皇極歷》的上元起于甲子日,由上式得到冬至大余為29,可知公元600年歷取冬至日的日名為“癸巳”。這一年理論冬至?xí)r刻是600年12月18日21時36分,日名也是“癸巳”[27],將“21時36分”化為以日為單位是0.9日。由日名可知恰與歷取冬至在同一日,這一日的儒略日序號是1940560。上式算得歷取冬至不足一日的小余為0.8118日。因此,這一年的歷取冬至?xí)r刻是600年12月18日19時29分。

    γ2不足整日數(shù)的部分為0.4274日,它與歷取冬至小余0.8118日兩者之和大于1日,則土星定見時刻與其前一個冬至所在日夜半的時距約是339.2391日。由此推算出土星定見日期是601年11月22日。

    第4步:求公元600年冬至?xí)r刻之后的土星定見所在日夜半的黃經(jīng)。

    (1)求公元600年歷取冬至?xí)r刻的太陽黃經(jīng)

    方法是在理論冬至?xí)r刻的太陽黃經(jīng)270°的基礎(chǔ)上,增減在歷取與理論冬至?xí)r刻差期間太陽運行的度數(shù)?!痘蕵O歷》依據(jù)日行盈縮將回歸年以“二分”為界分為兩部分,“秋分后春分前為盈汎,春分后秋分前為虧總”([21],頁1936)。據(jù)術(shù)文得知冬至每氣日數(shù)是:

    根據(jù)步日躔術(shù),得到冬至這一日太陽的中心差值為:

    因此,公元600年歷取冬至?xí)r刻的太陽黃經(jīng)是:

    (2)求公元600年冬至?xí)r刻之后的土星定見時刻的黃經(jīng)

    依據(jù)“求星見所在度”,得知關(guān)鍵步驟是求太陽在γ2內(nèi)實際運行度數(shù)d。定見時刻入冬至后第22氣約第7.5183日,根據(jù)步日躔術(shù),得到在γ2內(nèi)太陽實行度與平行度之差約為0.7518度。由于太陽平行度為每日1度,所以太陽從冬至?xí)r刻至定見時刻運行的度數(shù)為:

    d=γ2-0.7518=338.4274-0.7518=337.6755(度)。

    將其結(jié)果代入公式(4),即

    由此得知土星定見時刻的黃經(jīng)是226.4620°。

    (3)求公元600年冬至?xí)r刻之后的土星定見所在日夜半的黃經(jīng)

    (7)

    根據(jù)第3步所得土星定見時刻與其前一個冬至所在日夜半時距約是339.2391日,可知土星定見小余r′約為0.2391日,且由表1得到土星的初日速度v0為4364/46644(度/日)。將其代入(7)式得到:

    由此得知土星定見所在日夜半的黃經(jīng)是226.4400°。

    第5步:求公元600年冬至?xí)r刻之后的土星定見所在日至夕始伏每日夜半的黃經(jīng)。

    在第4步計算結(jié)果基礎(chǔ)上,累加表1土星動態(tài)表中的自晨始見至夕始伏土星的每日行度,得到定見之后每一日夜半的土星黃經(jīng)值,其結(jié)果如下文表2中的“歷取黃經(jīng)”一欄所示。

    3 《皇極歷》外行星計算精度及分析

    3.1 《皇極歷》外行星位置計算精度結(jié)果

    上一節(jié)中我們以土星為例,利用《皇極歷》五星推步方法已得到公元600年冬至之后的土星一個會合周期自晨始見至夕始伏每日夜半的歷取極黃經(jīng),所以將其與理論黃經(jīng)求差值的做法(歷取黃經(jīng)-理論黃經(jīng)),即可算得土星這個周期的黃經(jīng)絕對誤差(8)行星視位置的極黃經(jīng)和真黃經(jīng)的確有一些差別,但由于外行星運行軌道與黃道的夾角較小,因此這兩者的差別就非常小,且與行星理論的系統(tǒng)誤差比較,外行星的這些誤差目前可以暫時忽略。參見參考文獻[1]。。

    土星的理論黃經(jīng)結(jié)果可以通過天文軟件Skymap提取并計算得出,寧曉玉的研究表明Skymap軟件滿足古天文研究數(shù)據(jù)的精確度和穩(wěn)定性,是準確可靠的天文軟件[28]。所提取的理論數(shù)據(jù)是土星自601年11月22日定見之后341日的每日夜半時刻赤經(jīng)和赤緯值,需要將其換算到對應(yīng)的黃經(jīng)值,其結(jié)果如表2中的“理論黃經(jīng)”一欄所示。

    表2 《皇極歷》以601年11月22日為土星定見的會合周期晨始見至夕始伏位置精度

    上一節(jié)已得到公元600年冬至之后的土星定見日期是601年11月22日,其儒略日序號數(shù)是1940899。這里限于篇幅僅給出該會合周期自晨始見至夕始伏每間隔20日的土星黃經(jīng)數(shù)值。將《皇極歷》五星推步算法得到的歷取黃經(jīng)與理論黃經(jīng)求差值,以推算土星這個周期的黃經(jīng)絕對誤差,如表2所示。

    根據(jù)表2最后一欄“黃經(jīng)誤差”,以每日結(jié)果繪制土星黃經(jīng)絕對誤差散點圖,如圖4所示。經(jīng)統(tǒng)計分析,土星數(shù)據(jù)量有341個,黃經(jīng)誤差絕對值的平均值為1.46°,誤差絕對值的最大值為2.51°。

    圖4 601年11月22日定見至始伏土星黃經(jīng)絕對誤差散點圖

    以同樣方式分別計算出木星開皇二十年冬至后定見發(fā)生于601年6月2日,火星定見發(fā)生于601年11月29日。分別繪制出木星和火星在這個會合周期內(nèi)自晨始見至夕始伏期間的黃經(jīng)絕對誤差散點圖,如圖5和圖6所示。

    圖5 601年6月2日定見至始伏木星黃經(jīng)絕對誤差散點圖

    圖6 601年11月29日定見至始伏火星黃經(jīng)絕對誤差散點圖

    經(jīng)統(tǒng)計,木星數(shù)據(jù)量有363個,黃經(jīng)誤差絕對值的平均值為1.83°,誤差絕對值的最大值為3.12°;火星數(shù)據(jù)量有624個,黃經(jīng)誤差絕對值的平均值為4.21°,誤差絕對值的最大值為8.67°。

    為了更深入全面探討《皇極歷》外行星計算精度,接下來將推算《皇極歷》制成之后約30年長時段內(nèi)木星、火星和土星在每日夜半的黃經(jīng),并將其與行星理論黃經(jīng)比較,以得到絕對誤差結(jié)果。對于木星,計算了開皇二十年冬至后的28個完整會合周期自晨始見至夕始伏的每日夜半位置,數(shù)據(jù)量為10165個;對于火星,計算了15個完整會合周期的每日夜半位置,數(shù)據(jù)量為9204個;對于土星,計算了29個完整會合周期的每日夜半位置,數(shù)據(jù)量為10231個。

    根據(jù)計算結(jié)果分別繪制出木星、火星和土星的黃經(jīng)絕對誤差散點圖,如圖7—9所示。圖中的縱坐標表示行星黃經(jīng)絕對誤差,單位為“°”;橫坐標表示從公元600年冬至之后外行星首次定見所在日的夜半起算的日數(shù),單位為“日”。圖7中的橫坐標始于601年6月2日夜半,圖8始于601年11月29日夜半,圖9始于601年11月22日夜半。

    圖7 《皇極歷》木星黃經(jīng)絕對誤差散點圖

    圖8 《皇極歷》火星黃經(jīng)絕對誤差散點圖

    圖9 《皇極歷》土星黃經(jīng)絕對誤差散點圖

    從圖7—9可以看出,外行星的黃經(jīng)絕對誤差在約30年內(nèi)呈現(xiàn)出周期性的變化。然而,在此之中火星有兩個會合周期的誤差比較反常。通過查驗,這兩個周期的晨始見分別處于白露、秋分兩氣,處于這兩氣的修正方案是“白露至寒露,初日行半度,四十日行二十度”。已在第1章中論述火星動態(tài)描述時,指出此處“半度之行”是對火星晨始見位于“白露至寒露”段內(nèi)的二次修正,日率減去40日,度率減去20度。若不考慮術(shù)文中的“初日行半度”的修正,可繪制出火星黃經(jīng)絕對誤差散點圖,如圖10所示。

    圖10 《皇極歷》火星黃經(jīng)絕對誤差散點圖(不修正“初日行半度”)

    如圖10,若不考慮火星晨始見位于“白露至寒露”段內(nèi)“初日行半度”的修正,這兩個反常周期(圖8)的誤差會大幅度減小,且與前后會合周期的誤差銜接度更好。

    因此,根據(jù)以上分析,《皇極歷》外行星計算精度可由表3中的相關(guān)數(shù)據(jù)反映。

    表3 《皇極歷》601—630年外行星計算精度(9)表3中火星計算精度結(jié)果均為不考慮“初日行半度”情況下所得。

    印度天文學(xué)著作《五大歷書匯編》成書于6世紀中葉,其中的《寶莉莎歷數(shù)書》與中國傳統(tǒng)歷法處理行星運動的方式相同,即以代數(shù)方法計算行星的位置。唐泉對《寶莉莎歷數(shù)書》中的行星運動作了深入研究并給出外行星及金星的計算精度[11],外行星計算精度數(shù)據(jù)摘錄于表4。通過比較,《皇極歷》中推算外行星黃經(jīng)誤差精度相比《寶莉莎歷數(shù)書》要稍微高一些。

    表4 《寶莉莎歷數(shù)書》行星計算精度

    3.2 影響《皇極歷》外行星位置精度的原因分析

    根據(jù)上述對《皇極歷》行星推步算法的解讀分析,不難看出,定見時刻被視為行星在一個會合周期上的關(guān)節(jié)點,再由太陽位置及見去日推算行星定見時刻的位置,最后以行星每日的行度累加,即可得到行星任意一日的位置??梢?,行星定見時刻的推算,以及行星在一個會合周期內(nèi)自晨始見至夕始伏的每日行度都有可能影響到《皇極歷》外行星位置計算精度。接下來,從這兩方面展開探討。

    一方面,有必要對行星在約30年時段內(nèi)每次定見時刻的精度進行分析。依據(jù)《皇極歷》行星推步算法得到開皇二十年冬至后約30年內(nèi)木星、火星和土星各完整會合周期的歷取定見日期,以及由中國黃道星空天文軟件所得理論定見日期(10)“中國黃道星空”是劉次沅專門為研究中國古代天象記錄而開發(fā)的一款天文軟件,它可以計算并演示日、月、五星在黃道天區(qū)的位置和運動。此處的“理論定見日期”就是通過《皇極歷》所給外行星“見去日”常數(shù),進而借助該軟件推算這一天象發(fā)生的日期。,如表5所示。

    表5 《皇極歷》601—630年木星、火星和土星定見日期

    表5中結(jié)果表明木星歷取定見日期與理論定見日期的最大誤差是6日,平均誤差約為2日;火星的最大誤差是41日,平均誤差約為21日;土星的最大誤差是5日,平均誤差約為3日??梢?,《皇極歷》推算外行星定見時刻的計算誤差中火星最大,木星和土星基本相當。

    另一方面,我們對外行星在一個會合周期內(nèi)自晨始見至夕始伏的速度進行分析。根據(jù)表2的計算結(jié)果,可繪制以600年11月22日為土星定見的一個會合周期內(nèi)的理論速度與歷取速度,即呈現(xiàn)出開皇二十年冬至之后土星在一個會合周期內(nèi)的運動狀態(tài)。如圖11所示。

    圖11 土星理論速度與歷取速度圖示

    《皇極歷》中的土星在一個周期內(nèi)的運動狀態(tài)被分為前順行、前留、逆行、后留和后順行這五個段目,且順行與逆行的歷取速度均以勻速運動給出,如圖11所示。對比理論速度曲線可以看出,土星的理論速度實際在每時每刻都是改變的。歷取速度劃分的段目較少,很難達到與理論速度較好擬合。此外,歷取速度在前留、后留這兩段目各占39日,而實際理論速度中“留”卻僅是存有一瞬間的狀態(tài)。因此,《皇極歷》對土星一個周期內(nèi)各段目的運行狀態(tài)描述也是影響土星計算精度誤差的一個重要因素。

    以同樣的方式,繪制出開皇二十年冬至之后木星和火星在一個會合周期內(nèi)的理論速度與歷取速度,如圖12和圖13所示。

    圖12 木星理論速度與歷取速度圖示

    圖13 火星理論速度與歷取速度圖示

    值得注意的是,與木星和土星不同,圖13中所呈現(xiàn)的火星在一個會合周期內(nèi)的理論速度與歷取速度并不對稱,換言之,其晨始見的速度與該會合周期夕始伏的速度不相等。從理論速度曲線上同樣反映出這一現(xiàn)象,且歷取速度與之擬合較好。如圖14所示,我們給出更大的時間范圍內(nèi)速度擬合曲線,即約公元600—630年火星15個會合周期內(nèi)晨始見至夕始伏期間運行理論速度與歷取速度圖示。

    圖14 《皇極歷》火星公元601年11月29日夜半后每日理論速度與歷取速度圖示

    《皇極歷》的火星歷取速度之所以能達到這樣好的效果,其原因在于它設(shè)計火星前疾段和后疾段運動狀態(tài)時的巧妙算法,以一個會合周期中前疾初日與后疾初日所值節(jié)氣,來分別規(guī)定火星在這兩個段目運行的日率和度率。

    4 結(jié)論

    本文通過對《皇極歷》外行星入氣加減差、動態(tài)描述及五星推步算法構(gòu)建思路的解讀,并結(jié)合上述分析討論,得出下面幾點結(jié)論:

    (1)根據(jù)表3中的數(shù)據(jù),得知對《皇極歷》外行星而言,木星和土星黃經(jīng)誤差絕對值的最大值分別是5.27°和5.10°,誤差絕對值的平均值分別是1.94°和2.61°。這表明《皇極歷》中木星和土星的計算精度大體相當,而火星黃經(jīng)計算誤差要比木星和土星大得多?;鹦屈S經(jīng)誤差絕對值的最大值是17.67°,約是木星和土星黃經(jīng)最大誤差的3倍還要多,誤差絕對值的平均值也達到4.80°。此外,通過對外行星約30年中的每個定見日期的誤差分析,如表5,同樣得到火星的定見日期誤差最大,木星和土星誤差相當且較火星小很多。這是因為火星在三個外行星中,它的軌道偏心率最大、周期長、速度快,因而火星的計算相對木星和土星要更困難一些。

    (2)《皇極歷》討論火星“前疾”段的日率和度率時,附加修正法則所提到的“半度之行”,這一項修正對火星黃經(jīng)計算精度的提高似乎沒有發(fā)揮正面作用,如圖8和圖10所示。引入“半度之行”修正的火星黃經(jīng)誤差絕對值的最大值高達31.49°,誤差絕對值的平均值是6.65°;而未考慮“半度之行”修正的火星黃經(jīng)誤差絕對值的最大值是17.67°,誤差絕對值的平均值是4.80°??煽闯鰪木确矫婵疾欤@個算法顯得多余。但是從初唐《戊寅元歷》《麟德歷》這兩部歷法對“半度之行”的論述分析而言,其涵義就是對原有日率和度率的二次修正,理解上無誤。因而關(guān)于“半度之行”的合理性問題還需要作進一步探討。

    (3)《皇極歷》外行星黃經(jīng)絕對誤差散點圖呈現(xiàn)周期性變化的規(guī)律,與其公轉(zhuǎn)周期基本吻合,如圖7—10所示。文章中所討論自公元600年冬至之后首次始見的當日夜半起算約12000日時間段內(nèi),大約包含了2.7個木星公轉(zhuǎn)周期、1.1個土星公轉(zhuǎn)周期、17.5個火星公轉(zhuǎn)周期。木星和土星黃經(jīng)絕對誤差的變化規(guī)律與它們的公轉(zhuǎn)周期比較吻合,而火星由于其公轉(zhuǎn)周期較短,其黃經(jīng)絕對誤差變化不太容易與公轉(zhuǎn)周期對應(yīng)。然而,火星呈現(xiàn)出七個會合周期循環(huán)規(guī)律性的波動趨勢,這與火星的大沖周期更為匹配。

    (4)由于目前隋代之前的歷法行星計算精度的研究甚少,本人所在項目組的其他成員提供了一些未發(fā)表的魏晉南北朝時期歷法中外行星誤差計算精度的初步結(jié)果,得知這時期木星和土星的黃經(jīng)誤差絕對值的最大值達到十幾二十多度,火星更是有高達四五十度的誤差??梢姡痘蕵O歷》在引入太陽和行星不均勻運動修正算法之后,確實在很大程度上改善了其行星計算水平。雖然隋代及唐初時期歷法中的“入氣加減差”和“五星動態(tài)表”均屬于基于經(jīng)驗性觀測的數(shù)值算法系統(tǒng),但是這項研究也進一步表明古代歷法家在探索行星運動理論、設(shè)計數(shù)值算法逐漸逼近真值的過程中,不懈努力的精神是值得肯定和學(xué)習(xí)的。

    致 謝感謝兩位匿名審稿專家提出寶貴的修改建議。論文寫作過程中還得到了北京天文館古觀象臺楊帆老師提供的材料和幫助,在此表示衷心感謝!

    猜你喜歡
    黃經(jīng)外行星土星
    外行星
    土星之旅
    春分
    土星環(huán),你要去哪里
    大寒
    廈門航空(2019年1期)2019-12-17 14:02:06
    系外行星那些事——“呼啦圈”法
    土星
    系外行星探索與發(fā)現(xiàn)
    太空探索(2016年4期)2016-07-12 15:17:52
    系外行星的探測
    太空探索(2016年3期)2016-07-12 09:58:42
    按陽歷算的清明節(jié)
    久久97久久精品| 十分钟在线观看高清视频www | 一区二区av电影网| 高清av免费在线| 国产亚洲最大av| 好男人视频免费观看在线| 中文乱码字字幕精品一区二区三区| 插阴视频在线观看视频| 美女内射精品一级片tv| 中文字幕精品免费在线观看视频 | 欧美日韩精品成人综合77777| 中国三级夫妇交换| 天美传媒精品一区二区| 曰老女人黄片| 亚洲图色成人| 国产欧美日韩一区二区三区在线 | 丝瓜视频免费看黄片| 丝瓜视频免费看黄片| 九九爱精品视频在线观看| 国产日韩欧美亚洲二区| 午夜影院在线不卡| 精品亚洲成a人片在线观看| 两个人免费观看高清视频 | 七月丁香在线播放| 激情五月婷婷亚洲| 国产精品一区二区在线不卡| 卡戴珊不雅视频在线播放| 乱系列少妇在线播放| 国产精品伦人一区二区| 国产精品伦人一区二区| 97在线人人人人妻| 国产探花极品一区二区| 精品国产一区二区三区久久久樱花| 插逼视频在线观看| 深夜a级毛片| 精品人妻偷拍中文字幕| 亚洲第一区二区三区不卡| 日本午夜av视频| 色哟哟·www| 我要看黄色一级片免费的| 18禁动态无遮挡网站| 亚洲av成人精品一区久久| 精品人妻熟女毛片av久久网站| 全区人妻精品视频| 在线天堂最新版资源| 两个人的视频大全免费| 蜜臀久久99精品久久宅男| 欧美精品一区二区免费开放| 最新中文字幕久久久久| 久久99精品国语久久久| 精品久久久久久久久亚洲| 七月丁香在线播放| 亚洲欧洲精品一区二区精品久久久 | 日本午夜av视频| av免费观看日本| 久久精品国产自在天天线| 视频区图区小说| 男的添女的下面高潮视频| 久久久欧美国产精品| 99久久精品一区二区三区| 99久久人妻综合| 麻豆精品久久久久久蜜桃| 亚洲va在线va天堂va国产| 国产免费一区二区三区四区乱码| av国产久精品久网站免费入址| 日本爱情动作片www.在线观看| 国产精品女同一区二区软件| 免费久久久久久久精品成人欧美视频 | 日本vs欧美在线观看视频 | 亚洲综合精品二区| 亚洲va在线va天堂va国产| 99九九在线精品视频 | 五月开心婷婷网| 国产av精品麻豆| av有码第一页| 美女脱内裤让男人舔精品视频| 成人毛片60女人毛片免费| 久久人妻熟女aⅴ| 精品熟女少妇av免费看| 午夜福利网站1000一区二区三区| 亚洲av日韩在线播放| 寂寞人妻少妇视频99o| 交换朋友夫妻互换小说| 免费人妻精品一区二区三区视频| 色婷婷av一区二区三区视频| 亚洲精品乱久久久久久| 色哟哟·www| 久久97久久精品| 99久久精品国产国产毛片| 亚洲精品第二区| 51国产日韩欧美| 国精品久久久久久国模美| 国产亚洲欧美精品永久| 国产极品粉嫩免费观看在线 | 精品亚洲成国产av| av线在线观看网站| 观看免费一级毛片| 在线观看人妻少妇| 久久久久久久亚洲中文字幕| 久久久久久久亚洲中文字幕| 国产午夜精品一二区理论片| 人妻一区二区av| 午夜福利,免费看| 久久国产精品大桥未久av | 日本爱情动作片www.在线观看| 一级毛片久久久久久久久女| 下体分泌物呈黄色| 极品少妇高潮喷水抽搐| 大又大粗又爽又黄少妇毛片口| 大又大粗又爽又黄少妇毛片口| 久久精品国产鲁丝片午夜精品| 免费黄网站久久成人精品| 美女视频免费永久观看网站| 亚洲精品中文字幕在线视频 | 九九爱精品视频在线观看| 欧美成人午夜免费资源| 久久久久久久亚洲中文字幕| 蜜桃在线观看..| 欧美老熟妇乱子伦牲交| 午夜免费男女啪啪视频观看| 日韩电影二区| 欧美精品高潮呻吟av久久| 99热这里只有精品一区| 国产乱来视频区| 最近2019中文字幕mv第一页| 春色校园在线视频观看| 一区二区三区免费毛片| 国产黄片视频在线免费观看| xxx大片免费视频| 国产高清国产精品国产三级| 亚洲欧洲精品一区二区精品久久久 | 日韩av不卡免费在线播放| 天堂中文最新版在线下载| 国产免费一级a男人的天堂| 亚洲av欧美aⅴ国产| 久久99蜜桃精品久久| 最黄视频免费看| 久久精品国产自在天天线| 欧美日韩精品成人综合77777| 亚洲欧美精品专区久久| freevideosex欧美| 久久热精品热| 亚洲精品aⅴ在线观看| 国产免费福利视频在线观看| 欧美区成人在线视频| 黄色日韩在线| av有码第一页| h日本视频在线播放| 免费大片黄手机在线观看| 午夜福利,免费看| 十八禁网站网址无遮挡 | 一级毛片久久久久久久久女| 久久99热这里只频精品6学生| 久久久久网色| 丝袜在线中文字幕| 国产精品蜜桃在线观看| 香蕉精品网在线| 成年av动漫网址| 老女人水多毛片| 亚洲国产毛片av蜜桃av| 青春草视频在线免费观看| 日韩一区二区视频免费看| 国产av精品麻豆| 成人毛片a级毛片在线播放| 天堂中文最新版在线下载| 欧美精品人与动牲交sv欧美| 亚洲精品自拍成人| 99热网站在线观看| 又粗又硬又长又爽又黄的视频| 青春草亚洲视频在线观看| 国产黄片视频在线免费观看| 国产美女午夜福利| 亚洲国产成人一精品久久久| av专区在线播放| 亚洲av二区三区四区| 国产免费一级a男人的天堂| 如日韩欧美国产精品一区二区三区 | 狂野欧美白嫩少妇大欣赏| 欧美日韩亚洲高清精品| 日韩一区二区视频免费看| 国产午夜精品一二区理论片| 欧美激情国产日韩精品一区| 国产精品一区二区在线观看99| 午夜福利网站1000一区二区三区| 国产精品偷伦视频观看了| 日韩制服骚丝袜av| 一边亲一边摸免费视频| 亚洲欧洲日产国产| 久久综合国产亚洲精品| 观看美女的网站| 日本av手机在线免费观看| 久久韩国三级中文字幕| av不卡在线播放| 国产成人aa在线观看| 极品教师在线视频| 成人美女网站在线观看视频| 少妇的逼水好多| 熟女电影av网| 日韩人妻高清精品专区| 精品国产乱码久久久久久小说| 99久久中文字幕三级久久日本| 有码 亚洲区| 啦啦啦中文免费视频观看日本| 婷婷色综合www| 妹子高潮喷水视频| 国产日韩欧美亚洲二区| 夫妻性生交免费视频一级片| 街头女战士在线观看网站| 亚洲自偷自拍三级| 久久久久人妻精品一区果冻| 在线观看美女被高潮喷水网站| 毛片一级片免费看久久久久| 麻豆乱淫一区二区| 久久国产精品男人的天堂亚洲| 欧美激情极品国产一区二区三区| 又大又爽又粗| 亚洲综合色网址| 天天躁日日躁夜夜躁夜夜| 黑人巨大精品欧美一区二区mp4| 伊人久久大香线蕉亚洲五| 久久久欧美国产精品| 国产精品av久久久久免费| 日本av免费视频播放| 国产亚洲午夜精品一区二区久久| 99久久99久久久精品蜜桃| 丝袜在线中文字幕| 久久天堂一区二区三区四区| 日韩欧美免费精品| 如日韩欧美国产精品一区二区三区| 精品第一国产精品| 老司机影院成人| 无遮挡黄片免费观看| 美女中出高潮动态图| 久久中文字幕一级| 国产精品成人在线| 中文字幕人妻熟女乱码| 免费久久久久久久精品成人欧美视频| 国产成人免费无遮挡视频| 啦啦啦在线免费观看视频4| 9191精品国产免费久久| 黄色视频在线播放观看不卡| 黑人巨大精品欧美一区二区蜜桃| 十八禁网站免费在线| 一区二区三区精品91| 正在播放国产对白刺激| 亚洲欧美精品综合一区二区三区| 国产淫语在线视频| 91精品三级在线观看| 免费女性裸体啪啪无遮挡网站| 视频区图区小说| 国产日韩一区二区三区精品不卡| 亚洲精品av麻豆狂野| 一本久久精品| 国产在线观看jvid| 十八禁网站免费在线| 亚洲精品日韩在线中文字幕| 国产激情久久老熟女| 日韩免费高清中文字幕av| 欧美日韩亚洲高清精品| 日本黄色日本黄色录像| 欧美日韩视频精品一区| 高潮久久久久久久久久久不卡| 亚洲精品国产区一区二| 国产免费av片在线观看野外av| 性高湖久久久久久久久免费观看| 久久天躁狠狠躁夜夜2o2o| 色婷婷av一区二区三区视频| 高清在线国产一区| 咕卡用的链子| 精品国产超薄肉色丝袜足j| 午夜影院在线不卡| 亚洲中文av在线| 亚洲国产欧美一区二区综合| 亚洲伊人久久精品综合| 热re99久久国产66热| 激情视频va一区二区三区| 国产精品久久久av美女十八| 欧美日本中文国产一区发布| 最新的欧美精品一区二区| 国产男人的电影天堂91| 美国免费a级毛片| 午夜两性在线视频| avwww免费| 日本黄色日本黄色录像| 久久久久国内视频| 国产亚洲午夜精品一区二区久久| 老司机亚洲免费影院| 国产精品影院久久| 国产成人精品在线电影| 18禁黄网站禁片午夜丰满| 人人澡人人妻人| 亚洲精品第二区| 精品人妻在线不人妻| 国产伦人伦偷精品视频| 99久久国产精品久久久| 黄色怎么调成土黄色| 国产精品国产三级国产专区5o| 无遮挡黄片免费观看| 欧美日韩中文字幕国产精品一区二区三区 | 老鸭窝网址在线观看| 久久综合国产亚洲精品| 亚洲少妇的诱惑av| 丝袜美腿诱惑在线| 免费不卡黄色视频| 亚洲国产欧美网| 建设人人有责人人尽责人人享有的| 日韩视频在线欧美| 亚洲九九香蕉| 色婷婷av一区二区三区视频| 亚洲国产欧美网| 两个人免费观看高清视频| 最近最新免费中文字幕在线| av一本久久久久| 97人妻天天添夜夜摸| 美女中出高潮动态图| 久久精品人人爽人人爽视色| 一区二区av电影网| 脱女人内裤的视频| 亚洲全国av大片| 欧美激情 高清一区二区三区| 久久免费观看电影| 久久久精品94久久精品| 一级a爱视频在线免费观看| 少妇被粗大的猛进出69影院| 男女下面插进去视频免费观看| 国产亚洲欧美精品永久| 手机成人av网站| 欧美97在线视频| 欧美精品一区二区大全| 亚洲欧洲精品一区二区精品久久久| 蜜桃国产av成人99| 狠狠婷婷综合久久久久久88av| 亚洲免费av在线视频| 亚洲欧美一区二区三区久久| 777久久人妻少妇嫩草av网站| 不卡一级毛片| 欧美+亚洲+日韩+国产| 精品一品国产午夜福利视频| 成年美女黄网站色视频大全免费| 最新在线观看一区二区三区| 午夜日韩欧美国产| 成年av动漫网址| 精品人妻熟女毛片av久久网站| 欧美日韩一级在线毛片| 免费一级毛片在线播放高清视频 | av网站在线播放免费| 妹子高潮喷水视频| 肉色欧美久久久久久久蜜桃| 国产一区二区 视频在线| 国产福利在线免费观看视频| 啦啦啦 在线观看视频| 精品久久蜜臀av无| 中文字幕制服av| 1024视频免费在线观看| 亚洲欧美一区二区三区久久| 亚洲成av片中文字幕在线观看| 大码成人一级视频| 又紧又爽又黄一区二区| 久久中文字幕一级| 人人澡人人妻人| 国产淫语在线视频| 欧美日韩福利视频一区二区| 男人爽女人下面视频在线观看| 十八禁网站免费在线| 2018国产大陆天天弄谢| 婷婷成人精品国产| 国产成人精品无人区| 日日夜夜操网爽| 九色亚洲精品在线播放| 超碰97精品在线观看| www.999成人在线观看| 后天国语完整版免费观看| 欧美 亚洲 国产 日韩一| 在线永久观看黄色视频| 少妇人妻久久综合中文| 久久久久国产精品人妻一区二区| 久久久国产精品麻豆| 中文字幕最新亚洲高清| 精品亚洲成a人片在线观看| 中亚洲国语对白在线视频| 9色porny在线观看| 国产高清视频在线播放一区 | 午夜精品久久久久久毛片777| av国产精品久久久久影院| 成年av动漫网址| 日韩视频一区二区在线观看| 99热网站在线观看| 亚洲专区字幕在线| 欧美在线黄色| 亚洲精品自拍成人| 久久久久久久久免费视频了| 12—13女人毛片做爰片一| 精品久久久久久久毛片微露脸 | 国产99久久九九免费精品| 老熟妇仑乱视频hdxx| 俄罗斯特黄特色一大片| 91麻豆av在线| 狂野欧美激情性bbbbbb| 欧美黄色片欧美黄色片| 国产成人影院久久av| 亚洲av电影在线观看一区二区三区| 国产视频一区二区在线看| 在线精品无人区一区二区三| 亚洲一卡2卡3卡4卡5卡精品中文| 美女福利国产在线| 国产野战对白在线观看| 亚洲欧美精品自产自拍| 多毛熟女@视频| 在线观看www视频免费| av有码第一页| 一本色道久久久久久精品综合| 老司机在亚洲福利影院| 久久中文看片网| 一进一出抽搐动态| 国产亚洲精品一区二区www | 日韩制服骚丝袜av| 我要看黄色一级片免费的| 热99久久久久精品小说推荐| 国产成人av教育| 老司机亚洲免费影院| 国产主播在线观看一区二区| 欧美成狂野欧美在线观看| 免费一级毛片在线播放高清视频 | 黄色片一级片一级黄色片| 夫妻午夜视频| 国产亚洲av高清不卡| 宅男免费午夜| 99香蕉大伊视频| 一本综合久久免费| 一区二区日韩欧美中文字幕| 黄片大片在线免费观看| 窝窝影院91人妻| 国产高清视频在线播放一区 | 99国产精品一区二区三区| 中文字幕人妻丝袜一区二区| 美女扒开内裤让男人捅视频| 国产精品久久久久久精品古装| 精品国产一区二区三区久久久樱花| 久久久久国产一级毛片高清牌| 亚洲精品第二区| 午夜福利免费观看在线| 午夜激情久久久久久久| 亚洲熟女毛片儿| 色94色欧美一区二区| 久久久久国产一级毛片高清牌| 97精品久久久久久久久久精品| 国产欧美日韩精品亚洲av| av片东京热男人的天堂| 国产老妇伦熟女老妇高清| 亚洲精品美女久久久久99蜜臀| 王馨瑶露胸无遮挡在线观看| 乱人伦中国视频| cao死你这个sao货| 精品国内亚洲2022精品成人 | 狠狠狠狠99中文字幕| 久久综合国产亚洲精品| 精品亚洲成国产av| 免费看十八禁软件| 99久久99久久久精品蜜桃| av视频免费观看在线观看| 国产极品粉嫩免费观看在线| 手机成人av网站| 成年人免费黄色播放视频| 777米奇影视久久| 18禁观看日本| 亚洲国产欧美一区二区综合| 国产日韩一区二区三区精品不卡| 99国产极品粉嫩在线观看| 亚洲激情五月婷婷啪啪| 美女福利国产在线| 久久精品aⅴ一区二区三区四区| 精品国产一区二区久久| 水蜜桃什么品种好| 91成年电影在线观看| 黑人猛操日本美女一级片| 极品少妇高潮喷水抽搐| 黄片大片在线免费观看| 我要看黄色一级片免费的| 亚洲久久久国产精品| 亚洲精品国产一区二区精华液| 三上悠亚av全集在线观看| 久久精品成人免费网站| netflix在线观看网站| 欧美黄色片欧美黄色片| 国产精品亚洲av一区麻豆| 久久久久久免费高清国产稀缺| 久久精品国产亚洲av高清一级| 亚洲色图综合在线观看| 欧美黑人精品巨大| 99久久精品国产亚洲精品| 婷婷色av中文字幕| 天天操日日干夜夜撸| 真人做人爱边吃奶动态| 蜜桃在线观看..| 一本色道久久久久久精品综合| 亚洲精品中文字幕在线视频| a级毛片在线看网站| 美国免费a级毛片| 亚洲国产日韩一区二区| 精品一区在线观看国产| av在线老鸭窝| 一区二区av电影网| 激情视频va一区二区三区| 满18在线观看网站| 王馨瑶露胸无遮挡在线观看| 深夜精品福利| 法律面前人人平等表现在哪些方面 | 国产在线观看jvid| 操出白浆在线播放| 日韩一区二区三区影片| cao死你这个sao货| 国产亚洲精品一区二区www | 久久久久久久久免费视频了| 精品人妻在线不人妻| 狂野欧美激情性xxxx| 成人国产一区最新在线观看| 欧美精品一区二区免费开放| 天堂俺去俺来也www色官网| 欧美亚洲 丝袜 人妻 在线| 亚洲五月色婷婷综合| 麻豆av在线久日| 欧美日韩av久久| 成年美女黄网站色视频大全免费| 午夜日韩欧美国产| 一二三四社区在线视频社区8| 国产精品 欧美亚洲| 国产亚洲av片在线观看秒播厂| 亚洲人成电影免费在线| 天堂俺去俺来也www色官网| 亚洲精华国产精华精| 免费日韩欧美在线观看| 狠狠婷婷综合久久久久久88av| 大片免费播放器 马上看| 一二三四社区在线视频社区8| 国产高清国产精品国产三级| 欧美黑人精品巨大| 午夜久久久在线观看| 少妇被粗大的猛进出69影院| 高清av免费在线| 啦啦啦免费观看视频1| 中文精品一卡2卡3卡4更新| 国产成人系列免费观看| 桃花免费在线播放| 国产精品影院久久| 三级毛片av免费| 久久久久久久国产电影| 亚洲五月婷婷丁香| 午夜两性在线视频| 亚洲一区二区三区欧美精品| 97人妻天天添夜夜摸| 亚洲七黄色美女视频| 一本久久精品| 欧美日韩亚洲综合一区二区三区_| 久久国产亚洲av麻豆专区| 国产三级黄色录像| 啦啦啦在线免费观看视频4| 在线精品无人区一区二区三| 国产高清国产精品国产三级| 亚洲视频免费观看视频| √禁漫天堂资源中文www| 亚洲国产看品久久| 大码成人一级视频| 国产成人系列免费观看| 亚洲成国产人片在线观看| 亚洲国产毛片av蜜桃av| 欧美日韩亚洲国产一区二区在线观看 | 欧美+亚洲+日韩+国产| 国产国语露脸激情在线看| 天天躁夜夜躁狠狠躁躁| 一进一出抽搐动态| 十八禁高潮呻吟视频| 性少妇av在线| 国产精品久久久av美女十八| 首页视频小说图片口味搜索| 亚洲成人免费av在线播放| 最新在线观看一区二区三区| 丁香六月天网| 欧美黄色淫秽网站| 国产xxxxx性猛交| 一本—道久久a久久精品蜜桃钙片| 国产精品久久久av美女十八| 美女扒开内裤让男人捅视频| 国产主播在线观看一区二区| www.熟女人妻精品国产| 他把我摸到了高潮在线观看 | 少妇的丰满在线观看| 国产日韩欧美在线精品| 可以免费在线观看a视频的电影网站| 新久久久久国产一级毛片| 一进一出抽搐动态| av天堂在线播放| 狂野欧美激情性xxxx| 日韩 欧美 亚洲 中文字幕| 亚洲国产中文字幕在线视频| 久久精品aⅴ一区二区三区四区| 男女国产视频网站| a级片在线免费高清观看视频| 91大片在线观看| 国产精品国产av在线观看| 正在播放国产对白刺激| 久久天堂一区二区三区四区| 国产日韩一区二区三区精品不卡| 国产伦理片在线播放av一区| 天堂8中文在线网| 欧美日韩中文字幕国产精品一区二区三区 | 国产亚洲精品一区二区www | 国产精品久久久久久精品电影小说| 精品国产乱码久久久久久男人| 9色porny在线观看| 久久毛片免费看一区二区三区| 国产成人av激情在线播放| 欧美在线黄色| 亚洲伊人久久精品综合| 久久久国产欧美日韩av| 午夜激情久久久久久久|