📄 eph.js
字号:
//=========日月升降物件=============
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;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -