function [LBH2] = diff_LBH(LBH1,type1,type2)

% diff_LBH 能對不同的大地坐標系進行互相轉換（GDZ80,BJ54）
%  LBH1(1:n,1:3)為輸入，LBH2(1:n,1:3)為輸出，LBH1(1:n,2)為緯度B，
% type1是輸入大地系的類型（GDZ80,BJ54），type2是輸出大地系的類型
 
 
n = size(LBH1,1);

if strcmp(type1,'BJ54')
    a1 = 6378245; alpha1 = 1/298.3 ;
    e21 = 2*alpha1-alpha1^2; 
    
    XYZ01 = [3 4 5];    % 原點的坐標，用于計算橢球平移參數
elseif strcmp(type1,'GDZ80') 
        a1 = 6378140;  alpha1 = 1/298.257;
        e21 = 2*alpha1-alpha1^2;
        XYZ01 = [6 7 8];    % 原點的坐標，用于計算橢球平移參數
 end
   
  if strcmp(type2,'BJ54')   
      a2 = 6378245; alpha2 = 1/298.3 ;
      e22 = 2*alpha2-alpha2^2; 
  elseif strcmp(type2,'GDZ80')
        a2 = 6378140;  alpha2 = 1/298.257;
        e22 = 2*alpha2-alpha2^2;

  end
   
   if strcmp(type1,'BJ54')&strcmp(type2,'GDZ80')
      dxyz = [111 -95 -75];  % GDZ80-BJ54
  elseif strcmp(type1,'GDZ80')&strcmp(type2,'GDZ80')
      dxyz = [-111 95 75];
  end
  

   da = a2-a1; dalpha = alpha2-alpha1;
      
   N = a1./sqrt(1-e21.*sin(LBH1(1:n,2)).^2);   % 卯酉圈曲率半徑
   M = a1*(1-e21)./sqrt((1-e21.*sin(LBH1(1:n,2)).^2).^3);  % 子午圈曲率半徑
   
  for k = 1:n                   
      A(1:3,1:3) = [ -sin(LBH1(k,1))/(N(k)+LBH1(k,2))*cos(LBH1(k,2))      cos(LBH1(k,1))/((N(k)+LBH1(k,3))*cos(LBH1(k,2)))     0;
                     -sin(LBH1(k,2))*cos(LBH1(k,1))/(M(k)+LBH1(k,3))        sin(LBH1(k,2))*sin(LBH1(k,1))/(M(k)+LBH1(k,3))    cos(LBH1(k,2))/(M(k)+LBH1(k,3));
                     cos(LBH1(k,2))*cos(LBH1(k,1))           cos(LBH1(k,2))*sin(LBH1(k,1))        sin(LBH1(k,2));  ];
      B(1:3,1:2) = [ 0                         0;
                     N(k)*e21*sin(LBH1(k,2))*cos(LBH1(k,2))/(a1*(N(k)+LBH1(k,3)))      M(k)*(2-e21*sin(LBH1(k,2))^2)*sin(LBH1(k,2))*cos(LBH1(k,2))/((1-a1)*(M(k)+LBH1(k,3)));
                     -N(k)*(1-e21*sin(LBH1(k,2))^2)/a1            M(k)*(1-e21*sin(LBH1(k,2))^2)*sin(LBH1(k,2))^2/(1-a1);  ];
      dLBH(k,1:3) = ( A*dxyz'+B*[da dalpha]')' ;  %  把dLBH由列向量化為行向量。dLBH為大地坐標修正量。
      dLBH(k,1:2) = 206265.*dLBH(k,1:2);    % 將dL,dB化為以弧度秒″為單位
            
      LBH2(k,1:3) = LBH1(k,1:3)+dLBH(k,1:3);   % 將原坐標值加上坐標修正值，即為新的坐標值。
  end

  
  %　測試數據   LBH1(1:4,1)=[20;21;22;23];LBH1(1:4,2)=[20;21;22;23];LBH1(1:4,3)=[100;120;140;160];
  %　type1='GDZ80',type2 = 'BJ54', [LBH2]=diff_LBH(LBH1,type1,type2)
  
      
                 
%%%%%%%%%%%%%%%%%%%%    *** last line of diff_LBH.m ***   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%