| 以文本方式查看主题 - 中华农历论坛 (http://bbs.nongli.net/index.asp) -- 历法知识 (http://bbs.nongli.net/list.asp?boardid=2) ---- 就JPL星历表文件结构再度请教xjw01老师 (http://bbs.nongli.net/dispbbs.asp?boardid=2&id=18951) |
| -- 作者：蔡越 -- 发布时间：2009/5/10 18:41:00 -- 就JPL星历表文件结构再度请教xjw01老师 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老师可否完善您的代码, 在最后加上日心坐标呢? 致谢! |
| -- 作者：xjw01 -- 发布时间：2009/5/10 20:18:00 -- 我没想到你有兴趣调试。还是给你完整的工程包吧，这样调试起来比较方便。JPL的ascii星历文件你应该有了，我就不上传了，二进制星历文件可以用我的这个工程包生成。
使用C++Builder6.0调试
http://www.fjptsz.com/xxjs/xjw/myjpl.rar
[此贴子已经被作者于2009-5-10 21:45:44编辑过]
|
| -- 作者：xjw01 -- 发布时间：2009/5/10 20:21:00 -- 在JPL中本来就可以计算日心坐标。我的代码中提取坐标时getCoordOne()，中心用DE_sun就可以了 或在calc(double jd, int center)时也可以直接设置中心天体。
其实calc()之后才能调用getCoordOne() |
| -- 作者：xjw01 -- 发布时间：2009/5/10 21:50:00 -- 以下是引用蔡越在2009-5-10 18:41:00的发言：
// 如果必要,以下将太阳系质心坐标改为其它中心坐标(如日心或地心) 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老师可否完善您的代码, 在最后加上日心坐标呢? 致谢! 如果center 不等于 DE_SSBARY，即center等于太阳或地球或其它行星，也就是在这时候要进行坐标平移变换。
最后一行，中心天体的位移和速度应置0。 |
| -- 作者：蔡越 -- 发布时间：2009/5/10 23:40:00 -- 以下是引用xjw01在2009-5-10 20:21:00的发言： 在JPL中本来就可以计算日心坐标。我的代码中提取坐标时getCoordOne()，中心用DE_sun就可以了 或在calc(double jd, int center)时也可以直接设置中心天体。
其实calc()之后才能调用getCoordOne() 谢谢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我又不怎么熟手, 不过这些不能再麻烦老师了, 等我自己再消化后再来请教. 感谢! |
| -- 作者：xjw01 -- 发布时间：2009/5/11 6:35:00 -- 是的，先call()，得到太阳系质心坐标 再调用getCoordOne()得到日心坐标 |
| -- 作者：xjw01 -- 发布时间：2009/5/11 6:37:00 -- -0.128233975334908292D+08 会用到的
你用testpro测试一下精度，误差应在10负十几次方。-0.128233975334908292D+08不用，误差就达到10的正7次方 |
| -- 作者：蔡越 -- 发布时间：2009/5/11 12:41:00 -- 以下是引用xjw01在2009-5-11 6:37:00的发言： -0.128233975334908292D+08 会用到的
你用testpro测试一下精度，误差应在10负十几次方。-0.128233975334908292D+08不用，误差就达到10的正7次方 十分感谢xjw01老师指点! |
| -- 作者：蔡越 -- 发布时间：2009/5/11 13:11:00 -- 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 老师有时间就看看吧, 实在打扰您, 不好意思. |
| -- 作者：xjw01 -- 发布时间：2009/5/11 23:07:00 -- cheby(tx, block+pD, nC, h.p3[i], pv[i]); //插值计算 |