ver90togeo90.m

来自「本程序设计了常用的几个坐标系的转换」· M 代码 · 共 31 行

% PE-90直角坐标转换成PE-90地理坐标
clear;
input('\n\ninput vertical_90:');                                                % 输入PE-90坐标系下的直角坐标值
X_90=input('\nX_90=');                                                          % X轴分量（米）
Y_90=input('\nY_90=');                                                          % Y轴分量（米）
Z_90=input('\nZ_90=');                                                          % Z轴分量（米）
a=6378137;
f=1/298.257223563;
b=(1-f)*a;
e=2*f-f^2;
lamda_90=atan(Y_90/X_90)*180/pi;
if lamda_90<0
    lamda_90=180+lamda_90;
end
N(1)=a;
H(1)=sqrt(X_90^2+Y_90^2+Z_90^2)-sqrt(a*b);
phi(1)=atan(Z_90/sqrt(X_90^2+Y_90^2)*(1-e*N(1)/(N(1)+H(1)))^(-1))*180/pi;
N(2)=a/sqrt(1-e*(sin(phi(1)*pi/180))^2);
H(2)=sqrt(X_90^2+Y_90^2)/cos(phi(1)*pi/180)-N(2);
phi(2)=atan(Z_90/sqrt(X_90^2+Y_90^2)*(1-e*N(2)/(N(2)+H(2)))^(-1))*180/pi;
i=2;
while((abs(phi(i)-phi(i-1))>(0.00001/(60*60)))|(abs(H(i)-H(i-1))>0.001))
    i=i+1;
    N(i)=a/sqrt(1-e*(sin(phi(i-1)*pi/180))^2);
    H(i)=sqrt(X_90^2+Y_90^2)/cos(phi(i-1)*pi/180)-N(i);
    phi(i)=atan(Z_90/sqrt(X_90^2+Y_90^2)*(1-e*N(i)/(N(i)+H(i)))^(-1))*180/pi;
end
fprintf('\ni=%d \n',i);
fprintf('\nphi_90=%f \n',phi(i));
fprintf('\nlamda_90=%f \n',lamda_90);
fprintf('\nH_90=%f \n',H(i));