东北天坐标到大地坐标转化matlab代码

以下是将东北天坐标转换为大地坐标的 MATLAB 代码：

function [lat,lon,h] = ned2lla(n,e,d,lat0,lon0,h0)
% NED to LLA coordinate conversion
% Input:
% n,e,d - North, East, Down (meters)
% lat0, lon0, h0 - Reference point in degrees and meters
% Output:
% lat,lon,h - Latitude, Longitude (degrees), Height (meters)

% WGS84 ellipsoid parameters
a = 6378137;
f = 1/298.257223563;
b = a*(1-f);

% Convert reference point to ECEF
[rx,ry,rz] = lla2ecef(lat0,lon0,h0);

% Rotation matrix from NED to ECEF
R = [-sin(lat0)*cos(lon0) -sin(lat0)*sin(lon0) cos(lat0);
     -sin(lon0) cos(lon0) 0;
     -cos(lat0)*cos(lon0) -cos(lat0)*sin(lon0) -sin(lat0)];

% Convert NED to ECEF
p = [n;e;-d];
p_ecef = R*p + [rx;ry;rz];

% Convert ECEF to LLA
[lat,lon,h] = ecef2lla(p_ecef(1),p_ecef(2),p_ecef(3));

end

function [x,y,z] = lla2ecef(lat,lon,h)
% Convert from geodetic coordinates to ECEF coordinates
% Input:
% lat,lon,h - Latitude and longitude in degrees, altitude in meters
% Output:
% x,y,z - ECEF coordinates in meters

% WGS84 ellipsoid parameters
a = 6378137;
f = 1/298.257223563;
b = a*(1-f);
e2 = 1 - (b/a)^2;

% Convert to radians
lat = lat*pi/180;
lon = lon*pi/180;

% Radius of curvature in prime vertical
N = a./sqrt(1-e2*sin(lat).^2);

% ECEF coordinates
x = (N+h).*cos(lat).*cos(lon);
y = (N+h).*cos(lat).*sin(lon);
z = (N.*(1-e2)+h).*sin(lat);

end

function [lat,lon,h] = ecef2lla(x,y,z)
% Convert from ECEF coordinates to geodetic coordinates
% Input:
% x,y,z - ECEF coordinates in meters
% Output:
% lat,lon,h - Latitude and longitude in degrees, altitude in meters

% WGS84 ellipsoid parameters
a = 6378137;
f = 1/298.257223563;
b = a*(1-f);
e2 = 1 - (b/a)^2;

% Distance from z axis
r = sqrt(x.^2 + y.^2);

% Longitude
lon = atan2(y,x);

% Bowring's formula for initial estimate of latitude
p = sqrt(x.^2 + y.^2)/a;
lat = atan(z./(p.*(1-e2)));

% Iterate until convergence
for i=1:10
    N = a./sqrt(1-e2*sin(lat).^2);
    h = r./cos(lat) - N;
    lat = atan(z./(p.*(1-e2*N./(N+h))));
end

% Convert to degrees
lat = lat*180/pi;
lon = lon*180/pi;

end

其中，`n`，`e` 和 `d` 分别是东北天坐标系下的三个坐标值，`lat0`，`lon0` 和 `h0` 是参考点的大地坐标系下的纬度、经度和高度。函数的输出是大地坐标系下的纬度、经度和高度。