ECEF(X,Y,Z) to Longitude and Latitude.

clc;

% !!! This code following WGS84 Model !!!;
% Define X, Y, Z samples that have Latitude 50; Longitude 100; Height 200M;
% http://www.apsalin.com/convert-geodetic-to-cartesian.aspx
X = -713345.437320888;
Y = +4045583.0097852;
Z = +4862942.24652593;

% Estimate radius of curvature;
a = 6378137; % Meter; % Semi-major axis;
f = 1/298.257223563; % Ellipsoid flattening;
b = a*(1-f); % Define Semi-minor axis;

% Estimate auxiliary values;
P = sqrt(X^2 + Y^2);
Theta = atan(Z*a/P*b);
e = sqrt(((a^2) - (b^2))/a^2); % First eccentricity of The Earth;
e2 = sqrt((a^2 - b^2)/b^2); % Second eccentricity of The Earth;

% Initial value of Latitude;
ppi = atan2(Z,(P*(1-e^2)));

% Iteration loop for estimate Latitude(ppi);
for i = 1:10000
N = a / sqrt(1-e^2*sin(ppi).^2); % Prime vertical;
h = (P/cos(ppi)) - N; % Height;
ppi = atan2(Z,(P*(1-e^2*(N/(N+h))))); % Latitude
end

lamda = (atan2(Y,X))*180/pi;
ppi = ppi*180/pi;

% Display results
disp('Longitutde is :');
disp(lamda);
disp('Latitude is :');
disp(ppi);
disp('Height(Meter) is :');
disp(h);