中华农历论坛知识讨论区历法知识 → 日食搜索器


  共有27540人关注过本帖树形打印

主题:日食搜索器

帅哥哟,离线,有人找我吗?
xjw01
  1楼 个性首页 | 信息 | 搜索 | 邮箱 | 主页 | UC


加好友 发短信
等级:蝙蝠侠 帖子:866 积分:3804 威望:3 精华:20 注册:2008/3/20 22:14:00
日食搜索器  发帖心情 Post By:2010/1/15 22:05:00

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

[此贴子已经被作者于2010-1-16 20:04:03编辑过]

支持(0中立(0反对(0单帖管理 | 引用 | 回复 回到顶部
帅哥哟,离线,有人找我吗?
浪-淘-沙
  2楼 个性首页 | 信息 | 搜索 | 邮箱 | 主页 | UC


加好友 发短信
等级:版主 帖子:2068 积分:4263 威望:5 精华:4 注册:2008/11/13 21:03:00
  发帖心情 Post By:2010/1/15 22:39:00

运行结果:

2008-08-01 12:00:00 T ac=1
2009-01-25 16:24:23 A ac=1
2009-07-21 20:48:46 T ac=1
2010-01-15 01:13:09 A ac=1
2010-07-11 05:37:32 T ac=1
2011-01-04 10:01:55 P ac=1
2011-06-01 01:42:14 P ac=1
2011-06-30 14:26:18 P ac=1
2011-11-25 06:06:37 P ac=1
2012-05-20 10:31:00 A ac=1
2012-11-13 14:55:24 T ac=1
2013-05-09 19:19:47 A ac=1
2013-11-02 23:44:10 H3 ac=1
2014-04-29 04:08:33 A0 ac=0
2014-10-23 08:32:56 P ac=1
2015-03-20 00:13:15 T ac=1
2015-09-13 04:37:38 P ac=1
2016-03-08 09:02:01 T ac=1
ac=0不很准确1个,ac=1精确18


支持(0中立(0反对(0单帖管理 | 引用 | 回复 回到顶部
帅哥哟,离线,有人找我吗?
浪-淘-沙
  3楼 个性首页 | 信息 | 搜索 | 邮箱 | 主页 | UC


加好友 发短信
等级:版主 帖子:2068 积分:4263 威望:5 精华:4 注册:2008/11/13 21:03:00
  发帖心情 Post By:2010/1/15 22:42:00

用寿星计算结果:

中心点特征
影轴地心距 γ = 0.4002
中心地标 (经,纬) = -69.29,1.63
中心时刻 tm = 2010-01-15 07:07:40
太阳方位 (经,纬) = 345,66
日食类型 LX = A 环食
食分=0.9190, 食延=11m8s, 食带=335km

 

**************

用一楼的程序得到:

2010-01-15 01:13:09 A ac=1

时间上相差很多啊?


支持(0中立(0反对(0单帖管理 | 引用 | 回复 回到顶部
帅哥哟,离线,有人找我吗?
xjw01
  4楼 个性首页 | 信息 | 搜索 | 邮箱 | 主页 | UC


加好友 发短信
等级:蝙蝠侠 帖子:866 积分:3804 威望:3 精华:20 注册:2008/3/20 22:14:00
  发帖心情 Post By:2010/1/16 8:59:00

这是个输出显示的问题
测试代码中的
 document.write( JD.JD2str(jd+J2000)+' '+r.lx+' ac='+r.ac+'<br>');
改为
 document.write( JD.JD2str(r.jdSuo+J2000)+' '+r.lx+' ac='+r.ac+'<br>');
这样误差就只有几分钟。

此算法的平均误差为30公里左右,所以的“误算”估计不会超30公里*2/地球半径=1%,为了保险,程序输出一个ac参数,表示结果的可信度,如果为1,表示计算结果准确,如果ac=0,则表示计算结果有危险。判断依据100公里为容差指标,它远小于30公里。即,如果日食类型处于某些临界,可判为环食,也可判为全食,设环食与全食的临界条件为d公里,计算值如果与d靠近到不足100公里,则输出ac=0

计算过程中,各种误差最大值一般控制在100公里以内,平均误差小于50公里,误差叠加后,会有相消效果,实际误差一般更小。

 

影轴地心距误差约为0.005

 

还有,寿星输出的是最大食时间,这里输出的是朔时间。

 

 

 

此算法比一般的简易算法要精确得多。

 

算法原理描述如下。

 

1、坐标中心(O点)设在地心,将太阳与地心的连线坐为z坐标,过地心沿黄道逆方向为x轴,y轴指向黄极

2、计算日月合朔时间jd,那么就有:jd时刻,这个坐标中太阳经度为0,月亮经度也为0,也就是说日月同时在yz平面。太阳黄为不用计算,它的值为0,再计算太阳距离sR,再计算月亮黄纬B和月亮距离R。过太阳和月亮作一直线,交于y轴上的D点,那么OD=R*sin(B),这个约等式中R为地月距,B为月亮黄纬

3、求变差,即单位时间内月亮黄纬与黄纬与距离的变化速度vL,vB,vR,这里说的速度,是相对于上述坐标的。精度不要求很高,两位有效数据足够。

4、D点以及vL,vB求出后,可确定影轴在xy平面扫过的直线,并求出该直线与地球交线半弦长h以及影轴和地心距离gm,进而判断出有无中心线。

5、利用vR和R和h,vL,等参数,可确定影子移动过程中本影锥顶点位置的移动,准确判断是全食还是环食还是全环食等。

 

 

对精度要求不高的情况下,食计算是比较简的。


支持(0中立(0反对(0单帖管理 | 引用 | 回复 回到顶部
帅哥哟,离线,有人找我吗?
浪-淘-沙
  5楼 个性首页 | 信息 | 搜索 | 邮箱 | 主页 | UC


加好友 发短信
等级:版主 帖子:2068 积分:4263 威望:5 精华:4 注册:2008/11/13 21:03:00
  发帖心情 Post By:2010/1/16 9:44:00

以下是引用xjw01在2010-1-16 8:59:00的发言:

这是个输出显示的问题
测试代码中的
 document.write( JD.JD2str(jd+J2000)+' '+r.lx+' ac='+r.ac+'<br>');
改为
 document.write( JD.JD2str(r.jdSuo+J2000)+' '+r.lx+' ac='+r.ac+'<br>');
这样误差就只有几分钟。

我测试了一下,果然如此。

多谢讲解。


支持(0中立(0反对(0单帖管理 | 引用 | 回复 回到顶部
帅哥哟,离线,有人找我吗?
xjw01
  6楼 个性首页 | 信息 | 搜索 | 邮箱 | 主页 | UC


加好友 发短信
等级:蝙蝠侠 帖子:866 积分:3804 威望:3 精华:20 注册:2008/3/20 22:14:00
  发帖心情 Post By:2010/1/16 10:06:00

以上计算没有采用复杂的算法,大多采用直线与圆的关系来解决问题,精度虽然不高,用于食带位置计算意义不大,但仍总体准确把握日食情况。

我不太清楚古人具体计算日食的方法,我觉得以下几个问题没有解决,日食常常会失算,如果这些问题解决了,一般不会失算:

1、地球是一个球形,月地距离是球半径的60倍左右,这个概念是比较重要的。月亮黄经与黄纬。如果朔日前几天,天气良好,是可以测得月亮位置的。关键还是视差问题要解决。月亮的纬度受视差影响,它关系到“食与不食”,黄经误差(合朔误差)首先关系到日食发生时间,误差多了也会影响到“食与不食”。如果许可食带误差300公里,那么允计黄经误差1500角秒(约1小时)。黄纬的允计误差为黄经误差的1/10。在ELP方法计算月球星历,为达到以上精度,所需的黄经项数与黄纬项数差不多,黄纬略多一些。

在古代,如果能提前12小时至24小时测到月亮黄经,再得日月黄经差,再用这个差值除以平均日月黄经差变化速度,就可使合朔时刻误差控制在1小时以内。当然,古代后来也已经知道了月亮的黄经变化速度的规律,所以他们还可以把合朔时间控制在0.5小时以内,甚至10分钟以内。

如果认为天圆地方,那么正确的视差改正就无从谈起,所做的视差改正都将变成经验性的,一点也不可靠。把一个圆的东西看成方的,误差自然会很大。

2、太阳离地球很远,看成无穷远也无妨。但太阳黄经要算准,估计古天文这方面的计算准确度应该没有太在问题。

3、坐标变换问题要解决。

 

当然,古代要准确计算全食位置不大可能,没有精密星历解决不了这个问题。

 

以上是我昨天设计算法时进行误差测算得到的一些心得。


支持(0中立(0反对(0单帖管理 | 引用 | 回复 回到顶部
帅哥哟,离线,有人找我吗?
浪-淘-沙
  7楼 个性首页 | 信息 | 搜索 | 邮箱 | 主页 | UC


加好友 发短信
等级:版主 帖子:2068 积分:4263 威望:5 精华:4 注册:2008/11/13 21:03:00
  发帖心情 Post By:2010/1/16 11:02:00

以下是引用xjw01在2010-1-16 10:06:00的发言:

 

当然,古代要准确计算全食位置不大可能,没有精密星历解决不了这个问题。

许兄所断甚是。

 

从《明史》有关历法的记录(卷三十一开始),徐光启等人谈到日食推算问题,提到,就算是郭守敬等人编的《授时历》,对于日食的预测也往往不准的。郭当年在世时,预测日食也同样有偏差的。

而采用西洋天文学能准确预测日食。当时也进行过校验。这使得《崇祯历书》(其实就是仿西洋历法)能够得以编撰的主要原因之一。


支持(0中立(0反对(0单帖管理 | 引用 | 回复 回到顶部
帅哥哟,离线,有人找我吗?
anita315
  8楼 个性首页 | 信息 | 搜索 | 邮箱 | 主页 | UC


加好友 发短信
等级:论坛游民 帖子:62 积分:1214 威望:0 精华:0 注册:2009/4/23 20:37:00
  发帖心情 Post By:2010/1/19 17:33:00

請教本程式如何新增對某區域(或某經緯度) "是否可視的" 功能 ?
感謝樓主的文章 , 使我獲益良多 !

支持(0中立(0反对(0单帖管理 | 引用 | 回复 回到顶部
帅哥哟,离线,有人找我吗?
xjw01
  9楼 个性首页 | 信息 | 搜索 | 邮箱 | 主页 | UC


加好友 发短信
等级:蝙蝠侠 帖子:866 积分:3804 威望:3 精华:20 注册:2008/3/20 22:14:00
  发帖心情 Post By:2010/1/19 22:08:00

你是不是说的增加地方日食情况的计算?

以上算法是低精度的,误差较大,误差大约30公里左右,不太适合计算地方食。如果精度再提高3倍,就可以了。


支持(0中立(0反对(0单帖管理 | 引用 | 回复 回到顶部
帅哥哟,离线,有人找我吗?
anita315
  10楼 个性首页 | 信息 | 搜索 | 邮箱 | 主页 | UC


加好友 发短信
等级:论坛游民 帖子:62 积分:1214 威望:0 精华:0 注册:2009/4/23 20:37:00
  发帖心情 Post By:2010/1/20 9:04:00

我的想法是某區域(或某經緯度) "是否可視的" 功能 ,
而此區域是以 "省" , 或 "縣" 為單位 , 所以 "误差大约30公里左右" ,
應該是可以接受的 . 感謝樓主的回應與幫忙 ! ^_^...

支持(0中立(0反对(0单帖管理 | 引用 | 回复 回到顶部
总数 15 1 2 下一页

返回版面帖子列表

日食搜索器








签名