xjw01老师您好!

  我想请教有关您写的程序 DE405/406星历表算法 其中的2个问题, 如果您有时间可否百忙指点一二, 感谢感谢!

1. 在调试您的程序的时候, 发现:

比如在计算水星轨道的时候
double testjd=2444208.55
然后
cd.getCoordOne(DE_mercury, DE_SSBARY, r2);
这时候发现

    1  1018                                                                   
  0.244420850000000000D+07  0.244424050000000000D+07 -0.399530505986556709D+08 
-0.128233975334908292D+08  0.120336360588081973D+07  0.250943620144140441D+05 
-0.531482357268504802D+04  0.249470403150736502D+03  0.560640239484045910D+01 
-0.173195792463098552D+01  0.109235362514551099D+00  0.108598532966162008D-02 
-0.762930481634479249D-03  0.588936556536435332D-04 -0.147038968266065450D-06

只用了  -0.399530505986556709D+08 这一个, 后面的-0.128233975334908292D+08  0.120336360588081973D+07  ....... 根本不需要用到.

我想请教, 究竟是我调试出错呢还是 计算水星轨道的xyz坐标根本用不到后面这些-0.128233975334908292D+08  0.120336360588081973D+07? 究竟第4个开始的数据 -0.128233975334908292D+08 是什么呢?

第二个问题就是:

在 int DE_coord::calc(double jd, int center) 这个函数中的最后有这么一段我有点不明白:

// 如果必要,以下将太阳系质心坐标改为其它中心坐标(如日心或地心)
if (center!=DE_SSBARY) {
  for (i=0; i<13; i++) {
    if="if"  (i == center) i++; // 中心天体跳过
    if="if"  (i >= 13) continue;
    for (j=0; j<6; j++)
      pv[i][j] -= pv[center][j];
  }
  for(i=0;i<6;i++) pv[center][i] = 0; // 中心天体位置及速度置0
}


请问 如果 center 不等于 DE_SSBARY , 那么这最后一段执行下去究竟得到的是什么? 是地心的坐标吗?
因为在debugPro()过程中您是这么使用的:

  cd.getCoordOne(DE_sun,  DE_earth, r1); //取太阳地心坐标
  cd.getCoordOne(DE_venus, DE_earth, r2); //取金星地心坐标

因为
const DE_earth  = 2;
const DE_SSBARY  = 11; //太阳系质心

所以这段程序的意思就是说, 如果center 不等于 DE_SSBARY, 那么就继续执行, 不管center等于2还是其他.
所以我实在有点理会不了, 是否就是指地心坐标呢?

最后有个请求,  xjw01老师可否完善您的代码, 在最后加上日心坐标呢?

致谢!

我没想到你有兴趣调试。还是给你完整的工程包吧，这样调试起来比较方便。JPL的ascii星历文件你应该有了，我就不上传了，二进制星历文件可以用我的这个工程包生成。

 

使用C++Builder6.0调试

 

http://www.fjptsz.com/xxjs/xjw/myjpl.rar

 

在JPL中本来就可以计算日心坐标。我的代码中提取坐标时getCoordOne()，中心用DE_sun就可以了

或在calc(double jd, int center)时也可以直接设置中心天体。

 

其实calc()之后才能调用getCoordOne()

如果center 不等于 DE_SSBARY，即center等于太阳或地球或其它行星，也就是在这时候要进行坐标平移变换。

 

最后一行，中心天体的位移和速度应置0。

谢谢xjw01老师的代码!

等我消化消化. 刚回来看到xjw01老师的回复太感动了!

如此说来 calc()之后

cd.getCoordOne(DE_venus, DE_sun, r2); //取金星日心坐标
cd.getCoordOne(DE_mercury, DE_sun, r2); //取水星日心坐标
......
其他类推

而

calc()之后

cd.getCoordOne(DE_venus, DE_earth, r2); //取金星地心坐标
cd.getCoordOne(DE_mercury, DE_earth, r2); //取水星地心坐标
......
其他类推

这样可以吗? 谢谢!

还有就是
1 1018
0.244420850000000000D+07 0.244424050000000000D+07 -0.399530505986556709D+08
-0.128233975334908292D+08 0.120336360588081973D+07 0.250943620144140441D+05

我用VC++ 反复调试 -0.128233975334908292D+08 确实用不到, C++Builder 6我又不怎么熟手, 不过这些不能再麻烦老师了, 等我自己再消化后再来请教.

感谢!

是的，先call()，得到太阳系质心坐标

再调用getCoordOne()得到日心坐标

-0.128233975334908292D+08 会用到的

 

你用testpro测试一下精度，误差应在10负十几次方。-0.128233975334908292D+08不用，误差就达到10的正7次方

十分感谢xjw01老师指点!

  for (i=0; i<13; i++) {

    cheby(tx, block+pD, nC, h.p3[i], pv[i]); //插值计算

}

xjw01老师请看:

int blockN = (jd-h.JD1)/h.Ta; //JD所在块号-----------------------固定的
int    pD = h.p1[i]-1 + nD*nC*h.p3[i]; // 得到段的起点位置--------------------固定的

Chebyshev多项式系数序列的变化是cheby函数的第2个参数(double* y)决定的, 在这里即 block+pD 的值, 而block+pD 的值是由2个固定的值相加, i 循环13次下来, 只能计算到每13个(水金木火...的)第一个开始的位置.

所以从程序上理解, 似乎没有计算到后面的 0.128233975334908292D+08 0.120336360588081973D+07 0.250943620144140441D+05


这是调试的结果:

jd  =  2444210.550000
nC = 3
i=0  tx="-0.487500"  (block+pD)=-39953050.598656  h.p3[i]=14  pv[i]=-34305708.291438
nC = 3
i=1  tx="-0.743750"  (block+pD)=84198597.276258  h.p3[i]=10  pv[i]=72823247.226682
nC = 3
i=2  tx="-0.743750"  (block+pD)=35178045.725858  h.p3[i]=13  pv[i]=50274731.924044
nC = 3
i=3  tx="-0.871875"  (block+pD)=-139386516.619474  h.p3[i]=11  pv[i]=-116342599.155236
nC = 3
i=4  tx="-0.871875"  (block+pD)=-696135702.255049  h.p3[i]=8  pv[i]=-688068950.966355
nC = 3
i=5  tx="-0.871875"  (block+pD)=-1388573120.335987  h.p3[i]=7  pv[i]=-1386089390.893083
nC = 3
i=6  tx="-0.871875"  (block+pD)=-1723348709.712084  h.p3[i]=6  pv[i]=-1729754435.309103
nC = 3
i=7  tx="-0.871875"  (block+pD)=-744469516.001180  h.p3[i]=6  pv[i]=-750889294.191384
nC = 3
i=8  tx="-0.871875"  (block+pD)=-4057743033.962875  h.p3[i]=6  pv[i]=-4060750304.658238
nC = 3
i=9  tx="0.025000"  (block+pD)=172217.152528  h.p3[i]=13  pv[i]=178521.636864
nC = 3
i=10  tx="-0.743750"  (block+pD)=1170575.186604  h.p3[i]=11  pv[i]=1167174.622964
nC = 2
i=11  tx="-0.487500"  (block+pD)=-0.000042  h.p3[i]=10  pv[i]=-0.000043
nC = 3
i=12  tx="-0.487500"  (block+pD)=-0.029365  h.p3[i]=10  pv[i]=-0.029673

老师有时间就看看吧, 实在打扰您, 不好意思.

cheby(tx, block+pD, nC, h.p3[i], pv[i]); //插值计算
这里block+pD传入的是一个指针，简单的说，传入的是一个数组，数组的第1个元素对应block[pD]，在cheby中用y接收数组的地址，y[0]与block[pD]读取的是同一个内存地址，y[1]与block[pD+1]读取的是同一个内存地址，这是C语言的基本语法。也许你你前使用VB或其它语言，那么要理解这种语法，可能会给你造成障碍，C++中大量使用指针