用matlab代码实现霍曼轨道转移

霍曼轨道转移是一种将卫星从一圆形轨道转移到另一个圆形轨道的方法。以下是使用MATLAB实现霍曼轨道转移的示例代码：

% 霍曼轨道转移
% 地球轨道参数
muE = 398600; % 地球引力常数，单位 km^3/s^2
rE = 6378; % 地球半径，单位 km
h1 = 200; % 第一个圆形轨道高度，单位 km
r1 = rE + h1; % 第一个圆形轨道半径，单位 km
v1 = sqrt(muE/r1); % 第一个圆形轨道速度，单位 km/s

% 火星轨道参数
muM = 42828; % 火星引力常数，单位 km^3/s^2
rM = 3397; % 火星半径，单位 km
h2 = 500; % 第二个圆形轨道高度，单位 km
r2 = rM + h2; % 第二个圆形轨道半径，单位 km
v2 = sqrt(muM/r2); % 第二个圆形轨道速度，单位 km/s

% 转移椭圆参数
a = (r1 + r2)/2; % 椭圆半长轴长度，单位 km
e = (r2 - r1)/(r1 + r2); % 椭圆离心率
p = a*(1 - e^2); % 椭圆焦距距离，单位 km
delta_v = sqrt(muE/r1)*(sqrt((2*r2)/(r1+r2))-1); % 转移所需速度增量，单位 km/s

% 计算出发点和到达点的位置向量
r1_vec = [r1 0 0]; % 第一个卫星位置向量，单位 km
r2_vec = [0 r2 0]; % 第二个卫星位置向量，单位 km

% 计算转移轨道参数
theta = acos(dot(r1_vec, r2_vec)/(norm(r1_vec)*norm(r2_vec))); % 转移角度，单位 rad
T = pi*sqrt(a^3/muE); % 周期，单位 s
t = 0:100:1.5*T; % 时间序列，单位 s
M = sqrt(muE/a^3)*t'; % 平近点角，单位 rad
E = zeros(size(M)); % 初值
for i = 1:length(M)
    E(i) = fzero(@(x) x - e*sin(x) - M(i), M(i));
end
phi = 2*atan(sqrt((1+e)/(1-e))*tan(E/2)); % 真近点角，单位 rad
r = p./(1+e*cos(phi)); % 距离，单位 km
h = r - r1; % 高度，单位 km
theta2 = theta - phi; % 第二个圆形轨道上的真近点角，单位 rad

% 计算速度矢量
v_vec1 = [0 v1 0]; % 第一个卫星速度矢量，单位 km/s
v_vec2 = [0 v2 0]; % 第二个卫星速度矢量，单位 km/s
v_vec1_t = [-v1*sin(theta) v1*cos(theta) 0]; % 第一个卫星切向速度矢量，单位 km/s
v_vec1_n = [v1*cos(theta) v1*sin(theta) 0]; % 第一个卫星法向速度矢量，单位 km/s
v_vec2_t = [-v2*sin(theta2) v2*cos(theta2) 0]; % 第二个卫星切向速度矢量，单位 km/s
v_vec2_n = [v2*cos(theta2) v2*sin(theta2) 0]; % 第二个卫星法向速度矢量，单位 km/s
delta_v_vec1 = delta_v*(v_vec2_t/norm(v_vec2_t) - v_vec1_t/norm(v_vec1_t)); % 第一个卫星速度增量矢量，单位 km/s
delta_v_vec2 = delta_v*(v_vec2_n/norm(v_vec2_n) - v_vec1_n/norm(v_vec1_n)); % 第二个卫星速度增量矢量，单位 km/s
v_vec1_new = v_vec1_t + delta_v_vec1 + delta_v_vec2; % 第一个卫星新速度矢量，单位 km/s
v_vec2_new = v_vec2_t + delta_v_vec1 + delta_v_vec2; % 第二个卫星新速度矢量，单位 km/s

% 绘制轨道
figure;
hold on;
grid on;
axis equal;
title('Hohmann Transfer Orbit');
xlabel('x (km)');
ylabel('y (km)');
zlabel('z (km)');
plot3(0, 0, 0, 'ko', 'MarkerSize', 10);
plot3(r1_vec(1), r1_vec(2), r1_vec(3), 'bo', 'MarkerSize', 5);
plot3(r2_vec(1), r2_vec(2), r2_vec(3), 'ro', 'MarkerSize', 5);
plot3(r.*cos(phi), r.*sin(phi), zeros(size(phi)), 'k');
plot3(r1*cos(theta), r1*sin(theta), 0, 'bo', 'MarkerSize', 5);
plot3(r2*cos(theta2), r2*sin(theta2), 0, 'ro', 'MarkerSize', 5);
quiver3(r1_vec(1), r1_vec(2), r1_vec(3), v_vec1_t(1), v_vec1_t(2), v_vec1_t(3), 'b');
quiver3(r2_vec(1), r2_vec(2), r2_vec(3), v_vec2_t(1), v_vec2_t(2), v_vec2_t(3), 'r');
quiver3(r1_vec(1), r1_vec(2), r1_vec(3), delta_v_vec1(1), delta_v_vec1(2), delta_v_vec1(3), 'g');
quiver3(r2_vec(1), r2_vec(2), r2_vec(3), delta_v_vec2(1), delta_v_vec2(2), delta_v_vec2(3), 'g');
quiver3(r1_vec(1), r1_vec(2), r1_vec(3), v_vec1_new(1), v_vec1_new(2), v_vec1_new(3), 'c');
quiver3(r2_vec(1), r2_vec(2), r2_vec(3), v_vec2_new(1), v_vec2_new(2), v_vec2_new(3), 'm');
legend('Sun', 'Earth', 'Mars', 'Transfer Orbit', 'First Circular Orbit', 'Second Circular Orbit', 'First Orbit Velocity', 'Second Orbit Velocity', 'Delta V1', 'Delta V2', 'New First Orbit Velocity', 'New Second Orbit Velocity');
view(-30, 30);

这段代码将计算地球和火星两个卫星的圆形轨道参数，然后计算转移椭圆的参数，最后绘制出轨道图。在图中，蓝色点代表地球，红色点代表火星，黑色轨道代表转移轨道，绿色箭头代表第一个卫星的速度增量，紫色箭头代表第二个卫星的速度增量，蓝色箭头代表第一个卫星的新速度矢量，红色箭头代表第二个卫星的新速度矢量。