function [phi,lamdac,h]=e2g(x,y,z)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% filename e2g.m
% E2G   Earth-centered cartesian coordinates to geodetic coordinates conversion.
% [phi,lamdac,h]=e2g(x,y,z) gives 
% phi : geodetic latitude (rad),
% lamdac : celestial longitude (rad) and 
% h : altitude normal to ellipsoid (ft),
% corresponding to
% x, y and z : ECI cartesian position coordinates (ft).
% User can change the unit of length replacing the value of "a", such as
% a=6378.137;         %(km)
% a=6378137;          %(m)
% a=20925646.3254593; %(ft)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a=20925646.3254593; %(ft)
e2=0.00669437999;
b=a*sqrt(1-e2);
w=sqrt(x^2+y^2);
if w ~= 0
l=e2/2;
m=(w/a)^2;
n=((1-e2)*z/b)^2;
c=-(2*l^2+m+n)/2;
k=l^2*(l^2-m-n);
q=(m+n-4*l^2)^3/216+m*n*l^2;
D=sqrt((2*q-m*n*l^2)*m*n*l^2);
B=c/3-(q+D)^(1/3)-(q-D)^(1/3);
t=sqrt(sqrt(B^2-k)-(B+c)/2)-sign(m-n)*sqrt((B-c)/2);
w1=w/(t+l);
z1=(1-e2)*z/(t-l);
phi=atan2(z1,(1-e2)*w1);
lamdac=2*atan2(w-x,y);
h=sign(t-1+l)*sqrt((w-w1)^2+(z-z1)^2);
else
lamdac=0;
h=sign(z)*z-b;
phi=sign(z)*pi/2;
end
%end of file
