2010-11-13

如何计算日食的时间和类型（附R源代码）？

日食的原理在小学课本中就有介绍，也就是月亮的影子遮挡住了太阳。日食每次必定发生在农历初一，因为农历初一是朔，此时从地心看起来月亮距离太阳最近。为什么要强调地心呢？因为在地球表面的不同位置，有视差的存在。天文学家在计算天象的时候，很多时候会先计算地心坐标。具体到某一个地点，再归算到地平坐标。

日食的原理说来简单，但是计算起来并不容易。在我国古代，对日食的时间和类型的计算，是检验历法是否准确的重要标准。很多部历法由于不能准确的预测日月食，而不得不被新的历法所替代。如《新唐书.历志一》记载：“唐终始二百九十馀年，而历八改。初曰戊寅元历，曰麟德甲子元历，曰开元大衍历。”《新唐书.历志三上》又记载：”开元九年，麟德历署日蚀不效，诏僧一行作新历，推大衍数立术以应之…”。被誉为唐代历法首屈一指的《大衍历》，最后也难免被更为精确的历法所代替。

这里采用的算法，是基于比利时天文学家Jean Meeus提出的现代天文算法。详情请参见许剑伟老师翻译的《天文算法》第52章。

全部R程序改变自剑伟老师的《寿星万年历》javascript源代码。

### 快速日食搜索,jd为朔时间

### 改编自许剑伟老师的《寿星万年历》javascript源代码

### 输入儒略日2000，计算最近的朔的时间，并判断该次朔是否会发生日食。

DateSolarEclipse <-

function (jd){

rad = 180***3600/**pi; #每弧度的角秒数

radd = 180**/**pi; #每弧度的度数

pi2 = pi*****2;

pi_2 = pi**/**2;

J2000 = 2451545;

re**=list()**;

W = floor**((jd+8)/29.5306)*pi***2; #合朔时的日月黄经差

#合朔时间计算,2000前+-4000年误差1小时以内，+-2000年小于10分钟

t = ( W + 1.08472 **)/**7771.37714500204; #平朔时间

re**$jd = re$jdNewMoon = t***36525;

t2**=t***t

t3**=t2***t

t4**=t3***t;

L = ( 93.2720993**+483202.0175273***t**-0.0034029***t2**-t3/3526000+t4/**863310000 )/180*pi;

re**$ac=**1

re**$type=**’N’;

if**(abs(sin(L))>0.4)** return**(re)**; #一般大于21度已不可能

t = t - ( -0.0000331*tt + 0.10976 cos( 0.785 + 8328.6914*t)** **)/**7771;

t2**=** t*****t;

L = **(-1.084719 +7771.377145013*t **-0.0000331*t2 +

(22640 ***** cos(0.785+ 8328.6914*t +0.000152*t2)

+4586 ***** cos(0.19 + 7214.063*t -0.000218*t2**)**

+2370 ***** cos(2.54 + 15542.754*t -0.000070*t2**)**

+ 769 ***** cos**(3.1 + 16657.383***t**)**

+ 666 ***** cos**(1.5 + 628.302***t**)**

+ 412 ***** cos**(4.8 + 16866.93***t**)**

+ 212 ***** cos**(4.1 -1114.63*t)**

+ 205 ***** cos**(0.2 + 6585.76***t**)**

+ 192 ***** cos**(4.9 + 23871.45***t**)**

+ 165 ***** cos**(2.6 + 14914.45***t**)**

+ 147 ***** cos**(5.5 -7700.39*t)**

+ 125 ***** cos**(0.5 + 7771.38***t**)**

+ 109 ***** cos**(3.9 + 8956.99***t**)**

+ 55 ***** cos**(5.6 -1324.18*t)**

+ 45 ***** cos**(0.9 + 25195.62***t**)**

+ 40 ***** cos**(3.8 -8538.24*t)**

+ 38 ***** cos**(4.3 + 22756.82***t**)**

+ 36 ***** cos**(5.5 + 24986.07***t**)**

-6893 ***** cos(4.669257+628.3076*t**)**

- 72 ***** cos**(4.6261 +1256.62*t)**

- 43 ***** cos**(2.67823 +628.31*t)***t

+ 21**)** / rad**)**;

t = t + ( W - L ) / ( 7771.38

- 914 ***** sin**(** 0.7848 + 8328.691425t + 0.0001523t2 )

- 179 ***** sin**(** 2.543 +15542.7543*t )

- 160 ***** sin**(** 0.1874 + 7214.0629*****t ) );

re**$jd = re$jdNewMoon = jd = t***36525; #朔时刻

#纬 52,15 (角秒)

t2**=t***t**/**10000;

t3**=t2***t**/**10000;

mB**=** ( 18461cos(0.0571+ 8433.46616****t -0.640*t2 -1*t3**)**

+ 1010***cos(**2.413 + 16762.1576 **t + 0.88 t2 + 25t3)*

+ 1000***cos(**5.440 **-**104.7747 **t + 2.16 t2 + 26t3)*

+ 624***cos(**0.915 + 7109.2881 **t + 0 t2 + 7t3)*

+ 199***cos(**1.82 + 15647.529 *****t **-**2.8 *t2 -19*t3)

+ 167***cos(**4.84 **-**1219.403 *****t **-**1.5 *t2 -18*t3)

+ 117***cos(**4.17 + 23976.220 *****t -1.3 t2 + 6t3)

+ 62***cos(**4.8 + 25090.849 **t + 2 t2 + 50t3)*

+ 33***cos(**3.3 + 15437.980 **t + 2 t2 + 32t3)*

+ 32***cos(**1.5 + 8223.917 **t + 4 t2 + 51t3)*

+ 30***cos(**1.0 + 6480.986 **t + 0 t2 + 7t3)*

+ 16***cos(**2.5 **-**9548.095 *****t **-**3 *t2 -43*t3)

+ 15***cos(**0.2 + 32304.912 **t + 0 t2 + 31t3)*

+ 12***cos(**4.0 + 7737.590 *t)

+ 9***cos(**1.9 + 15019.227 *t)

+ 8***cos(**5.4 + 8399.709 *t)

+ 8***cos(**4.2 + 23347.918 *t)

+ 7***cos(**4.9 **-**1847.705 *t)

+ 7***cos(**3.8 **-**16133.856 *t)

+ 7***cos(**2.7 + 14323.351 **t))*;

mB = mB / rad;

#月地距离 106, 23 (千米)

mR = **(**385001

+20905*cos**(5.4971+** 8328.691425t*+ 1.52 t2 + 25t3)**

+ 3699cos(4.900 + 7214.06287t **-**2.18 *t2 -19*t3)

+ 2956cos(0.972 + 15542.75429t -0.66 t2 + 6t3)

+ 570***cos(**1.57 + 16657.3828 **t + 3.0 t2 + 50t3)*

+ 246***cos(**5.69 **-**1114.6286 *****t **-**3.7 *t2 -44*t3)

+ 205***cos(**1.02 + 14914.4523 *****t -1 t2 + 6t3)

+ 171***cos(**3.33 + 23871.4457 **t + 1 t2 + 31t3)*

+ 152***cos(**4.94 + 6585.761 *****t **-**2 *t2 -19*t3)

+ 130***cos(**0.74 **-**7700.389 *****t **-**2 *t2 -25*t3)

+ 109***cos(**5.20 + 7771.377 *t)

+ 105***cos(**2.31 + 8956.993 **t + 1 t2 + 25t3)*

+ 80***cos(**5.38 **-**8538.241 **t + 2.8 t2 + 26t3)*

+ 49***cos(**6.24 + 628.302 *t)

+ 35***cos(**2.7 + 22756.817 *****t **-**3 *t2 -13*t3)

+ 31***cos(**4.1 + 16171.056 *****t -1 t2 + 6t3)

+ 24***cos(**1.7 + 7842.365 *****t **-**2 *t2 -19*t3)

+ 23***cos(**3.9 + 24986.074 **t + 5 t2 + 75t3)*

+ 22***cos(**0.4 + 14428.126 *****t **-**4 *t2 -38*t3)

+ 17***cos(**2.0 + 8399.679 **t))*;

mR = mR**/** 6378.1366;

t**=jd/**365250;

t2**=t***t;

t3**=t2***t;

#误0.0002AU

sR = **(**10001399 #日地距离

+167070*cos**(3.098464 + 6283.07585***t**)**

+ 1396***cos(**3.0552 + 12566.1517 *t)

+ 10302cos(1.10749 + 6283.07585t**)***t

+ 172***cos(**1.064 + 12566.152 ***t)***t

+ 436***cos(**5.785 + 6283.076 ***t)***t2

+ 14***cos(**4.27 + 6283.08 *t)*t3)

sR = sR ***** 1.49597870691**/6378.1366***10;

#经纬速度

t**=jd/**36525;

vL = **(**7771 #月日黄经差速度

-914*sin**(0.785 + 8328.6914***t**)**

-179*sin**(2.543 +15542.7543*t)**

-160*sin**(0.187 + 7214.0629***t**))**;

vB = (-755*sin**(0.057 + 8433.4662***t**)** #月亮黄纬速度

- 82sin*(**2.413 +16762.1576*t));

vR =(-27299*sin**(5.497 + 8328.691425***t**)**

- 4184sin(4.900 + 7214.06287t**)**

- 7204sin*(**0.972 +15542.75429*t));

vL = vL**/**36525

vB = vB**/**36525

vR = vR**/**36525; #每日速度

gm = mRsin(mB)*vL/sqrt(vBvB**+vL***vL**)**;

smR**=** sR**-**mR; #gm伽马值,smR日月距

mk = 0.2725076;

sk = 109.1222;

f1 = (sk+mk)/smR; r1 = mk+f1*mR; #tanf1半影锥角, r1半影半径

f2 = (sk-mk)/smR; r2 = mk-f2*mR; #tanf2本影锥角, r2本影半径

b = 0.9972;

Agm = abs**(gm)**;

Ar2 = abs**(r2)**;

fh2 = mR**-mk/**f2;

if**(Agm < 1){**

h = sqrt**(1-gm***gm**)**

}else{

h = 0

}

## h = Agm<1 ? sqrt(1-gm*gm) : 0; #fh2本影顶点的z坐标

## ls1,ls2,ls3,ls4;

if**(fh2<h)** {

re**$**type = ‘T’; #全食

}else{

re**$**type = ‘A’; #环食

}

ls1 = Agm**-(b+**r1 );

if**(abs(ls1)<0.016)** re**$ac=**0; #无食分界

ls2 = Agm**-(b+Ar2)**;

if**(abs(ls2)<0.016)** re**$ac=**0; #偏食分界

ls3 = Agm**-(**b );

if**(abs(ls3)<0.016)** re**$ac=**0; #无中心食分界

ls4 = Agm**-(b-Ar2)**;

if**(abs(ls4)<0.016)** re**$ac=**0; #有中心食分界(但本影未全部进入)

if (ls1>0) {

re**$**type = ‘N’; #无日食

} else {

if**(ls2>0)** {

re**$**type = ‘P’; #偏食

} else {

if**(ls3>0)** {

re**$type = paste(re$type, ‘0’, sep = “”)**; #无中心

} else {

if**(ls4>0)** {

re**$type = paste(re$type, ‘1’, sep = “”)**; #有中心(本影未全部进入)

} else { #本影全进入

if**(abs(fh2-h)<0.019)** re**$ac=**0;

if**(** abs**(fh2)<**h ){

dr = vR***h/vL/**mR;

H1 = mR**-dr-mk/**f2; #入点影锥z坐标

H2 = mR**+dr-mk/**f2; #出点影锥z坐标

if**(H1>0)** re**$type=**’H3’; #环全全

if**(H2>0)** re**$type=**’H2’; #全全环

if**(H1>0&H2>0)** re**$type=**’H’; #环全环

if**(abs(H1)<0.019)** re**$ac=**0;

if**(abs(H2)<0.019)** re**$ac=**0;

}

return (re);

}

#取整数部分

int2 <-

function (v) {

return (floor(v));

}

### 公历转儒略日

JD <-

function**(y,m,d){**

n**=0;G=**0;

if**(y***372**+m***31**+int2(d)>=588829)**

G = 1; #判断是否为格里高利历日1582372+1031+15

if**(m<=2)** {

m = m **+**12;

y = y - 1;

}

if**(G > 0)** {

n**=int2(y/100)**;

n**=2-n+int2(n/4)**; #加百年闰

}

res <- int2**(365.25*(y+4716))** + int2**(30.6001*(m+1))+d+**n - 1524.5;

class**(res)** <- “JD”

return**(res)**

}

### 将儒略日转换为儒略日2000

JD2000 <-

function**(jd){**

res <- jd - JD**(2000,1,1)**

return**(res)**

}

### 将日期转换为儒略日2000

date2JD2000 <-

function**(y,m,d){**

juliand <- JD**(y,m,d)**

res <- JD2000**(juliand)**

return**(res)**

}

### 打印JD类型对象的方法

print.JD <-

function**(JD){**

print**(sprintf(“%f”,JD))**

}

##儒略日转为公历

Julian2Date <-

function**(jd){**

r = list**()**;

D = int2**(jd+0.5)**;

F**=** jd**+0.5-**D #取得日数的整数部份A及小数部分F

if**(D >= 2299161)** {

c = int2**((D-1867216.25)/36524.25)**;

D = D + 1 + c - int2**(c/4)**;

}

D = D + 1524;

r**$Y = int2((D - 122.1)/365.25)**;#年数

D = D - int2**(365.25***r**$Y)**;

r**$M = int2(D/30.601)**; #月数

D = D - int2**(30.601***r**$M)**;

r**$**D = D; #日数

if**(r$M>13)** {

r**$M = r$**M - 13;

r**$Y = r$**Y - 4715;

} else {

r**$M = r$**M - 1;

r**$Y = r$**Y - 4716;

}

#日的小数转为时分秒

F = F ***** 24;

r**$h=int2(F)**;

F = F - r**$**h;

F = F ***** 60;

r**$m=int2(F)**;

F = F **-r$**m;

F = F ***** 60;

r**$s=**F;

return**(r)**;

}

#日期转为字符串

DD2str <-

function**(r){**

Y = r**$**Y

M = r**$**M

D = r**$**D;

h = r**$**h;

m = r**$**m;

s = int2**(r$s+0.5)**;

if**(s >= 60)** {

s = s - 60;

m = m + 1;

}

if**(m >= 60)** {

m = m - 60

h = h + 1;

}

h = paste**(“0”, h, sep = “”)**;

m = paste**(“0”, m, sep = “”)**;

s = paste**(“0”, s, sep = “”)**;

Y = substr**(Y, nchar(Y)-** 4, nchar**(Y)** );

M = substr**(M, nchar(M)-** 1, nchar**(M)** );

D = substr**(D, nchar(D)-** 1, nchar**(D)** );

h = substr**(h, nchar(h)-** 1, nchar**(h)** );

m = substr**(m, nchar(m)-** 1, nchar**(m)** );

s = substr**(s, nchar(s)-** 1, nchar**(s)** );

res <- paste**(Y,”-“,M,”-“,D,” “,h,”:”,m,”:”,s, sep = “”)**;

return**(res)**

}

JD2str <-

function**(jd){** #JD转为串

r = DD**(** jd );

res <- DD2str**(** r );

return**(res)**

}

### 判断距离输入日期最近的朔，有没有日食发生。

## “P”为偏食

## “T”为全食

## “A”为日环食

## “N”没有日食发生

find.solar.eclipse <-

function**(y,m,d){**

res <- list**()**

dd <- DateSolarEclipse**(date2JD2000(y,m,d))**

date.solar.eclipse <- DD2str**(Julian2Date(dd$jdNewMoon + JD(2000,1,1)** + 0.5**))**

type <- dd**$**type

res**$**NewMoonTime <- date.solar.eclipse

res**$**EclipseType <- type

return**(res)**

}

##距离1987年9月23日最近的朔，以及日食类型

find.solar.eclipse**(1987,9,23)**