给定地心地固坐标系下三个速度分量，得到东北天坐标系下的三个速度分量的matlab代码

以下是给定地心地固坐标系下三个速度分量，得到东北天坐标系下的三个速度分量的 Matlab 代码：

function [vn, ve, vd] = ecef2ned_vel(x, v, lat, lon, h)

    % WGS-84 ellipsoid constants
    a = 6378137; % Semi-major axis
    b = 6356752.3142; % Semi-minor axis
    e = sqrt(1-(b/a)^2); % Eccentricity
    ep = sqrt((a/b)^2-1); % Second Eccentricity

    % Compute the necessary trigonometric functions
    sinLat = sin(lat);
    cosLat = cos(lat);
    sinLon = sin(lon);
    cosLon = cos(lon);

    % Compute the local radius of curvature in the prime vertical
    N = a^2./sqrt(a^2*cosLat.^2+b^2*sinLat.^2);

    % Compute the local radius of curvature in the meridian
    M = a*(1-e^2)./((1-(e^2)*sinLat.^2).^1.5);

    % Compute the rotation matrix from the ECEF to the local tangent plane
    R = [-sinLat*cosLon          -sinLat*sinLon           cosLat     ;
         -sinLon                 cosLon                  0          ;
         -cosLat*cosLon          -cosLat*sinLon          -sinLat     ];

    % Convert the ECEF velocity vector to the local tangent plane
    Ve = v(1);
    Vn = v(2);
    Vd = v(3);
    vn_dot = R(1,1)*Ve + R(1,2)*Vn + R(1,3)*Vd;
    ve_dot = R(2,1)*Ve + R(2,2)*Vn + R(2,3)*Vd;
    vd_dot = R(3,1)*Ve + R(3,2)*Vn + R(3,3)*Vd;

    % Compute the velocity in the NED frame
    vn = vn_dot;
    ve = ve_dot;
    vd = -vd_dot;

    % Return the NED velocity vector
    vn = vn/N;
    ve = ve/(M*cosLat);
    vd = vd;
end

其中，输入参数 `x` 表示 ECEF 坐标系下的位置向量，`v` 表示 ECEF 坐标系下的速度向量，`lat` 表示纬度，`lon` 表示经度，`h` 表示高度。函数的输出是 NED 坐标系下的速度向量，分别为 `vn`、`ve` 和 `vd`。