function elaz=sum_elaz(xobs,yobs,zobs,xsat,ysat,zsat)

robs=[xobs yobs zobs];

%make Up unit column vector: rhat
rhat = robs/sqrt(sum(robs.^2));
if find(size(rhat) == 3) == 2; rhat = rhat.'; end

%
%make North unit column vector: nhat
zhat = [0 0 1].';
north = zhat - sum(rhat.*zhat)*rhat;
nhat = north/sqrt(sum(north.^2));

%
%make East unit column vector: ehat
z_cross_I = [0 -1 0; 1 0 0; 0 0 0];
east = z_cross_I * rhat;
ehat = east/sqrt(sum(east.^2));

%
%compute vector from observer to satellite
xsatobs = xsat - robs(1);
ysatobs = ysat - robs(2);
zsatobs = zsat - robs(3);
xyzsatobs = [xsatobs, ysatobs, zsatobs];

%
%compute east, north and up components: ecomp, ncomp and rcomp
ecomp = xyzsatobs * ehat;
ncomp = xyzsatobs * nhat;
rcomp = xyzsatobs * rhat;

%
%compute zenith angle and elevation
mag = sqrt((xyzsatobs.^2)*[1;1;1]);
zenith = acos(rcomp./mag) * 180/pi;
el = 90 - zenith;

%
%compute azimuth
az = unwrap(atan2(ecomp, ncomp));
if (az<0)
az = az + 2 * pi;
end
azdeg = az * 180/pi;
elaz=[el azdeg];

return