xyz2blh.m

这是国外关于卫星导航方面一书的源代码

function [bb,l,h]=xyz2blh(x,y,z,ell)
% conversion: cartesian to ellipsoidal coordinates
%
% how:     [bb,l,h]=xyz2blh(x,y,z)
% how:     [bb,l,h]=xyz2blh(x,y,z,ellips)
% Input:   X,Y,Z in Meter
%          ELLIPS is a text string: 'wgs84' (default), 'grs80', or 'bessel'
% Output:  B,L in Degree.Decimals, H in Meter
%
% Formulae from Hofmann-Wellenhof et al. (1994):"GPS in der Praxis"
% Do not use negative heights !
%
% (c) iapg - 1996 - 2002

if nargin < 4, ell = 'wgs84'; end	% WGS84 default
ell = lower(deblank(ell));
if strcmp(ell,'bessel')			% Bessel-Ellipsoid 
 a=6377397.155;
 f=1/299.1528128;
elseif strcmp(ell,'wgs84') 		% WGS84-Ellipsoid
 a=6378137;
 f=1/298.257223563;
elseif strcmp(ell,'grs80') 		% GRS80-Ellipsoid
 a=6378137;
 f=1/298.257222101;
else
 error('undefined Ellipsoid')
end


b=(1-f)*a;
e1q=((a^2)-(b^2))/(a^2);
e2q=((a^2)-(b^2))/(b^2);

p=sqrt(x.^2+y.^2);
t=atan2(z.*a,(p.*b));
bb=atan2((z+e2q.*b.*(sin(t)).^3),(p-e1q.*a.*(cos(t)).^3));
l=atan2(y,x);
n=a./sqrt(1-e1q.*(sin(bb)).^2);
h=p./cos(bb)-n;

rho=180/pi;
bb=bb*rho;
l=l*rho;

% --------------------------------------------------------------------