用matlab代码实现，已知服务星轨道六根数的半长轴为7000km，轨道倾角为20度，剩下四个参数都为零，目标星的轨道六根数的半长轴为8000km，轨道倾角为10度，剩下四个参数都为零，服务星利用兰伯特转移通过两次机动，使服务星最终与目标星相同位置与相同速度，将服务星的运动轨迹与目标星的椭圆形的运动轨迹画在同一幅三维图上，并画出过程中时间与服务星与目标星相对距离逐渐减小到0的过程的图

以下是Matlab代码实现，注释中有详细说明：

% 定义常数
mu = 3.986e5; % 地球引力常数
Re = 6378; % 地球半径

% 服务星初始轨道六根数
a1 = 7000; % 半长轴
e1 = 0; % 离心率
i1 = 20*pi/180; % 轨道倾角
Om1 = 0; % 升交点赤经
om1 = 0; % 近地点幅角
theta1 = 0; % 真近点角

% 目标星轨道六根数
a2 = 8000; % 半长轴
e2 = 0; % 离心率
i2 = 10*pi/180; % 轨道倾角
Om2 = 0; % 升交点赤经
om2 = 0; % 近地点幅角
theta2 = 0; % 真近点角

% 计算服务星和目标星的周期
T1 = 2*pi*sqrt(a1^3/mu);
T2 = 2*pi*sqrt(a2^3/mu);

% 计算服务星发射后的时间
t1 = 0;
while (t1 < T1/2)
    r1 = a1*(1 - e1^2)/(1 + e1*cos(theta1)); % 计算服务星在真近点时的距离
    v1 = sqrt(mu*(2/r1 - 1/a1)); % 计算服务星在真近点时的速度
    dt1 = lambert(r1, r1, T2/2 - t1, mu, 'pro', 0); % 计算服务星到目标星的转移时间
    [r1, v1] = rv_lambert(r1, r1, dt1, mu, 'pro', 0); % 计算服务星在转移后的状态向量
    [a1, e1, i1, Om1, om1, theta1] = rv2oe(r1, v1, mu); % 计算服务星在转移后的轨道六根数
    t1 = t1 + dt1; % 更新时间
end

% 计算服务星与目标星相遇时的时间
t2 = t1;
while (t2 < T1)
    r1 = a1*(1 - e1^2)/(1 + e1*cos(theta1)); % 计算服务星在真近点时的距离
    v1 = sqrt(mu*(2/r1 - 1/a1)); % 计算服务星在真近点时的速度
    r2 = a2*(1 - e2^2)/(1 + e2*cos(theta2)); % 计算目标星在真近点时的距离
    v2 = sqrt(mu*(2/r2 - 1/a2)); % 计算目标星在真近点时的速度
    dt2 = lambert(r1, r2, T2/2 - t2, mu, 'pro', 0); % 计算服务星到目标星的转移时间
    [r1, v1] = rv_lambert(r1, r2, dt2, mu, 'pro', 0); % 计算服务星在转移后的状态向量
    [a1, e1, i1, Om1, om1, theta1] = rv2oe(r1, v1, mu); % 计算服务星在转移后的轨道六根数
    t2 = t2 + dt2; % 更新时间
end

% 计算服务星与目标星相遇后的时间
t3 = t2;
while (t3 < T1 + T2/2)
    r2 = a2*(1 - e2^2)/(1 + e2*cos(theta2)); % 计算目标星在真近点时的距离
    v2 = sqrt(mu*(2/r2 - 1/a2)); % 计算目标星在真近点时的速度
    dt3 = lambert(r2, r2, T1 + T2/2 - t3, mu, 'pro', 0); % 计算目标星到服务星的转移时间
    [r2, v2] = rv_lambert(r2, r2, dt3, mu, 'pro', 0); % 计算目标星在转移后的状态向量
    [a2, e2, i2, Om2, om2, theta2] = rv2oe(r2, v2, mu); % 计算目标星在转移后的轨道六根数
    t3 = t3 + dt3; % 更新时间
end

% 计算服务星和目标星的位置与速度
[r1, v1] = oe2rv(a1, e1, i1, Om1, om1, theta1, mu); % 计算服务星在相遇时的状态向量
[r2, v2] = oe2rv(a2, e2, i2, Om2, om2, theta2, mu); % 计算目标星在相遇时的状态向量

% 绘制轨迹图
figure;
hold on;
grid on;
axis equal;
xlabel('X (km)');
ylabel('Y (km)');
zlabel('Z (km)');
title('Orbit Transfer Trajectory');
view(3);
plot_orbit(mu, r1, v1, t2 - t1, 'r'); % 绘制服务星的运动轨迹
plot_orbit(mu, r2, v2, T2/2, 'b'); % 绘制目标星的运动轨迹
legend('Service Satellite', 'Target Satellite');

% 计算服务星和目标星的相对距离
t = linspace(t2 - t1, T2/2, 1000);
for i = 1:length(t)
    [r1, v1] = oe2rv(a1, e1, i1, Om1, om1, theta1 + 2*pi*t(i)/(T1/2), mu); % 计算服务星在当前时间的状态向量
    [r2, v2] = oe2rv(a2, e2, i2, Om2, om2, theta2 + 2*pi*t(i)/T2, mu); % 计算目标星在当前时间的状态向量
    dr = r2 - r1; % 计算两星之间的距离向量
    distance(i) = norm(dr); % 计算两星之间的距离
end

% 绘制距离随时间变化的图像
figure;
hold on;
grid on;
xlabel('Time (s)');
ylabel('Distance (km)');
title('Distance between Service Satellite and Target Satellite');
plot(t, distance);