<script language=javascript>
//========================
//今天日环食，时长11分8秒，特写本程序用于快速搜索日食，以此做个纪念
//设计本算法费时1整天
//在E5300/2G/G41几秒钟内可搜索数1万年的日食
//许剑伟,于家里,2010.1.15
//========================
function int2(x){return Math.floor(x);}
var JD={ //日期元件
JD:function(y,m,d){ //公历转儒略日
var n=0, G=0;
if(y*372+m*31+int2(d)>=588829) G=1; //判断是否为格里高利历日1582*372+10*31+15
if(m<=2) m+=12, y--;
if(G) n=int2(y/100), n=2-n+int2(n/4); //加百年闰
return int2(365.25*(y+4716)) + int2(30.6001*(m+1))+d+n - 1524.5;
},
DD:function(jd){ //儒略日数转公历
var r=new Object();
var D=int2(jd+0.5), F=jd+0.5-D, c; //取得日数的整数部份A及小数部分F
if(D>=2299161) c=int2((D-1867216.25)/36524.25),D+=1+c-int2(c/4);
D += 1524; r.Y = int2((D-122.1)/365.25);//年数
D -= int2(365.25*r.Y); r.M = int2(D/30.601); //月数
D -= int2(30.601*r.M); r.D = D; //日数
if(r.M>13) r.M -= 13, r.Y -= 4715;
else r.M -= 1, r.Y -= 4716;
//日的小数转为时分秒
F*=24; r.h=int2(F); F-=r.h;
F*=60; r.m=int2(F); F-=r.m;
F*=60; r.s=F;
return r;
},
DD2str:function(r){ //日期转为串
var Y=" "+r.Y, M="0"+r.M, D="0"+r.D;
var h=r.h, m=r.m, s=int2(r.s+.5);
if(s>=60) s-=60,m++;
if(m>=60) m-=60,h++;
h="0"+h; m="0"+m; s="0"+s;
Y=Y.substr(Y.length-5,5); M=M.substr(M.length-2,2); D=D.substr(D.length-2,2);
h=h.substr(h.length-2,2); m=m.substr(m.length-2,2); s=s.substr(s.length-2,2);
return Y+"-"+M+"-"+D+" "+h+":"+m+":"+s;
},
JD2str:function(jd){ //JD转为串
var r=this.DD( jd );
return this.DD2str( r );
}
};
/*****
ecFast()函数返回参数说明
r.jdSuo 朔时刻
r.lx 日食类型
*****/
function ecFast(jd){ //快速日食搜索,jd为朔时间(J2000起算的儒略日数不必很精确)
var re=new Object();
var rad=180*3600/Math.PI;
var t,t2,t3,t4;
var L,mB,mR,sR, vL,vB,vR;
var W = Math.floor((jd+8)/29.5306)*Math.PI*2; //合朔时的日月黄经差
//合朔时间计算,2000前+-40000年误差1小时以内，+-2000年小于10分钟
t = ( W + 1.08472 )/7771.37714500204; //平朔时间
re.jdSuo = 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*Math.PI;
re.ac=1, re.lx=\'N\';
if(Math.abs(Math.sin(L))>0.4) return re; //一般大于21度已不可能
t -= ( -0.0000331*t*t + 0.10976 *Math.cos( 0.785 + 8328.6914*t) )/7771;
t2=t*t;
L = -1.084719 +7771.377145013*t -0.0000331*t2 +
(22640 * Math.cos(0.785+ 8328.6914*t +0.000152*t2)
+4586 * Math.cos(0.19 + 7214.063*t -0.000218*t2)
+2370 * Math.cos(2.54 + 15542.754*t -0.000070*t2)
+ 769 * Math.cos(3.1 + 16657.383*t)
+ 666 * Math.cos(1.5 + 628.302*t)
+ 412 * Math.cos(4.8 + 16866.93*t)
+ 212 * Math.cos(4.1 -1114.63*t)
+ 205 * Math.cos(0.2 + 6585.76*t)
+ 192 * Math.cos(4.9 + 23871.45*t)
+ 165 * Math.cos(2.6 + 14914.45*t)
+ 147 * Math.cos(5.5 -7700.39*t)
+ 125 * Math.cos(0.5 + 7771.38*t)
+ 109 * Math.cos(3.9 + 8956.99*t)
+ 55 * Math.cos(5.6 -1324.18*t)
+ 45 * Math.cos(0.9 + 25195.62*t)
+ 40 * Math.cos(3.8 -8538.24*t)
+ 38 * Math.cos(4.3 + 22756.82*t)
+ 36 * Math.cos(5.5 + 24986.07*t)
-6893 * Math.cos(4.669257+628.3076*t)
- 72 * Math.cos(4.6261 +1256.62*t)
- 43 * Math.cos(2.67823 +628.31*t)*t
+ 21) / rad;
t += ( W - L ) / ( 7771.38
- 914 * Math.sin( 0.7848 + 8328.691425*t + 0.0001523*t2 )
- 179 * Math.sin( 2.543 +15542.7543*t )
- 160 * Math.sin( 0.1874 + 7214.0629*t ) );
re.jdSuo = jd = t*36525; //朔时刻
//纬 52,15 (角秒)
t2=t*t/10000,t3=t2*t/10000;
mB=
18461*Math.cos(0.0571+ 8433.46616*t -0.640*t2 -1*t3)
+ 1010*Math.cos(2.413 + 16762.1576 *t + 0.88 *t2 + 25*t3)
+ 1000*Math.cos(5.440 -104.7747 *t + 2.16 *t2 + 26*t3)
+ 624*Math.cos(0.915 + 7109.2881 *t + 0 *t2 + 7*t3)
+ 199*Math.cos(1.82 + 15647.529 *t -2.8 *t2 -19*t3)
+ 167*Math.cos(4.84 -1219.403 *t -1.5 *t2 -18*t3)
+ 117*Math.cos(4.17 + 23976.220 *t -1.3 *t2 + 6*t3)
+ 62*Math.cos(4.8 + 25090.849 *t + 2 *t2 + 50*t3)
+ 33*Math.cos(3.3 + 15437.980 *t + 2 *t2 + 32*t3)
+ 32*Math.cos(1.5 + 8223.917 *t + 4 *t2 + 51*t3)
+ 30*Math.cos(1.0 + 6480.986 *t + 0 *t2 + 7*t3)
+ 16*Math.cos(2.5 -9548.095 *t -3 *t2 -43*t3)
+ 15*Math.cos(0.2 + 32304.912 *t + 0 *t2 + 31*t3)
+ 12*Math.cos(4.0 + 7737.590 *t)
+ 9*Math.cos(1.9 + 15019.227 *t)
+ 8*Math.cos(5.4 + 8399.709 *t)
+ 8*Math.cos(4.2 + 23347.918 *t)
+ 7*Math.cos(4.9 -1847.705 *t)
+ 7*Math.cos(3.8 -16133.856 *t)
+ 7*Math.cos(2.7 + 14323.351 *t);
mB/=rad;
//距 106, 23 (千米)
mR = 385001
+20905*Math.cos(5.4971+ 8328.691425*t+ 1.52 *t2 + 25*t3)
+ 3699*Math.cos(4.900 + 7214.06287*t -2.18 *t2 -19*t3)
+ 2956*Math.cos(0.972 + 15542.75429*t -0.66 *t2 + 6*t3)
+ 570*Math.cos(1.57 + 16657.3828 *t + 3.0 *t2 + 50*t3)
+ 246*Math.cos(5.69 -1114.6286 *t -3.7 *t2 -44*t3)
+ 205*Math.cos(1.02 + 14914.4523 *t -1 *t2 + 6*t3)
+ 171*Math.cos(3.33 + 23871.4457 *t + 1 *t2 + 31*t3)
+ 152*Math.cos(4.94 + 6585.761 *t -2 *t2 -19*t3)
+ 130*Math.cos(0.74 -7700.389 *t -2 *t2 -25*t3)
+ 109*Math.cos(5.20 + 7771.377 *t)
+ 105*Math.cos(2.31 + 8956.993 *t + 1 *t2 + 25*t3)
+ 80*Math.cos(5.38 -8538.241 *t + 2.8 *t2 + 26*t3)
+ 49*Math.cos(6.24 + 628.302 *t)
+ 35*Math.cos(2.7 + 22756.817 *t -3 *t2 -13*t3)
+ 31*Math.cos(4.1 + 16171.056 *t -1 *t2 + 6*t3)
+ 24*Math.cos(1.7 + 7842.365 *t -2 *t2 -19*t3)
+ 23*Math.cos(3.9 + 24986.074 *t + 5 *t2 + 75*t3)
+ 22*Math.cos(0.4 + 14428.126 *t -4 *t2 -38*t3)
+ 17*Math.cos(2.0 + 8399.679 *t);
mR/=6378.1366;
t=jd/365250, t2=t*t, t3=t2*t;
//误0.0002AU
sR = 10001399 //日地距离
+167070*Math.cos(3.098464 + 6283.07585*t)
+ 1396*Math.cos(3.0552 + 12566.1517 *t)
+ 10302*Math.cos(1.10749 + 6283.07585*t)*t
+ 172*Math.cos(1.064 + 12566.152 *t)*t
+ 436*Math.cos(5.785 + 6283.076 *t)*t2
+ 14*Math.cos(4.27 + 6283.08 *t)*t3;
sR*=1.49597870691/6378.1366*10;
//经纬速度
t=jd/36525;
vL = 7771 //月日黄经差速度
-914*Math.sin(0.785 + 8328.6914*t)
-179*Math.sin(2.543 +15542.7543*t)
-160*Math.sin(0.187 + 7214.0629*t);
vB =-755*Math.sin(0.057 + 8433.4662*t) //月亮黄纬速度
- 82*Math.sin(2.413 +16762.1576*t);
vR =-27299*Math.sin(5.497 + 8328.691425*t)
- 4184*Math.sin(4.900 + 7214.06287*t)
- 7204*Math.sin(0.972 +15542.75429*t);
vL/=36525, vB/=36525, vR/=36525; //每日速度
var gm = mR*Math.sin(mB)*vL/Math.sqrt(vB*vB+vL*vL), smR=sR-mR; //gm伽马值,smR日月距
var mk = 0.2725076, sk = 109.1222;
var f1 = (sk+mk)/smR, r1 = mk+f1*mR; //tanf1半影锥角, r1半影半径
var f2 = (sk-mk)/smR, r2 = mk-f2*mR; //tanf2本影锥角, r2本影半径
var b = 0.9972, Agm = Math.abs(gm), Ar2 = Math.abs(r2);
var fh2 = mR-mk/f2, h = Agm<1 ? Math.sqrt(1-gm*gm) : 0; //fh2本影顶点的z坐标
var ls1,ls2,ls3,ls4;
if(fh2<h) re.lx = \'T\';
else re.lx = \'A\';
ls1 = Agm-(b+r1 ); if(Math.abs(ls1)<0.016) re.ac=0; //无食界
ls2 = Agm-(b+Ar2); if(Math.abs(ls2)<0.016) re.ac=0; //偏食界
ls3 = Agm-(b ); if(Math.abs(ls3)<0.016) re.ac=0; //无中心食界
ls4 = Agm-(b-Ar2); if(Math.abs(ls4)<0.016) re.ac=0; //有中心食界(但本影未全部进入)
if (ls1>0) re.lx = \'N\'; //无日食
else if(ls2>0) re.lx = \'P\'; //偏食
else if(ls3>0) re.lx += \'0\'; //无中心
else if(ls4>0) re.lx += \'1\'; //有中心(本影未全部进入)
else{ //本影全进入
if(Math.abs(fh2-h)<0.019) re.ac=0;
if( Math.abs(fh2)<h ){
var dr = vR*h/vL/mR;
var H1 = mR-dr-mk/f2; //入点影锥z坐标
var H2 = mR+dr-mk/f2; //出点影锥z坐标
if(H1>0) re.lx=\'H3\'; //环全全
if(H2>0) re.lx=\'H2\'; //全全环
if(H1>0&&H2>0) re.lx=\'H\'; //环全环
if(Math.abs(H1)<0.019) re.ac=0;
if(Math.abs(H2)<0.019) re.ac=0;
}
}
return re;
}
//==============代码测试====================
var J2000=2451545;
var i,jd,n1=0,n2=0;
for(i=0;i<100;i++){
jd=2454680-J2000+29.5306*i;
r=ecFast(jd); //A0
if(!r.ac) n1++;
if(r.lx==\'N\') continue;
n2++;
document.write( JD.JD2str(r.jdSuo+J2000)+\' \'+r.lx+\' ac=\'+r.ac+\'<br>\');
}
document.write(\'ac=0不很准确\'+n1+\'个,ac=1精确\'+n2);
</script>