地月二体轨道模型

3. 2地月二体轨道模型
假定卫星在地球中心引力场运动,忽略各种摄动力因素影响的卫星轨道称为二体轨道。

此时，卫星的运动轨迹是一个不变的椭圆轨道，分析这种轨道的特性称为二体问题。

其相应的运动方程为：
223d dt r
μ
=-r r
(1)
称之为开普勒方程[4]，其中μ（32398601.19/km s μ=）为地心引力常数。

为仿真计算的需要，下面把各摄动力在地心惯性坐标系下用力的分量形式表示出来，并用下标,,x y z 分别表示力在i F 坐标系中的三个坐标方向上的三个分量。

由(1)式可以直接得到0F 的分量表示：
0003
3
3
,,x y z F x F y F z r r
r
μ
μ
μ
=-
=-
=-
(2)
3. 3理想模型
假定月球在地球中心引力场运动,忽略各种摄动力因素影响，地球为圆形球体，那么月球在地心惯性坐标系中的位置和速度分别为[4] ：
20cos cos()sin sin()cos (1)sin cos()cos sin()cos 1cos sin()sin x f f i a e y f f i e f z f i ωωωωωΩ+-Ω+⎡⎤
⎡⎤
-⎢⎥⎢⎥==
Ω++Ω+⎢⎥⎢⎥+⎢⎥⎢⎥+⎣⎦
⎣⎦
r
(3)
]0sin (cos )(1f e f a μ⎡⎤
⎢⎥==-++⎢⎥-⎢⎥⎣⎦
x v y z E Q
(4)
cos cos sin sin cos cos sin sin cos cos sin sin i i i ωωωωωΩ-Ω⎡⎤
⎢⎥=Ω+Ω⎢⎥
⎢⎥⎣⎦
E
(5)
sin cos cos sin cos sin sin cos cos cos cos sin i i i ωωωωω-Ω-Ω⎡⎤
⎢⎥=-Ω+Ω⎢⎥
⎢⎥⎣⎦
Q
(6)
由上面四个公式可得月球在地心惯性坐标系中假设初始点的位置r 和速度v （参看附录1）分别为
（
-173450,72740,-311210）
(7)
和
（
0.2011,-0.9981,-0.3454）
(8)
记0000(,,)x y z =r 和0000(,,)x y z v v v =v 分别为月球的位置矢量和速度矢量。

由于(1)式等价于[5]： 00x dx v dt = 00y dy v dt = 00z dz
v dt
= 003x dv x dt r μ=- 003y dv y dt r μ=- 003z dv z
dt r
μ=- 记300000(,),(,)(,/)f t r μ==-r v v r X X ，则上面6个方程构成的方程组等价于：
()
00/,d dt f t =X X
(9)
取采样周期（采样时间间隔）为N ，把(10)式离散化得：
001011[]
(,)
k k k k f t N ----=X X X
(10)
即：
0010101(,)k k k k f t N
---=+X X X
(11)
其中00()k k t =X X 为k t 时刻位置速度向量。

利用(11)式和初始时刻的位置和速度就可以求出以后各时刻的位置和速度，从而确定月球轨道。

已知月球的轨道平均周期为29.53天，使用MATLAB 软件编程[6]计算出任一周期的各个时刻的月球在地心坐标系下的位置矢量，自探月卫星发射时刻开始的以1秒为间隔的各时刻月球位置坐标参看表1.
表1 各时刻月球位置坐标
图2是自探月卫星发射时刻开始的以1秒为间隔模拟的一个周期的月球轨道示意图：
附录
附录2 计算月球初始位置和速度
f=0;
图2 月球绕地进动轨道示意图
5
u=398601.19;
r=a*(1-e^2)/(1+e*cos(f))*[cos(ou)*cos(w+f)-sin(ou)*sin(w+f)*cos(i);si n(ou)*cos(w+f)+cos(ou)*sin(w+f)*cos(i);sin(w+f)*sin(i)];
E=[cos(w)*cos(ou)-sin(w)*sin(ou)*cos(i);cos(w)*sin(ou)+sin(w)*cos(ou) *cos(i);sin(w)*sin(i)];
Q=[-sin(w)*cos(ou)-cos(w)*sin(ou)*cos(i);-sin(w)*sin(ou)+cos(w)*cos(o u)*cos(i);cos(w)*sin(i)];
v=sqrt(u/a/(1-e^2))*((-sin(f)*E)+(e+cos(f))*Q);
r
v
附录3（画出月球轨道并求出各个时刻月球位置与速度）
clear
clc
u=398601.19;
u2=u/81;
rs2=rs*9/10;
vs2=[-0.8563,4.5805,-1.7957];
x0=[rs2,vs2];
r2=sqrt((rs2(1)+rs(1))^2+(rs2(2)++rs(2))^2+(rs2(3)++rs(3))^2);
N=60;
n=29.53*24*60;
for k=1:n-1;
a(k)=(X(k+1)-X(k))/N; b(k)=(Y(k+1)-Y(k))/N;c(k)=(Z(k+1)-Z(k))/N; vs(k,:)=[a(k),b(k),c(k)];
end
for k=1:n-2;
a1(k)=(a(k+1)-a(k))/N;
b1(k)=(b(k+1)-b(k))/N;c1(k)=(c(k+1)-c(k))/N;
vs1(k,:)=[a1(k),b1(k),c1(k)];
end
for k=1:n-2;
x=x0+[vs2,(-u/r2^3-u2/r2^3)/2*rs2-u2/r2^3/2*rs-vs1(k,:)]*N; X2(k)=x(1); Y2(k)=x(2); Z2(k)=x(3);
x0=x;
r2=sqrt((rs2(1)+rs(1))^2+(rs2(2)++rs(2))^2+(rs2(3)++rs(3))^2); rs2=x0(1:3);
vs2=x0(4:6);
end
plot3(X2,Y2,Z2)
hold on
comet3(X2,Y2,Z2)。