我们看到一般的行星轨道有7个参数包括6个轨道参数和一个历元EPOCH.doc

我们看到一般的行星轨道有7个参数，包括6个轨道参数和一个历元epoch。根据这7个参数，在忽略近日点近动和其他轨道变动的情况下，就可以确定确定任意时刻行星的位置。这里，我们假定行星转动方向和地球公转方向一致。这6个参数是半长轴a偏心率e倾角i升交点黄经近日点幅角平近点角M其中，前5个为和时间无关的参数，可以确定轨道，平近点角是和时间有关的参数，可以确定某个时刻行星在轨道的什么位置。为什么用这6个参数就可以确定轨道，我后面在计算过程中会逐步给出。这6个参数后3个不好理解，我们根据图来理解。下载 (12.97 KB)2011-4-29 13:56先来看近日点幅角。近日点幅角是从升交点沿行星运动轨道逆时针量到近日点的角度。就是上面这个图的。我们这样定义一个坐标系，以太阳为中心，近日点方向为x轴，轨道平面为xy平面，如上图。在这个坐标系中，我们根据a和e两个参数可以确定轨道。这个轨道是（参数方程）下载 (3.24 KB)2011-4-29 13:56其中，是从椭圆参数方程里来的，和时间有关，意思是这个时刻行星位置由近日点逆时针转过的角度。这个参数，和平近点角M有关。我们来看平近点角M：天体从近点起假想地以平均角速度运动时其向径扫过的角度。百科链接：/view/1862822.htm什么意思呢？我们可以这样理解这个概念，根据M和行星的周期，我们可以算出一个时间，在我们前面提到的7个参数之一的那个历元epoch倒退这个时间后，行星就正好在近日点。我们取上面提到的倒退后的时间为时间零点，以日为时间单位，以天文单位为长度单位，由角动量守恒： 其中，是以天文单位和儒略日为单位得出的太阳常数。将带入上面的式子，可以得到关于t的微分方程。附录里有更详细的计算过程。可以解除，再代回这个r的式子求出r，从而算出坐标。下载 (28.29 KB)2011-4-29 13:56这样，我们根据a、e、M三个常数，确定了以历元epoch为时间原点，在上面那个坐标系下的坐标随时间变化的关系。但这个坐标是行星自己的轨道坐标系下的，实际应用中，我们要转为黄道坐标系下使用。这里说的黄道坐标系和通常说的不太一样，通常说的黄道坐标系是地心黄道坐标系，这里指的是日心黄道坐标系，我们这样取：取太阳为原点，地球公转平面为xy平面，春分点方向为x方向，夏至点方向为y方向，北黄极为z方向。如图：下载 (13.22 KB)2011-4-29 13:56要完成这两个坐标系间的转换，要用到一个知识：欧拉角。 百科链接：/view/90087.htm大家可以参照链接中给的图体会倾角i，升交点黄经和近日点幅角的关系，其中i就是欧拉角中的，就是欧拉角中的，就是欧拉角中的。那么，根据欧拉角，日心黄道坐标系到行星运动坐标系的变换矩阵为：如果是从行星运动坐标系转换到日心黄道坐标系，就是上面矩阵的逆矩阵，这个矩阵是正交矩阵，所以逆矩阵就是转置矩阵：我们根据a、e、M、epoch得到的行星运动坐标系下坐标(x,y,z)，那么，用这个坐标乘以上面的逆矩阵，就得到日心黄道坐标系下坐标(x,y,z)。实际上，我们观测的结果的不是这个坐标，因为这个坐标是相对于太阳的位置，我们观测的是从地球上看的位置。那么就需要我们求出地球的坐标，然后相减得到的就是地心黄道坐标。地球也可以用同样的方法求得位置，地球的a=1，e=0.0167，i=0，地球没有升交点的概念，不妨认为春分点是升交点，在这个意义下，=0，就是地球的近日点黄经，为102.3度。M跟历元有关，实际上，不用求地球的平近点角，这是因为在上面那个计算行星运动坐标随时间变化的关系中，不需要用t=0确定边界条件，应该用春分点确定边界条件。求得地心黄道坐标后，就可以根据坐标变换求出地心赤道坐标。这些都是后续问题了。附录：详细计算参数方程的解：带入得到：，其中为周期，为质量修正因子。可以积分，积分得到：其中，边界条件为(t=0)=0 。解这个超越方程就得到关于t的表达式。但是这个式子有问题，因为左边的积分是在一定积分区间（-,）内取得的，出了这个区间积分不对。我们可以对其做这样的改正：令，那么正确的积分是：，其中表示取整。实际操作时，这样做方程不好解，换一种思路。假设f()在（-,）上的反函数为g，求得反函数对应区间为(-1/2,1/2)，令，则有：。附录：完整mathematica程序(*参数输入部分*)(*填入轨道参数：a半长轴，e椭率，M平近点角，i倾角，W升交点黄经，w近日点幅角，后面加1表示角度制，填入参数历元year0、month0、day0以及你所需要的日期year1、month1、day1、hour1。如果有周期数据，请输入单位为日的周期数据T1，不知道周期请输入0。如果不知道周期数据且行星质量较大时，要输入行星质量，单位亿亿亿千克*)a:=5.20336301;e:=0.04839266;M1:=19.65053;i1:=1.30530;W1:=100.55615;w1:=274.1977;year0:=2000;month0:=1;day0:=1;T1:=4332.6;mass:=1898.6;(*填入要计算位置的时间*)year1=2008;month1=1;day1=1;hour1=0;(*一些函数*)(*日期转儒略日函数*)JDyeaR_,montH_,daY_:=(a1=Floor(14-montH)/12;y=yeaR+4800-a1;m=montH+12a1-3;daY+Floor(153m+2)/5+365y+Floory/4-32083);(*二分法解方程函数*)Dichotomyfunc_,xmin_,xmax_,m_,nmax_:=(Forcount=1;x1=xmin;x2=xmax,count0&funcx*funcxmax0,x1=x,Iffuncx1*funcx0,x2=x; ;x);(*欧拉逆矩阵函数，注意是逆矩阵，欧拉矩阵应为此矩阵的转置*)ETii_,WW_,ww_:=(E11:=CosWWCosww-CosiiSinWWSinww;E12:=-CosWWSinww-CosiiSinWWCosww;E13:=SiniiCosWW;E21:=SinWWCosww+CosiiCosWWSinww;E22:=-SinWWSinww+CosiiCosWWCosww;E23:=-SiniiCosWW;E31:=SiniiSinww;E32:=SiniiCosww;E33:=Cosii;E11,E12,E13,E21,E22,E23,E31,E32,E33);(*直角坐标转球坐标函数*)CartToSphexx_:=ArcTanxx2/xx1*57.29578+Whichxx10&xx20,0,xx10&xx2=0,ArcCoscosa,2Pi-ArcCoscosa;af*57.29578,df*57.29578);(*计算部分*)(*基本参数计算*)i:=i1*0.01745329;W:=W1*0.01745329;w:=w1*0.01745329;M:=M1*0.01745329;massS:=1988910;T:=IfT10,T1,2Pi*(massS/(massS+mass)0.5*a1.5/0.0172021;jd1:=JDyear1,month1,day1+hour1/24-1/3;jd0:=JDyear0,month0,day0;t:=jd1-jd0+M/2/Pi*T;t1:=t-T*Floort/T+0.5;ae:=1;ee:=0.01671022;Me:=-102.94719*0.01745329;ie:=0;We:=0;we:=102.94719*0.01745329;Te:=365.2564;yeare:=2010;monthe:=9;daye:=23.465;(*地球参数*)jde:=JDyeare,monthe,daye;te:=jd1-jde+Me/2/Pi*Te;t1e:=te-Te*Floorte/Te+0.5;(*计算位置角j的值,得到行星运动平面坐标系下坐标,经欧拉矩阵变换后，变换到日心黄道坐标系的坐标，并求出日心黄经黄纬*)ffai_:=ArcTanSqrt(1-e)/(1+e)*Tanfai/2/Pi-e*(1-e2)0.5*Sinfai/(1+e*Cosfai)/2/Pi-t1/T;j:=Dichotomyf,-3.14159,3.14159,10-6,100;r:=a*(1-e2)/(1+e*Cosj);Print轨道坐标Coordinate1=r*Cosj,r*Sinj,0Print日心黄道空间坐标Coordinate=ETi,W,w.Coordinate1Print日心黄道球面坐标ElipHelioCent=CartToSpheCoordinate(*地球坐标以及地心黄道坐标、地心赤道坐标J2000*)fefai_:=ArcTanSqrt(1-ee)/(1+ee)*Tanfai/2/Pi-ee*(1-ee2)0.5*Sinfai/(1+ee*Cosfai)/2/Pi-t1e/Te;je:=Dichotomyfe,-3.14159,3.14159,10-6,100;re:=ae*(1-ee2)/(1+ee*Cosje);Print地球轨道坐标x1e=re*Cosje,re*Sinje,0Print地球黄道空间坐标xe=ETie,We,we.x1ePrint地心黄道空间坐标Elip=Coordinate-xePrint地心黄道球面坐标J2000ElipGeoCent=CartToSpheElipPrint地心赤道坐标J2000Equa=ElipToEquaElipGeoCent(*木星数据a=5.20336301e=0.0483