wgs84tope90.m

本程序设计了常用的几个坐标系的转换

% WGS-84坐标转换成PE-90坐标
clear;
input('\n\ninput geography_84:');                                                   % 输入WGS-84坐标系下的地理坐标值
i=1;
while i<=6 
fprintf('\ni=%d \n',i);
phi_84(i)=input('\nphi_84(i)= ');                                                   % WGS-84下的纬度（度）
lamda_84(i)=input('\nlamda_84(i)= ');                                               % WGS-84下的经度（度）
H_84(i)=input('\nH_84(i)= ');                                                       % WGS-84下的高度（米）
a=6378137;                                                                          % WGS-84椭球体长半轴
f=1/298.257223563;                                                                  % WGS-84椭球体扁率
e=sqrt(0.00669437999013);                                                           % WGS-84椭球体第一偏心率
N=a/sqrt(1-e^2*(sin(phi_84(i)/180*pi))^2);                                          % N为卯酉圈曲率半径
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%WGS-84地理坐标转换成WGS-84直角坐标
X_84(i)=(N+H_84(i))*cos(phi_84(i)/180*pi)*cos(lamda_84(i)/180*pi);
Y_84(i)=(N+H_84(i))*cos(phi_84(i)/180*pi)*sin(lamda_84(i)/180*pi);
Z_84(i)=[N*(1-e^2)+H_84(i)]*sin(phi_84(i)/180*pi);
fprintf('\nX_84(i)=%.3f \n',X_84(i));                                               % 以小数形式输出，保留3位小数，且输出换行
fprintf('\nY_84(i)=%.3f \n',Y_84(i));
fprintf('\nZ_84(i)=%.3f \n',Z_84(i));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%WGS-84直角坐标转换成PE-90直角坐标
X_90(i)=1/(1+1.9e-6*1.9e-6)*[X_84(i)+1.9e-6*(Y_84(i)-2.5)];
Y_90(i)=1/(1+1.9e-6*1.9e-6)*[-1.9e-6*X_84(i)+(Y_84(i)-2.5)];
Z_90(i)=1/(1+1.9e-6*1.9e-6)*[(1+1.9e-6*1.9e-6)*Z_84(i)];
fprintf('\nX_90(i)=%.3f \n',X_90(i));                                           
fprintf('\nY_90(i)=%.3f \n',Y_90(i));
fprintf('\nZ_90(i)=%.3f \n',Z_90(i));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%PE-90直角坐标转换成PE-90地理坐标
lamda_90(i)=atan(Y_90(i)/X_90(i))*180/pi;                                            %保证lamda_90在0～180度内
if lamda_90(i)<0
    lamda_90(i)=180+lamda_90(i);
end
b=(1-f)*a;
N(1)=a;
H(1)=sqrt(X_90(i)^2+Y_90(i)^2+Z_90(i)^2)-sqrt(a*b);
phi(1)=atan(Z_90(i)/sqrt(X_90(i)^2+Y_90(i)^2)*(1-e^2*N(1)/(N(1)+H(1)))^(-1))*180/pi;
N(2)=a/sqrt(1-e^2*(sin(phi(1)*pi/180))^2);
H(2)=sqrt(X_90(i)^2+Y_90(i)^2)/cos(phi(1)*pi/180)-N(2);
phi(2)=atan(Z_90(i)/sqrt(X_90(i)^2+Y_90(i)^2)*(1-e^2*N(2)/(N(2)+H(2)))^(-1))*180/pi;
k=2;
while((abs(phi(k)-phi(k-1))>=(0.00001/(60*60)))|(abs(H(k)-H(k-1))>=0.001))
    k=k+1;
    N(k)=a/sqrt(1-e^2*(sin(phi(k-1)*pi/180))^2);
    H(k)=sqrt(X_90(i)^2+Y_90(i)^2)/cos(phi(k-1)*pi/180)-N(k);
    phi(k)=atan(Z_90(i)/sqrt(X_90(i)^2+Y_90(i)^2)*(1-e^2*N(k)/(N(k)+H(k)))^(-1))*180/pi;
end
phi_90(i)=phi(k);
H_90(i)=H(k);
fprintf('\nphi_90(i)=%f \n',phi_90(i));
fprintf('\nlamda_90(i)=%f \n',lamda_90(i));
fprintf('\nH_90(i)=%f \n',H_90(i));  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 输出两坐标系下的误差值
Delta_X(i)=abs(X_84(i)-X_90(i));
Delta_Y(i)=abs(Y_84(i)-Y_90(i));
Delta_Z(i)=abs(Z_84(i)-Z_90(i));
Delta_phi(i)=abs(phi_84(i)-phi_90(i));
Delta_lamda(i)=abs(lamda_84(i)-lamda_90(i));
Delta_H(i)=abs(H_84(i)-H_90(i));
fprintf('\nDelta_X(i)=%f \n',X_84(i)-X_90(i));                                                                                    
fprintf('\nDelta_Y(i)=%f \n',Y_84(i)-Y_90(i));
fprintf('\nDelta_Z(i)=%f \n',Z_84(i)-Z_90(i));
fprintf('\nDelta_phi(i)=%f \n',phi_84(i)-phi_90(i));                                       
fprintf('\nDelta_lamda(i)=%f \n',lamda_84(i)-lamda_90(i));
fprintf('\nDelta_H(i)=%f \n',H_84(i)-H_90(i));
i=i+1;
end