/* --------------------------------------------------------------------------------------
出自： 编程中国 http://www.bc-cn.net
作者： yu_hua
时间： 2007-11-8 编程论坛首发
声明： 尊重作者劳动，转载请保留本段文字

“嫦娥一号”发射这么久了，已经变成了绕月运行的卫星。
下面给出几个简单的不考虑摄动影响下的轨道参数计算公式。

1。轨道周期T的公式：T = 3.3小时 / sqrt(中心天体的平均密度，克/立方厘米)

⑴卫星绕地球运行时：因为地球平均密度为5.517克/立方厘米，所以不考虑空气阻力
近地圆轨道的运行周期
T₀ = 3.3×60分钟 / sqrt(5.517) = 84.3分钟 = 5058秒，从而第一宇宙速度为
V₀ = 2π×6371千米 / T₀ = 7.914千米/秒，地球的第二宇宙速度(逃逸速度)则为
V_∞ = V₀×sqrt(2) = 11.19千米/秒

⑵如果圆轨道的高度H(千米)不能忽略不计，则绕地圆轨道的运行周期
T_H = T₀×(1 + H/R)^1.5 = 84.3分钟×(1 + H千米/6371)^1.5
式中，中心天体即地球的平均半径R = 6371千米。

⑶通常卫星绕地球运行时的轨道不是正圆，而是椭圆，则上式中轨道高度应理解为
平均高度 H = (近地点高度 ＋ 远地点高度) / 2
应当指出，如果把平均高度Ｈ加上中心天体半径Ｒ，恰好就是椭圆轨道的半长径a。
可见这样做理论上是精确的，因为符合开普勒行星运动第三定律。

⑷针对嫦娥一号给出两个算例：
①已知：近地点高度600千米，远地点高度51000千米，试求环绕地球一圈的周期。
平均高度 H = (600 ＋ 51000) / 2 = 25800 千米
绕行周期 T_H = 84.3分钟×(1 + 25800/6371)^1.5 = 956.56分钟 ≈ 16小时
②已知：近月点高度200千米，远月点高度8600千米，试求环绕月球一圈的周期。
平均高度 H = (200 ＋ 8600) / 2 = 4400 千米
绕行周期 T_H = 108.2分钟×(1 + 4400/1738)^1.5 = 718.1分钟 ≈ 12小时
式中，1738千米为月球的半径。108.2分钟为近月轨道周期，推导如下：
T₀ = 3.3×60分钟 / sqrt(3.348)=108.2分，式中：3.348为月球平均密度。
所以，月球的“第一宇宙速度”为
V₀ = 2π×1738千米 / T₀ = 1.682千米/秒
相应地，月球的“逃逸速度”为
V_∞ = V₀×sqrt(2) = 2.38千米/秒

由于重力场是个保守力场，所以卫星的机械能守恒，从而随着与中心天体间的
距离的增大，重力势能上升，卫星动能下降，速度变慢。于是，过近地点时卫
星的速度最快，而飞过远地点时卫星的速度最慢。故卫星绕中心天体运行时的
2。瞬间速率V的公式：V米/秒 = sqrt[Ｇ×Ｍ吨×(2/r千米 - 1/a千米)]
式中，引力常数Ｇ＝6.672×10^-11，Ｍ为中心天体的质量。
已知：地球质量Ｍ_地＝5.976×10²¹吨，月球质量Ｍ_月＝7.35×10¹⁹吨。
还有：r 为卫星到地心或月心的距离(千米)，a 为椭圆轨道的半长径(千米)

仍针对嫦娥一号给出两个算例：
①已知：近地点高度600千米，远地点高度51000千米，试求嫦娥一号经过这两个
点时的瞬间速率V_近和V_远。
平均高度 H = (600 ＋ 51000) / 2 = 25800 千米
轨道半长径 a = R + H = 6371 + 25800 = 32171 千米
与地心距离 r近 = R + 600 = 6971 千米
r远 = R + 51000 = 57371 千米
从而 V_近 = sqrt[6.672×10^-11×5.976×10²¹×(2/6971 - 1/32171)]
= 10099.5 米/秒 ≈ 10.1 千米/秒
同理 V_远 = sqrt[6.672×10^_-11×5.976×10²¹×(2/57371 - 1/32171)]
= 1227 米/秒 = 1.227 千米/秒
②已知：近月点高度200千米，远月点高度8600千米，试求嫦娥一号经过这两个
点时的瞬间速率V_近和V_远。
平均高度 H = (200 ＋ 8600) / 2 = 4400 千米
轨道半长径 a = R_月 + H = 1738 + 4400 = 6138 千米
与月心距离 r_近 = R_月 + 200 = 1938 千米
r_远 = R_月 + 8600 = 10338 千米
从而 V_近 = sqrt[6.672×10^-11×7.35×10¹⁹×(2/1938 - 1/6138)]
= 2064 米/秒 = 2.064 千米/秒
同理 V_远 = sqrt[6.672×10^-11×7.35×10¹⁹×(2/10338 - 1/6138)]
= 387 米/秒 = 0.387 千米/秒

3。尝试用C语言编写一个程序，对以上两项结果进行动力学验证。换句话说，
就是用数值积分的方法做下列定积分：
_远地点
∫ dt
^近地点
看看这个定积分的值是否非常接近于周期Ｔ的二分之一。
[提示]速度的瞬间方向总是沿着椭圆的切线方向。此外，本题不能
　　　用开普勒第二定律(面积定律)算，否则将发生循环论证错误。

[参考答案]
物理学告诉我们：时间＝距离/速度，即 dt＝ds/V
微积分指出，极坐标系的线元 ds＝sqrt[(dr)²＋(rdθ)²]＝sqrt[(r')²＋r²]dθ
所以前面给出的那个定积分就演变成

远地点 π
∫ dt ＝∫sqrt[(r')²＋r²]dθ/V
近地点 0

式中 r'＝dr/dθ 即r对θ的一阶导数。下面的工作就是导出 r(θ)的函数关系。
----------------------------------------------------------------------
由椭圆的参数方程
x = a·cosφ
y = b·sinφ
可以证明：椭圆的曲率半径 R_曲＝(a²sin²φ＋b²cos²φ)^1.5/(ab)
当参数φ＝0度，得到曲率半径的最小值：R_min＝b²/a
当参数φ＝90度,得到曲率半径的最大值：R_max＝a²/b
椭圆轨道的极坐标方程是：r = R_min/(1＋e·cosθ) ＝ b²/(a＋c·cosθ)
式中，偏心率e＝c/a，R_min为曲率半径最小值。
设中心天体的平均半径为R，卫星(相对于中心天体表面)的高度为H，另外
用下标max和min以区分最大高度和最小高度，用H_ave表示平均高度。于是
轨道半长径 a ＝ R＋H_ave＝ R＋(H_max＋H_min)/2
此外 焦距 c ＝ (H_max－H_min)/2
从而半短径 b ＝ sqrt(a²－c²) ＝ sqrt((R＋H_max)×(R－H_min))
椭圆偏心率 e ＝ (H_max－H_min)/(2R＋H_max＋H_min)
曲率半径 R_min ＝ (R＋H_max)(R＋H_min)/(R＋H_ave)
轨道的极坐标方程还可写成
r(θ)＝(R＋H_max)(R＋H_min)/[R＋H_ave＋0.5×(H_max－H_min)cosθ]
----------------------------------------------------------------------
经推导，r(θ)的一阶导数 r'(θ)＝er²sinθ/R_min
于是 sqrt[(r')²＋r²]＝r·sqrt[1+(e·r·sinθ/R_min)²]
＝sqrt(1＋e²＋2e·cosθ)·R_min/(1＋e·cosθ)²
又因 卫星速率 V＝sqrt[ＧＭ(2/r－1/a)]
所以 被积函数
f(θ)＝sqrt(1＋e²＋2e·cosθ)·R_min/(1＋e·cosθ)²/sqrt[ＧＭ(2/r－1/a)]
其中 r＝R_min/(1＋e·cosθ)
π
要完成的积分就是 ∫f(θ)dθ ，并验证积分的结果是否非常接近Ｔ/2
0
第16楼的C程序完成了上述任务。这个程序很“土”，相当于梯形法定积分。
将整个积分区间 0°≤θ≤180°切割为10800等分,子区间宽度△θ为１角分。
运行结果举例：
远地点51000km，近地点600km，轨道长度166133km，卫星周期15小时57.0分钟(嫦娥1号)
远地点71150km，近地点600km，轨道长度210103km，卫星周期24小时0.0分钟 (嫦娥1号)
远地点2384km，近地点439km，轨道长度48707km，卫星周期113.9分钟(中国第1颗人造卫星)
远地点964km，近地点228km，轨道长度43744km，卫星周期96.4分钟(前苏联第1颗人造卫星)

4。月球公转轨道(要同时考虑地球和太阳的影响)计算机模拟。

5。嫦娥一号的地月转移轨道的计算机模拟。这种情况下,必须考虑月球引力场
对嫦娥一号的影响。不过，为了不要把事情弄得过于复杂，我们假设：月球的
公转轨道平面(即白道)与嫦娥一号的地月转移轨道是同一个平面；并且，白道
是一个正圆(半径约为384400千米)、月球的平均公转速度可取1.023千米/秒。

朋友，您愿意接受这一挑战吗？俺暂时不公布答案！！！

-------------------------------------------------------------------------------------- */

[此贴子已经被作者于2007-11-12 19:18:57编辑过]

bai chi yi ge, except big words.

斑竹也爱灌水？有本事来点干货。

//楼主问题3的C语言程序公布如下
#include <stdio.h>
#include <math.h>
double GM=6.672e-11 * 5.976e+21; //如果把 5.976e+21 改成 7.35e+19，
double R=6371; //地球的平均半径 //并把6371改成1738，就可算月球卫星。
int main( )
{ double T=0; //秒数累加
double S=0; //里程累加
long i, Nmax=180*60L;
double dth=3.1415926535898/Nmax; //theta增量
double ds; //极坐标下的弧长元
double dt; //时间增量
double Hmax,Hmin,Have,Rmin;
double aa,ee,hh,rr,the,V,temp;
printf("绕地球飞行时的最大、最小高度(km)：");
scanf("%lf%*c%lf",&Hmax,&Hmin);
Have = (Hmax+Hmin)/2; //平均高度
aa = R+Have; //椭圆轨道的半长径
ee = (Hmax-Hmin)/2/aa; //偏心率，也称离心率
Rmin = (R+Hmax)*(R+Hmin)/aa; //最小曲率半径
for(i=0; i<Nmax; i++)
{
the = (i+0.5)*dth; //确定theta值
temp=1+ee*cos(the); //中间变量
rr = Rmin/temp; //卫星与地心之间的距离(千米)
hh = rr-R; //卫星相对于地面的高度(千米)
V=sqrt(GM*(2/rr - 1/aa)); //卫星的瞬时速率(米/秒)
ds=dth*Rmin*sqrt(2*temp-1+ee*ee)/(temp*temp)*1000; //以米为单位
dt=ds/V;
T+=dt;
S+=ds; //以米为单位
}
S+=S;
printf("轨道长度的理论值是%0.1lf公里\n",S/1e3);
T+=T; //加上另外半圈所需要的绕行时间
T/=60; //以分钟为单位
if(T<180)
printf("卫星周期的理论值是%0.1lf分钟\n",T);
else
{
long Hour=(long)(T/60);
T-=60*Hour;
printf("卫星周期的理论值是%ld小时%.1lf分钟\n",Hour,T);
}
return 0;
}

[此贴子已经被作者于2007-11-10 21:46:21编辑过]

　　网络上看到北京网友一个帖子，大意说太阳对月球的引力大于地球对月球的引力，为什么
月球只围绕地球转圈子呢？所以近期我如有空闲将优先模拟日、月、地三天体的绕行模拟程序。
-----------------------------------------------------------------------------------
首先实地计算一下，月球受到的两种引力分别有多大？
已知太阳质量Ｍ_日＝1.989×10³⁰千克、地球质量Ｍ_地＝5.976×10²⁴千克、Ｍ_月＝7.35×10²²千克
月日平均距离 r_日＝日地平均距离＝1天文单位＝1.496×10¹¹米
月地平均距离 r_地＝3.844×10⁸米
引力常数Ｇ＝6.672×10^-11
根据万有引力定律：Ｆ＝Ｇm₁m₂/r²
太阳对月球的引力: Ｇ×Ｍ_日×Ｍ_月/r_日²＝4.358×10²⁰牛顿
地球对月球的引力: Ｇ×Ｍ_地×Ｍ_月/r_地²＝1.983×10²⁰牛顿
可见太阳对月球的引力约为地月引力的2.2倍，北京网友所言非虚。然而北京这位网友却荒唐地
认为日、月之间万有引力定律失效。实际上当人们说月球的公转轨道是个近地点363300千米
远地点405500千米的椭圆时，参照系是地球,不是太阳。如果站在太阳的立场上看月球运动，
则看到的是两种运动的叠加:⑴在陪着地球一块绕日旋转的基础上追加⑵相对于地球的旋转运动。
显然与太阳相比，地球算不上好的惯性参照系。推而广之,太阳与银心(银河系中心)相比，也算
不上好的惯性参照系呢:据天文观测太阳系每时每刻以每秒300公里的速度绕着银心作近似的匀速
圆周运动，周期约为两亿年！这表明
太阳系与银心之间的距离 r_银＝300千米/秒×(3600秒×24)×365×2×10⁸/2π＝30×10¹⁶千米
哇噻，三十亿亿公里！一个天文数字。光线在真空中每秒走三十万公里，两者除一下，得到：
一万亿秒(人的寿命仅为20～30亿秒)约合三万二千年，故太阳系到银心的距离是3.2万光年。
根据向心力与万有引力相等的物理学原理，可以求得
银心质量Ｍ_银＝r_银Ｖ²/Ｇ＝40×10⁴⁰千克，约合太阳质量的两千亿倍。
当然并不像太阳质量那样地高度聚集，这里的所谓银心质量分布在半径为3.2万光年的球内。
所以天文学家有理由宣称银河系大约由两三千亿颗恒星组成。具有讽刺意味的是，银河系的
平均密度却小得可怜，不妨计算一下：Ｍ_银/[4π/3×(r_银)³]＝3.5×10^-21千克/米³
请注意：标准状况下空气密度为1.29千克/米³。星空、星空，原来真的很“空”。
------------------------------------------------------------------------------------
<下转23楼>

[此贴子已经被作者于2007-11-14 10:06:24编辑过]

<续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编辑过]

/*----------------------------------------------------------------
　　卫星绕地球运动的模拟程序(地心坐标系并且设卫星的质量极小)
----------------------------------------------------------------*/
#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米/秒

/*----------------------------------------------------------------
卫星绕地球运动的模拟程序(地心坐标系并且设卫星的质量极小)
<改进29楼程序的时分秒处理>
----------------------------------------------------------------*/
#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 hms(double tim,int*hh,int*mm,double*ss)
// 将tim中的秒数化成 hh小时 mm分 ss秒
{
*hh = (int)(tim/3600);
tim = tim - *hh*3600 ;
*mm = (int)(tim / 60);
*ss = tim - *mm * 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; //小时、分
double 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ａ＝%.0f千米，ｂ＝%.0f千米，ｃ＝%.0f千米，ｅ＝%.3f\n\n",
aa/1e3,bb/1e3,cc/1e3,ee);
printf(" %.3f \n",bb*bb/aa/1e3);
printf("极坐标方程为ｒ＝────────\n");
printf(" １＋%0.3fcosθ\n\n",ee);

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

T=TPI/sqrt(GM)*pow(aa,1.5); //理论周期(秒)
hms(T,&hh,&mm,&ss);
printf("Ｔ＝%02d:%02d:%.3f\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("\n高度(千米) 速度(米/秒) 仰角(度) θ角(度) hh:mm:ss.cc 里程(千米)\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; //刷新计时器的秒数
hms(tim,&hh,&mm,&ss);
if(fabs(ss-second)>=1){ second=(int)ss;
printf("%-11.2f%-12.3f%-9.3f%-9.3f%02d:%02d:%5.2f %-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("轨道周期＝%.3f秒\n",tim);
printf("平均速度＝%.3f米/秒\n",S/tim);
return 0;
}

基于固定的直角坐标系的小质量卫星轨道模拟程序设计

[题]地球半径R=6378千米，引力常数与地球质量乘积GM = 398603×10⁶(千米³/千秒²)。用固定的
直角坐标系模拟小质量卫星的运行轨道，仅考虑地球对卫星的吸引力。卫星轨道的近地点高Hmin和
远地点高Hmax作为输入参数。
[解]就本题而言，轨道形状是已知的：椭圆。像解析几何那样，把二维直角坐标系的原点定义在
椭圆的对称中心，从而让轨道的方程为
(x/a)²＋(y/b)²＝1
式中，半长径a＝R＋(Hmax＋Hmin)/2
半短径b＝(a²－c²)^1/2
焦距c＝(Hmax－Hmin)/2
我们将地球放到右焦点上，即令地心的直角坐标为(c ,0)。这样一来，近地点直角坐标为(a ,0)，
远地点则为(-a ,0)。下面制订ｔ→ｔ+△ｔ时段内，卫星的飞行轨迹以便建立迭代方案。
设ｔ时刻卫星的位置矢量是r、速度矢量是v，ｔ+△t时刻到达r＋△r。
最终，确立下列循环求解方案

⑴输入(或由理论确定)卫星的初始位置r与初速度v

-GM·(r－re)
⑵计算加速度a₁(r)＝─────── 式中re表示地心的位置矢量
|r－re|³

⑶粗算ｔ+△ｔ时刻卫星位置 r＋＝v△t (暂不考虑含△t平方的项)

⑷再算此位置的加速度a₂(r)，公式不变，但卫星的位置r有变化

⑸求出速度增量△v＝(a₁＋a₂)·(△t/2)

⑹ｔ+△t时刻的位置 r＋＝△v·(△t/2)

⑺ｔ+△t时刻的速度 v＋＝△v

⑻如果还没有绕行完一圈，则回到第⑵步

/*----------------------------------------------------------------
　　 小卫星绕地球运动的模拟程序(基于二维直角坐标系)

程序流程：输入近地点和远地点高度，由解析几何确定椭圆轨道的
　　方程，再用“活力公式”确定卫星过这两个点时的理论速度。模拟
　　程序将地心设定在右焦点(c,0)把近地点(a,0)作为模拟运动的初始
　　位置，把近地点的理论速度Vmax作为初速度(仰角等于零)。卫星在
　　XOY平面内沿椭圆轨道逆时针运动。

　　　　　 以下用的是“千米•吨•千秒”制
-----------------------------------------------------------------*/
#include <stdio.h>
#include <math.h>

const double dt = 10e-6; //步长10毫秒
const double TPI=2*3.1415926535897932;
const double GM=398603e6;//6.672e-11*5.976e+21吨
const double Re=6378; //地球平均半径6371(千米)

struct Vector
{ double x,y; } re; //地心的直角坐标

double distance(struct Vector rp,struct Vector rs)
// 从源点rs到场点rp的距离
{
return sqrt(pow(rp.x-rs.x,2)+pow(rp.y-rs.y,2));
}

struct Vector acc( struct Vector r )
// 由地心引力产生的加速度(矢量)
{ static
struct Vector acc;
double temp = -GM/pow(distance(r,re),3);
acc.x = temp*(r.x-re.x);
acc.y = temp*(r.y-re.y);
return acc;
}

void hms(double sec,int*hh,int*mm,double*ss)
// 将 sec 中的秒数化成 hh 小时 mm 分 ss 秒
{
*hh = (int)(sec/3600);
sec = sec - *hh*3600 ;
*mm = (int)(sec / 60);
*ss = sec - *mm * 60 ;
}

double Atan2( double y,double x )
// 带累计功能的反正切(可大于2π)
{
static double old=0,turns=0;
double now=atan2(y,x);
if(now<old)turns+=TPI;
old=now;
return now+turns;
}

int main( )
{
struct Vector rr,vv,dv,ac1,ac2;
double H; //卫星距地球表面的高度(千米)
double V; //卫星相对于地心的速率(米/秒)
double alpha; //卫星运动的仰角
double tim=0; //计时器(秒)
double the=0; //卫星的真近点角
double S=0.0; //卫星的当前里程(千米)
double r; //卫星的地心距离(千米)
int second=0; //存放前一秒钟
int hh,mm; //小时、分
double ss; //秒(含小数部)
double dt2=dt/2; //步长之半
double Hmax,Hmin,Hav,aa,bb,cc,ee,Vmax,Vmin,T;

printf("近地点高度(千米)、远地点高度(千米)＝");
scanf("%lf%*c%lf",&Hmin,&Hmax);
Hav=(Hmax+Hmin)/2;
aa=Re+Hav;
cc=(Hmax-Hmin)/2;
re.x=cc;re.y=0;
bb=sqrt(aa*aa-cc*cc);
ee=cc/aa;
printf("\nａ＝%.0f千米，ｂ＝%.0f千米，ｃ＝%.0f千米，ｅ＝%.4f\n\n",aa,bb,cc,ee);
printf(" %.3f千米 \n",bb*bb/aa);
printf("极坐标方程为ｒ＝────────\n");
printf(" 1＋%0.4fcosθ\n\n",ee);
Vmax=sqrt(GM*(2/(aa-cc)-1/aa)); //用活力公式计算近地点速度
Vmin=sqrt(GM*(2/(aa+cc)-1/aa)); //用活力公式计算远地点速度
printf("Ｖmax＝%.3f(米/秒)，Ｖmin＝%.3f(米/秒)\n\n",Vmax,Vmin);
T=TPI/sqrt(GM)*pow(aa,1.5); //理论周期(千秒)
hms(T*1e3,&hh,&mm,&ss);
printf("Ｔ＝%02d:%02d:%.3f\n\n",hh,mm,ss);
printf("press Enter key to continue....");
getchar();getchar();

alpha=0; //近地点处仰角等于零
r=Re+Hmin; //近地点到地心(千米)
rr.x=cc+r; rr.y=0; //得到位置矢量的初值
vv.x=0; vv.y=Vmax; //得到速度矢量的初值
printf("\n高度(千米) 速度(米/秒) 仰角(度) θ角(度) hh:mm:ss.cc 里程(千米)\n");
do
{
double product;
struct Vector r0=rr;
ac1 = acc(rr); //起点加速度矢量
rr.x+=vv.x*dt;
rr.y+=vv.y*dt;
ac2 = acc(rr); //终点加速度矢量(有误差)
dv.x=(ac1.x+ac2.x)*dt2;
dv.y=(ac1.y+ac2.y)*dt2;
rr.x+=dv.x*dt2;
rr.y+=dv.y*dt2;
vv.x+=dv.x;
vv.y+=dv.y;
product=vv.x*(rr.x-re.x)+vv.y*(rr.y);
r=distance(rr,re); //t＋dt时刻地心距
H=r-Re; //t＋dt时刻的高度
V=sqrt(vv.x*vv.x+vv.y*vv.y);//t＋dt时刻的速度
alpha=asin(product/(r*V)); //t＋dt时刻的仰角
the=Atan2(rr.y,rr.x-re.x); //t＋dt时刻近点角
S+=distance(rr,r0); //t＋dt时刻的里程
tim+=dt; //t＋dt时刻千秒数
hms(tim*1e3,&hh,&mm,&ss);
if(fabs(ss-second)>=1){ second=(int)ss;
printf("%-11.2f%-12.3f%-9.3f%-9.3f%02d:%02d:%5.2f %-9.2f\n",
H,V,alpha/TPI*360,the/TPI*360,hh,mm,ss,S); }
}
while(the<TPI);
S*=TPI/the;
tim*=TPI/the;
printf("轨道周长＝%.3f千米\n",S);
printf("轨道周期＝%.3f秒\n",tim*1e3);
printf("平均速度＝%.3f米/秒\n",S/tim);
return 0;
}

/*
本程序的准确度高于31楼的，在轨道形状、轨道周期
和能量守恒诸方面都与理论值吻合得非常好。这就为
多天体相互绕行模拟程序的设计奠定了一个良好基础
*/

[此贴子已经被作者于2007-11-21 20:16:05编辑过]