日食预测代码之R代码

 这个是科学网的大牛改编的R代码。我计划将它改写为C或Matlab代码。后者简单些但是会牺牲运行效率。

### 快速日食搜索,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*t*t + 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.691425*t + 0.0001523*t2 )

- 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= ( 18461*cos(0.0571+ 8433.46616*t -0.640*t2 -1*t3)

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

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

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

+ 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 + 6*t3)

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

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

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

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

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

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

+ 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.691425*t+ 1.52 *t2 + 25*t3)

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

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

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

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

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

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

+ 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 + 25*t3)

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

+ 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 + 6*t3)

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

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

+ 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)

+ 10302*cos(1.10749 + 6283.07585*t)*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) #月亮黄纬速度

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

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

- 4184*sin(4.900 + 7214.06287*t)

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

vL = vL/36525

vB = vB/36525

vR = vR/36525; #每日速度

gm = mR*sin(mB)*vL/sqrt(vB*vB+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)

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; #判断是否为格里高利历日1582*372+10*31+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

return(res)

}

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

date2JD2000 <-

function(y,m,d){

juliand

res

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

return(res)

}

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

## "P"为偏食

## "T"为全食

## "A"为日环食

## "N"没有日食发生

find.solar.eclipse <-

function(y,m,d){

res

dd

date.solar.eclipse

type

res$NewMoonTime

res$EclipseType

return(res)

}

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

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

代码来自科学网

By 高达数字实验室 http://gndrive.org/

原文链接：日食预测代码之R代码，转载请注明来源！