历法知识-回复帖子 [中华农历论坛]

主题最新回顾(发布时间：2009/5/12 1:09:00)
-- 作者：蔡越 -- *以下是引用xjw01在2009-5-11 23:07:00的发言：* 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++中大量使用指针谢谢xjw01老师! 我的编程完全是自学的, 老底让老师看出来了. 感谢! 既学天文知识, 又学了编程知识!

主题最新回顾(发布时间：2009/5/12 1:09:00)

-- 作者：蔡越
--

以下是引用xjw01在2009-5-11 23:07:00的发言：

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++中大量使用指针

谢谢xjw01老师! 我的编程完全是自学的, 老底让老师看出来了.
感谢! 既学天文知识, 又学了编程知识!

主题最新回顾(发布时间：2009/5/11 23:07:00)
-- 作者：xjw01 -- 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++中大量使用指针

主题最新回顾(发布时间：2009/5/11 23:07:00)

-- 作者：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 + nDnCh.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 老师有时间就看看吧, 实在打扰您, 不好意思.

主题最新回顾(发布时间：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

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

主题最新回顾(发布时间：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 12:41:00)

-- 作者：蔡越
--

以下是引用xjw01在2009-5-11 6:37:00的发言：

-0.128233975334908292D+08 会用到的

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

十分感谢xjw01老师指点!

主题最新回顾(发布时间：2009/5/11 6:37:00)
-- 作者：xjw01 -- -0.128233975334908292D+08 会用到的你用testpro测试一下精度，误差应在10负十几次方。-0.128233975334908292D+08不用，误差就达到10的正7次方

主题最新回顾(发布时间：2009/5/11 6:37:00)

-- 作者：xjw01
--

-0.128233975334908292D+08 会用到的

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

主题最新回顾(发布时间：2009/5/11 6:35:00)
-- 作者：xjw01 -- 是的，先call()，得到太阳系质心坐标再调用getCoordOne()得到日心坐标

主题最新回顾(发布时间：2009/5/11 6:35:00)

-- 作者：xjw01
--

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

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

主题最新回顾(发布时间：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我又不怎么熟手, 不过这些不能再麻烦老师了, 等我自己再消化后再来请教. 感谢!

主题最新回顾(发布时间：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我又不怎么熟手, 不过这些不能再麻烦老师了, 等我自己再消化后再来请教.

感谢!

主题最新回顾(发布时间：2009/5/10 21:50:00)
-- 作者：xjw01 -- 以下是引用蔡越在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 21:50:00)

-- 作者：xjw01
--

以下是引用蔡越在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 20:21:00)
-- 作者：xjw01 -- 在JPL中本来就可以计算日心坐标。我的代码中提取坐标时getCoordOne()，中心用DE_sun就可以了或在calc(double jd, int center)时也可以直接设置中心天体。其实calc()之后才能调用getCoordOne()

主题最新回顾(发布时间：2009/5/10 20:21:00)

-- 作者：xjw01
--

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

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

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

主题最新回顾(发布时间：2009/5/10 20:18:00)
-- 作者：xjw01 -- 我没想到你有兴趣调试。还是给你完整的工程包吧，这样调试起来比较方便。JPL的ascii星历文件你应该有了，我就不上传了，二进制星历文件可以用我的这个工程包生成。使用C++Builder6.0调试 http://www.fjptsz.com/xxjs/xjw/myjpl.rar [此贴子已经被作者于2009-5-10 21:45:44编辑过]

主题最新回顾(发布时间：2009/5/10 20:18:00)

-- 作者：xjw01
--

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

使用C++Builder6.0调试

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

[此贴子已经被作者于2009-5-10 21:45:44编辑过]