大地坐标与空间直角坐标互相转换

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