zhijiao_dadi.m

zhijiao_dadi将直角坐标XYZ( BJ54 ,GDZ80)转换为大地坐标下的BLH.请多多指教！

function [BLH]=zhijiao_dadi(XYZ,type,L)    %  直角系'WGS' or 'BJ54' or 'GDZ80' or 'DXZ' －－－－－大地系
                                                 
%  zhijiao_dadi将直角坐标XYZ( BJ54 ,GDZ80)转换为大地坐标下的BLH.
%  数据格式: XYZ(1:n,1:3)  ,BLH(1:n,1:3)  , L 为迭代次数
% XYZ( BJ54 ,GDZ80 )--------BLH      迭代计算法                
% type = 'WGS' or 'BJ54' or 'GDZ80' or 'DXZ' 
%  B 纬度，L 经度，H  高度
% BJ54 (XYZ) ------------ (BLH)

n=size(XYZ,1); 

if strcmp(type,'BJ54')              % 测试 X=3694472.468; Y=3694472.468;  Z=5194534.424;  L=45 B=45 H=1e+6
        a_BJ54 = 6378245; alpha_BJ54 = 1/298.3 ;
        e2_BJ54 = 2*alpha_BJ54-alpha_BJ54^2;      % BJ54的长半轴和扁率
        
        t1(1:n,1:L)=0;                        % t1 存放迭代计算的tanB的值。           
        BLH(1:n,2) = atan(XYZ(1:n,2)/XYZ(1:n,1))*180/pi;      % L
        t1(1:n,1) = XYZ(1:n,3).*(1+e2_BJ54)./sqrt(XYZ(1:n,1).^2+XYZ(1:n,2).^2);  % tanB的初值
                      
        
        for k=2:L               % L为迭代次数    
            t1(1:n,k) = (1./sqrt(XYZ(1:n,1).^2+XYZ(1:n,2).^2)).*(XYZ(1:n,3)+a_BJ54*e2_BJ54.*t1(1:n,k-1)./sqrt( 1+(1-e2_BJ54).*t1(1:n,k-1).^2 ) );
        end
        BLH(1:n,1) = atan(t1(1:n,L));                % B
         N_BJ54(1:n,1) = a_BJ54./sqrt(1-e2_BJ54.*sin(BLH(1:n,1)).^2);    % N 卯酉圈曲率半径
      
        BLH(1:n,3) = sqrt(XYZ(1:n,1).^2+XYZ(1:n,2).^2)./cos(BLH(1:n,1))-N_BJ54;   % H
        BLH(1:n,1) = BLH(1:n,1)*180/pi;  

% GDZ80  (X,Y,Z) ------------ (L,B,H)
elseif strcmp(type,'GDZ80')
        a_GDZ80 = 6378140;  alpha_GDZ80 = 1/298.257;
        e2_GDZ80 = 2*alpha_GDZ80-alpha_GDZ80^2;
          % GDZ80的长半轴和扁率
          
        t1(1:n,1:L)=0;                                  
        BLH(1:n,2) = atan(XYZ(1:n,2)./XYZ(1:n,1));      % L
        
        t1(1:n,1) = XYZ(1:n,3).*(1+e2_GDZ80)./sqrt(XYZ(1:n,1).^2+XYZ(1:n,2).^2);  % tanB的初值
      
        
        for k=2:L               % L为迭代次数    
            t1(1:n,k) = (1./sqrt(XYZ(1:n,1).^2+XYZ(1:n,2).^2)).*(XYZ(1:n,3)+a_GDZ80*e2_GDZ80.*t1(1:n,k-1)./sqrt( 1+(1-e2_GDZ80).*t1(1:n,k-1).^2 ) );
        end
        BLH(1:n,1) = atan(t1(1:n,L));                % B
        N_GDZ80(1:n,1) = a_GDZ80./sqrt(1-e2_GDZ80.*sin(BLH(1:n,1)).^2);    % N
          
        BLH(1:n,3) = sqrt(XYZ(1:n,1).^2+XYZ(1:n,2).^2)./cos(BLH(1:n,1))-N_GDZ80;   % H    
        BLH(1:n,1) = BLH(1:n,1)*180/pi;
end    
%%%%%%%%%%%%%% *** last line of zhijiao_dadi.m%%%%%%%%%%%%%%%%%%%%%%%%%%%v





% 测试　 XYZ=[3694472.468  3694472.468  5194534.424];type='BJ54',BLH=zhijiao_dadi(XYZ,type,3)
% 测试结果　 BLH = [0.785398  45  1e+6]
%  t1(1,1)=1.00087  t1(1,2)=1 t1(1,2)=1