function [vn,svn,ve,sve,vu,svu,corr_ne,corr_eu,corr_nu]=vxyz2vneu(vx,vy,vz,covvxyz,lat,lon)
%function [vn,svn,ve,sve,vu,svu,corr_ne,corr_eu,corr_nu]=vxyz2vneu(vx,vy,vz,covvxyz,lat,lon)
%
% make sure lat, lon are in degrees;
%
% note change to vertical sign 3/2/05

lat=lat*pi/180;
lon=lon*pi/180;

Tnx = -sin(lat).*cos(lon);
Tex = -sin(lon);
Tdx = -cos(lat).*cos(lon);
Tny = -sin(lat).*sin(lon);
Tey =  cos(lon);
Tdy = -cos(lat).*sin(lon);
Tnz =  cos(lat);
Tez =  0;
Tdz = -sin(lat);

T=[Tnx Tny Tnz; Tex Tey Tez; Tdx Tdy Tdz];

% vn = Tnx.*vx + Tny.*vy + Tnz.*vz;
% ve = Tex.*vx + Tey.*vy + Tez.*vz;
% vu = -Tdx.*vx - Tdy.*vy - Tdz.*vz;

Vxyz=[vx;vy;vz];

Vneu = T*Vxyz;
vn=Vneu(1);
ve=Vneu(2);
vu=Vneu(3);

% to make up velocity positive (as per p. 155 "Plate Tectonics: How it Works" Cox & Hart)
vu = -vu;

covvneu = T*covvxyz*(T');

% svx=sqrt(covvxyz(1,1));
% svy=sqrt(covvxyz(2,2));
% svz=sqrt(covvxyz(3,3));

% svn = sqrt((Tnx.*svx).^2 + (Tny.*svy).^2 + (Tnz.*svz).^2);
% sve = sqrt((Tex.*svx).^2 + (Tey.*svy).^2 + (Tez.*svz).^2);
% svu = sqrt((Tdx.*svx).^2 + (Tdy.*svy).^2 + (Tdz.*svz).^2);
% 
% cov_ne = 0;
% corr_ne =  cov_ne./sqrt(svn.*sve);

svn = sqrt(covvneu(1,1));
sve = sqrt(covvneu(2,2));
svu = sqrt(covvneu(3,3));

cov_ne = covvneu(1,2);
corr_ne =  cov_ne./(svn.*sve);

cov_eu = covvneu(2,3);
corr_eu =  cov_eu./(sve.*svu);

cov_nu = covvneu(1,3);
corr_nu =  cov_nu./(svn.*svu);
