请用matlab代码示例一个简单的异面轨道转移问题的lambert转移，并画出轨迹图像

异面轨道转移问题是指将一个物体从一个轨道转移到另一个不同的轨道，这两个轨道可能不在同一平面内。Lambert问题是解决这类问题的一种方法，它通过计算两个位置矢量和时间差来确定转移轨道。

以下是一个简单的Matlab代码示例，用于计算在两个位置矢量r1和r2之间的Lambert转移轨道，并绘制出轨迹图像。

% 定义初始位置矢量和结束位置矢量
r1 = [10000, 0, 0];
r2 = [0, 15000, 0];

% 定义初始时间和结束时间
t1 = 0;
t2 = 3600;

% 定义中间时间差
dt = t2 - t1;

% 计算初始位置矢量和结束位置矢量之间的距离和方向
d = norm(r2 - r1);
c = cross(r1, r2);
if c(3) >= 0
    theta = acos(dot(r1, r2)/(norm(r1)*norm(r2)));
else
    theta = 2*pi - acos(dot(r1, r2)/(norm(r1)*norm(r2)));
end

% 使用lambert方法计算转移轨道的速度矢量
[v1, v2] = lambert(r1, r2, dt, 0, 0, 0, 0);

% 计算轨道
state1 = [r1, v1];
state2 = [r2, v2];
[~, orbit1] = ode45(@rates, [t1 t2], state1);
[~, orbit2] = ode45(@rates, [t1 t2], state2);

% 绘制轨迹图像
figure;
plot3(orbit1(:,1), orbit1(:,2), orbit1(:,3), 'b');
hold on;
plot3(orbit2(:,1), orbit2(:,2), orbit2(:,3), 'r');
grid on;
axis equal;
legend('初始轨道', '结束轨道');
title('Lambert转移轨道');
xlabel('x');
ylabel('y');
zlabel('z');

% 定义lambert方法
function [v1, v2] = lambert(r1, r2, dt, mu, dir, tol, maxiter)

if nargin < 7
    maxiter = 100;
end

% 计算初始位置矢量和结束位置矢量之间的距离和方向
d = norm(r2 - r1);
c = cross(r1, r2);
if c(3) >= 0
    theta = acos(dot(r1, r2)/(norm(r1)*norm(r2)));
else
    theta = 2*pi - acos(dot(r1, r2)/(norm(r1)*norm(r2)));
end

% 计算半焦距和delta
a = sin(theta)/sqrt((1 + cos(theta)));
if dir == 0
    delta = 0;
else
    delta = acos(1/d + 1/a);
end

% 迭代求解
for i = 1:maxiter
    psi = a*(delta - sin(delta));
    if psi > 0 && dir == 1
        psi = 2*pi - psi;
    elseif psi < 0 && dir == 0
        psi = 2*pi + psi;
    end
    y = sqrt(mu*a)*dt;
    f = 1 - psi^2/a*d + (psi*y)/d;
    g = dt - (1/sqrt(mu*a^3))*psi^3 + (psi^2*y)/d;
    fdot = sqrt(mu/(a*y))*(psi*d - 1);
    gdot = 1 - psi^2/(2*a)*d - (psi*y)/d;
    v1 = (r2 - f*r1)/g;
    v2 = fdot*r1 + gdot*v1;
    error = norm(r2 - lambert_propagate(r1, v1, mu, dt)') + norm(v2 - lambert_propagate_v(r1, v1, mu, dt)');
    if error < tol
        break;
    end
    delta = delta + error*10^(-5);
end

end

% 定义轨道微分方程
function dydt = rates(~, state, mu)

r = state(1:3);
v = state(4:6);
r_mag = norm(r);
dydt = [v; -mu*r/r_mag^3];
end

% 定义lambert方法中的位置矢量更新
function r = lambert_propagate(r1, v1, mu, dt)

r_mag = norm(r1);
v_mag = norm(v1);
psi = acos(dot(r1, v1)/(r_mag*v_mag));
if dot(r1, v1) < 0
    psi = 2*pi - psi;
end
a = (r_mag*v_mag^2/mu)/(1 - r_mag*v_mag^2/mu*cos(psi));
if a < 0
    a = -a;
end
f = 1 - r_mag/a*(1 - cos(psi));
g = dt + sqrt(a^3/mu)*(psi - sin(psi));
r = f*r1 + g*v1;

end

% 定义lambert方法中的速度矢量更新
function v = lambert_propagate_v(r1, v1, mu, dt)

r_mag = norm(r1);
v_mag = norm(v1);
psi = acos(dot(r1, v1)/(r_mag*v_mag));
if dot(r1, v1) < 0
    psi = 2*pi - psi;
end
a = (r_mag*v_mag^2/mu)/(1 - r_mag*v_mag^2/mu*cos(psi));
if a < 0
    a = -a;
end
fdot = sqrt(mu*a)/(r_mag*v_mag)*sin(psi);
gdot = 1 - r_mag/a*(1 - cos(psi));
v = fdot*r1 + gdot*v1;

end

运行上述代码，将会得到一个包含两条轨道的3D图形，其中蓝色轨道表示初始轨道，红色轨道表示结束轨道，它们之间的Lambert转移轨道是通过迭代计算得出的。