% ************************************************************************
% 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);
