【GNSS】【Matlab】 XYZ_2_BLH BLH_2_XYZ
function [X,Y,Z]=BLH2xyz(B,L,H) %transform the (B L H) into (X Y Z) under the WGS84 coordinate system a=6378140.00; %a=6378137.00; % ??????????ο?WGS-84????????? b=6356752.3142; % ????? AE_WGS84=6378137.0; % Major axis of the WGS84 ellipsoid Alfa_WGS84=298.257223563; %1 over flattening of the WGS 84 ellipsoid e2=0.00669437999013; % ???????????? N=AE_WGS84/sqrt(1.0-e2*sin(B)*sin(B)); TmpV=N+H; X=TmpV*cos(B)*cos(L); Y=TmpV*cos(B)*sin(L); Z=(N*(1-e2)+H)*sin(B);
www.pudn.com > xyz2blh.rar > xyz2BLH.m, change:2012-08-08,size:1742b
Search codes
function BLH=xyz2BLH(X,Y,Z) % Transform the (X Y Z) into (B L H) under the WGS84 coordinate system % ?????????????? %%%%%%%%%%%%%%%%%%%%%???????ó??????????????%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%???%%% %{ a=6378137.00; % ??????????ο?WGS-84????????? b=6356752.3142; % ????? alfa=1/298.257223563; % ???? e2=0.00669437999013; % ???????????? L=atan2(Y,X); %L=azimuth(X,Y); % ?????????????λ?? N0=a; H0=sqrt(X^2+Y^2+Z^2)-sqrt(a*b); B0=atan(Z/sqrt(X^2+Y^2)/(1.0-e2*N0/(N0+H0))); while(true) N1=a/sqrt(1.0-e2*sin(B0)^2); H1=sqrt(X^2+Y^2)/cos(B0)-N1; B1=atan(Z/sqrt(X^2+Y^2)/(1.0-e2*N1/(N1+H1))); if ((abs(H1-H0)>0.001)&&(abs(B1-B0)>4.848132257047973e-11)) % B??????0.00001" N0=N1; H0=H1; B0=B1; else break; end end %} %%%???????ppt?????????д%%% a=6378137.00; % ??????????ο?WGS-84????????? b=6356752.3142; % ????? alfa=1/298.257223563; % ???? e2=0.00669437999013; % ???????????? AE_WGS84=6378137.0; % Major axis of the WGS84 ellipsoid Alfa_WGS84=298.257223563; %1 over flattening of the WGS 84 ellipsoid L=atan2(Y,X); %%????????B %N=AE_WGS84/sqrt(1.0-e2*sin(B)*sin(B)); B0=atan(Z/sqrt(X^2+Y^2)); N0=AE_WGS84/sqrt(1.0-e2*sin(B0)*sin(B0)); while(true) B1=atan((Z+N0*e2*sin(B0))/sqrt(X^2+Y^2)); N1=AE_WGS84/sqrt(1.0-e2*sin(B1)*sin(B1)); if (abs(B1-B0)>4.848132257047973e-11) B0=B1; N0=N1; else break; end end H1=sqrt(X^2+Y^2)*sec(B1)-N1; % if H1<0 H1=0; else H1=H1; end % if B1==NaN B1=pi/2; else B1=B1; end B1=B1*180/pi; L=L*180/pi; BLH.B=B1; BLH.L=L ; BLH.H=H1;