eph.js

寿星万年历 基于天文算法的万年历 javascript 年代范围大 精度高

  v = (g.x-G.x)/dt;
  dt= -(G.y*u+G.x*v)/(u*u+v*v); jd += dt; //极值时间

  //求直线到影子中心的最小值
  var x=G.x+dt*v, y=G.y+dt*u, rmin=Math.sqrt(x*x+y*y);
  //注意,以上计算得到了极值及最小距rmin,但没有再次计算极值时刻的半径,对以下的判断造成一定的风险,必要的话可以再算一次。不过必要性不很大，因为第一次极值计算已经很准确了,误差只有几秒
  //求月球与影子的位置关系
  if(rmin<=G.mr+G.er){ //食计算
   this.lT[0] = jd; //食甚

   this.lT[1] = this.lineT(G,v,u, G.mr+G.er, 0); //初亏
   this.lecXY(this.lT[1],g);
   this.lT[1] = this.lineT(g,v,u, g.mr+g.er, 0); //初亏再算一次

   this.lT[2] = this.lineT(G,v,u, G.mr+G.er, 1); //复圆
   this.lecXY(this.lT[2],g);
   this.lT[2] = this.lineT(g,v,u, g.mr+g.er, 1); //复圆再算一次
  }
  if(rmin<=G.mr+G.Er){ //半影食计算
   this.lT[3] = this.lineT(G,v,u, G.mr+G.Er, 0); //半影食始
   this.lecXY(this.lT[3],g);
   this.lT[3] = this.lineT(g,v,u, g.mr+g.Er, 0); //半影食始再算一次

   this.lT[4] = this.lineT(G,v,u, G.mr+G.Er, 1); //半影食终
   this.lecXY(this.lT[4],g);
   this.lT[4] = this.lineT(g,v,u, g.mr+g.Er, 1); //半影食终再算一次
  }
  if(rmin<=G.er-G.mr){ //全食计算
   this.lT[5] = this.lineT(G,v,u, G.er-G.mr, 0); //食既
   this.lecXY(this.lT[5],g);
   this.lT[5] = this.lineT(g,v,u, g.er-g.mr, 0); //食既再算一次

   this.lT[6] = this.lineT(G,v,u, G.er-G.mr, 1); //生光
   this.lecXY(this.lT[6],g);
   this.lT[6] = this.lineT(g,v,u, g.er-g.mr, 1); //生光再算一次
  }
 }
};

//====================================

var rsPL={ //日食批量快速计算器
 Zs:new Array(), Zm:new Array(), Zjd:0, Zdt:1/24, Zn:5,  //日月赤道坐标插值表
 sT:new Array(), //地方日食时间表
 MScoord_init:function(jd) { //准备根数表
  var kn1=Math.floor((jd+8)/29.5306);
  var kn2=Math.floor((this.Zjd+8)/29.5306);
  if(this.Zjd && kn1==kn2) return;

  this.Zjd = jd = XL.MS_aLon_t2( Math.floor((jd+8)/29.5306)*Math.PI*2 )*36525; //低精度的朔(误差10分钟),与食甚的误差1到2小时

  this.Zs.length=0, this.Zm.length=0;
  var i,z=new Array();
  var E,T, a=this.Zs, b=this.Zm;
  for(i=-this.Zn;i<=this.Zn;i++){ //插值点范围不要超过360度(约1个月)
   T=(this.Zjd+i*this.Zdt)/36525;
   ZB.nutation(T);
   E = ZB.hcjj(T) + ZB.dE; //真黄赤交角

   XL.E_coord(T,z,-1,-1,-1);   //地球坐标
   z[0] = rad2mrad(z[0]+Math.PI+ZB.gxc_sunLon(T)+ZB.dL);  z[1] = -z[1] + ZB.gxc_sunLat(T); //补上太阳光行差及章动
   ZB.llrConv( z, E ); //转为赤道坐标
   a[a.length]=z[0], a[a.length]=z[1], a[a.length]=z[2]; //太阳视经,月球赤纬

   XL.M_coord(T,z,-1,-1,-1); //月球坐标
   z[0] = rad2mrad( z[0]+ZB.gxc_moonLon(T)+ZB.dL );  z[1] += ZB.gxc_moonLat(T);  //补上月球光行差及章动
   ZB.llrConv( z, E ); //转为赤道坐标
   b[b.length]=z[0], b[b.length]=z[1], b[b.length]=z[2];
  }
 },
 MScoord:function(jd,z,xt){//日月坐标快速计算(贝赛尔插值法)，xt='sun'时计算太阳
  var i,B, t,p,n;
  if(xt=='sun') B=this.Zs; else B=this.Zm; //日月选择
  t=(jd-this.Zjd)/this.Zdt; //相对于总中心的距离
  p=Math.floor(t+0.5), n=t-p, p=3*(p+this.Zn); //p插值中心,n为插值因子,p再转为插值中心在数据中的位置
  if(p<=0 || p>=B.length) return 'err';
  for(i=0; i<3; i++,p++)
   z[i] = B[p] + ( B[p+3]-B[p-3] + (B[p+3]+B[p-3]-B[p]*2)*n ) * n/2;
  z[0]=rad2mrad(z[0]);
  return 'ok';
 },
 secXY:function(jd,L,fa,high,re){ //日月xy坐标计算。参数：jd是力学时,站点经纬L,fa,海拔high(千米)
  //基本参数计算
  var deltat = JD.deltatT2(jd); //TD-UT
  ZB.nutation(jd/36525);
  var gst= ZB.gst(jd-deltat,deltat) + ZB.dL*Math.cos(ZB.hcjj(jd/36525) + ZB.dE); //真恒星时(不考虑非多项式部分)

  var z=new Array();
  //=======月亮========
  this.MScoord(jd,z,'moon'); re.mCJ=z[0]; re.mCW=z[1]; re.mR=z[2]; //月亮视赤经,月球赤纬
  var mShiJ = rad2rrad(gst - L - z[0]); //得到此刻月亮时角
  ZB.parallax(z, mShiJ,fa, high); re.mCJ2=z[0], re.mCW2=z[1], re.mR2=z[2]; //修正了视差的赤道坐标

  //=======太阳========
  this.MScoord(jd,z,'sun'); re.sCJ=z[0]; re.sCW=z[1]; re.sR=z[2]; //太阳视赤经,月球赤纬
  var sShiJ = rad2rrad(gst - L - z[0]); //得到此刻太阳时角
  ZB.parallax(z,sShiJ,fa,high); re.sCJ2=z[0], re.sCW2=z[1], re.sR2=z[2]; //修正了视差的赤道坐标

  //=======视半径========
  re.mr = 358473400/re.mR2/rad;
  re.sr = 959.63/re.sR2/rad;
  //=======日月赤经纬差转为日面中心直角坐标(用于日食)==============
  re.x = rad2rrad(re.mCJ2-re.sCJ2) * Math.cos((re.mCW2+re.sCW2)/2);
  re.y = re.mCW2-re.sCW2;
  re.t = jd;
 },
 lineT:function(G, v,u, r, n){//已知t1时刻星体位置、速度，求x*x+y*y=r*r时,t的值
  var b=G.y*v-G.x*u, A=u*u+v*v, B=u*b, C=b*b-r*r*v*v, D=B*B-A*C;
  if(D<0) return 0;
  D=Math.sqrt(D); if(!n) D=-D;
  return G.t+((-B+D)/A-G.x)/v;
 },
 secMax:function(jd,L,fa,high){ //日食的食甚计算(jd为近朔的力学时,误差几天不要紧)
  var i;
  for(i=0;i<5;i++) this.sT[i]=0; //分别是:食甚,初亏,复圆,食既,生光
  this.LX='';

  this.MScoord_init(jd);
  jd=this.Zjd; //食甚初始估值为插值表中心时刻(粗朔)

  var G=new Object(), g=new Object();
  this.secXY(jd,L,fa,high,G);
  jd -= G.x/0.2128; //与食甚的误差在20分钟以内

  var i,u,v,dt=60/86400,dt2;
  for(i=0;i<2;i++){
   if( this.secXY(jd,L,fa,high,G)   =='err') return;
   if( this.secXY(jd+dt,L,fa,high,g)=='err') return;
   u = (g.y-G.y)/dt;
   v = (g.x-G.x)/dt;
   dt2 = -(G.y*u+G.x*v)/(u*u+v*v);
   jd += dt2; //极值时间
  }

  //求直线到太阳中心的最小值
  var x=G.x+dt2*v, y=G.y+dt2*u, rmin=Math.sqrt(x*x+y*y);

  if(rmin<=G.mr+G.sr){ //食计算
   this.sT[1] = jd; //食甚
   this.LX='偏';

   this.sT[0] = this.lineT(G,v,u, G.mr+G.sr, 0); //初亏
   for(i=0;i<2;i++) { //初亏再算2次
    this.secXY(this.sT[0],L,fa,high,g);
    this.sT[0] = this.lineT(g,v,u, g.mr+g.sr, 0);
   }
   this.sT[2] = this.lineT(G,v,u, G.mr+G.sr, 1); //复圆
   for(i=0;i<2;i++) { //复圆再算2次
    this.secXY(this.sT[2],L,fa,high,g);
    this.sT[2] = this.lineT(g,v,u, g.mr+g.sr, 1);
   }
  }
  if(rmin<=G.mr-G.sr){ //全食计算
   this.LX='全';
   this.sT[3] = this.lineT(G,v,u, G.mr-G.sr, 0); //初亏
   this.secXY(this.sT[3],L,fa,high,g);
   this.sT[3] = this.lineT(g,v,u, g.mr-g.sr, 0); //初亏再算1次

   this.sT[4] = this.lineT(G,v,u, G.mr-G.sr, 1); //复圆
   this.secXY(this.sT[4],L,fa,high,g);
   this.sT[4] = this.lineT(g,v,u, g.mr-g.sr, 1); //复圆再算1次
  }
  if(rmin<=G.sr-G.mr){ //环食计算
   this.LX='环';
   this.sT[3] = this.lineT(G,v,u, G.sr-G.mr, 0); //初亏
   this.secXY(this.sT[3],L,fa,high,g);
   this.sT[3] = this.lineT(g,v,u, g.sr-g.mr, 0); //初亏再算1次

   this.sT[4] = this.lineT(G,v,u, G.sr-G.mr, 1); //复圆
   this.secXY(this.sT[4],L,fa,high,g);
   this.sT[4] = this.lineT(g,v,u, g.sr-g.mr, 1); //复圆再算1次
  }
 },

 //以下涉及南北界计算
 A:new Array(), B:new Array(), //本半影锥顶点坐标
 P : {S:new Array(), M:new Array(), g:0},//t1时刻的日月坐标,g为恒星时
 Q : {S:new Array(), M:new Array(), g:0},//t2时刻的日月坐标
 V : new Array(), //食界表
 Vc: '', Vb: '', //食中心类型,本影南北距离

 zb0:function(jd){
  //基本参数计算
  var deltat = JD.deltatT2(jd); //TD-UT
  var E=ZB.hcjj(jd/36525);
  ZB.nutation(jd/36525);

  this.P.g = ZB.gst(jd-deltat, deltat) + ZB.dL*Math.cos(E+ZB.dE); //真恒星时(不考虑非多项式部分)
  this.MScoord(jd,this.P.S, 'sun');
  this.MScoord(jd,this.P.M, 'moon');

  var t2=jd+60/86400;
  this.Q.g = ZB.gst(t2-deltat,deltat) + ZB.dL*Math.cos(E+ZB.dE);
  this.MScoord(t2,this.Q.S, 'sun');
  this.MScoord(t2,this.Q.M, 'moon');

  //转为直角坐标
  var z1=new Array(), z2=new Array();
  ZB.llr2xyz(this.P.S[0],this.P.S[1],this.P.S[2]*cs_AU, z1);
  ZB.llr2xyz(this.P.M[0],this.P.M[1],this.P.M[2], z2);

  var k=959.63/358473400*cs_AU, x,y,z; //k为日月半径比
  //本影锥顶点坐标计算
  x = (z1[0]-z2[0])/(1-k)+z2[0];
  y = (z1[1]-z2[1])/(1-k)+z2[1];
  z = (z1[2]-z2[2])/(1-k)+z2[2];
  ZB.xyz2llr(x,y,z,this.A);
  //半影锥顶点坐标计算
  x = (z1[0]-z2[0])/(1+k)+z2[0];
  y = (z1[1]-z2[1])/(1+k)+z2[1];
  z = (z1[2]-z2[2])/(1+k)+z2[2];
  ZB.xyz2llr(x,y,z,this.B);
 },

 zbXY:function(p,L,fa){
  var s=new Array(p.S[0],p.S[1],p.S[2]);
  var m=new Array(p.M[0],p.M[1],p.M[2]);
  ZB.parallax(s, p.g-L-p.S[0],fa, 0); //修正了视差的赤道坐标
  ZB.parallax(m, p.g-L-p.M[0],fa, 0); //修正了视差的赤道坐标
  //=======视半径========
  p.mr = 358473400/m[2]/rad;
  p.sr = 959.63/s[2]/rad;
  //=======日月赤经纬差转为日面中心直角坐标(用于日食)==============
  p.x = rad2rrad(m[0]-s[0]) * Math.cos((m[1]+s[1])/2);
  p.y = m[1]-s[1];
 },
 p2p:function(L,fa,re,fAB,f){ //f取+-1
  var p=this.P, q=this.Q;
  this.zbXY(this.P,L,fa);
  this.zbXY(this.Q,L,fa);

  var u=q.y-p.y, v=q.x-p.x, a=Math.sqrt(u*u+v*v),r=959.63/p.S[2]/rad;

  var W=p.S[1]+f*r*v/a, J=p.S[0]-f*r*u/a/Math.cos((W+p.S[1])/2), R=p.S[2]*cs_AU;

  var A = fAB ? this.A : this.B;

  ZB.line_earth( J,W,R, A[0],A[1],A[2] );
  re.J = rad2rrad(p.g-ZB.le_J);
  re.W = ZB.le_W;
 },
 pp0:function(re){ //食中心点计算
  var p=this.P;
  ZB.line_earth( p.M[0],p.M[1],p.M[2], p.S[0],p.S[1],p.S[2]*cs_AU );
  re.J = rad2rrad(p.g-ZB.le_J);
  re.W = ZB.le_W; //无解返回值是100
  
  if(re.W==100) { re.c = ''; return; }
  re.c='全';
  this.zbXY(p,re.J,re.W);
  if(p.sr>p.mr) re.c='环';
 },
 nbj:function(jd){ //南北界计算
  this.MScoord_init(jd);
  var i, G=new Object(), V=this.V;
  for(i=0;i<10;i++) V[i]=100; this.Vc='',this.Vb=''; //返回初始化,纬度值为100表示无角,经度100也是无解,但在以下程序中经度会被转为-PI到+PI
  if( Math.abs(jd-this.Zjd)>this.Zdt*this.Zn ) return; //时刻值无效

  this.zb0(jd);
  this.pp0(G); V[0]=G.J, V[1]=G.W, this.Vc=G.c; //食中心

  G.J=G.W=0; for(i=0;i<2;i++) this.p2p(G.J,G.W,G,1, 1); V[2]=G.J, V[3]=G.W; //本影北界,环食为南界(本影区之内,变差u,v基本不变,所以计算两次足够)
  G.J=G.W=0; for(i=0;i<2;i++) this.p2p(G.J,G.W,G,1,-1); V[4]=G.J, V[5]=G.W; //本影南界,环食为北界
  G.J=G.W=0; for(i=0;i<3;i++) this.p2p(G.J,G.W,G,0,-1); V[6]=G.J, V[7]=G.W; //半影北界
  G.J=G.W=0; for(i=0;i<3;i++) this.p2p(G.J,G.W,G,0, 1); V[8]=G.J, V[9]=G.W; //半影南界

  if(V[3]!=100&&V[5]!=100){ //粗算本影南北距离
    var x=(V[2]-V[4])*Math.cos((V[3]+V[5])/2), y=V[3]-V[5];
    this.Vb = (cs_rEarA*Math.sqrt(x*x+y*y)).toFixed(0)+'千米';
  }
  //Cal_zdzb.innerHTML=(V[8]*180/Math.PI).toFixed(5)+' '+(V[9]*180/Math.PI).toFixed(5);
 }
};