geo84tover84.asv

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

% WGS-84下地理坐标转换成PE-90下
clear;
input('input geography_84:');                                                   % 输入WGS-84坐标系下的地理坐标值
phi_84=input('phi_84=');                                                        % WGS-84下的纬度（度）
lamda_84=input('lamda_84=');                                                    % WGS-84下的经度（度）
H_84=input('H_84=');                                                            % WGS-84下的高度（米）
a=6378137;                                                                      % WGS-84椭球体长半轴
f=1/298.257223563;                                                              % WGS-84椭球体扁率
e=2*f-f^2;                                                                      % WGS-84椭球体第一偏心率
N=a/sqrt(1-e^2*(sin(phi_84/180*pi))^2);                                         % N为卯酉圈曲率半径
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X_84=(N+H_84)*cos(phi_84/180*pi)*cos(lamda_84/180*pi);
Y_84=(N+H_84)*cos(phi_84/180*pi)*sin(lamda_84/180*pi);
Z_84=[N*(1-e^2)+H_84]*sin(phi_84/180*pi);
fprintf('X_84=%.3f \n',X_84);                                                   % 以小数形式输出，保留3位小数，且输出换行
fprintf('Y_84=%.3f \n',Y_84);
fprintf('Z_84=%.3f \n',Z_84);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X_90=1/(1+1.9e-6*1.9e-6)*[X_84+1.9e-6*(Y_84-2.5)];
Y_90=1/(1+1.9e-6*1.9e-6)*[-1.9e-6*X_84+(Y_84-2.5)];
Z_90=1/(1+1.9e-6*1.9e-6)*[(1+1.9e-6*1.9e-6)*Z_84];
fprintf('X_90=%.3f \n',X_90);                                           
fprintf('Y_90=%.3f \n',Y_90);
fprintf('Z_90=%.3f \n',Z_90);
fprintf('DeltaX=%f \n',X_84-X_90);
fprintf('DeltaY=%f \n',Y_84-Y_90);
fprintf('DeltaZ=%f \n',Z_84-Z_90);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lamda_90=180+atan(Y_90/X_90)*180/pi;
b=(1-f)*a;
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^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^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^2*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^2*(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^2*N(i)/(N(i)+H(i)))^(-1))*180/pi;
end
fprintf('phi_90=%f \n',phi(i));
fprintf('lamda_90=%f \n',lamda_90);
fprintf('H_90=%f \n',H(i));