function [NEU,CNEU,LLH] = xyz2neu(O,V,SV,COV);

% XYZ2NEU  Convert ECEF coordinates into local (topocentric) frame.
%
%     Input:
%       O = origin vector in ECEF frame (m)
%       V = position or velocity vector in ECEF frame (m or m/yr)
%       SV = stdev in ECEF frame (m or m/yr)
%       COV = covariances, XY XZ YZ (m^2 or m^2/yr^2)
%       (NOTE: O, V, SV, COV can be n x 3 matrices, n = # of sites)
%
%     Output:
%       NEU  = output on NEU frame (m or m/yr)
%       CNEU = associated covariance (m or m^2/yr^2), format is:
%                    Cnn Cne Cnu Cee Ceu Cuu
%       (NOTE: NEU and CNEU will be matrices with n rows)
%       LLH  = coordinates of origin vector in WGS84
%              time lon lat elevation
%
%     Call: [NEU,CNEU,LLH] = xyz2neu(O,V,SV,COV);

% if O is a single point, make it the same size as V
if (size(O,1) == 1)
  XR = ones(size(V,1),1) .* O(1);  
  YR = ones(size(V,1),1) .* O(2);  
  ZR = ones(size(V,1),1) .* O(3);  
else
  XR = O(:,1); YR = O(:,2); ZR = O(:,3);
end

% read rest of input
vx = V(:,1); vy = V(:,2); vz = V(:,3);
svx = SV(:,1); svy = SV(:,2); svz = SV(:,3);
cxy = COV(:,1); cxz = COV(:,2); cyz = COV(:,3);

% convert origin vector to ellipsoidal coordinates
T = zeros(size(XR,1),1);
LLH = xyz2wgs([T XR YR ZR]);

% compute sines and cosines
cp = cos(LLH(:,2).*pi/180); sp = sin(LLH(:,2).*pi/180); % longitude
cl = cos(LLH(:,3).*pi/180); sl = sin(LLH(:,3).*pi/180); % latitude

% for each site
NEU = [];
CNEU = [];
for i=1:size(V,1)
  % build the rotation matrix
  R = [ -sl(i)*cp(i)   -sl(i)*sp(i)    cl(i);
         -sp(i)            cp(i)           0;
         cl(i)*cp(i)    cl(i)*sp(i)    sl(i)];

  % apply the rotation
  NEUi = R * [vx(i);vy(i);vz(i)];

  % build covariance for that site
  CVi = [svx(i)^2 cxy(i)   cxz(i);
         cxy(i)   svy(i)^2 cyz(i);
         cxz(i)   cyz(i)   svz(i)^2];

  % propagate covariance
  CNEUi = R * CVi * R';

  % increment result matrices
  NEU = [NEU;
         NEUi(1) NEUi(2) NEUi(3)];
  CNEU = [CNEU;
          CNEUi(1,1) CNEUi(1,2) CNEUi(1,3) CNEUi(2,2) CNEUi(2,3) CNEUi(3,3)];
end

