用开普勒方程求解行星运动 [中华农历论坛]

//说明：本程序列是我早期利用开普勒定律编写的，参考了国外的一篇文章，精度一般。其中地球坐标的精度约为2角秒

double poly2(double *d,double x,int n){ //多项式计算
//计算:d0 + d1*x + d2*x^2 + d3*x^3 +...+ dn*x^n
int i; double c=0,xn=1;
for(i=0;i<n;i++,xn*=x)
c+=d[i]*xn;
return c;
}
double rad2mrad(double v){ //对超过0-2PI的角度转为0-2PI
v=fmod(v,2*M_PI);
if(v<0) return v+2*M_PI;
return v;
}
void llrConv(double *JWsrc, double *JW,double E){ //球面坐标旋转
//黄道赤道坐标变换,赤到黄E取负
double HJ=rad2mrad(JWsrc[0]), HW=JWsrc[1];
double sinE =sin(E),cosE =cos(E);
double sinW=cosE*sin(HW)+sinE*cos(HW)*sin(HJ);
double J=atan2( sin(HJ)*cosE-tan(HW)*sinE, cos(HJ) );
JW[0]=rad2mrad(J);
JW[1]=asin(sinW);
}

struct Orb_elem{//行星轨道要素表
//主要参数
double N; //升交点的黄经,J2000坐标中
double i; //行星轨道与地球轨道的夹角(倾角)
double w; //轨道近日点到升交点的夹角
double M; //轨道平近点角(行星当前位置到近日点的角度),角度按时间积累
double a; //轨道半长径
double e; //轨道离心率(e=0是圆,e=0-1是椭贺,e=1的抛物线)
double pM; //平近点角速度(用于速度计算)
//其它参数(以下参数可由主参数算出)
double L; //平黄经
double P; //P=a^1.5,轨道周期,a的单位是天文单au时,P的单位是地球公转年
double q,Q;//q=a*(1-e)近日距离,Q=a*(1+e);远日距离
double v,E;//v=真近点角(0至360度),E=偏近点角(正圆上的仿射角)
double llr[3];//在黄道上的坐标
double Vr,Vh; //Vr径向速度,Vh横向速度
};

//xn=8月球,地心黄道分点坐标,a的单位千米。xn=其它：日心Date黄道分点坐标,a为AU
//xn值:0水,1金,2地,3火,4木,5土,6天,7海(日心坐标)
void getPos_xx(double jd,int xn,Orb_elem &ob,int jing=0);

//=========以下是利用轨道要素计算行星坐标的函数===============
//这些计算不含章动岁差等,粗算Date平黄道坐标下的行星坐标等参数
//============================================================
void getPosFromOrb(Orb_elem &ob){//转入轨道要素返回经纬等
//求偏近点角E,它与平近点角M的关系是:M=E-e*sin(E),式中e是离心率
//该式可由开普勒周期定律或动力学方程导出
double c, E=ob.M, M=ob.M, e=ob.e; int i;
for(i=0;i<10;i++){ //牛顿求根法计算偏近点角E
    c=(E-e*sin(E)-M)/(1-e*cos(E)); //迭代修正量
    if(fabs(c)<1e-13) break;
    E-=c;
}
ob.E=E;
double rx=ob.a*(cos(E)-e);         //分解在长轴上的分量
double ry=ob.a*sqrt(1-e*e)*sin(E); //分解短轴上的分量
ob.v=atan2(ry,rx);                 //真近点角
ob.llr[0]=ob.w+ob.v;ob.llr[1]=0;   //在行星天道中的经度及纬度
ob.llr[2]=sqrt(rx*rx+ry*ry);       //该时刻对坐标中心距离
llrConv(ob.llr,ob.llr,ob.i);       //从行星天道坐标转到Date黄道坐标
ob.llr[0]+=ob.N;                   //ob.N是升交点在黄道中的经度
//其它参数计算
double a=ob.a, pE=ob.pM*a*a/ob.llr[2]; //M'=C[1]平近点角速度,(单位:au/每日)
ob.Vh=pE*sqrt(1-e*e);//横向速度Vh=r*v'
ob.Vr=pE*(E-M);      //径向速度
ob.L =ob.w+ob.N+ob.M;//平经度
ob.P =a*sqrt(a);     //周期
ob.q =a*(1-ob.e);    //近点距
ob.Q =a*(1+ob.e);    //远点距
}
void getPos_xx(double jd,int xn,Orb_elem &ob,int jing){
//计算行星位置(Date黄道坐标),含岁差,不含章动光行差等修正
//本函数计算误差较大,适用于快速计算,xn是行星的编号
//jing是精度控制(计算地球时jing=1采用较高精度)
if(xn==9) {ob.llr[0]=ob.llr[1]=ob.llr[2]=0; return;} //太阳
static double cs[9][12]={
   //xn=8月球,地心黄道分点坐标,a的单位千米。xn=其它：日心Date黄道分点坐标,a为AU
   //xn值:0水,1金,2地,3火,4木,5土,6天,7海(日心坐标), 每行为N(0,1),i(0,1),w(0,1),M(0,1),a(0,1),e(0,1)
   { 48.3313, 3.24587E-5, 7.0047, 5.000E-8, 29.1241, 1.01444E-5, 174.7947,4.0923344368, 0.387098,0, 0.205635, 5.59E-10},//0水
   { 76.6799, 2.46590E-5, 3.3946, 2.750E-8, 54.891, 1.38374E-5, 50.4084,1.6021302244, 0.723330,0, 0.006773,-1.302E-9},//1金
   { 0,      0.00000   , 0,      0.000,    102.9405, 4.70935E-5, 357.5254,0.9856002585, 1.0,     0, 0.016709,-1.151E-9},//2地
   { 49.5574, 2.11081E-5, 1.8497,-1.780E-8, 286.5016, 2.92961E-5, 19.3881,0.5240207766, 1.523688,0, 0.093405, 2.516E-9},//3火
   {100.4542, 2.76854E-5, 1.3030,-1.557E-7, 273.8777, 1.64505E-5, 20.0196,0.0830853001, 5.20256, 0, 0.048498, 4.469E-9},//4木
   {113.6634, 2.38980E-5, 2.4886,-1.081E-7, 339.3939, 2.97661E-5, 317.0172,0.0334442282, 9.55475, 0, 0.055546,-9.499E-9},//5土
   { 74.0005, 1.39780E-5, 0.7733, 1.900E-8, 96.6612, 3.05650E-5, 142.6081,0.0117258060, 19.18171,-1.550E-8, 0.047318,7.45E-9},//6天
   {131.7806, 3.01730E-5, 1.7700,-2.550E-7, 272.8461,-6.02700E-6, 260.2561,0.0059951470, 30.05826, 3.313E-8, 0.008606,2.15E-9},//7海
   {125.0434,-0.0529538083, 5.1454,0, 318.3099,0.1643573223, 134.9629,13.0649929509, 384402, 0, 0.054900, 0       } //8月
};
static int init=0; int i,j;
if(!init){ init=1;
    for(i=0;i<9;i++) //全部转为弧度
    for(j=0;j<8;j++) cs[i][j]=D2R(cs[i][j]);
}
double *p=cs[xn];
ob.N=poly2(p, jd,2); ob.i=poly2(p+2,jd,2);
ob.w=poly2(p+4,jd,2); ob.M=poly2(p+6,jd,2);
ob.a=poly2(p+8,jd,2); ob.e=poly2(p+10,jd,2);
ob.pM=p[1];         //平近点角速度

//由已上要素计算出经纬及速度等
getPosFromOrb(ob);

//处理主要摄动问题
double d0=0,d1=0,d2=0, tt=jd/365250;
if(xn==0) d0+=(3.0+29*tt+108.4*tt*tt)/rad; //平均误差8秒
if(xn==1) d0+=(2.8+12*tt+108.2*tt*tt)/rad; //平均误差11秒
if(xn==2){ //对地球受到的摄动进行修正
    //地球的千年内精度为2角秒,最大误正负6角秒
    if(jing){
     d0+=3.417571E-5 * cos( 2.829+     3.523*tt); d0+=3.497056E-5 * cos( 2.744+ 5753.385*tt);
     d0+=3.135896E-5 * cos( 3.628+ 77713.771*tt); d0+=2.676218E-5 * cos( 4.418+ 7860.419*tt);
     d0+=2.342687E-5 * cos( 6.135+ 3930.210*tt); d0+=1.273166E-5 * cos( 2.037+   529.690*tt);
     d0+=1.324292E-5 * cos( 0.742+ 11506.770*tt); d0+=0.901855E-5 * cos( 2.045+    26.298*tt);
     d0+=1.199167E-5 * cos( 1.110+ 1577.344*tt); d0+=0.857223E-5 * cos( 3.508+   398.149*tt);
     d0+=0.779786E-5 * cos( 1.179+ 5223.694*tt); d0+=0.990250E-5 * cos( 5.233+ 5884.927*tt);
     d0+=0.753141E-5 * cos( 2.533+ 5507.553*tt); d0+=0.505264E-5 * cos( 4.583+ 18849.228*tt);
     d0+=0.492379E-5 * cos( 4.205+   775.523*tt); d0+=0.356655E-5 * cos( 2.920+     0.067*tt);
     d0+=0.284125E-5 * cos( 1.899+   796.298*tt); d0+=0.242810E-5 * cos( 0.345+ 5486.778*tt);
     d0+=0.317087E-5 * cos( 5.849+ 11790.629*tt); d0+=0.271039E-5 * cos( 0.315+ 10977.079*tt);
     d0+=0.206160E-5 * cos( 4.806+ 2544.314*tt); d0+=0.205385E-5 * cos( 1.869+ 5573.143*tt);
     d0+=0.202261E-5 * cos( 2.458+ 6069.777*tt);
     d0+=0.425264E-5 * cos( 1.590+     3.523*tt)*tt; d0+=0.108977E-5 * cos( 2.966+1577.344*tt)*tt;
     d0+=0.093478E-5 * cos( 2.592+ 18849.228*tt)*tt; d0+=0.119261E-5 * cos( 5.796+ 26.298*tt)*tt;
     d0+=8.719837E-5 * cos( 1.572+ 6283.076*tt)*tt*tt/3.3;
     //纬度摄动修正
     d1+=0.279620E-5 * cos( 3.20+ 84334.662*tt); d1+=0.101643E-5 * cos( 5.42+   5507.553*tt);
     d1+=0.080445E-5 * cos( 3.88+   5223.694*tt); d1+=0.043806E-5 * cos( 3.70+   2352.866*tt);
     d1+=0.031933E-5 * cos( 4.00+   1577.344*tt); d1+=0.022724E-5 * cos( 3.98+   1047.747*tt);
    }
    d0+=(2.264089+10.667444*tt+108.680755*tt*tt)/rad; //长期项
}
if(xn==3){ //误差20秒(J2000前后1000年)
    d0+=2.7745E-4 *cos( 5.97+     3.523*tt); d0+=1.0610E-4 *cos( 2.94+ 2281.230*tt);
    d0+=1.2316E-4 *cos( 0.85+ 2810.921*tt); d0+=0.8927E-4 *cos( 4.16+     0.017*tt);
    d0+=0.8716E-4 *cos( 6.11+ 13362.450*tt); d0+=0.6798E-4 *cos( 0.36+   398.149*tt);
    d0+=0.7775E-4 *cos( 3.34+ 5621.843*tt); d0+=(-43.1-122.6*tt+111.2*tt*tt)/rad;
}
if(xn==4){//处理木星摄动,误差1分
    d0+=5.7361E-3 *cos( 1.44+   7.114*tt); d0+=0.9718E-3 *cos( 4.14+ 632.784*tt);
    d0+=0.7290E-3 *cos( 3.64+ 522.577*tt); d0+=0.6426E-3 *cos( 3.41+ 103.093*tt);
    d0+=0.3981E-3 *cos( 2.29+ 419.485*tt); d0+=0.3886E-3 *cos( 1.27+ 316.392*tt);
    d0+=0.2796E-3 *cos( 1.78+ 536.805*tt); d0+=2.2892E-3 *cos( 6.03+   7.114*tt)*tt;
}
if(xn==5){//处理土星摄动,误差3分
    d0+=0.014142 *cos( 4.586+   7.114*tt); d0+=0.003984 *cos( 0.521+ 206.186*tt);
    d0+=0.002068 *cos( 0.247+ 103.093*tt); d0+=0.005643 *cos( 2.885+   7.114*tt)*tt;
    d0+=0.001077 *cos( 2.278+ 206.186*tt)*tt;
}
if(xn==6) d0+=(40.1-125.7*tt+2903*tt*tt)/rad; //处理天王星摄动,误差5分
if(xn==8){ //处理月球摄动
    double ws=D2R(282.9405+4.70935E-5*jd); //取太阳平近点角
    double Ms=D2R(357.5254+0.9856002585*jd);//太阳近地点到升交点夹角
    double Ls = Ms + ws;        //太阳升交点平黄经Ns=0
    double Mm = ob.M;           //取月球平近点角
    double Lm = ob.L;           //月球平黄经Mm + wm + Nm
    double D = Lm - Ls;        //日月对地角距(夹角)
    double F = Mm + ob.w;      //月球到升角点角度
    d0+=-0.02224*sin(Mm-2*D);    d0+= 0.01148*sin(2*D);      //经度摄动值
    d0+=-0.00325*sin(Ms);        d0+=-0.00103*sin(2*Mm-2*D);
    d0+=-0.00099*sin(Mm-2*D+Ms); d0+= 0.00093*sin(Mm+2*D);
    d0+= 0.00080*sin(2*D-Ms);    d0+= 0.00072*sin(Mm-Ms);
    d0+=-0.00061*sin(D);         d0+=-0.00054*sin(Mm+Ms);
    d1 =-0.00302*sin(F-2*D);     d1+=-0.00096*sin(Mm-F-2*D); //纬度摄动值
    d1+=-0.00080*sin(Mm+F-2*D); d1+= 0.00058*sin(F+2*D);
    d1+= 0.00030*sin(2*Mm+F);
    d2 =-3700*cos(Mm-2*D);       d2+=-2934*cos(2*D); //距离摄动值
}
ob.llr[0] =rad2mrad(ob.llr[0]+d0); //加上黄纬摄动量
ob.llr[1]+=d1; //加上黄纬摄动
ob.llr[2]+=d2; //加上距离摄动
}

回复标题：
上传附件：

	签名：不显示显示

主题：用开普勒方程求解行星运动

用开普勒方程求解行星运动