空间直角坐标XYZ换算为经纬度BL的matlab编程

空间直角坐标XYZ换算为经纬度BL程序
（matlab编程）
度分秒转弧度函数代码：
function hd=dzh(a)
hh=sign(a);
a=abs(a);
hd=hh*(fix(a)+fix((a-fix(a))*100)/60+((a-fix(a))*100-...
fix((a-fix(a))*100))*100/3600)*pi/180;%度分秒转化为弧度
end
弧度转度分秒函数代码：
function jd=hzd(x)
jd=fix(x*206264.8062470964/3600)+fix((x*206264.8062470964/3600-...
fix(x*206264.8062470964/3600))*60)/100+((x*206264.8062470964/3600-...
fix(x*206264.8062470964/3600))*60-fix((x*206264.8062470964/3600-...
fix(x*206264.8062470964/3600))*60))*60/10000;%弧度转化为度分秒
end
主程序代码：
fprintf('-----克拉索夫斯基椭球体请输入1；1975年国际椭球体请输入2；WGS 84椭球请输入3-----')
kk=input('请输入：');
if kk==1
a=6378245;%长半轴克拉索夫斯基椭球体
b=6356863.019;%短半轴
elseif kk==2
a=6378140; %长半轴 1975年国际椭球体
b=6356755.288; %短半轴
else
a=6378137;%长半轴 1975年国际椭球体
b=6356752.314; %短半轴
end
e1=sqrt(a^2-b^2)/a; %第一偏心率
c=a^2/b;
X=input('请输入X:');
Y=input('请输入Y:');
Z=input('请输入Z:');
L=atan(Y/X);%求大地经度L
B1=atan(Z/sqrt(X^2+Y^2));%求纬度初值
for i=1:+inf %设定循环，进行迭代
yita=e1*cos(B1);
V=sqrt(1+yita^2);
N1=c/V ;%%卯酉圈曲率半径
B2=atan((Z+N1*e1^2*sin(B1))/sqrt(X^2+Y^2));%迭代后的B
if abs(B2-B1)<0.00000001%如果两次相邻的B相差小于0.002秒则终止循环 break
end
B1=B2;%用新的B替换旧B
end
B=hzd(B2)%将结果转化为度分秒
L=hzd(L)。