xyz转blh matlab函数

以下是一个将 XYZ 坐标转换为 BLH 坐标的 MATLAB 函数：

function [lat, lon, alt] = xyz2blh(x, y, z, a, b)
% XYZ2BLH Convert XYZ coordinates to geodetic coordinates.
%   [LAT, LON, ALT] = XYZ2BLH(X, Y, Z, A, B) converts the input Cartesian
%   coordinates X, Y, and Z (in meters) to geodetic coordinates LAT and LON
%   (in degrees) and ALT (in meters) using the specified reference ellipsoid
%   defined by semi-major axis A and semi-minor axis B (in meters).

% Convert Cartesian coordinates to spherical coordinates
[theta, phi, r] = cart2sph(x, y, z);

% Compute ellipsoid parameters
e2 = (a^2 - b^2) / a^2;
ep2 = (a^2 - b^2) / b^2;

% Compute geodetic latitude
lat = atan2((z / r),((1 - e2) * sqrt(x^2 + y^2) / r));
old_lat = 0;
while (abs(lat - old_lat) > 1e-10)
    old_lat = lat;
    N = a / sqrt(1 - e2 * sin(old_lat)^2);
    lat = atan2((z / r + ep2 * N * sin(old_lat)), sqrt(x^2 + y^2));
end

% Compute geodetic longitude
lon = atan2(y, x);

% Compute geodetic altitude
N = a / sqrt(1 - e2 * sin(lat)^2);
alt = sqrt(x^2 + y^2) / cos(lat) - N;

% Convert latitude and longitude to degrees
lat = rad2deg(lat);
lon = rad2deg(lon);
end

要使用此函数，请提供 XYZ 坐标和 reference ellipsoid 的参数 a 和 b。例如，要将 XYZ 坐标 `[100, 50, 200]` 转换为 BLH 坐标，使用 WGS84 ellipsoid 的参数，可以执行以下操作：

x = 100;
y = 50;
z = 200;
a = 6378137;    % Semi-major axis of WGS84 ellipsoid
b = 6356752.3142;    % Semi-minor axis of WGS84 ellipsoid
[lat, lon, alt] = xyz2blh(x, y, z, a, b);