function [xyz,VCxyz] = lph2xyz(lph,a,f_inv,VClph)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function: lph2xyz.m
% purpose: Convert longitude, latitude, & ellipsoidal height (lph)
%          to X, Y, & Z coord's
% author: David J. VandenBerg
% status: developed
% version: Matlab 5.2
% last date of development: May '99
% variables:
%   Input:
%     lph=1x3 column mtx of coordinates
%     a=semi-major axis
%     f_inv=inverse flattening
%     VClph=3x3 variance-covariance mtx of coordinates
%   Output:
%     xyz=1x3 column mtx of coordinates
%     VCxyz=3x3 variance-covariance mtx of coordinates
% example usage: [vec,mtx]=lph2xyz(vec,a,f_inv,mtx)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
eSQ=2*1/f_inv-(1/f_inv)^2;        % first eccentricity
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

xyz=[(N+lph(:,3)).*cos(lph(:,2)).*cos(lph(:,1)), ...      % calc. X
	(N+lph(:,3)).*cos(lph(:,2)).*sin(lph(:,1)), ...        % calc. Y
   (N*(1-eSQ)+lph(:,3)).*sin(lph(:,2))];          % calc. Z

% check for std dev's
if nargin==3
   % no error propagation
elseif nargin==4
   % jacobian: xyz WRT lph
   J=[-(N+lph(3))*sin(lph(1))*cos(lph(2))  -(M+lph(3))*cos(lph(1))*sin(lph(2))  cos(lph(1))*cos(lph(2))
      (N+lph(3))*cos(lph(1))*cos(lph(2))   -(M+lph(3))*sin(lph(1))*sin(lph(2))  sin(lph(1))*cos(lph(2))
      0                                    (M+lph(3))*cos(lph(2))               sin(lph(2))];
   % compile solution mtx
   VCxyz=J*VClph*J';            % solution mtx for loop
else
   disp('error: incorrect number of inputs for "lph2xyz"');
end;
return;