这个是科学网的大牛改编的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/
代码好长啊~~