天文算法——第1集
xjw01 2009.1.30
//==================基本常数================
J2000 = 2451545;
rad = 180/Math.PI*3600; //每弧度的角秒数
//==================数学部分================
Number.prototype.toFixed=function(m){ //IE6.0的toFixed()有Bug,所以重写
var n=this, f='', p=Math.pow(10,m); //p为10进制移位量;
if(n<0) n = -n, f = '-'; //把负数转为正数
var a=Math.floor(n),b=n-a; //分离整数与小数
b = Math.round(b*p ); //移位并四舍五入
if(b>=p) a++,b-=p; //进位
if(m) b = '.'+(p+b+'').substr(1,m); //小数部分左边补0
else b = '';
return f+a+b;
};
function fix(f,v){ //格式化输出数字,如fix("6.1",5.36)
if(!f.length) return v+'';
var i,p=f.indexOf('.'), a,b,c=' ';
if(p<0) a=f-0, b=7;
else a=f.substr(0,p)-0, b=f.substr(p+1,f.length-p-1);
c+=(v-0).toFixed(b);
return c.substr(c.length-a,a);
}
function printf(f){ //仿c语言printf(format,n1,n2……)格式化输出,格式化方式：串对应%s,数字对应%d或%8.4d或%4e
var i,j,k=1,c,c2,c3,s='';
for(i=0;i<f.length;i++){
c=f.substr(i,1);
if(c!='%') { s+=c; continue; }
if(f.substr(i+1,1)=='%') { s+=c; i++; continue; }
c2=c3='';
if(k<arguments.length) c2=arguments[k++];
for(i++,j=i;i<f.length;i++){
c=f.substr(i,1);
if(c=='s'||c=='d'||c=='e') { c3=f.substr(j,i-j); break; }
}
if(c=='s') s+=c2;
if(c=='d') s+=fix(c3,c2);
if(c=='e') s+=(c2-0).toExponential(c3);
}
return s;
}
function sqrt(x) { return Math.sqrt(x); }
function floor(x){ return Math.floor(x); }
function abs(x) { return Math.abs(x); }
function sin(x) { return Math.sin(x); }
function cos(x) { return Math.cos(x); }
function tan(x) { return Math.tan(x); }
function asin(x) { return Math.asin(x); }
function atan(x) { return Math.atan(x); }
function atan2(y,x){ return Math.atan2(y,x); }
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 llrConv(JW,E){ //球面坐标旋转
//黄道赤道坐标变换,赤到黄E取负
var r=new Array(),J=JW[0],W=JW[1];
r[0]=atan2( sin(J)*cos(E) - tan(W)*sin(E), cos(J) );
r[1]=asin ( cos(E)*sin(W) + sin(E)*cos(W)*sin(J) );
r[2]=JW[2];
r[0]=rad2mrad(r[0]);
return r;
}
function xyzConv(z,A){ //直角坐标旋转
//A[9]表示新坐标三轴的单位矢
var r = new Array();
r[0] = A[0]*z[0] + A[1]*z[1] + A[2]*z[2];
r[1] = A[3]*z[0] + A[4]*z[1] + A[5]*z[2];
r[2] = A[6]*z[0] + A[7]*z[1] + A[8]*z[2];
return r;
}
function xyzConvPQ(z,P,Q){ //小角量坐标旋转(直角坐标)
//如果z为天体的J2000黄道从标,P和Q为此刻Date黄极的J2000坐标(即岁差表中的P和Q),那么本函数把z旋转到Date黄道(注意,不函黄纬总岁差)
//如果P和Q为岁差表中的-P和-Q，本函数则把z转到J2000坐标
P/=2,Q/=2;
var P2=P*P, Q2=Q*Q, R=sqrt(1-P2-Q2);
var A=new Array(
1-2*P2, 2*P*Q, -2*P*R,
2*P*Q, 1-2*Q2, 2*Q*R,
2*P*R, -2*Q*R, 1-2*(P2+Q2) );
return xyzConv(z,A);
}
function llrConvPQ(z,P,Q){ //小角量坐标旋转(球面坐标)
var r=llr2xyz(z);
r=xyzConvPQ(r,P,Q);
return xyz2llr(r);
}
function llr2xyz(JW){ //球面转直角坐标
var r=new Array(),J=JW[0],W=JW[1],R=JW[2];
r[0]=R*cos(W)*cos(J);
r[1]=R*cos(W)*sin(J);
r[2]=R*sin(W);
return r;
}
function xyz2llr(xyz){ //直角坐标转球
var r=new Array(), x=xyz[0],y=xyz[1],z=xyz[2];
r[2] = sqrt(x*x+y*y+z*z);
r[1] = asin(z/r[2]);
r[0] = rad2mrad( atan2(y,x) );
return r;
}
function poly2(d,x,n){ //多项式计算
//计算:d0 + d1*x + d2*x^2 + d3*x^3 +...+ dn*x^n
var i, c=0,xn=1;
for(i=0;i<n;i++,xn*=x)
c+=d[i]*xn;
return c;
}
function dxs_c2j(a){ //(t+a1)(t+a2)(t+a3)……多项式_乘式转加式(即展开多项式),a是数组(a1,a2,a3,……)
var i,j;
var r=new Array(a[0],1);
for(i=1;i<a.length;i++){//从第二项开始作乘法,共length-1次
for(j=i+1,r[j]=0;j>0;j--)
r[j]=r[j-1]+r[j]*a[i];
r[0]*=a[i];
}
return r;
}
function jfcz(a){ //最大主元高斯法解方程组,a数组表示方程组(a11,a12,a13,……,b1, a21,a22,a23,……,b2, a31……)
var hn=floor((-1+sqrt(4*a.length+1)+1e-5)/2); //hn是方程行数
var i,j,k,hw=hn+1,len=hn*hw; //hw为行宽,len数组元素总个数
var p,pm,c,max; //p,pm元素定位变量等临时变量
var r=new Array();
for(i=0;i<len;i++) r[i]=a[i];
for(i=0;i<hn;i++){
p=i*hw+i; //p待定主元位置
for(j=p,pm=p,max=0;j<len;j+=hw){ //选最大主元,位置在pm
c=abs(r[j]);
if(c>max) max=c,pm=j;
}
for(j=0;j<hw-i;j++) c=r[p+j],r[p+j]=r[pm+j],r[pm+j]=c; //行交换
for(j=i;j<len;j+=hw){ //高斯消去法
if(j==p) continue; //当前行不能消元
c=r[j]/r[p]; //行相减系数
for(k=0;k<hw-i;k++) r[j+k]-=r[p+k]*c; //行相减消元
}
}
for(i=0;i<hn;i++) r[i]=r[i*hw+hn]/r[i*hw+i]; //结果
r.length=hn; return r;
}
function dxs_zxec(f,x,n){//多项式拟合(使用最小二乘法),序列的点数同x数组确定,f数组是x的函数,n是多项式的项数
var fc=new Array(), Xk=new Array();
var i,j,k,p,len=n*(n+1),xn=x.length;
for(i=0;i<len;i++) fc[i]=0;
for(k=0;k<xn;k++){
for(i=1,Xk[0]=1;i<n;i++) Xk[i]=Xk[i-1]*x[k]; //算出x[k]的各次方
for(i=0;i<n;i++){
p=i*(n+1);
for(j=0;j<n;j++) fc[p+j]+=Xk[j]*Xk[i];
fc[p+n]+=f[k]*Xk[i];
}
}
return jfcz(fc);
}
function dxsXS_lglr(f,x){//拉格朗日插值求多项式系数
var i,j,k,c,n=x.length;
var r=new Array(),b=new Array(),g;
for(i=0;i<n;i++) r[i]=0;
for(i=0;i<n;i++){
for(j=0,k=0,c=f[i];j<n;j++)
if(i!=j) c/=x[i]-x[j],b[k++]=-x[j];
g=dxs_c2j(b);
for(j=0;j<n;j++) r[j]+=g[j]*c;
}
return r;
}
function dxsXS_jfcz(f,x){ //解方程组法求多项式系数
var i,j,p,c,n=x.length;
var a=new Array();
for(i=0;i<n;i++){
p=i*(n+1),c=1;
for(j=0;j<n;j++) a[p+j]=c,c*=x[i];
a[p+n]=f[i];
}
return jfcz(a);
}
//==============岁差部分===============
var preceTab_IAU1976=new Array( //IAU1976岁差表
0, 5038.7784, -1.07259, -0.001147, //fi
84381.448, 0, +0.05127, -0.007726, //w
0, +4.1976, +0.19447, -0.000179, //P
0, -46.8150, +0.05059, +0.000344, //Q
84381.448,-46.8150, -0.00059, +0.001813, //E
0, +10.5526, -2.38064, -0.001125, //x
0, 47.0028, -0.03301, +0.000057, //pi(自导出),根据P=sin(pi)*sin(II)
629554.886,-869.8192, +0.03666, -0.001504, //II(自导出),根据Q=sin(pi)*cos(II)
0, 5029.0966, +1.11113, +0.000006, //p
0, 2004.3109, -0.42665, -0.041833, //th
0, 2306.2181, +0.30188, +0.017998, //Z
0, 2306.2181, +1.09468, +0.018203 //z
);
var preceTab_IAU2000=new Array( //IAU2000岁差表
0, 5038.478750, -1.07259, -0.001147,0,0, //fi
84381.448, -0.025240, +0.05127, -0.007726,0,0, //w,
0, +4.1976, +0.19447, -0.000179,0,0, //P,
0, -46.8150, +0.05059, +0.000344,0,0, //Q,
84381.448, -46.84024, -0.00059, +0.001813,0,0, //E,
0, +10.5526, -2.38064, -0.001125,0,0, //x
0, 47.0028, -0.03301, +0.000057,0,0, //pi(自导出)
629554.886,-869.8192, +0.03666, -0.001504,0,0, //II(自导出)
0, 5028.79695, +1.11113, +0.000006,0,0, //p
0, 2004.1917476, -0.4269353, -0.0418251, -0.0000601, -0.0000001, //th
+2.5976176,2306.0809506, +0.3019015, +0.0179663, -0.0000327, -0.0000002, //Z
-2.5976176,2306.0803226, +1.0947790, +0.0182273, +0.0000470, -0.0000003 //z
);
var preceTab_P03=new Array(
0, 5038.481507, -1.0790069, -0.00114045, +0.000132851, -9.51e-8, //fi
84381.406000, -0.025754, +0.0512623, -0.00772503, -4.67e-7, +3.337e-7, //w
0, 4.199094, +0.1939873, -0.00022466, -9.12e-7, +1.20e-8, //P
0, -46.811015, +0.0510283, +0.00052413, -6.46e-7, -1.72e-8, //Q
84381.406000, -46.836769, -0.0001831, +0.00200340, -5.76e-7, -4.34e-8, //E
0, 10.556403, -2.3814292, -0.00121197, +0.000170663, -5.60e-8, //x
0, 46.998973, -0.0334926, -0.00012559, +1.13e-7, -2.2e-9, //pi
629546.7936, -867.95758, +0.157992, -0.0005371, -0.00004797, +7.2e-8, //II
0, 5028.796195, +1.1054348, +0.00007964, -0.000023857, +3.83e-8, //p
0, 2004.191903, -0.4294934, -0.04182264, -7.089e-6, -1.274e-7, //th
2.650545, 2306.083227, +0.2988499, +0.01801828, -5.971e-6, -3.173e-7, //Z
-2.650545, 2306.077181, +1.0927348, +0.01826837, -0.000028596, -2.904e-7 //z
);
var preceTab_B03=new Array( //B03岁差表
0.0, 5038.478750, -1.0719530, -0.00114366, 0.000132832, -9.40e-8, -3.50e-9, 1.7e-10, //fi
84381.40880, -0.026501, 0.0512769, -0.00772723, -0.000000492, 3.329e-7, -3.1e-10,-6.0e-11, //w
0.0, 4.199604, 0.1939715, -0.00022350, -0.000001035, 1.9e-9, 0.0, 0.0, //P
0.0, -46.809550, 0.0510421, 0.00052228, -0.000000569, -1.4e-9, 1.0e-11, 0.0, //Q
84381.40880, -46.836051, -0.0001667, 0.00199911, -0.000000523, -2.48e-8, -3.0e-11, 0.0, //E
0.0, 10.557686, -2.3813769, -0.00121258, 0.000170238, -7.70e-8, -3.99e-9, 1.6e-10, //x
0.0, 46.997560, -0.0335050, -0.00012370, 0.000000030, 0.0, 0.0, 0.0, //pi
629543.988, -867.9218, 0.15342, 0.000005, -0.0000037, -1.0e-8, 0.0, 0.0, //II
0.0, 5028.792262, 1.1124406, 0.00007699, -0.000023479, -1.78e-8, 1.8e-10, 1.0e-11, //p
0.0, 2004.190936, -0.4266980, -0.04182364, -0.000007291, -1.127e-7, 3.6e-10, 9.0e-11, //th
2.72767, 2306.080472, 0.3023262, 0.01801752, -0.000005708, -3.040e-7, -1.3e-10, 0.0, //Z
-2.72767, 2306.076070, 1.0956768, 0.01826676, -0.000028276, -2.486e-7, -5.0e-11, 0.0 //z
);
function prece(t,sc,mx){ //t是儒略世纪数(T力学时),sc是岁差量名称,mx是模型
var i,tn=1,c=0, p,n;
if(mx=='IAU1976') n=4, p=preceTab_IAU1976;
if(mx=='IAU2000') n=6, p=preceTab_IAU2000;
if(mx=='P03') n=6, p=preceTab_P03;
if(mx=='B03') n=8, p=preceTab_B03;
sc=String("fi w P Q E x pi II p th Z z ").indexOf(sc+' ')/3;
for(i=0;i<n;i++,tn*=t) c+=p[sc*n+i]*tn;
return c/rad;
}
function HDllr_J2D(t,llr,mx){//黄道球面坐标_J2000转Date分点,t为儒略世纪数
//J2000黄道旋转到Date黄道(球面对球面),也可直接由利用球面旋转函数计算,但交角接近为0时精度很低
var r=new Array( llr[0],llr[1],llr[2] );
r[0]+=prece(t,'fi',mx); r=llrConv(r, prece(t,'w',mx));
r[0]-=prece(t,'x', mx); r=llrConv(r,-prece(t,'E',mx));
return r;
}
function HDllr_D2J(t,llr,mx){//黄道球面坐标_Date分点转J2000,t为儒略世纪数
var r=new Array( llr[0],llr[1],llr[2] );
r=llrConv(r, prece(t,'E',mx)); r[0]+=prece(t,'x', mx);
r=llrConv(r,-prece(t,'w',mx)); r[0]-=prece(t,'fi',mx);
r[0]=rad2mrad(r[0]);
return r;
}
加好友 
发短信
Post By：2008/10/22 18:30:00
