matlab坐标转换

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 < 1.e-20
   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