//=========日月升降物件=============

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*cs_rEar/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 = rad2rrad(ZB.gst(jd,this.dt) - this.L - z[0]); //得到此刻天體時角
   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;
 }
};



//========太陽月亮計算類=============

function sun_moon(){
 this.calc=function(T,L,fa,high){ //sun_moon類的成員函數。參數：T是力學時,站點經緯L,fa,海拔high(千米)
  //基本參數計算
  this.T=T, this.L=L, this.fa=fa;
  this.dt = JD.deltatT2(T); //TD-UT
  this.jd = T - this.dt;    //UT
  T/=36525; ZB.nutation(T);
  this.dL = ZB.dL;   //黃經章
  this.dE = ZB.dE;   //交角章動
  this.E  = ZB.hcjj(T) + this.dE; //真黃赤交角
  this.gst= ZB.gst(this.jd,this.dt) + this.dL*Math.cos(this.E); //真恒星時(不考慮非多項式部分)
  var z=new Array();

  //=======月亮========
  //月亮黃道坐標
  XL.M_coord(T,z,-1,-1,-1); //月球坐標
  z[0]  = rad2mrad( z[0]+ZB.gxc_moonLon(T)+this.dL );  z[1] += ZB.gxc_moonLat(T);  //補上月球光行差及章動
  this.mHJ = z[0]; this.mHW = z[1]; this.mR  = z[2]; //月球視黃經,視黃緯,地月質心距

  //月球赤道坐標
  ZB.llrConv( z, this.E ); //轉為赤道坐標
  this.mCJ = z[0]; this.mCW = z[1]; //月球視赤經,月球赤緯

  //月亮時角計算
  this.mShiJ = rad2mrad(this.gst - L - z[0]); //得到此刻天體時角
  if( this.mShiJ>Math.PI ) this.mShiJ -= pi2;

  //修正了視差的赤道坐標
  ZB.parallax(z, this.mShiJ,fa, high); //視差修正
  this.mCJ2 = z[0], this.mCW2 = z[1], this.mR2=z[2];

  //月亮時角坐標
  z[0] += Math.PI/2-this.gst+L;  //轉到相對于地平赤道分點的赤道坐標(時角坐標)

  //月亮地平坐標
  ZB.llrConv (z, Math.PI/2-fa );    //轉到地平坐標(只改經緯度)
  z[0] = rad2mrad( Math.PI/2-z[0] );
  this.mDJ = z[0]; this.mDW = z[1]; //方位角,高度角
  if(z[1]>0) z[1] += ZB.AR2(z[1]); //大氣折射修正
  this.mPJ = z[0]; this.mPW = z[1]; //方位角,高度角

  //=======太陽========
  //太陽黃道坐標
  XL.E_coord(T,z,-1,-1,-1);   //地球坐標
  z[0]  = rad2mrad(z[0]+Math.PI+ZB.gxc_sunLon(T)+this.dL);  //補上太陽光行差及章動
  z[1]  =-z[1] + ZB.gxc_sunLat(T); //z數組為太陽地心黃道視坐標
  this.sHJ = z[0]; this.sHW = z[1]; this.sR  = z[2]; //太陽視黃經,視黃緯,日地質心距

  //太陽赤道坐標
  ZB.llrConv( z, this.E ); //轉為赤道坐標
  this.sCJ = z[0]; this.sCW = z[1]; //太陽視赤經,視赤緯

  //太陽時角計算
  this.sShiJ = rad2mrad(this.gst - L - z[0]); //得到此刻天體時角
  if( this.sShiJ>Math.PI ) this.sShiJ -= pi2;

  //修正了視差的赤道坐標
  ZB.parallax(z,this.sShiJ,fa,high); //視差修正
  this.sCJ2=z[0], this.sCW2=z[1], this.sR2=z[2];

  //太陽時角坐標
  z[0] += Math.PI/2-this.gst+L;  //轉到相對于地平赤道分點的赤道坐標

  //太陽地平坐標
  ZB.llrConv( z, Math.PI/2-fa );
  z[0] = rad2mrad( Math.PI/2-z[0] );
  //z[1] -= 8.794/rad/z[2]*Math.cos(z[1]); //直接在地平坐標中視差修正(這里把地球看為球形,精度比ZB.parallax()稍差一些)
  this.sDJ = z[0]; this.sDW = z[1]; //方位角,高度角

  if(z[1]>0) z[1] += ZB.AR2(z[1]); //大氣折射修正
  this.sPJ = z[0]; this.sPW = z[1]; //方位角,高度角

  //=======其它========
  //時差計算
  t=T/10; var t2=t*t,t3=t2*t,t4=t3*t,t5=t4*t;
  var Lon = ( 1753469512 + 6283319653318*t + 529674*t2 + 432*t3 - 1124*t4 - 9*t5 + 630 * Math.cos(6+3*t) )/1000000000 + Math.PI - 20.5/rad; //修正了光行差的太陽平黃經
  Lon = rad2mrad( Lon - (this.sCJ-this.dL*Math.cos(this.E)) ); //(修正了光行差的平黃經)-(不含dL*cos(E)的視赤經)
  if(Lon>Math.PI) Lon-=pi2; //得到時差,單位是弧度
  this.sc  = Lon/pi2;   //時差(單位:日)

  //真太陽與平太陽
  this.pty = this.jd-L/pi2;  //平太陽時
  this.zty = this.jd-L/pi2+this.sc; //真太陽時

  //視半徑
  //this.mRad = XL.moonRad(this.mR,this.mDW);  //月亮視半徑(角秒)
  this.mRad = 358473400/this.mR2; //月亮視半徑(角秒)
  this.sRad = 959.63/this.sR2; //太陽視半徑(角秒)
  this.e_mRad = 358473400/this.mR; //月亮地心視半徑(角秒)
  this.eShadow = (cs_rEarA/this.mR*rad-(959.63-8.794)/this.sR )*51/50; //地本影在月球向徑處的半徑(角秒),式中51/50是大氣厚度補償
  this.eShadow2= (cs_rEarA/this.mR*rad+(959.63+8.794)/this.sR )*51/50; //地半影在月球向徑處的半徑(角秒),式中51/50是大氣厚度補償
  this.mIll = XL.moonIll(T); //月亮被照面比例

  //中心食計算
  if( Math.abs(rad2rrad(this.mCJ-this.sCJ))<50/180*Math.PI ){
   ZB.line_earth( this.mCJ,this.mCW,this.mR, this.sCJ,this.sCW,this.sR*cs_AU );
   this.zx_J = rad2rrad(this.gst-ZB.le_J);
   this.zx_W = ZB.le_W; //無解返回值是100
  } else this.zx_J=this.zx_W=100;
 };
 this.toHTML=function(fs){
  var s = '<table width="100%" cellspacing=1 cellpadding=0 bgcolor="#FFC0C0">';

  s += '<tr><td bgcolor=white align=center>';
  s += '平太陽 ' + JD.timeStr(msc.pty) + ' 真太陽 <font color=red>' + JD.timeStr(msc.zty) + '</font><br>';
  s += '時差 ' + m2fm(msc.sc*86400,2,1) + " 月亮被照亮 " + (msc.mIll*100).toFixed(2)+'% ';
  s += '</td></tr>';

  s += '<tr><td bgcolor=white><center><pre style="margin-top: 0; margin-bottom: 0"><font color=blue><b>表一       月亮            太陽</b></font>\r\n';
  s += '視黃經 ' + rad2str(msc.mHJ,0) +'  '+ rad2str(msc.sHJ,0) + '\r\n';
  s += '視黃緯 ' + rad2str(msc.mHW,0) +'  '+ rad2str(msc.sHW,0) + '\r\n';
  s += '視赤經 ' + rad2str(msc.mCJ,1) +'  '+ rad2str(msc.sCJ,1) + '\r\n';
  s += '視赤緯 ' + rad2str(msc.mCW,0) +'  '+ rad2str(msc.sCW,0) + '\r\n';
  s += '距離     ' + msc.mR.toFixed(0)  +'千米          '+ msc.sR.toFixed(6)+'AU'+'\r\n';
  s += '</pre></center></td></tr>';

  s += '<tr><td bgcolor=white><center><pre style="margin-top: 0; margin-bottom: 0"><font color=blue><b>表二       月亮            太陽</b></font>\r\n';
  s += '方位角 ' + rad2str(msc.mPJ,0) +'  '+ rad2str(msc.sPJ,0) + '\r\n';
  s += '高度角 ' + rad2str(msc.mPW,0) +'  '+ rad2str(msc.sPW,0) + '\r\n';
  s += '時角   ' + rad2str(msc.mShiJ,0)+'  '+rad2str(msc.sShiJ,0)+'\r\n';
  s += '視半徑(觀測點) ' + m2fm(msc.mRad,2,0) +'     '+ m2fm(msc.sRad,2,0)+'\r\n';
  s += '</pre></center></td></tr>';

  if(fs){
   s += '<tr><td bgcolor=white align=center>';
   s += '力學時 ' +JD.setFromJD_str(msc.T+J2000);
   s += ' ΔT=' + (msc.dt*86400).toFixed(1) +'秒<br>';
   s += '黃經章 '+(msc.dL/pi2*360*3600).toFixed(2) +'" ';
   s += '交角章 '+(msc.dE/pi2*360*3600).toFixed(2) +'" ';
   s += 'ε='+trim(rad2str(msc.E,0));
   s += '</td></tr>';
  }
  s += '</table>';
  return s;
 };
}


//====================================
var ysPL={ //月食快速計算器
 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;
 },
 lecXY:function(jd,re){//日月黃經緯差轉為日面中心直角坐標(用于月食)
  var T=jd/36525, zm=new Array(), zs=new Array();

  //=======太陽月亮黃道坐標========
  XL.E_coord(T,zs,-1,-1,-1);   //地球坐標
  zs[0]  = rad2mrad(zs[0]+Math.PI+ZB.gxc_sunLon(T));  zs[1]  =-zs[1] + ZB.gxc_sunLat(T); //補上太陽光行差
  XL.M_coord(T,zm,-1,-1,-1); //月球坐標
  zm[0]  = rad2mrad( zm[0]+ZB.gxc_moonLon(T) );  zm[1] += ZB.gxc_moonLat(T);  //補上月球光行差就可以了

  //=======視半徑=======
  re.e_mRad = 358473400/zm[2]; //月亮地心視半徑(角秒)
  re.eShadow = (cs_rEarA/zm[2]*rad-(959.63-8.794)/zs[2] )*51/50; //地本影在月球向徑處的半徑(角秒),式中51/50是大氣厚度補償
  re.eShadow2= (cs_rEarA/zm[2]*rad+(959.63+8.794)/zs[2] )*51/50; //地半影在月球向徑處的半徑(角秒),式中51/50是大氣厚度補償

  re.x = rad2rrad(zm[0]+Math.PI-zs[0]) * Math.cos((zm[1]-zs[1])/2);
  re.y = zm[1]+zs[1];
  re.mr= re.e_mRad/rad,  re.er=re.eShadow/rad, re.Er=re.eShadow2/rad;
  re.t = jd;
 },
 lecMax:function(jd){ //月食的食甚計算(jd為近朔的力學時,誤差幾天不要緊)
  this.lT=new Array();
  for(var i=0;i<7;i++) this.lT[i]=0; //分別是:食甚,初虧,復圓,半影食始,半影食終,食既,生光

  jd = XL.MS_aLon_t2( Math.floor((jd-4)/29.5306)*Math.PI*2 +Math.PI)*36525; //低精度的朔(誤差10分鐘),與食甚相差10分鐘左右

  var g=new Object(), G=new Object(), u,v;

  //求極值(平均誤差數秒)
  u = -18461 * Math.sin(0.057109+0.23089571958*jd)*0.23090/rad; //月日黃緯速度差
  v = (XL.M_v(jd/36525)-XL.E_v(jd/36525))/36525; //月日黃經速度差
  this.lecXY(jd,G);
  jd -= (G.y*u+G.x*v)/(u*u+v*v); //極值時間

  //精密求極值
  var dt=60/86400;
  this.lecXY(jd,G); this.lecXY(jd+dt,g); //精密差分得速度,再求食甚
  u = (g.y-G.y)/dt;