    %% Bowring iterative method to determine User Geodetic Coordinates %%
%Input =  user coordinates, semimajor (a) and semiminor (b) aixs lenghts of
%Ellipsoid
%Output = geodetic latitude (phi), longitude (lambda) & height (h)
% [rad] [rad] [meters]
function [phi,lambda,h]=Bowring(xu,yu,zu,a,b)
e=sqrt(1-(b/a)^2);  %Ellipsoid 
e2=(a/b)*e;         %second eccentricity

p=sqrt(xu^2+yu^2);
u1=atan((zu/p)*(a/b)); %parametric latitude (first approximation)
if isnan(u1)
    u1=0;
end
phi=atan((zu+(e2^2)*b*sin(u1)^3)/(p-(e^2)*a*cos(u1)^3));
u=(b/a)*tan(phi);
while u-u1>10e-2
    u1=u;
    %Geodetic longitude
    phi=atan((zu+(e2^2)*b*sin(u1)^3)/(p-(e^2)*a*cos(u1)^3));
    u=(b/a)*tan(phi);
end
N=a/(sqrt(1-(e^2)*sin(phi)^2));
%Geodetic height
if phi~=pi/2 && phi~=-pi/2
    h=abs(p/cos(phi)-N);
elseif phi~=0
    h=abs(zu/sin(phi)+N*(e^2-1));
end
%Geodetic longitude
if xu>=0
    lambda=atan(yu/xu);
elseif xu<0 && yu>=0
    lambda=pi+atan(yu/xu);
elseif xu<0 && yu<0
    lambda=-pi+atan(yu/xu);
end
if isnan(lambda)
    lambda=0;
end
end