用matlab代码实现利用拉格朗日系数法求解兰伯特转移所需的两个速度增量

兰伯特转移是一种航天器在两个不同轨道之间进行转移的方法，其中需要计算两个速度增量。下面是用 MATLAB 实现拉格朗日系数法求解兰伯特转移所需的两个速度增量的代码：

% 兰伯特转移求解速度增量
% 输入参数：r1和r2为起始和终止位置向量，dt为转移时间
% 输出参数：v1和v2为起始和终止速度向量
function [v1, v2] = lambert(r1, r2, dt)
    mu = 398600.44; % 地球引力常数
    r = norm(r1); % 起始位置向量的模长
    r_ = norm(r2); % 终止位置向量的模长
    c = cross(r1, r2); % 叉积
    sin_theta = norm(c) / (r * r_); % sin(theta)的值
    cos_theta = dot(r1, r2) / (r * r_); % cos(theta)的值
    if c(3) >= 0
        theta = atan2(sin_theta, cos_theta); % theta的值
    else
        theta = 2*pi - atan2(sin_theta, cos_theta); % theta的值
    end

    A = sin(theta) * sqrt(r * r_ / (1 - cos_theta)); % A的值
    if A == 0
        A = 1e-6;
    end

    z = 0;
    while z == 0
        for i = 1:2
            y = r + r_ + A * ((-1)^i) * sqrt(2 * (r * r_ + r * A + r_ * A) - A^2 * (1 + cos_theta));
            f = (1 - y / r) * cos(sqrt(y^2 - r^2 - 2 * r * y * cos(theta)) / A) ...
                + (y / r_) * cos(sqrt(y^2 - r_^2 - 2 * r_ * y * cos(theta)) / A) - 1;
            fp = sqrt(2) * sin(sqrt(y^2 - r^2 - 2 * r * y * cos(theta)) / A) * ...
                (sqrt(y^2 - r^2 - 2 * r * y * cos(theta)) - A * cos(sqrt(y^2 - r^2 - 2 * r * y * cos(theta)) / A)) / (sqrt(y^2 - r^2 - 2 * r * y * cos(theta)) * A^2) ...
                - (y / r^2) * cos(sqrt(y^2 - r^2 - 2 * r * y * cos(theta)) / A) ...
                + (1 - y / r_) * sqrt((y^2 - r_^2 - 2 * r_ * y * cos(theta))) / (r_ * A) .* ...
                sin(sqrt(y^2 - r_^2 - 2 * r_ * y * cos(theta)) / A);
            if abs(f) < 1e-8
                z = 1;
                break;
            else
                A = A - f / fp;
            end
        end
    end

    psi = sqrt(r * r_ * (1 - cos_theta) / A^2); % psi的值
    chi = (r_ - r) / (A * sqrt(sin_theta)); % chi的值

    f = 1 - psi^2 / r;
    g = dt - sqrt(r^3 / mu) * (psi * sin(sqrt(r * mu / r^3) * dt) - sqrt(r) * chi * (1 - cos(sqrt(r * mu / r^3) * dt)));
    gdot = 1 - psi^2 / r_ * (1 - cos(sqrt(r_ * mu / r_^3) * dt)) / sqrt(r_ / mu);
    v1 = (r2 - f * r1) / g; % 起始速度向量
    v2 = (gdot * r2 - r1) / g; % 终止速度向量
end

在上面的代码中，`r1` 和 `r2` 分别为起始和终止位置向量，`dt` 为转移时间。函数返回值 `v1` 和 `v2` 分别为起始和终止速度向量。