请用MATLAB语言写一个将位置速度矢量转换为轨道六根数的程序

以下是一个简单的MATLAB程序，它将位置速度矢量转换为轨道六根数：

% 输入位置速度矢量
r = [7000 0 0]; % km
v = [0 7.5 0]; % km/s

% 计算轨道参数
mu = 398600; % km^3/s^2，地球的标准引力常数
r_mag = norm(r);
v_mag = norm(v);
h_vec = cross(r,v);
h_mag = norm(h_vec);
inc = acosd(h_vec(3)/h_mag);
if inc == 0 || inc == 180
    raan = 0;
else
    raan = acosd(h_vec(1)/(h_mag*sind(inc)));
    if h_vec(2) < 0
        raan = 360 - raan;
    end
end
ecc_vec = ((v_mag^2 - mu/r_mag)*r - dot(r,v)*v)/mu;
ecc = norm(ecc_vec);
if ecc == 0
    argp = 0;
    nu = 0;
else
    if inc == 0 || inc == 180
        argp = 0;
    else
        argp = acosd(dot(h_vec,ecc_vec)/(h_mag*ecc));
        if ecc_vec(3) < 0
            argp = 360 - argp;
        end
    end
    nu = acosd(dot(ecc_vec,r)/(ecc*r_mag));
    if dot(r,v) < 0
        nu = 360 - nu;
    end
end
semi_major_axis = 1/(2/r_mag - v_mag^2/mu);

% 输出轨道参数
fprintf('semi-major axis: %.2f km\n', semi_major_axis);
fprintf('eccentricity: %.4f\n', ecc);
fprintf('inclination: %.2f deg\n', inc);
fprintf('right ascension of ascending node: %.2f deg\n', raan);
fprintf('argument of periapsis: %.2f deg\n', argp);
fprintf('true anomaly: %.2f deg\n', nu);

该程序的输入为位置矢量 `r` 和速度矢量 `v`，以公里和公里/秒为单位。输出为轨道的半长轴 `semi_major_axis`（以公里为单位）、离心率 `ecc`、倾角 `inc`、升交点赤经 `raan`、近地点幅角 `argp` 和真近点角 `nu`，以度为单位。