function?[dphidlambdah]?=?togeod（afinvXYZ）
%TOGEOD???Subroutine?to?calculate?geodetic?coordinates
%?????????latitude?longitude?height?given?Cartesian
%?????????coordinates?XYZ?and?reference?ellipsoid
%?????????values?semi-major?axis?（a）?and?the?inverse
%?????????of?flattening?（finv）.

%??The?units?of?linear?parameters?XYZa?must?all?agree?（mkmmift..etc）
%??The?output?units?of?angular?quantities?will?be?in?decimal?degrees
%??（15.5?degrees?not?15?deg?30?min）.??The?output?units?of?h?will?be?the
%??same?as?the?units?of?XYZa.

%??Copyright?（C）?1987?C.?Goad?Columbus?Ohio
%??Reprinted?with?permission?of?author?1996
%??Fortran?code?translated?into?MATLAB
%??Kai?Borre?03-30-96
%?Changed?according?to?Matlab?ver.?6.5.1?9?February?2004

h?=?0;
tolsq?=?1.e-10;
maxit?=?10;
%?compute?radians-to-degree?factor
rtd?=?180/pi;
%?compute?square?of?eccentricity
if?finv????esq?=?0;
else??
???esq?=?（2-1/finv）/finv;?
end
oneesq?=?1-esq;
%?first?guess
%?P?is?distance?from?spin?axix
P?=?sqrt（X^2+Y^2）;
%?direct?calculation?of?longitude
if?P?>?1.e-20
???dlambda?=?atan2（YX）*rtd;
el