<续19楼> 日、月、地 三天体运动模拟程序设计
-----------------------------------------------------------------------
　　这个任务仍然挺大的，所以将按先后顺序进行下列程序设计
一。卫星围绕地球的椭圆轨道模拟(地心坐标系即假设地心不动)
二。月亮与地球相互绕行的轨道模拟(质心坐标系即假设二者的质心不动)
三。日、月、地 三天体运动模拟(日心坐标系也即假设太阳的中心不动)
-----------------------------------------------------------------------

　　卫星绕地球运动的模拟(地心坐标系并且设卫星的质量极小)

这里的卫星是指质量远小于中心天体地球的那样的卫星，如人造卫星；而月亮就
不满足这一要求，因为月球的质量已经达到地球的81.3分之一，对地球运动已经
造成不算小的影响，例如海洋的潮汐现象等。
　　[题]已知地球质量Ｍ＝5.976×10²⁴千克，半径Ｒ＝6371千米。有1颗卫星以
Ｖ₀(米/秒)的初速度、仰角为α(度)从离地Ｈ(千米)的高度发射出去，试求:此后
卫星的运动轨迹。(只考虑地心引力、忽略大气阻力等)
　　[思路]卫星的每一个瞬间，都可看成是在作“斜抛运动”。具体地说就是在
t₀→t₀+△t 时段内,质点以初速Ｖ₀(米/秒)和仰角α(度)在地心引力下作抛物线运
动。只要步长△t足够短，这样做误差就非常小。所谓“仰角”应当怎样理解呢？
首先你要建立“水平面”。卫星高高在上，固然离海平面十分遥远，但是你可以
设想将海平面沿着地心与卫星的连线方向平移Ｈ(千米)使与卫星重合，这就得到
了“水平面”。前述“仰角α(度)”就是卫星瞬间速度与水平面的夹角。在α大于
零的时候，是真正的仰角，此时卫星离地球越来越远；反之，在α小于零的时候，
实际是“俯冲”，此时卫星离地球越来越近。理论上，只有在α恒等于零的状态下
卫星才会围绕地球作精确的匀速圆周运动。应当指出，水平面并不是1个几何平面
所以假设卫星沿着与地心连线相垂直的切线方向笔直笔直的运动，则其仰角α并不
恒等于零，而是不断增大的。换句话说，每经历一个△t时间间隔，都要重新定位
“水平面”。
　　有人可能疑惑: 无数个α＝0的平抛运动怎么会形成圆轨道呢？事实上,每经过
△t 卫星沿水平方向移动了：△x＝Ｖ₀△t，竖直方向受地心引力作用反而下降了
△y＝0.5×ｇ△t²的距离,式中ｇ是卫星所在高度Ｈ处的重力加速度。用地心坐标
系描述，此时卫星的横坐标Ｘ＝△x＝Ｖ₀△t，纵坐标Ｙ＝Ｒ＋Ｈ－0.5×ｇ△t²
因而,此时卫星距地心 r＝sqrt(Ｘ²＋Ｙ²)。强令 r＝Ｒ＋Ｈ，即可解出圆周轨道
所需要的速度Ｖ(米/秒)。过程如下
　　　　Ｖ²△t²＋(Ｒ＋Ｈ－0.5×ｇ△t²)²＝(Ｒ＋Ｈ)²
　　　　Ｖ²△t²－2×(Ｒ＋Ｈ)×0.5×ｇ△t²＝０ (高阶小量已略去)
　　　　Ｖ²＝2×(Ｒ＋Ｈ)×0.5×ｇ
　　　　Ｖ＝sqrt[(Ｒ＋Ｈ)×ｇ]
[解]根据上述思路，兹建立解题步骤如下
　⑴输入初始参数：距地面高度Ｈ(千米)、速度Ｖ(米/秒)、仰角α(度)
　　必要时还可输入初始时刻ｔ＝ｔ0，让模拟更加逼真，否则ｔ＝０。
　　卫星初始位置：Ｘ＝０　Ｙ＝Ｒ＋Ｈ　 近点角θ＝０
　⑵确定ｔ时刻卫星速度的水平分量Ｖx和垂直分量Ｖy
　　Ｖx＝Ｖcosα　　Ｖy＝Ｖsinα
　⑶确定ｔ＋△ｔ时刻卫星的位置
　Ｘ(ｔ＋△ｔ)＝Ｘ(ｔ)＋Ｖx△ｔ
　Ｙ(ｔ＋△ｔ)＝Ｙ(ｔ)＋Ｖy△ｔ－ 0.5ｇ△t²
　　其中ｇ＝GM/Ｙ²，G为引力常数，M为中心天体地球的质量。
　⑷确定ｔ＋△ｔ时刻卫星的速度
　　Ｖx(ｔ＋△ｔ)＝Ｖx(ｔ)
　　Ｖy(ｔ＋△ｔ)＝Ｖy(ｔ)－ｇ△t
　　Ｖ(ｔ＋△ｔ)＝sqrt(Ｖx²＋Ｖy²)
　　α (ｔ＋△ｔ)＝arctan(Ｖy/Ｖx)
　　此外用机械能守恒原理确定Ｖ(ｔ＋△ｔ)也是可以的。
⑸确定ｔ＋△ｔ时刻卫星的近点角增量
　　△θ＝arctan[Ｘ(ｔ＋△ｔ)/Ｙ(ｔ＋△ｔ)]
θ(ｔ＋△ｔ)＝θ(ｔ)＋△θ
　⑹坐标旋转△θ角，即更新“水平面”。具体如下
　　新的Ｙ＝sqrt(Ｘ²＋Ｙ²)
　　新的Ｈ＝Ｙ－Ｒ
　　新的Ｘ＝０
　　新的α ＝α ＋△θ
　　注意，速率Ｖ不因坐标旋转而变化
　⑺新的ｔ＝ｔ＋△ｔ
　　输出高度Ｈ、速率Ｖ、仰角α 、时刻ｔ、近点角θ
　　如果 θ<2π 则转第⑵步；
　　否则,检验卫星旋转1圈后轨道参数是否等于初值。
显然，接近于初值的程度反映此算法是否有使用价值。

<下转29楼>

[此贴子已经被作者于2007-11-20 8:58:22编辑过]

匀速圆周运动...很多年没做过了...公式都忘了...G是多少来着?

初中你就学积分了,你强

/*----------------------------------------------------------------
　　卫星绕地球运动的模拟程序(地心坐标系并且设卫星的质量极小)
----------------------------------------------------------------*/
#include <stdio.h>
#include <math.h>

const double dt = 10e-3; //步长10毫秒
const double TPI=2*3.1415926535897932;
const double GM=6.672e-11 * 5.976e+24;
const double Re=6371000; //地球平均半径(米);

double g(double r)//距地心r处的重力加速度
{
return GM/(r*r);
}

void time(double tim,int*hh,int*mm,int*ss)
{
long sec=(long)(tim+0.5);
*hh = sec/3600;
*mm = sec%3600/60;
*ss = sec%60;
}

int main( )
{
double H; //卫星距地球表面的高度
double V; //卫星相对于地心的速率
double alpha; //卫星运动的仰角
double tim=0; //计时器(秒)
double the=0; //卫星的真近点角
double S=0.0; //卫星的当前里程
double r; //卫星的地心距离
double dth; //真近点角的增量
int second=0; //存放前一秒钟
int hh,mm,ss; //小时、分、秒
alpha=0; //仰角的缺省值
{ //这些红色代码可验证本模拟程序,删除掉也可以
double Hmax,Hmin,Hav,aa,bb,cc,ee,Vmax,Vmin,T;
printf("近地点高度(千米)、远地点高度(千米)＝");
scanf("%lf%*c%lf",&Hmin,&Hmax);
Hmin*=1e3;Hmax*=1e3;
Hav=(Hmax+Hmin)/2;
aa=Re+Hav;
cc=(Hmax-Hmin)/2;
bb=sqrt(aa*aa-cc*cc);
ee=cc/aa;
printf("\nａ＝%.0lf千米，ｂ＝%.0lf千米，ｃ＝%.0lf千米，ｅ＝%.3lf\n\n",
aa/1e3,bb/1e3,cc/1e3,ee);
printf(" %.3lf \n",bb*bb/aa/1e3);
printf("极坐标方程为ｒ＝────────\n");
printf(" １＋%0.3lfcosθ\n\n",ee);

Vmax=sqrt(GM*(2/(aa-cc)-1/aa)); //用活力公式计算近地点速度
Vmin=sqrt(GM*(2/(aa+cc)-1/aa)); //用活力公式计算远地点速度
printf("Ｖmax＝%.3lf(米/秒)，Ｖmin＝%.3lf(米/秒)\n\n",Vmax,Vmin);

T=TPI/sqrt(GM)*pow(aa,1.5); //理论周期(秒)
time(T,&hh,&mm,&ss);
printf("Ｔ＝%02d:%02d:%02d\n\n",hh,mm,ss);
printf("press Enter key to continue....");
getchar();getchar();
}
printf("\n高度(千米) 速度(米/秒) 仰角(度) ＝");
scanf("%lf%*c%lf%*c%lf",&H,&V,&alpha);
alpha*=TPI/360;
H*=1000;
r=Re+H;
printf("高度(千米) 速度(米/秒) 仰角(度) θ角(度) hh:mm:ss 里程(千米)\n");
do
{
double ds=V*dt;
double sa=sin(alpha),ca=cos(alpha),gt=g(r)*dt;
double r0=r;
r=sqrt(r*r+ds*ds+2*r*ds*sa)-0.5*gt*dt;
dth=V*ca*dt/r;//坐标系随卫星旋转dθ角
S+=sqrt(pow(r-r0,2)+r*r0*dth*dth);
V-=gt*sa; //t＋dt时刻的速度
alpha+=dth-gt*ca/V; //刷新仰角
the+=dth; //刷新卫星的近点角
tim+=dt; //刷新计时器的秒数
time(tim,&hh,&mm,&ss);
if(ss != second){ second=ss;
printf("%-11.2f%-12.3f%-9.3f%-9.3f%02d:%02d:%02d %-9.2f\n",
(r-Re)/1e3,V,alpha/TPI*360,the/TPI*360,hh,mm,ss, S/1e3);}
}
while(the<TPI);
S*=TPI/the;
tim*=TPI/the;
printf("轨道周长＝%.3f千米\n",S/1e3);
printf("轨道周期＝%.2f秒\n",tim);
printf("平均速度＝%.1f米/秒\n",S/tim);
return 0;
}

[运行实录，红字是解释]

近地点高度(千米)、远地点高度(千米)＝439,2384 下划线表示输入。这是我国第一颗人造卫星的数据。

半长径ａ＝7783千米，半短径ｂ＝7721千米，焦距ｃ＝973千米，偏心率ｅ＝0.125

7660.977
极坐标方程为ｒ＝──────── r表示从卫星到地心的距离、θ为近点角。
１＋0.125cosθ

Ｖmax＝8115.743(米/秒)，Ｖmin＝6312.759(米/秒) 分别为近地点理论速度和远地点理论速度。

Ｔ＝01:53:52 这是根据开普勒定律求出的理论周期1小时53分52秒，用于校验模拟运行的准确程度。

press Enter key to continue....

高度(千米) 速度(米/秒) 仰角(度) ＝ 439 8115.743 0 将近地点高度和近地点理论速度以及
零仰角作为模拟程序的输入数据。

高度(千米) 速度(米/秒) 仰角(度) θ角(度) hh:mm:ss 里程(千米)
439.00 8115.743 0.004 0.034 00:00:01 4.06
439.00 8115.742 0.011 0.102 00:00:02 12.17
439.00 8115.739 0.019 0.171 00:00:03 20.37
439.01 8115.736 0.027 0.240 00:00:04 28.49
439.01 8115.731 0.034 0.308 00:00:05 36.60
439.02 8115.726 0.042 0.376 00:00:06 44.72
439.02 8115.719 0.049 0.445 00:00:07 52.83
439.03 8115.711 0.057 0.513 00:00:08 60.95
439.04 8115.702 0.065 0.581 00:00:09 69.06
439.05 8115.692 0.072 0.649 00:00:10 77.18
。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。
2383.99 6312.765 0.014 179.904 00:56:54 24339.08
2384.00 6312.764 0.008 179.946 00:56:55 24345.39
2384.00 6312.763 0.002 179.987 00:56:56 24351.70 模拟结果
远地点2384千米，速度6312.763米/秒，仰角为0，56分56.3秒跑24353.7千米
2384.00 6312.763 -0.004 180.028 00:56:57 24358.02
2384.00 6312.764 -0.010 180.069 00:56:58 24364.33
2383.99 6312.765 -0.016 180.111 00:56:59 24370.64
。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。
439.04 8115.705 -0.047 359.580 01:53:46 48657.56
439.03 8115.711 -0.039 359.649 01:53:47 48665.68
439.03 8115.717 -0.031 359.717 01:53:48 48673.79
439.02 8115.721 -0.024 359.785 01:53:49 48681.91
439.02 8115.724 -0.016 359.853 01:53:50 48690.03
439.02 8115.726 -0.009 359.922 01:53:51 48698.14
439.02 8115.726 -0.001 359.990 01:53:52 48706.26
模拟结果 近地点439.02千米，速度8115.726米/秒，仰角为0
而初始数据 近地点439.00千米，速度8115.743米/秒，仰角为0
由此可见 模拟1圈下来，近地点虚高20米，速度丢失0.017米/秒
轨道周长＝48707.455千米 跑完全程48707.455千米
轨道周期＝6831.65秒 历时:1小时53分51.65秒 与理论周期一致
平均速度＝7129.7米/秒