function RAH = NEU2POLAR ( NEU )
%YSNTAX:
% RAH = NEU2POLAR (NEU)
% 站心坐标系转换成站心极坐标
% transform topocentric coordinate (NEU) into topocentric polar coordinate
% (RAH)
%
%INPUTS:
% NEU: (3x1) vector of topocentric coordinates
%
%OUTPUT:
% RAH: (3x1) vector the topocentric polar coordinates
% (1): the polar range in meter
% (2): the azimuth in radian
% (3): the elevation in radian
%
%e-mail:Bofeng_Li163.com
%%===================BEGIN PROGRAM========================%%
if (nargin < 1)
error('not enough input arguments');
end
[m, n] = size(NEU);
if (m~=3)
error('incorrect dimension of input');
end
RAH = zeros(3, n);
RAH(1,:) = sqrt(NEU(1,:).^2 + NEU(2,:).^2 + NEU(3,:).^2);
RAH(2,:) = get_azim(NEU(1,:), NEU(2,:));
RAH(3,:) = atan(NEU(3,:)./sqrt(NEU(1,:).^2 + NEU(2,:).^2));
return
%%===========================END PROGRAM=================================%%
function alpha = get_azim(dx, dy)
%SYNTAX:
% alpha = get_azim(dx,dy)
%
% compute the azimuth based on the coordinate differences
%
%INPUTS:
% dx: coordinate difference of north component between two points
% dy: coordinate difference of east component between two points
%
%OUTPUT:
% alpha: calculated azimuth in radian degree
%
%%=============================BEGIN PROGRAM=========================%%
if (nargin < 2)
error('not enough input arguments');
end
if(size(dx)~=size(dy))
error('incorrect dimensions of inputs');
end
n = length(dx);
alpha = zeros(n,1);
for i = 1 : n
if ( dx(i) == 0)
if (dy(i)>0)
alpha(i) = pi/2;
else
alpha(i) = 3*pi/2;
end
else
R = atan(dy(i)/dx(i));
if (dx(i)>0)
if (dy(i)>=0)
alpha(i) = R;
else
alpha(i) = R + 2*pi;
end
else
alpha(i) = R + pi;
end
end
end
return
