xyz2ell.m

function [lat, lon, h] = xyz2ell(X, Y, Z, a, e2)
% XYZ2ELL  Converts cartesian coordinates to ellipsoidal coodinates.

% We can add some default values here, i.e. a and e2.

% Latitude and height convergence criteria
eps_lat = 1.e-12;
eps_height = 1.e-5;

% Initial values...

p = sqrt(X .* X + Y .* Y);
lat = atan2(Z, p .* (1 - e2));
h = 0;
dh = 1;
dlat = 1;

% Iterate until convergence, the difference between lat(i) and lat(i - 1) <= eps, and so for the height.

while sum(dlat > eps_lat) | sum(dh > eps_height)
  lat0 = lat;
  h0 = h;
  v = a ./ sqrt(1 - e2 .* sin(lat) .* sin(lat));
  h = p ./ cos(lat) - v;
  lat = atan2(Z, p .* (1 - e2 .* v ./ (v + h)));
  dlat = abs(lat-lat0);
  dh = abs(h - h0);
end
lon = atan2(Y, X);

% For reference for atan2() uses refer to cens.ioc.ee/local/matlab/techdoc/ref/atan2.html