1. 程式人生 > >【GNSS】【Matlab】 XYZ_2_BLH BLH_2_XYZ

【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;