有关“现代历法”——主要编者：songgz

地平线——天地的交线，它又是一个圆。圆中心位于观测点。考虑到天球很大，我们可以把圆(以及圆中心)平移到地心或日心，这样，黄道、赤道、地平线的中心就相同了。也许你会问，平移后，会不会影响日月在地平坐标中的位置？当然会有一点影响，不过我们可以通过视差修正的方法来解决问题。

p，γo'到γm的角距离就是Date黄道上的岁差ψ，γo到γ1的角距离是J2000黄道上的岁差χ，γ1到γm是Date赤道上的岁差εo，是历元2000.0的黄赤交角ε，是当日的黄赤交角ω，Date平赤道与J2000黄道的夹角π，是两个黄道之前的夹角П，是N到γo的角距离，它等于N到γo'的角距离另一组常用的岁差计算参数90°-ζ，γo到Q的角距离，是一个常用的岁差计算参数90°+z，γm到Q的角距离，是一个常用的岁差计算参数θ是两个赤道面的夹角

---------------------[问题]：程序中使用了“J2000.0起算的儒略世纪数”，你会计算它吗？[例题]：比如，2002年1月1日12:00:00TD的儒略世纪数T是多少？TD表示力学时。[解]：J2000.0对应2000年1月1日12:00:00TD由于2000年是闰年，含有366天，也就是说：2001年1月1日12:00:00加上366天后变为2001年1月1日12:00:00又因为2001年是平年，所以再加上365天就是2002年1月1日12:00:00所以，所求时刻相对于J2000.0相差366+365=731天那么，T=731/36525=0.020013689

行星光行差计算起来比较麻烦，也没有现成的简单公式可以应用，我们总是根据具体的行星问题展开计算。用下文所述的方法得到的光行差是非常准确的，最后得到的行星视位置与skymap的计算结果分毫不差。当然，日月计算也会涉及行星光行差问题，但这只以行星光行差中的最简单的一个特例。

（1）虽然星体B的运动改变了光的方向(伽俐略原理)，但对于观测者，只能接收到到来自SA路径的光线，所以观测者并没有觉察到运动方向的改变。

计算量：一次精确计算位置坐标的计算量为B，那么该算法的计算量小于1.3*B

这是一个超高速的算法，基于ELP/MPP02，最大理论误差7秒，理论平均误差1秒，程序大小8.5k。

请修改程序中的testT，以便计算其它时刻的情形。

本程序未经检验，可能有误，有空的时候我会与《中国天文年历》及“瑞士星历表”比对一下。到时再做报告。

xjw01于十中，2008年5月27日晚编写

M_L1=newArray(1.677,4.669,628.30196,-0.03,0,0,0.516,3.372,6585.7609,-2.16,-19,0.1,0.414,5.728,14914.4523,-0.64,6,0,0.371,3.969,7700.3895,1.55,25,0,0.276,0.74,8956.9934,1.5,25,0,0.246,4.23,-2.3012,1.5,25,0,0.071,0.14,7842.365,-2.2,-19,0,0.061,2.50,16171.056,-0.7,6,0,0.045,0.44,8399.679,-0.4,0,0,0.040,5.77,14286.150,-0.6,6,0,0.037,4.63,1256.604,-0.1,0,0,0.037,3.42,5957.459,-2.1,-19,0,0.036,1.80,23243.144,0.9,31,0,0.024,0,16029.081,3,50,0,0.022,1.0,-1742.931,-4,-44,0,0.019,3.1,17285.685,3,50,0,0.017,1.3,0.329,2,25,0,0.014,0.3,8326.390,3,50,0,0.013,4.0,7072.088,2,25,0,0.013,4.4,8330.993,0,0,0,0.013,0.1,8470.667,-2,-19,0,0.011,1.2,22128.515,-3,0,0,0.011,2.5,15542.754,-1,0,0,0.008,0.2,7214.06,-2,0,0,0.007,4.9,24499.75,1,0,0,0.007,5.1,13799.82,-4,0,0,0.006,0.9,-486.33,-4,0,0,0.006,0.7,9585.30,1,0,0,0.006,4.1,8328.34,2,0,0,0.006,0.6,8329.04,2,0,0,0.005,2.5,-1952.48,1,0,0,0.005,2.9,-0.71,0,0,0,0.005,3.6,30457.21,-1,0,0);

M_L2=newArray(0.0049,4.67,628.3020,0,0,0,0.0023,2.67,-2.301,1.5,25,0,0.0015,3.37,6585.761,-2.2,-19,0,0.0012,5.73,14914.452,-0.6,6,0,0.0011,3.97,7700.389,2,25,0,0.0008,0.7,8956.993,1,25,0);

functionMnn(F,t,t2,t3,t4,n){//计算M_L0或M_L1或M_L2n=Math.floor(n*F.length/M_L0.length+0.5)*6;//项数是以周期项为基准,n是周期项的计算项数if(n>F.length)n=F.length;vari,v=0;t2/=1e4,t3/=1e8,t4/=1e8;for(i=0;i

functionmoonV(t,jing){//月球速度计算,传入世经数varv=8399.70911033384;if(jing==0)returnv;v-=914*Math.sin(0.7848+8328.691425*t+1.523*t*t/10000);if(jing==1)returnv;//误差小于5%v-=179*Math.sin(2.543+15542.7543*t)//误差小于0.3%+160*Math.sin(0.1874+7214.0629*t)+62*Math.sin(3.14+16657.3828*t)+34*Math.sin(4.827+16866.9323*t)+22*Math.sin(4.9+23871.4457*t)+12*Math.sin(2.59+14914.4523*t)+7*Math.sin(0.23+6585.7609*t)+5*Math.sin(0.9+25195.624*t)+5*Math.sin(2.32-7700.3895*t)+5*Math.sin(3.88+8956.9934*t)+5*Math.sin(0.49+7771.3771*t)+4*Math.sin(5.5+24986.074*t)+4*Math.sin(4.3+22756.817*t)+2*Math.sin(1+32200.137*t)+2*Math.sin(1.9+14428.126*t)+2*Math.sin(0.4+31085.509*t)+2*Math.sin(1.528+628.302*t);returnv;}

testT=30;//世纪数L=moonLon(testT,-1);//正算t=moon_t(L);//反算dt=(testT-t)*35625*86400;

为了能够与《中国天文年历》比对，我们须把以上程序中的月球坐标转换为视坐标

这时我们应处理光行差及章动

月球光行差问题我们已讨论过精密的光行差计算，但是，在农历计算中根本无需精密计算，以下考虑简化计算。从ELP的中序列截取5%精度的距离及速速表达式：日月距离：r=385001+20905*sin(0.7848+8328.691425*t+1.523*t*t/10000);//最大误差2%黄经速度：v=8400-914*sin(0.7848+8328.691425*t+1.523*t*t/10000);//最大误差5%为了求证方便，令：r0=385001，v0=8400φ=0.7848+8328.691425*t+1.523*t*t/10000dr=20905*sin(φ)dv=914*sin(φ)光速：c=300000千米/秒=9.467E14千米/世纪那么有：r=r0+dr；v=v0+dv；τ=r/c；(式中τ光行时)光行期间黄经的变化量：dL=v*τ=v*r/c，代入得：dL=(v0+dv)*(r0+dr)/c=(v0*r0/c)*(1+dv/v+dr/r)(略去二阶小量)dL=(3.416E-6弧度)*(1+(20905/385001-914/8400)*sin(φ))=(0.7角秒)*(1-0.055*sin(φ))结论：月球的黄经光行差约为0.7角秒(最大误差16%)，如果考虑上式中的周期，最大误差5%

以下关于章动，取自IAU1980，最大误差不超过0.1角秒functionnutation(t){//章动计算,t为世纪数,返回值为角秒单位varc,dL=0,dE=0,t2=t*t,y0=-17.20-0.01742*t;c=2.1824-33.75705*t+36e-6*t2;dL=y0*Math.sin(c);dE=9.20*Math.cos(c);c=3.5069+1256.66393*t+11e-6*t2;dL+=-1.32*Math.sin(c);dE+=0.57*Math.cos(c);c=1.3375+16799.4182*t-51e-6*t2;dL+=-0.23*Math.sin(c);dE+=0.10*Math.cos(c);c=4.3649-67.5141*t+72e-6*t2;dL+=0.21*Math.sin(c);dE+=-0.09*Math.cos(c);c=0.04-628.302*t;dL+=-0.14*Math.sin(c);c=2.36+8328.691*t;dL+=0.07*Math.sin(c);c=3.46+1884.966*t;dL+=-0.05*Math.sin(c);dE+=0.02*Math.cos(c);c=5.44+16833.175*t;dL+=-0.04*Math.sin(c);dE+=0.02*Math.cos(c);c=3.69+25128.110*t;dL+=-0.03*Math.sin(c);c=3.55+628.362*t;dL+=0.02*Math.sin(c);this.dL=dL;//黄经章动this.dE=dE;//交角章动}

调用方法：nu=newnutation(testT);

那么：nu.dL为黄经章动，nu.dE为交角章动

/*

我们修改第八贴的测试程序的末尾部分，以便能够与《天文年历》比对

*/

rad=180*3600/Math.PI;functionrad2mrad(v){//对超过0-2PI的角度转为0-2PIv=v%(2*Math.PI);if(v<0)returnv+2*Math.PI;returnv;}functionrad2str(d,tim){//将弧度转为字串//tim=0输出格式示例:-23°59'48.23"//tim=1输出格式示例:18h29m44.52svars="+";varw1="°",w2="’",w3="”";if(d<0)d=-d,s='-';if(tim){d*=12/Math.PI;w1="h",w2="m",w3="s";}elsed*=180/Math.PI;vara=Math.floor(d);d=(d-a)*60;varb=Math.floor(d);d=(d-b)*60;varc=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;returns;}

testT=0.08+(1-1.5)/36525;//世纪数L=moonLon(testT,-1);//正算t=moon_t(L);//反算dt=(testT-t)*35625*86400;

//以下是视位置计算functionnutation(t){//章动计算,t为世纪数,返回值为角秒单位varc,dL=0,dE=0,t2=t*t,y0=-17.20-0.01742*t;c=2.1824-33.75705*t+36e-6*t2;dL=y0*Math.sin(c);dE=9.20*Math.cos(c);c=3.5069+1256.66393*t+11e-6*t2;dL+=-1.32*Math.sin(c);dE+=0.57*Math.cos(c);c=1.3375+16799.4182*t-51e-6*t2;dL+=-0.23*Math.sin(c);dE+=0.10*Math.cos(c);c=4.3649-67.5141*t+72e-6*t2;dL+=0.21*Math.sin(c);dE+=-0.09*Math.cos(c);c=0.04-628.302*t;dL+=-0.14*Math.sin(c);c=2.36+8328.691*t;dL+=0.07*Math.sin(c);c=3.46+1884.966*t;dL+=-0.05*Math.sin(c);dE+=0.02*Math.cos(c);c=5.44+16833.175*t;dL+=-0.04*Math.sin(c);dE+=0.02*Math.cos(c);c=3.69+25128.110*t;dL+=-0.03*Math.sin(c);c=3.55+628.362*t;dL+=0.02*Math.sin(c);this.dL=dL;//黄经章动this.dE=dE;//交角章动}

nu=newnutation(testT);L2=rad2mrad(L+nu.dL/rad-0.7/rad);//补章动与光行差Ls=rad2str(L2,0);document.write("视黄经："+Ls+"
");

测试结果如下：

随机插出几个月球视黄经数据进行比较

2008年01月01日0hTD+197°19'24.43"中国天文年历+197°19’24.91"本程序

2008年01月06日0hTD+256°54'36.32"中国天文年历+256°54'36.11"本程序

2008年01月18日0hTD+56°04'29.83"中国天文年历+56°04'29.68"本程序

2100年01月01日0hTD+157°24'01.183"瑞士星历表+157°24'01.96"本程序

2100年01月18日0hTD+22°14'39.400"瑞士星历表+22°14'40.47"本程序

2200年01月02日0hTD+108°26'45.916"瑞士星历表+108°26’46.12"本程序我们不难发现，误差与精密星历相差无几，黄经误差在0.5"左右，对应的时应误差为1秒左右。

如果用标准差表达误差，误差大约为3秒左右。

这已足已证明它是个精确的月历程序，程序不到10k

这是精简的超高速的算法

rad=180*3600/Math.PI;functionrad2mrad(v){//对超过0-2PI的角度转为0-2PIv=v%(2*Math.PI);if(v<0)returnv+2*Math.PI;returnv;}functionrad2str(d,tim){//将弧度转为字串//tim=0输出格式示例:-23°59'48.23"//tim=1输出格式示例:18h29m44.52svars="+";varw1="°",w2="’",w3="”";if(d<0)d=-d,s='-';if(tim){d*=12/Math.PI;w1="h",w2="m",w3="s";}elsed*=180/Math.PI;vara=Math.floor(d);d=(d-a)*60;varb=Math.floor(d);d=(d-b)*60;varc=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;returns;}

varE_L2=newArray(87198,1.0721,6283.07585,3091,0.8673,12566.1517,273,0.053,3.523,163,5.188,26.298,158,3.685,155.42,95,0.76,18849.23,89,2.06,77713.77,70,0.83,775.52,51,4.66,1577.34,41,1.03,7.11,38,3.44,5573.14,35,5.14,796.3,32,6.05,5507.55,30,1.19,242.73,29,6.12,529.69,27,0.31,398.15,25,2.28,553.57,24,4.38,5223.69,21,3.75,0.98,17,0.9,951.72,15,5.76,1349.87,14,4.36,1748.02,13,3.72,1194.45,13,2.95,6438.5,12,2.97,2146.17,11,1.27,161000.69,10,0.6,3154.69,10,5.99,6286.6,9,4.8,5088.63,9,5.23,7084.9,8,3.31,213.3,8,3.42,5486.78,7,6.19,4690.48,7,3.43,4694,6,1.6,2544.31,6,1.98,801.82,6,2.48,10977.08,5,1.44,6836.65,5,2.34,1592.6,5,1.31,4292.33,5,3.81,149854.4,4,0.04,7234.79,4,4.94,7632.94,4,1.57,71430.7,4,3.17,6309.37,3,0.99,6040.35,3,0.67,1059.38,3,3.18,2352.87,3,3.55,8031.09,3,1.92,10447.39,3,2.52,6127.66,3,4.42,9437.76,3,2.71,3894.18,3,0.67,25132.3,3,5.27,6812.77,3,0.55,6279.55);

varE_L3=newArray(2892,5.8438,6283.0758,168,5.488,12566.152,30,5.2,155.42,13,4.72,3.52,7,5.3,18849.23,6,5.97,242.73,4,3.79,553.57,1,4.3,6286.6,1,0.91,6127.66);

varE_L4=newArray(77,4.13,6283.08,8,3.84,12566.15,4,0.42,155.42);varE_L5=newArray(2,2.77,6283.08,1,2.01,155.42);

//地球黄纬数据,误差0.2"varE_B0=newArray(2796,3.1987,84334.6616,1016,5.4225,5507.5532,804,3.88,5223.694,438,3.704,2352.866,319,4,1577.344,227,3.985,1047.747);

varE_B1=newArray(90,3.9,5507.55,62,1.73,5223.69);varE_B2=newArray(17,1.63,84334.66);

//地球向径数据,误差0.00001AUvarE_R0=newArray(1000139888,0,0,16706996,3.09846351,6283.07584999,139560,3.055246,12566.1517,30837,5.19847,77713.77147,16285,1.17388,5753.38488,15756,2.84685,7860.41939,9248,5.4529,11506.7698,5424,4.5641,3930.2097,4721,3.661,5884.9268,3460,0.9637,5507.5532,3288,5.8998,5223.6939,3068,0.2987,5573.1428,2432,4.2735,11790.6291,2118,5.8471,1577.3435,1858,5.0219,10977.0788,1748,3.0119,18849.2275);

varE_R1=newArray(1030186,1.1074897,6283.07585,17212,1.06442,12566.1517,7022,3.1416,0);varE_R2=newArray(43594,5.78455,6283.07585,1236,5.5793,12566.1517,123,3.142,0,88,3.63,77713.77);varE_R3=newArray(1446,4.2732,6283.0758,67,3.92,12566.15);varE_R4=newArray(39,2.56,6283.08,3,2.27,12566.15);varE_R5=newArray(1,1.22,6283.08);

functionhcjj(t){//返回黄赤交角,t是千年数vart2=t*t,t3=t2*t,t4=t3*t;return(84381.4088-468.36051*t-0.01667*t2-1.99911*t3-0.00523*t4)/rad;}

z=newArray();da=1;testT=0.008+(da-1.5)/365250;//世纪数earthCoor(testT,z,-1);//地球黄经z[0]+=Math.PI,z[1]=-z[1];//太阳黄经黄纬

//L=z[0]-50287.9*(testT-0.008)/rad;//L=rad2mrad(L);//document.write(rad2str(L,0)+"
");

gx=gxc_sun(testT,z[0]);//太阳光行差nu=newnutation(testT*10);//章动计算z[0]+=(nu.dL+gx)/rad;//补章动与光行差E=hcjj(testT)+nu.dE/rad;//得到真黄赤交角HCconv(z,E);//转到赤道坐标

document.write("坐标精度测试
");document.write("testT："+testT+"
");document.write("视赤经："+rad2str(z[0],1)+"
");document.write("视赤纬："+rad2str(z[1],0)+"
");document.write("距离："+z[2]+"
");document.write("请与《2008中国天文年历》太阳视赤经赤纬比对。该年的第"+da+"天0hTD。注意：几乎分秒不差！
");document.write("本程序中VSOP87地球数据序列以做了长期项拟合到DE405/406
");