/****************************************
以下是天文計算部分,包含有：
常數(shù) rad   : 每弧度的角秒數(shù)
常數(shù) J2000 : 2000年1月1日 12:00:00 的儒略日數(shù)
常數(shù) pi2   : 圓周率的2倍,即2*3.14159...
rad2str()  : 弧度轉(zhuǎn)為"度分秒"或"時分秒"
rad2str2() : 弧度轉(zhuǎn)為"度分"
rad2mrad() : 將超出0到360度的角度轉(zhuǎn)為0到360度的解度
mod2(a,b)  : 臨界余數(shù)(a與最近的整倍數(shù)b相差的距離)
int2(a)    : 取小等于a的最大整數(shù)
物件 JD    : 公歷——儒略日——時間物件
物件 ZB    : 坐標變換物件，包含坐標旋轉(zhuǎn)、章動、日月光行差(時)、視差、大氣折射、恒星時的計算
物件 XL    : 日月黃道平分點坐標、視坐標、速度、已知經(jīng)度反求時間等方面的計算
物件 SZJ   : 用來計算日月的升起、中天、降落
注意，上述函數(shù)或物件是純天文學(xué)的，根據(jù)實際需要組合使用可以得到所需要的各種日月坐標，計算精度及計算速度也是可以根據(jù)需要有效控制的。
*****************************************/
Number.prototype.toFixed=function(len){ //IE的toFixed()有Bug,所以重寫
  var n=this,f='';
  if(n<0) n = -n, f = '-';
  var n1 = Math.floor(n); //整數(shù)部分
  var n2 = n-n1;          //小數(shù)部分
  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){ //將弧度轉(zhuǎn)為字串
 //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){ //將弧度轉(zhuǎn)為字串,精確到分
 //輸出格式示例: -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){//秒轉(zhuǎn)為分秒,fx為小數(shù)點位數(shù),fs為1轉(zhuǎn)為"分秒"格式否則轉(zhuǎn)為"角分秒"格式
 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的角度轉(zhuǎn)為0-2PI
  v=v % (2*Math.PI);
  if(v<0) return v+2*Math.PI;
  return v;
}
function rad2rrad(v){//對超過-PI到PI的角度轉(zhuǎn)為-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){ //臨界余數(shù)(a與最近的整倍數(shù)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); }  //取整數(shù)部分

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年之后幾年（一值到y(tǒng)1年）的速度估計。
     //sjd是y1年之后的加速度估計。瑞士星歷表jsd=31,NASA網(wǎng)站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){ //公歷轉(zhuǎn)儒略日
   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 ); }, //公歷轉(zhuǎn)儒略日
  setFromJD:function(jd){ //儒略日數(shù)轉(zhuǎn)公歷
   jd+=0.5;
   var A=int2(jd), F=jd-A, D;  //取得日數(shù)的整數(shù)部份A及小數(shù)部分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); //去除整年日數(shù)后余下日數(shù)
   this.M =int2(D/30.6001);       //月數(shù)
   this.D =D-int2(this.M*30.6001);//去除整月日數(shù)后余下日數(shù)
   this.Y-=4716; this.M--;
   if(this.M>12) this.M-=12;
   if(this.M<=2) this.Y++;
   //日的小數(shù)轉(zhuǎn)為時分秒
   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(){ //日期轉(zhuǎn)為串
   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);
  },
  //星期相關(guān)
  Weeks : new Array('日','一','二','三','四','五','六','七'),
  getWeek:function(jd){ return int2(jd+1.5)%7; }, //星期計算
  nnweek:function(y,m,n,w){ //求y年m月的第n個星期w的儒略日數(shù)
   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){ //球面坐標旋轉(zhuǎn)
  //黃道赤道坐標變換,赤到黃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){ //球面轉(zhuǎn)直角坐標
  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){ //直角坐標轉(zhuǎn)球
  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是世紀數(shù)
  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;  //黃經(jīng)章動
  this.dE=dE/100/rad;  //交角章動
 },
 nutationLon:function(t){ //只計算黃經(jīng)章動
  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是世紀數(shù)
  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年首起算的日數(shù)(UT),dt是deltatT(日),精度要求不高時dt可取值為0
  //返回格林尼治恒星時(不含赤經(jīng)章動及非多項式部分),即格林尼治子午圈的平春風(fēng)點起算的赤經(jīng)
  var t=(T+dt)/36525,t2=t*t, t3=t2*t, t4=t3*t;
  return pi2*(0.7790572732640 + 1.00273781191135448*T) //嚴格說這里的T是UT,下一行的t是力學(xué)時(世紀數(shù))
      + (0.014506 + 4612.15739966*t + 1.39667721*t2 - 0.00009344*t3 + 0.00001882*t4)/rad;
 },
 gxc_sunLon:function(t){ //太陽光行差,t是世紀數(shù)
  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; //黃經(jīng)光行差
 },
 gxc_sunLat:function(t) { return 0;       }, //黃緯光行差
 gxc_moonLon:function(t){ return -3.4E-6; }, //月球經(jīng)度光行差,誤差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){ //視差修正