function [ENUCov semiMinor semiMajor orientation] = xyzcov2enucov( Lat, Lon, xyzcov )

% First form the Jacobian Matrix:
%  J = | -sin(lat)   -sin(lat)*cos(long)  cos(lat)*cos(long) |
%      |  cos(long)  -sin(lat)*sin(long)  cos(lat)*sin(long) |
%      |     0              cos(lat)            sin(lat)     |

dRAD2DEG = 180 / pi;

Lat = Lat ./ dRAD2DEG;
Lon = Lon ./ dRAD2DEG;

sinLat = sin(Lat);
sinLon = sin(Lon);
cosLat = cos(Lat);
cosLon = cos(Lon);

for i=1:length(Lat)
  J(1,1) = -sinLon(i);
  J(1,2) =  cosLon(i);
  J(1,3) =  0;
  J(2,1) = -sinLat(i).*cosLon(i);
  J(2,2) = -sinLat(i).*sinLon(i);
  J(2,3) =  cosLat(i);
  J(3,1) =  cosLat(i).*cosLon(i);
  J(3,2) =  cosLat(i).*sinLon(i);
  J(3,3) =  sinLat(i);

  COV(1:3,1:3) = xyzcov(i,:,:);

  ENUCov(i,:,:) = J .* COV .* J';

end

% Now compute the Semi-major/Semi-minor and orientation in the
% horizontal plane

qEE = ENUCov(:,1,1);
qNN = ENUCov(:,2,2);
qEN = ENUCov(:,1,2);

% compute eigenvalues 
tmp1          = 0.5 .* (qEE + qNN) ;
qEE_minus_qNN = qEE - qNN ;
tmp2          = 0.5 .* sqrt( qEE_minus_qNN .* qEE_minus_qNN + 4.0 .* qEN .* qEN ) ;

lambda1 = tmp1 + tmp2 ;
lambda2 = tmp1 - tmp2 ;

% compute axis length
semiMajor = sqrt( lambda1 ) ;
semiMinor = sqrt( lambda2 ) ;

% compute semi-major eigenvector.  To avoid divide by 0, multiply
%   entire vector by -qEN:
x1 = -qEN ;
y1 = qEE - lambda1 ;

% compute semi-major orientation in degrees
orientation = atan2( x1, y1 ) ./ dRAD2DEG;


