空间直角坐标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)