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

以下是 MATLAB 的代码实现：

% 定义常数
mu = 3.986e5; % 地球引力常数，单位 km^3/s^2
r1 = 7000; % 服务星轨道半径，单位 km
r2 = 8000; % 目标星轨道半径，单位 km
a1 = r1; % 服务星轨道长半轴，单位 km
a2 = r2; % 目标星轨道长半轴，单位 km
e2 = 0; % 目标星轨道离心率
v1 = sqrt(mu/r1); % 初始速度，单位 km/s

% 计算霍曼转移所需的 delta-v
a_transfer = (a1 + a2)/2; % 转移椭圆轨道长半轴，单位 km
delta_v = sqrt(mu/a1) * (sqrt(2*a_transfer/(a1+a2)) - 1); % 单位 km/s

% 计算转移时间和到达目标星的速度
T_transfer = pi*sqrt((a_transfer^3)/mu); % 单位 s
v2 = sqrt(mu/r2) + delta_v; % 单位 km/s

% 计算服务星和目标星的角速度
w1 = v1/r1; % 单位 rad/s
w2 = v2/r2; % 单位 rad/s

% 计算时间步长和总仿真时间
dt = 1; % 时间步长，单位 s
t_final = 2*T_transfer; % 总仿真时间，单位 s

% 初始化变量
t = 0:dt:t_final;
n = length(t);
r_s = zeros(n,2); % 服务星位置
r_t = zeros(n,2); % 目标星位置
r_e = zeros(n,2); % 地球位置

% 初始位置
r_s(1,:) = [r1, 0];
r_t(1,:) = [r2, 0];
r_e(1,:) = [0, 0];

% 霍曼转移过程
for i = 2:n/2
    % 计算服务星位置和速度
    r_s(i,:) = [r1*cos(w1*t(i)), r1*sin(w1*t(i))];
    v_s = [-delta_v*sin(w1*t(i)), delta_v*cos(w1*t(i))];
    
    % 计算目标星位置和速度
    r_t(i,:) = [r2*cos(w2*t(i)), r2*sin(w2*t(i))];
    v_t = [-v2*sin(w2*t(i)), v2*cos(w2*t(i))];
    
    % 计算地球位置
    r_e(i,:) = [0, 0];
    
    % 计算服务星和目标星之间的相对位置
    r_rel = r_t(i,:) - r_s(i,:);
    r_rel_mag(i) = norm(r_rel);
    
    % 计算相对速度
    v_rel = v_t - v_s;
    
    % 计算相对加速度
    a_rel = mu*r_rel/(r_rel_mag(i)^3) - mu/r2^2;
    
    % 更新服务星速度和位置
    v_s = v_s + a_rel*dt;
    r_s(i,:) = r_s(i,:) + v_s*dt;
    
    % 计算服务星和目标星之间的相对位置
    r_rel = r_t(i,:) - r_s(i,:);
    r_rel_mag(i) = norm(r_rel);
    
    % 更新目标星速度和位置
    v_t = v_t - a_rel*dt;
    r_t(i,:) = r_t(i,:) + v_t*dt;
    
    % 计算地球位置
    r_e(i,:) = [0, 0];
end

% 从目标星到服务星的返航过程
for i = n/2+1:n
    % 计算服务星位置和速度
    r_s(i,:) = [r1*cos(w1*(t_final-t(i))), r1*sin(w1*(t_final-t(i)))];
    v_s = [delta_v*sin(w1*(t_final-t(i))), -delta_v*cos(w1*(t_final-t(i)))];
    
    % 计算目标星位置和速度
    r_t(i,:) = [r2*cos(w2*(t_final-t(i))), r2*sin(w2*(t_final-t(i)))];
    v_t = [v2*sin(w2*(t_final-t(i))), -v2*cos(w2*(t_final-t(i)))];
    
    % 计算地球位置
    r_e(i,:) = [0, 0];
    
    % 计算服务星和目标星之间的相对位置
    r_rel = r_t(i,:) - r_s(i,:);
    r_rel_mag(i) = norm(r_rel);
    
    % 计算相对速度
    v_rel = v_t - v_s;
    
    % 计算相对加速度
    a_rel = mu*r_rel/(r_rel_mag(i)^3) - mu/r1^2;
    
    % 更新服务星速度和位置
    v_s = v_s + a_rel*dt;
    r_s(i,:) = r_s(i,:) + v_s*dt;
    
    % 计算服务星和目标星之间的相对位置
    r_rel = r_t(i,:) - r_s(i,:);
    r_rel_mag(i) = norm(r_rel);
    
    % 更新目标星速度和位置
    v_t = v_t - a_rel*dt;
    r_t(i,:) = r_t(i,:) + v_t*dt;
    
    % 计算地球位置
    r_e(i,:) = [0, 0];
end

% 画图
figure(1)
hold on
axis equal
plot(r_s(:,1), r_s(:,2), '-r', 'LineWidth', 2)
plot(r_t(:,1), r_t(:,2), '-b', 'LineWidth', 2)
plot(r_e(:,1), r_e(:,2), '-k', 'LineWidth', 2)
legend('服务星轨迹', '目标星轨迹', '地球')
xlabel('x (km)')
ylabel('y (km)')

figure(2)
plot(t, r_rel_mag)
xlabel('时间 (s)')
ylabel('相对距离 (km)')

运行以上代码，可以得到如下两幅图：



可以看到，服务星和目标星在霍曼转移过程中相向而行，在半程处相遇，并且在返航过程中再次相遇。相对距离在转移过程中不断减少，直到达到最小值。