function [lph,VClph] = enu2lph(enu,lphOrig,a,f_inv,VCenu)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function: lph2enu.m
% purpose: Convert easting, northing, up coord's with errors
%          to lon, lat, ellipsoidal height with errors
% author: David J. VandenBerg
% status: developed
% version: Matlab 5.2
% last date of development: May, '99
% variables:
%   Input:
%     enu=1x3 column mtx of coordinates
%     lphOrig=1x3 column mtx of coordinates
%     a=semi-major axis
%     f_inv=inverse flattening
%     VCenu=3x3 variance-covariance mtx of coordinates
%   Output:
%     lph=1x3 column mtx of coordinates
%     VClph=3x3 variance-covariance mtx of coordinates
% example usage: [mtx,mtx]=enu2lph(mtx,vec,a,f_inv,mtx)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Rotation matrices for conversion to e,n,u
w=pi/2-lphOrig(2);                             % rotation about x-axis
k=pi/2+lphOrig(1);                             % rotation about z-axis
% rotation mtx
% R=[cos(k)  -cos(w)*sin(k) sin(w)*sin(k)
%    sin(k)  cos(w)*cos(k)  -sin(w)*cos(k)
%    0       sin(w)         cos(w)];
% Dxyz=R*enu';                               % rotate points

Dxyz=[cos(k)*enu(:,1)-cos(w)*sin(k)*enu(:,2)+sin(w)*sin(k)*enu(:,3), ...
   sin(k)*enu(:,1)+cos(w)*cos(k)*enu(:,2)-sin(w)*cos(k)*enu(:,3), ...
   sin(w)*enu(:,2)+cos(w)*enu(:,3)];

% translate coordinates
xyzOrig=lph2xyz(lphOrig,a,f_inv);         % convert origin to xyz
xyz=[xyzOrig(1)*ones(size(Dxyz,1),1)  xyzOrig(2)*ones(size(Dxyz,1),1)  xyzOrig(3)*ones(size(Dxyz,1),1)]+Dxyz;                         % translate points
%xyz=xyz';

% convert to geocentric to geodetic
lph=xyz2lph(xyz,a,f_inv);

% check for std dev's
if nargin==4
   % no error propagation
elseif nargin==5
   % point constants
   eSQ=2*1/f_inv-(1/f_inv)^2;      % first eccentricity squared
   W=sqrt(1-eSQ*sin(lph(2))^2);    % "W" constant for latitude
   M=a*(1-eSQ)/W^3;                % radius of curvature in the meridian
   N=a/W;                          % radius of curvature in the prime vertical
   
   % jacobian enu WRT lph
   H=[(N+lph(3))*cos(lph(2))  0           0
      0                       (M+lph(3))  0
      0                       0           1];
   
   % compile solution mtx
   H_inv=inv(H);                % inverse of metric mtx
   VClph=H_inv*VCenu*H_inv';    % solution mtx for loop
else
   disp('error: incorrect number of inputs for "enu2lph"');
end;
return;