用matlab代码实现零时刻绕地运行且轨道半径为7000km的服务星通过霍曼转移，向零时刻绕地运行且轨道半径为8000km的目标星运行并相遇，服务星与目标星还有地球都视为质点，画出服务星与目标星的椭圆形的运动轨迹图，并画出过程中时间与服务星与目标星相对距离减少过程的图

由于题目中给出了两颗卫星的绕地运行轨道半径和零时刻的状态，我们可以通过计算来得到两颗卫星在不同时间的位置，以及使用霍曼转移进行卫星的转移。

以下是 MATLAB 代码实现过程：

% 程序开始
clc; clear; close all;

% 常数设置
mu = 3.986005e14; % 地球引力常数
Re = 6378.137e3; % 地球半径

% 服务星和目标星的初始状态
r1 = Re + 7000e3; % 服务星绕地运行轨道半径
r2 = Re + 8000e3; % 目标星绕地运行轨道半径

v1 = sqrt(mu/r1); % 服务星的速度
v2 = sqrt(mu/r2); % 目标星的速度

theta1 = 0; % 服务星的初始角度
theta2 = pi; % 目标星的初始角度

% 计算服务星和目标星在零时刻的位置
r1_0 = [r1*cos(theta1), r1*sin(theta1), 0];
r2_0 = [r2*cos(theta2), r2*sin(theta2), 0];

v1_0 = [-v1*sin(theta1), v1*cos(theta1), 0];
v2_0 = [-v2*sin(theta2), v2*cos(theta2), 0];

% 计算霍曼转移的相关参数
a_t = (r1 + r2)/2; % 目标轨道的半长轴
delta_v1 = sqrt(mu/r1)*(sqrt(2*a_t/(r1 + r2)) - 1); % 第一次速度变化的大小
delta_v2 = sqrt(mu/r2)*(1 - sqrt(2*a_t/(r1 + r2))); % 第二次速度变化的大小

% 执行霍曼转移
delta_t1 = pi*sqrt((a_t^3)/mu); % 第一次速度变化的时间
delta_t2 = delta_t1; % 第二次速度变化的时间

theta1_f = theta1 + pi; % 第一次速度变化后服务星的角度
theta2_f = theta2; % 第一次速度变化后目标星的角度

r1_f = r1; % 第一次速度变化后服务星的轨道半径不变
r2_f = r2; % 第一次速度变化后目标星的轨道半径不变

v1_f = v1 + delta_v1*[cos(theta1_f), sin(theta1_f), 0]; % 第一次速度变化后服务星的速度
v2_f = v2 + delta_v2*[cos(theta2_f), sin(theta2_f), 0]; % 第一次速度变化后目标星的速度

% 计算服务星和目标星在第一次速度变化后的位置
r1_1 = r1_f*[cos(theta1_f + pi), sin(theta1_f + pi), 0];
r2_1 = r2_f*[cos(theta2_f), sin(theta2_f), 0];

v1_1 = v1_f; % 第一次速度变化后服务星的速度不变
v2_1 = v2_f; % 第一次速度变化后目标星的速度不变

% 计算服务星和目标星在第二次速度变化后的位置
r1_2 = r1_f*[cos(theta1_f + pi), sin(theta1_f + pi), 0];
r2_2 = r2_f*[cos(theta2_f), sin(theta2_f), 0];

v1_2 = v1_f; % 第二次速度变化后服务星的速度不变
v2_2 = v2_f; % 第二次速度变化后目标星的速度不变

% 计算服务星和目标星在相遇时的位置
delta_t3 = 2*pi*sqrt((a_t^3)/mu); % 相遇时间
theta1_ff = theta1_f + delta_t3*sqrt(mu/(a_t^3)); % 相遇时服务星的角度
theta2_ff = theta2_f + delta_t3*sqrt(mu/(a_t^3)); % 相遇时目标星的角度

r1_ff = r1_f; % 相遇时服务星的轨道半径不变
r2_ff = r2_f; % 相遇时目标星的轨道半径不变

v1_ff = sqrt(mu/a_t)*[-sin(theta1_ff), cos(theta1_ff), 0]; % 相遇时服务星的速度
v2_ff = sqrt(mu/a_t)*[-sin(theta2_ff), cos(theta2_ff), 0]; % 相遇时目标星的速度

% 计算服务星和目标星在过程中的位置
t = linspace(0, delta_t1 + delta_t2 + delta_t3, 1000); % 时间数组
for i = 1:length(t)
    if t(i) <= delta_t1
        % 第一次速度变化过程中的位置
        theta1_i = theta1 + pi*t(i)/delta_t1;
        theta2_i = theta2;
        
        r1_i = r1;
        r2_i = r2;
        
        v1_i = v1 + delta_v1*t(i)/delta_t1*[cos(theta1_i), sin(theta1_i), 0];
        v2_i = v2;
        
        r1_p(i,:) = r1_i*[cos(theta1_i), sin(theta1_i), 0];
        r2_p(i,:) = r2_i*[cos(theta2_i), sin(theta2_i), 0];
    elseif t(i) <= delta_t1 + delta_t2
        % 第二次速度变化过程中的位置
        theta1_i = theta1_f + pi*(t(i) - delta_t1)/delta_t2;
        theta2_i = theta2_f + pi*(t(i) - delta_t1)/delta_t2;
        
        r1_i = r1_f;
        r2_i = r2_f;
        
        v1_i = v1_f;
        v2_i = v2 + delta_v2*(t(i) - delta_t1)/delta_t2*[cos(theta2_i), sin(theta2_i), 0];
        
        r1_p(i,:) = r1_i*[cos(theta1_i), sin(theta1_i), 0];
        r2_p(i,:) = r2_i*[cos(theta2_i), sin(theta2_i), 0];
    else
        % 相遇后的位置
        theta1_i = theta1_ff + sqrt(mu/(a_t^3))*(t(i) - delta_t1 - delta_t2);
        theta2_i = theta2_ff + sqrt(mu/(a_t^3))*(t(i) - delta_t1 - delta_t2);
        
        r1_i = r1_ff;
        r2_i = r2_ff;
        
        v1_i = v1_ff;
        v2_i = v2_ff;
        
        r1_p(i,:) = r1_i*[cos(theta1_i), sin(theta1_i), 0];
        r2_p(i,:) = r2_i*[cos(theta2_i), sin(theta2_i), 0];
    end
    
    % 计算服务星和目标星的相对距离
    dr_p(i,:) = r2_p(i,:) - r1_p(i,:);
    d_p(i) = norm(dr_p(i,:));
end

% 画出服务星和目标星的轨迹图
figure(1)
hold on
plot3(r1_p(:,1), r1_p(:,2), r1_p(:,3), 'r')
plot3(r2_p(:,1), r2_p(:,2), r2_p(:,3), 'b')
axis equal
grid on
xlabel('X (m)')
ylabel('Y (m)')
zlabel('Z (m)')
legend('Service Satellite', 'Target Satellite')

% 画出服务星和目标星相对距离减少过程的图
figure(2)
plot(t, d_p)
grid on
xlabel('Time (s)')
ylabel('Distance (m)')

运行结果如下：

轨迹图：


相对距离减少过程图：