eph.js

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

/****************************************
以下是天文计算部分,包含有：
常数 rad   : 每弧度的角秒数
常数 J2000 : 2000年1月1日 12:00:00 的儒略日数
常数 pi2   : 圆周率的2倍,即2*3.14159...
rad2str()  : 弧度转为"度分秒"或"时分秒"
rad2str2() : 弧度转为"度分"
rad2mrad() : 将超出0到360度的角度转为0到360度的解度
mod2(a,b)  : 临界余数(a与最近的整倍数b相差的距离)
int2(a)    : 取小等于a的最大整数
物件 JD    : 公历——儒略日——时间物件
物件 ZB    : 坐标变换物件，包含坐标旋转、章动、日月光行差(时)、视差、大气折射、恒星时的计算
物件 XL    : 日月黄道平分点坐标、视坐标、速度、已知经度反求时间等方面的计算
物件 SZJ   : 用来计算日月的升起、中天、降落
注意，上述函数或物件是纯天文学的，根据实际需要组合使用可以得到所需要的各种日月坐标，计算精度及计算速度也是可以根据需要有效控制的。
*****************************************/
Number.prototype.toFixed=function(len){ //IE的toFixed()有Bug,所以重写
  var n=this,f='';
  if(n<0) n = -n, f = '-';
  var n1 = Math.floor(n); //整数部分
  var n2 = n-n1;          //小数部分
  var p = Math.pow(10,len);
  n2 = Math.round( n2*p ); //移位并四舍五入
  if(n2>=p) n1++,n2-=p; //进位
  n2 = (p+n2)+'';
  n2 = n2.substr(1,len);
  if(!len) return f+n1;
  return f+n1+'.'+n2;
};

cs_rEar = 6378.1366; //地球赤道半径(千米)
cs_rEarA= 0.99834*cs_rEar; //平均半径
cs_AU   = 1.49597870691e8; //天文单位长度(千米)
cs_sinP = cs_rEar/cs_AU;   //sin(太阳视差)
cs_PI   = Math.asin(cs_sinP); //太阳视差

rad  = 180*3600/Math.PI,
pi2  = Math.PI*2;
J2000 = 2451545;
function rad2str(d,tim){ //将弧度转为字串
 //tim=0输出格式示例: -23°59' 48.23"
 //tim=1输出格式示例:  18h 29m 44.52s
 var s=" ";
 var w1="°",w2="'",w3='"';
 if(d<0)  d=-d,s='-';
 if(tim){ d*=12/Math.PI; w1="h ",w2="m",w3="s"; }
 else     d*=180/Math.PI;
 var a=Math.floor(d); d=(d-a)*60;
 var b=Math.floor(d); d=(d-b)*60;
 var c=Math.floor(d); d=(d-c)*100;
 d=Math.floor(d+0.5);
 if(d>=100) d-=100, c++;
 if(c>=60)  c-=60,  b++;
 if(b>=60)  b-=60,  a++;
 a="   "+a, b="0"+b, c="0"+c, d="0"+d;
 s+=a.substr(a.length-3,3)+w1;
 s+=b.substr(b.length-2,2)+w2;
 s+=c.substr(c.length-2,2)+".";
 s+=d.substr(d.length-2,2)+w3;
 return s;
}
function rad2str2(d){ //将弧度转为字串,精确到分
 //输出格式示例: -23°59'
 var s="+";
 var w1="°",w2="'",w3='"';
 if(d<0)  d=-d,s='-';
 d*=180/Math.PI;
 var a=Math.floor(d);
 var b=Math.floor((d-a)*60+0.5);
 if(b>=60)  b-=60,  a++;
 a="   "+a, b="0"+b;
 s+=a.substr(a.length-3,3)+w1;
 s+=b.substr(b.length-2,2)+w2;
 return s;
}
function m2fm(v,fx,fs){//秒转为分秒,fx为小数点位数,fs为1转为"分秒"格式否则转为"角分秒"格式
 var gn='';
 if(v<0) v=-v, gn='-';
 var f=int2(v/60), m=v-f*60;
 if(fs) return gn + f + '分' +m.toFixed(fx) + '秒';
 else   return gn + f + "'"  +m.toFixed(fx) + '"' ;
}
function rad2mrad(v){ //对超过0-2PI的角度转为0-2PI
  v=v % (2*Math.PI);
  if(v<0) return v+2*Math.PI;
  return v;
}
function rad2rrad(v){//对超过-PI到PI的角度转为-PI到PI
  v=v % (2*Math.PI);
  if(v<=-Math.PI) return v+2*Math.PI;
  if(v>Math.PI) return v-2*Math.PI;
  return v;
}
function mod2(a,b){ //临界余数(a与最近的整倍数b相差的距离)
  var c=a/b;
  c -= Math.floor(c);
  if(c>0.5) c-=1;
  return c*b;
}
function int2(v){ return Math.floor(v); }  //取整数部分

var JD={ //日期元件
  Y:2000, M:1, D:1, h:12, m:0, s:0,
  dts:new Array( // TD - UT1 计算表
   -4000,108371.7,-13036.80,392.000, 0.0000, -500, 17201.0,  -627.82, 16.170,-0.3413,
    -150, 12200.6,  -346.41,  5.403,-0.1593,  150,  9113.8,  -328.13, -1.647, 0.0377,
     500,  5707.5,  -391.41,  0.915, 0.3145,  900,  2203.4,  -283.45, 13.034,-0.1778,
    1300,   490.1,   -57.35,  2.085,-0.0072, 1600,   120.0,    -9.81, -1.532, 0.1403,
    1700,    10.2,    -0.91,  0.510,-0.0370, 1800,    13.4,    -0.72,  0.202,-0.0193,
    1830,     7.8,    -1.81,  0.416,-0.0247, 1860,     8.3,    -0.13, -0.406, 0.0292,
    1880,    -5.4,     0.32, -0.183, 0.0173, 1900,    -2.3,     2.06,  0.169,-0.0135,
    1920,    21.2,     1.69, -0.304, 0.0167, 1940,    24.2,     1.22, -0.064, 0.0031,
    1960,    33.2,     0.51,  0.231,-0.0109, 1980,    51.0,     1.29, -0.026, 0.0032,
    2000,    63.87,    0.1,   0,     0,      2005),
  deltatExt:function(y,jsd){ var dy=(y-1820)/100; return -20+jsd*dy*dy; }, //二次曲线外推
  deltatT:function(y){ //计算世界时与原子时之差,传入年
   if(y>=2005){
     //sd是2005年之后几年（一值到y1年）的速度估计。
     //sjd是y1年之后的加速度估计。瑞士星历表jsd=31,NASA网站jsd=32,skmap的jsd=29
     var y1=2014, sd=0.4, jsd=31;
     if(y<=y1) return 64.7 + (y-2005) *sd; //直线外推
     var v = this.deltatExt(y,jsd);        //二次曲线外推
     var dv= this.deltatExt(y1,jsd) - ( 64.7+(y1-2005)*sd ); //y1年的二次外推与直线外推的差
     if(y<y1+100  ) v -= dv*(y1+100-y)/100;
     return v; 
   }
   var i,d=this.dts;
   for(i=0;i<d.length;i+=5) if(y<d[i+5]) break;
   var t1=(y-d[i])/(d[i+5]-d[i])*10, t2=t1*t1, t3=t2*t1;
   return d[i+1] +d[i+2]*t1 +d[i+3]*t2 +d[i+4]*t3;
  },
  deltatT2:function(t){ //传入儒略日(J2000起算),计算TD-UT(单位:日)
    return this.deltatT(t/365.2425+2000)/86400.0;
  },
  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)+0.01) + int2(30.6001*(m+1))+d+n - 1524.5;
  },
  toJD:function(){ return this.JD( this.Y, this.M, this.D+((this.s/60+this.m)/60+this.h)/24 ); }, //公历转儒略日
  setFromJD:function(jd){ //儒略日数转公历
   jd+=0.5;
   var A=int2(jd), F=jd-A, D;  //取得日数的整数部份A及小数部分F
   if(A>=2299161) D=int2((A-1867216.25)/36524.25),A+=1+D-int2(D/4);
   A     +=1524; //向前移4年零2个月
   this.Y =int2((A-122.1)/365.25);//年
   D      =A-int2(365.25*this.Y); //去除整年日数后余下日数
   this.M =int2(D/30.6001);       //月数
   this.D =D-int2(this.M*30.6001);//去除整月日数后余下日数
   this.Y-=4716; this.M--;
   if(this.M>12) this.M-=12;
   if(this.M<=2) this.Y++;
   //日的小数转为时分秒
   F*=24; this.h=int2(F); F-=this.h;
   F*=60; this.m=int2(F); F-=this.m;
   F*=60; this.s=F;
  },
  toStr:function(){ //日期转为串
   var Y="     "+this.Y,M="0"+this.M, D="0"+this.D;
   var h=this.h,m=this.m,s=int2(this.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;
  },
  setFromJD_str:function(jd){ this.setFromJD(jd); return this.toStr(); },
  timeStr:function(jd){ //提取jd中的时间(去除日期)
   var h,m,s;
   jd+=0.5; jd = (jd - int2(jd));
   jd*=24; h = int2(jd); jd-=h;
   jd*=60; m = int2(jd); jd-=m;
   jd*=60; s = int2(jd+0.5);
   if(s>=60) s-=60,m++;
   if(m>=60) m-=60,h++;
   h="0"+h; m="0"+m; s="0"+s;
   return h.substr(h.length-2,2)+':'+m.substr(m.length-2,2)+':'+s.substr(s.length-2,2);
  },
  //星期相关
  Weeks : new Array('日','一','二','三','四','五','六','七'),
  getWeek:function(jd){ return int2(jd+1.5)%7; }, //星期计算
  nnweek:function(y,m,n,w){ //求y年m月的第n个星期w的儒略日数
   var jd = JD.JD(y,m,1.5); //月首儒略日
   var w0 = (jd+1)%7;       //月首的星期
   var r  = jd-w0+7*n+w;    //jd-w0+7*n是和n个星期0,起算下本月第一行的星期日(可能落在上一月)。加w后为第n个星期w
   if(w>=w0) r-=7; //第1个星期w可能落在上个月,造成多算1周,所以考虑减1周
   if(n==5){
    m++; if(m>12) m=1, y++;   //下个月
    if(r>=JD.JD(y,m,1.5)) r-=7; //r跑到下个月则减1周
   }
   return r;
  }
};

var ZB={ //坐标类
 llrConv:function(JW,E){ //球面坐标旋转
  //黄道赤道坐标变换,赤到黄E取负
  var sinE =Math.sin(E),    cosE =Math.cos(E);
  var sinJ =Math.sin(JW[0]),cosJ =Math.cos(JW[0]);
  var sinW =Math.sin(JW[1]),cosW =Math.cos(JW[1]), tanW=Math.tan(JW[1]);
  JW[0]=Math.atan2( sinJ*cosE - tanW*sinE, cosJ );
  JW[1]=Math.asin ( cosE*sinW + sinE*cosW*sinJ  );
  JW[0]=rad2mrad(JW[0]);
 },
 llr2xyz:function(J,W,R,z){ //球面转直角坐标
  z[0]=R*Math.cos(W)*Math.cos(J);
  z[1]=R*Math.cos(W)*Math.sin(J);
  z[2]=R*Math.sin(W);
 },
 xyz2llr:function(x,y,z,q){ //直角坐标转球
  q[2] = Math.sqrt(x*x+y*y+z*z);
  q[1] = Math.asin(z/q[2]);
  q[0] = rad2mrad( Math.atan2(y,x) );
 },
 nutB:new Array(
   2.1824,  -33.75705, 36e-6,-1720,920,
   3.5069, 1256.66393, 11e-6,-132, 57,
   1.3375,16799.4182, -51e-6, -23, 10,
   4.3649,  -67.5141,  72e-6,  21, -9,
   0.04,   -628.302,   0,     -14,  0,
   2.36,   8328.691,   0,       7,  0,
   3.46,   1884.966,   0,      -5,  2,
   5.44,  16833.175,   0,      -4,  2,
   3.69,  25128.110,   0,      -3,  0,
   3.55,    628.362,   0,       2,  0),
 nutation:function(t){ //章动计算,t是世纪数
  var i,c,a, t2=t*t, B=this.nutB, dL=0,dE=0;
  for(i=0;i<B.length;i+=5){
   c = B[i]+B[i+1]*t+B[i+2]*t2;
   if(i==0) a=-1.742*t; else a=0;
   dL+=(B[i+3]+a)*Math.sin(c);
   dE+= B[i+4]   *Math.cos(c);
  }
  this.dL=dL/100/rad;  //黄经章动
  this.dE=dE/100/rad;  //交角章动
 },
 nutationLon:function(t){ //只计算黄经章动
  var i,a, t2=t*t, dL=0, B=this.nutB;
  for(i=0;i<B.length;i+=5){
   if(i==0) a=-1.742*t; else a=0;
   dL += (B[i+3]+a) * Math.sin( B[i]+B[i+1]*t+B[i+2]*t2 );
  }
  return dL/100/rad;
 },
 hcjj:function(t){ //返回黄赤交角,t是世纪数
  var t2=t*t, t3=t2*t,t4=t3*t;
  return (84381.4088 -46.836051*t -0.0001667*t2 -0.00199911*t3-0.000000523*t4)/rad;
 },
 j1_j2:function(J1,W1,J2,W2){ //求角度差(未测试)
   var dJ=rad2rrad(J1-J2), dW=W1-W2;
   if(Math.abs(dJ)<1/1000 && Math.abs(dW)<1/1000){
     dJ*=Math.cos((W1+W2)/2);
     return Math.sqrt(dJ*dJ+dW*dW);
   }
   return Math.acos( Math.sin(W1)*Math.sin(W2) + Math.cos(W1)*Math.cos(W2)*Math.cos(dJ) );
 },
 gst:function(T,dt){ //传入T是2000年首起算的日数(UT),dt是deltatT(日),精度要求不高时dt可取值为0
  //返回格林尼治恒星时(不含赤经章动及非多项式部分),即格林尼治子午圈的平春风点起算的赤经
  var t=(T+dt)/36525,t2=t*t, t3=t2*t, t4=t3*t;
  return pi2*(0.7790572732640 + 1.00273781191135448*T) //严格说这里的T是UT,下一行的t是力学时(世纪数)
      + (0.014506 + 4612.15739966*t + 1.39667721*t2 - 0.00009344*t3 + 0.00001882*t4)/rad;
 },
 gxc_sunLon:function(t){ //太阳光行差,t是世纪数
  var v =-0.043126+ 628.301955*t -0.000002732*t*t; //平近点角
  var e = 0.016708634-0.000042037*t-0.0000001267*t*t;
  return  ( -20.49552 * (1+e*Math.cos(v)) )/rad; //黄经光行差
 },
 gxc_sunLat:function(t) { return 0;       }, //黄纬光行差
 gxc_moonLon:function(t){ return -3.4E-6; }, //月球经度光行差,误差0.07"
 gxc_moonLat:function(t){ //月球纬度光行差,误差0.006"
  return 0.063*Math.sin(0.057+8433.4662*t+0.000064*t*t)/rad;
 },
 AR:function(ho){ return -0.0002909/Math.tan( ho+0.002227/(ho+0.07679) ); }, //大气折射,ho是视高度
 AR2:function(h){ return  0.0002967/Math.tan( h +0.003138/( h+0.08919) ); }, //大气折射,h是真高度
 parallax:function(z,H,fa,high){ //视差修正