Function (fi0,fi1,fi,lambda,N,h,p) = XYZnaGeo(x,y,z,e2,a)
p=sqrt(x*x+y*y);
fi0=atan((z/p)*((1-e2)^(-1)));
do
N=a/(sqrt(1-e2*((sin(fi0))^2)));
h=p/(cos(fi0))-N;
fi1=atan((z/p)/((1-e2*(N/(N+h)))^(-1)));
until (abs(fi1-fi0))<((0.0001*pi)/(180*3600));
fi=fi1;
lambda=atan(y/x);
N=a/(sqrt(1-e2*((sin(fi))^2)));
h=(p/cos(fi))-N;
endfunction;