##引言 座標轉換程式碼，由大地座標（經、緯度、高）與空間直角座標（XYZ）互相轉換。MATLAB程式碼。 ##BLH2XYZ

function [X,Y,Z] = BLh2XYZ(LAT,LON,h) 

% Example:
%(BLh)WGS84-(XYZ)WGS84
% LAT = 40.9987167395335;
% LON = 39.7652393428761;
% h = 51.403;

refell = 1;

switch refell
    case 1
        % IERS 2003 numerical standards
        % ellipsoid parameters for xyz2ellip.m
        a_tidefree = 6378136.6; %m      Equatorial radius of the Earth
        f_tidefree = 1/298.25642;     % Flattening factor of the Earth
        a = a_tidefree;  %m      Equatorial radius of the Earth
        f = f_tidefree;       % Flattening factor of the Earth
    case 2
        % GRS 80 (http://www.bkg.bund.de/nn_164850/geodIS/EVRS/EN/References/...
        %  Definitions/Def__GRS80-pdf,templateId=raw,property=publication...
        %  File.pdf/Def_GRS80-pdf.pdf)
        a_grs80    = 6378137;
        f_grs80    = 0.00335281068118;
        a = a_grs80;   %m      Equatorial radius of the Earth
        f = f_grs80;        % Flattening factor of the Earth
    case 3
        % WGS84
        a=6378137;
        f=1/298.25722356;
    case 4
        % Hayford
        a=6378388;
        f=1/297;
end

b = a-f*a;

lat=LAT*pi/180;
lon=LON*pi/180;
e2=(a^2-b^2)/a^2;
N=a/sqrt(1-e2*(sin(lat)^2));
X=(N+h)*cos(lat)*cos(lon);
Y=(N+h)*cos(lat)*sin(lon);
Z=(N*(1-e2)+h)*sin(lat);
 

##XYZ2BLH

% ************************************************************************
%   Description:
%   Transformation from Cartesian coordinates X,Y,Z to ellipsoidal 
%   coordinates lam,phi,elh. based on Occam subroutine transf.
%
%   Input:										
%      pos = [x,y,z]                 [m,m,m]
%               can be a matrix: number of rows = number of stations
%
%   Output:
%      coor_ell = [lat,lon,h]      [rad,rad,m]
% 
%   External calls: 	
%      global   a_...              Equatorial radius of the Earth [m]     
%               f_...              Flattening factor of the Earth
%
%   Coded for VieVS: 
%   17 Dec 2008 by Lucia Plank
%
%   Revision:
%   
% *************************************************************************
function [lat,lon,h]=xyz2ell(pos)

% choose reference ellipsoid
%       1 ...... tide free
%       2 ...... GRS80
refell =1;

switch refell
    case 1
        % IERS 2003 numerical standards (tide free crust)
        % ellipsoid parameters for xyz2ellip.m
        a = 6378136.6; %m      Equatorial radius of the Earth
        f = 1/298.25642;     % Flattening factor of the Earth
    case 2
        % GRS 80
        % (http://www.bkg.bund.de/nn_164850/geodIS/EVRS/EN/References/...
        %  Definitions/Def__GRS80-pdf,templateId=raw,property=publication...
        %  File.pdf/Def_GRS80-pdf.pdf)
        a = 6378137;  %m      Equatorial radius of the Earth
        f = 0.00335281068118;       % Flattening factor of the Earth
end

e2=2*f-f^2;

lon=angle(pos(:,1)+1i*pos(:,2));

lat=angle(sqrt(pos(:,1).^2+pos(:,2).^2)+1i*pos(:,3));
for j=1:10
  N=a./sqrt(1-e2*sin(lat).*sin(lat));
  h=sqrt(pos(:,1).^2+pos(:,2).^2)./cos(lat)-N;
  lat=angle(sqrt(pos(:,1).^2+pos(:,2).^2).*((1-e2)*N+h)+1i*pos(:,3).*(N+h));
end

%lat=cart2phigd(pos);