eph.js

返回 第一章 关于寿星万年历 第二章 公历与儒略日 第三章 回历计算 第四章 天文坐标 第五章 近现代农历的天文学性质及其规则 第六章 离散序列的直线拟合算法

new Array(//MR1
0.514,4.16,14914.4523,-0.6,6,0,0.382,1.80,6585.7609,-2.2,-19,0,
0.327,2.40,7700.3895,1.5,25,0,0.264,5.45,8956.9934,1.5,25,0,
0.123,3.10,628.302,0,0,0,0.078,0.93,16171.056,-0.7,6,0,
0.061,4.86,7842.365,-2.2,-19,0,0.050,4.2,14286.150,-1,6,0,
0.042,5.2,8399.679,0,0,0,0.032,0.2,23243.144,1,31,0,
0.025,2.6,-1742.931,-4,-44,0,0.025,1.8,5957.459,-2,-19,0,
0.018,4.8,16029.081,3,50,0,0.014,1.5,17285.68,3,50,0,
0.014,1.0,15542.75,-1,0,0,0.013,5.0,8326.39,3,50,0,
0.012,4.8,8470.67,-2,0,0,0.012,2.8,8330.99,0,0,0,
0.011,2.4,7072.09,2,0,0,0.010,5.9,22128.52,-3,0,0),
new Array(//MR2
0.0015,4.2,14914.452,-1,6,0,0.0011,1.8,6585.761,-2,-19,0,
0.0009,2.4,7700.389,2,25,0,0.0008,5.5,8956.993,1,25,0)
),

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

//星历函数(日月球面坐标计算)
Enn:function(ob,t,n){ //计算E_L0或E_L1或E_L2等
  var i,j,F,N,v=0,tn=1,c;
  if(ob==this.EL){
    var t2=t*t, t3=t2*t, t4=t3*t, t5=t4*t; //千年数的各次方
    v += 1753469512 + 6283319653318*t + 529674*t2 + 432*t3 - 1124*t4 - 9*t5 + 630 * Math.cos(6+3*t); //地球平黄经(已拟合DE406)
  }
  n*=3; if(n<0) n = ob[0].length;
  for(i=0;i<ob.length;i++,tn*=t){
   F = ob[i];
   N = int2( n*F.length/ob[0].length+0.5 );  if(i) N+=3;  if(N >= F.length) N=F.length;
   for(j=0,c=0;j<N;j+=3) c+=F[j]*Math.cos(F[j+1] +t*F[j+2]);
   v += c*tn;
  }
  return v/1000000000;
},
E_coord:function(t,re,n){ //返回地球坐标,t为世纪数
  t/=10;
  re[0]= this.Enn( this.EL, t, n);
  re[1]= this.Enn( this.EB, t,-1);
  re[2]= this.Enn( this.ER, t,-1);
},

E_Lon:function(t,n){ return this.Enn(this.EL,t/10,n); }, //地球经度计算,返回Date分点黄经,传入世纪数、取项数

Mnn:function(ob,t,n){ //计算ML0或ML1或ML2
  var i,j,F,N,v=0,tn=1,c;
  var t2=t*t,t3=t2*t,t4=t3*t,t5=t4*t,tx=t-10;
  if(ob==this.ML){
    v += (3.81034409 + 8399.684730072*t -3.319e-05*t2 + 3.11e-08*t3 - 2.033e-10*t4)*rad; //月球平黄经(弧度)
    v += 5028.792262*t + 1.1124406*t2 + 0.00007699*t3 - 0.000023479*t4 -0.0000000178*t5;  //岁差(角秒)
    if(tx>0) v += -0.866 +1.43*tx +0.054*tx*tx; //对公元3000年至公元5000年的拟合,最大误差小于10角秒
  }
  t2/=1e4,t3/=1e8,t4/=1e8;
  n*=6; if(n<0) n = ob[0].length;
  for(i=0;i<ob.length;i++,tn*=t){
   F=ob[i];
   N = int2( n*F.length/ob[0].length+0.5 );  if(i) N+=6;  if(N >= F.length) N=F.length;
   for(j=0,c=0;j<N;j+=6) c+=F[j]*Math.cos(F[j+1] +t*F[j+2] +t2*F[j+3] +t3*F[j+4] +t4*F[j+5]);
   v += c*tn;
  }
  if(ob!=this.MR) v/=rad;
  return v;
},
M_coord:function(t,re,n1,n2,n3){ //返回月球坐标,n1,n2,n3为各坐标所取的项数
  re[0] = this.Mnn( this.ML, t, n1 );
  re[1] = this.Mnn( this.MB, t, n2 );
  re[2] = this.Mnn( this.MR, t, n3 );
},
M_Lon:function(t,n){ return this.Mnn(this.ML,t,n); }, //月球经度计算,返回Date分点黄经,传入世纪数,n是项数比例

//=========================
E_v:function(t){ //地球速度,t是世纪数,误差小于万分3
 var f=628.307585*t;
 return 628.332 +21 *Math.sin(1.527+f) +0.44 *Math.sin(1.48+f*2) 
         +0.129*Math.sin(5.82+f)*t +0.00055*Math.sin(4.21+f)*t*t;
},
M_v:function(t){ //月球速度计算,传入世经数
  var v = 8399.71 - 914*Math.sin( 0.7848 + 8328.691425*t + 0.0001523*t*t ); //误差小于5%
  v-=179*Math.sin( 2.543  +15542.7543*t )  //误差小于0.3%
    +160*Math.sin( 0.1874 + 7214.0629*t )
    +62 *Math.sin( 3.14   +16657.3828*t )
    +34 *Math.sin( 4.827  +16866.9323*t )
    +22 *Math.sin( 4.9    +23871.4457*t )
    +12 *Math.sin( 2.59   +14914.4523*t )
    +7  *Math.sin( 0.23   + 6585.7609*t )
    +5  *Math.sin( 0.9    +25195.624 *t )
    +5  *Math.sin( 2.32   - 7700.3895*t )
    +5  *Math.sin( 3.88   + 8956.9934*t )
    +5  *Math.sin( 0.49   + 7771.3771*t );
   return v;
},

//=========================
MS_aLon:function(t,Mn,Sn){ //月日视黄经的差值
  return this.M_Lon(t,Mn) + ZB.gxc_moonLon(t) - ( this.E_Lon(t,Sn) + ZB.gxc_sunLon(t) + Math.PI );
},
S_aLon:function(t,n){  //太阳视黄经
 return this.E_Lon(t,n) + ZB.nutationLon(t) + ZB.gxc_sunLon(t) + Math.PI; //注意，这里的章动计算很耗时
},

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

E_Lon_t:function(W){ //已知地球真黄经求时间
  var t,v= 628.3319653318;
  t  = ( W - 1.75347          )/v; v=this.E_v(t);   //v的精度0.03%，详见原文
  t += ( W - this.E_Lon(t,10) )/v; v=this.E_v(t);   //再算一次v有助于提高精度,不算也可以
  t += ( W - this.E_Lon(t,-1) )/v;
  return t;
},
M_Lon_t:function(W){ //已知真月球黄经求时间
  var t,v= 8399.70911033384;
  t  = ( W - 3.81034          )/v;
  t += ( W - this.M_Lon(t,3 ) )/v; v=this.M_v(t);  //v的精度0.5%，详见原文
  t += ( W - this.M_Lon(t,20) )/v;
  t += ( W - this.M_Lon(t,-1) )/v;
  return t;
},
MS_aLon_t:function(W){ //已知月日视黄经差求时间
  var t,v= 7771.37714500204;
  t  = ( W + 1.08472               )/v;
  t += ( W - this.MS_aLon(t, 3, 3) )/v; v=this.M_v(t)-this.E_v(t);  //v的精度0.5%，详见原文
  t += ( W - this.MS_aLon(t,20,10) )/v;
  t += ( W - this.MS_aLon(t,-1,60) )/v;
  return t;
},
S_aLon_t:function(W){ //已知太阳视黄经反求时间
  var t,v= 628.3319653318;
  t  = ( W - 1.75347-Math.PI   )/v; v=this.E_v(t); //v的精度0.03%，详见原文
  t += ( W - this.S_aLon(t,10) )/v; v=this.E_v(t); //再算一次v有助于提高精度,不算也可以
  t += ( W - this.S_aLon(t,-1) )/v;
  return t;
},
/****
MS_aLon_t1:function(W){ //已知月日视黄经差求时间,高速低精度,误差不超过40秒
  var t,v = 7771.37714500204;
  t  = ( W + 1.08472               )/v;
  t += ( W - this.MS_aLon(t, 3, 3) )/v;  v=this.M_v(t)-this.E_v(t);  //v的精度0.5%，详见原文
  t += ( W - this.MS_aLon(t,50,20) )/v;
  return t;
},
S_aLon_t1:function(W){ //已知太阳视黄经反求时间,高速低精度,最大误差不超过50秒,平均误差15秒
  var t,v= 628.3319653318;
  t  = ( W - 1.75347-Math.PI   )/v; v = 628.332 + 21*Math.sin( 1.527+628.307585*t );
  t += ( W - this.S_aLon(t,3) )/v;
  t += ( W - this.S_aLon(t,40))/v;
  return t;
},
****/
MS_aLon_t2:function(W){ //已知月日视黄经差求时间,高速低精度,误差不超过600秒
  var t,v = 7771.37714500204;
  t  = ( W + 1.08472 )/v;
  var L,t2 = t*t;
  t -= ( -0.00003309*t2 + 0.10976*Math.cos( 0.784758 + 8328.6914246*t + 0.000152292*t2 ) + 0.02224 *Math.cos( 0.18740  + 7214.0628654*t - 0.00021848 *t2 ) - 0.03342 *Math.cos( 4.669257 + 628.307585*t ) )/v;
  L = this.M_Lon(t,20) - (4.8950632+ 628.3319653318*t + 0.000005297*t*t + 0.0334166*Math.cos(4.669257+628.307585*t) + 0.0002061*Math.cos(2.67823+628.307585*t)*t + 0.000349*Math.cos(4.6261+1256.61517*t) - 20.5/rad);
  v = 7771.38 - 914*Math.sin( 0.7848 + 8328.691425*t + 0.0001523*t*t ) - 179*Math.sin( 2.543 + 15542.7543*t ) - 160*Math.sin( 0.1874 + 7214.0629*t );
  t += ( W - L )/v;
  return t;
},
S_aLon_t2:function(W){ //已知太阳视黄经反求时间,高速低精度,最大误差不超过600秒
  var t,L,v= 628.3319653318;
  t =  ( W - 1.75347-Math.PI )/v;
  t -= ( 0.000005297*t*t + 0.0334166 * Math.cos( 4.669257 + 628.307585*t) + 0.0002061 * Math.cos( 2.67823  + 628.307585*t)*t )/v;
  t += ( W - this.E_Lon(t,8) - Math.PI + (20.5+17.2*Math.sin(2.1824-33.75705*t))/rad )/v;
  return t;
},
moonIll:function(t){ //月亮被照亮部分的比例
 var t2 = t*t, t3 =t2*t, t4 = t3*t;
 var D,M,m,a, dm=Math.PI/180;
 D = 297.8502042 + 445267.1115168*t - 0.0016300*t2 + t3/545868 - t4/113065000; //日月距角
 M = 357.5291092 +  35999.0502909*t - 0.0001536*t2 + t3/24490000;              //太阳平近点
 m = 134.9634114 + 477198.8676313*t + 0.0089970*t2 + t3/69699 - t4/14712000;   //月亮平近点
 D*=dm, M*=dm, m*=dm;
 a = Math.PI - D + (-6.289*Math.sin(m) +2.100*Math.sin(M) -1.274*Math.sin(D*2-m) -0.658*Math.sin(D*2) -0.214*Math.sin(m*2) -0.110*Math.sin(D))*dm;
 return (1+Math.cos(a))/2;
},
moonRad:function(r,h){ //转入地平纬度及地月质心距离,返回站心视半径(角秒)
 return 358473400/r*(1+Math.sin(h)*6378.14/r);
}
};

var SZJ={//日月的升中天降,不考虑气温和气压的影响
 L  : 0,  //站点地理经度,向东测量为负
 fa : 0,  //站点地理纬度
 dt : 0,  //TD-UT
 E  : 0.409092614, //黄赤交角
 getH:function(h,w){ //h地平纬度,w赤纬,把回时角
  var c = ( Math.sin(h) - Math.sin(this.fa)*Math.sin(w) ) / Math.cos(this.fa)/Math.cos(w);
  if(Math.abs(c)>1) return Math.PI;
  return Math.acos(c);
 },

 Mcoord:function(jd,H0,z){ //章动同时影响恒星时和天体坐标,所以不计算章动。返回时角及赤经纬
   XL.M_coord( (jd+this.dt)/36525, z, 30,20,8 ); //低精度月亮赤经纬
   ZB.llrConv( z, this.E ); //转为赤道坐标
   z.H = rad2mrad(ZB.gst(jd,this.dt) - this.L - z[0]);   if( z.H>Math.PI ) z.H -= pi2; //得到此刻天体时角
   if(H0) z.H0 = this.getH( 0.7275*6378.14/z[2]-34*60/rad, z[1] ); //升起对应的时角
 },
 Mt : function(jd){ //月亮到中升降时刻计算,传入jd含义与St()函数相同
  this.dt = JD.deltatT2(jd);
  this.E  = ZB.hcjj(jd/36525);
  jd -= mod2(0.1726222 + 0.966136808032357*jd - 0.0366*this.dt - this.L/pi2, 1); //查找最靠近当日中午的月上中天,mod2的第1参数为本地时角近似值

  var r = new Array(), sv = pi2*0.966;
  r.z = r.s = r.j = r.c = r.h = jd;
  this.Mcoord(jd,1,r); //月亮坐标
  r.s += (-r.H0 - r.H )/sv;
  r.j += ( r.H0 - r.H )/sv;
  r.z += (    0 - r.H )/sv;
  this.Mcoord(r.s,1,r);  r.s += ( -r.H0 - r.H )/sv;
  this.Mcoord(r.j,1,r);  r.j += ( +r.H0 - r.H )/sv;
  this.Mcoord(r.z,0,r);  r.z += (     0 - r.H )/sv;
  return r;
 },

 Scoord:function(jd,H0,H1,z){ //章动同时影响恒星时和天体坐标,所以不计算章动。返回时角及赤经纬
   z[0] = XL.E_Lon( (jd+this.dt)/36525, 5 ) + Math.PI - 20.5/rad;  //太阳坐标(修正了光行差)
   z[1] = 0; z[2]=1;
   ZB.llrConv( z, this.E ); //转为赤道坐标
   z.H = rad2mrad(ZB.gst(jd,this.dt) - this.L - z[0]);   if( z.H>Math.PI ) z.H -= pi2; //得到此刻天体时角
   if(H0) z.H0 = this.getH(-50*60/rad,z[1]);  //地平以下50分
   if(H1) z.H1 = this.getH(-Math.PI/30,z[1]); //地平以下6度
 },
 St : function(jd){ //太阳到中升降时刻计算,传入jd是当地中午12点时间对应的2000年首起算的格林尼治时间UT
  this.dt = JD.deltatT2(jd);
  this.E  = ZB.hcjj(jd/36525);
  jd -= mod2(jd - this.L/pi2, 1); //查找最靠近当日中午的日上中天,mod2的第1参数为本地时角近似值

  var r = new Array(), sv = pi2;
  r.z = r.s = r.j = r.c = r.h = jd;
  this.Scoord(jd,1,1,r); //太阳坐标
  r.s += (-r.H0 - r.H )/sv; //升起
  r.j += ( r.H0 - r.H )/sv; //降落
  r.c += (-r.H1 - r.H )/sv; //民用晨
  r.h += ( r.H1 - r.H )/sv; //民用昏
  r.z += (    0 - r.H )/sv; //中天
  this.Scoord(r.s,1,0,r);  r.s += ( -r.H0 - r.H )/sv;
  this.Scoord(r.j,1,0,r);  r.j += ( +r.H0 - r.H )/sv;
  this.Scoord(r.c,0,1,r);  r.c += ( -r.H1 - r.H )/sv;
  this.Scoord(r.h,0,1,r);  r.h += ( +r.H1 - r.H )/sv;
  this.Scoord(r.z,0,0,r);  r.z += (     0 - r.H )/sv;
  return r;
 },

 rts:new Array(),//多天的升中降
 calcRTS:function(jd,n,Jdl,Wdl,sq){ //多天升中降计算,jd是当地起始略日(中午时刻),sq是时区
  var i,c,r;
  if(!this.rts.length) { for(var i=0;i<31;i++) this.rts[i] = new Array(); }
  this.L = Jdl, this.fa = Wdl, sq/=24; //设置站点参数
  for(i=0;i<n;i++){ r=this.rts[i];  r.Ms=r.Mz=r.Mj=""; }
  for(i=-1;i<=n;i++){
   if(i>=0&&i<n){ //太阳
    r = SZJ.St(jd+i+sq);
    this.rts[i].s = JD.timeStr(r.s-sq); //升
    this.rts[i].z = JD.timeStr(r.z-sq); //中
    this.rts[i].j = JD.timeStr(r.j-sq); //降
    this.rts[i].c = JD.timeStr(r.c-sq); //晨
    this.rts[i].h = JD.timeStr(r.h-sq); //昏
    this.rts[i].ch = JD.timeStr(r.h-r.c-0.5); //光照时间,timeStr()内部+0.5,所以这里补上-0.5
    this.rts[i].sj = JD.timeStr(r.j-r.s-0.5); //昼长
   }
   r = SZJ.Mt(jd+i+sq); //月亮
   c=int2(r.s-sq+0.5)-jd;  if(c>=0&&c<n) this.rts[c].Ms = JD.timeStr(r.s-sq);
   c=int2(r.z-sq+0.5)-jd;  if(c>=0&&c<n) this.rts[c].Mz = JD.timeStr(r.z-sq);
   c=int2(r.j-sq+0.5)-jd;  if(c>=0&&c<n) this.rts[c].Mj = JD.timeStr(r.j-sq);
  }
  this.rts.dn = n;
 }
}