function [lph] = eas2lph(eas,lphOrig,a,f_inv)
%function [lph,VClph] = eas2lph(eas,lphOrig,a,f_inv,VCeas,VClphOrig)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function: eas2lph.m
% purpose: Convert Elev, Az, SR (& errors) to lon, lat, height (& errors)
% author: David J. VandenBerg
% reference: Leick, 1995, chap. 7
% status: developed
% version: Matlab 5.2
% last date of development: May '99
% variables:
%   Input:
%     eas=1x3 column mtx of coordinates
%     lphOrig=1x3 column mtx of enu origin coordinates
%     a=semi-major axis
%     f_inv=inverse flattening
%     VCeas=3x3 variance-covariance mtx of coordinates
%     VClphOrig=3x3 variance-covariance mtx of origin coordinates
%   Output:
%     lph=1x3 column mtx of coordinates
%     VClph=3x3 variance-covariance mtx of coordinates
% example usage: eas2lph(mtx,vec,a,f_inv,mtx,vec)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% convert E,N,U to El,Az,Sr
enu=[eas(:,3).*cos(eas(:,1)).*sin(eas(:,2)), ...     % calc. East
	eas(:,3).*cos(eas(:,1)).*cos(eas(:,2)), ...        % calc. North
    eas(:,3).*sin(eas(:,1))];                  % calc. Up

% conversion to lph
lph=enu2lph(enu,lphOrig,a,f_inv);         % convert to El,Az,Sr

% translate coordinates
% xyzOrig=lph2xyz(lphOrig,a,f_inv);          % convert origin to xyz
% xyz=lph2xyz(lph,a,f_inv);

% % check for std dev's
% if nargin==4
%    % no error propagation
% elseif nargin==6
%    K=zeros(3,6);                             % initialize solution mtx
%    
%    % origin station constants
%    eSQ=2*1/f_inv-(1/f_inv)^2;                % first eccentricity squared
%    W=sqrt(1-eSQ*sin(lphOrig(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 for origin station: xyz WRT lph
%    J1=[-(N+lphOrig(3))*sin(lphOrig(1))*cos(lphOrig(2))  -(M+lphOrig(3))*cos(lphOrig(1))*sin(lphOrig(2))  cos(lphOrig(1))*cos(lphOrig(2))
%       (N+lphOrig(3))*cos(lphOrig(1))*cos(lphOrig(2))    -(M+lphOrig(3))*sin(lphOrig(1))*sin(lphOrig(2))  sin(lphOrig(1))*cos(lphOrig(2))
%       0                                                 (M+lphOrig(3))*cos(lphOrig(2))                   sin(lphOrig(2))];
%    
%    % remote station constants
%    W2=sqrt(1-eSQ*sin(lph(2))^2);             % "W" constant for latitude
%    M2=a*(1-eSQ)/W2^3;                        % radius of curvature in the meridian
%    N2=a/W2;                                  % radius of curvature in the prime vertical
%    
%    % jacobian for remote station: xyz WRT lph
%    J2=[-(N2+lph(3))*sin(lph(1))*cos(lph(2))  -(M2+lph(3))*cos(lph(1))*sin(lph(2))  cos(lph(1))*cos(lph(2))
%       (N2+lph(3))*cos(lph(1))*cos(lph(2))    -(M2+lph(3))*sin(lph(1))*sin(lph(2))  sin(lph(1))*cos(lph(2))
%       0                                      (M2+lph(3))*cos(lph(2))               sin(lph(2))];
%    
%    % solution constants
%    Dx=xyz(1)-xyzOrig(1);
%    Dy=xyz(2)-xyzOrig(2);
%    Dz=xyz(3)-xyzOrig(3);
%    el=eas(1);
%    az=eas(2);
%    s=eas(3);
%    
%    % jacobian: eas WRT xyz
%    g11=(-s*cos(lphOrig(2))*cos(lphOrig(1))+sin(el)*Dx)/(s^2*cos(el));
%    g12=(-s*cos(lphOrig(2))*sin(lphOrig(1))+sin(el)*Dx)/(s^2*cos(el));
%    g13=(-s*sin(lphOrig(2))+sin(el)*Dz)/(s^2*cos(el));
%    g21=(-sin(lphOrig(2))*cos(lphOrig(1))*sin(az)+sin(lphOrig(1))*cos(az))/(s*cos(el));
%    g22=(-sin(lphOrig(2))*sin(lphOrig(1))*sin(az)-cos(lphOrig(1))*cos(az))/(s*cos(el));
%    g23=(cos(lphOrig(2))*sin(az))/(s*cos(el));
%    g31=-Dx/s;
%    g32=-Dy/s;
%    g33=-Dz/s;
%    G=[g11  g12  g13               % compile mtx
%       g21  g22  g23
%       g31  g32  g33
%    ];
%    
%    % compile solution mtx             % solution mtx for loop
%    VClph=inv(-G*J2)*(VCeas-G*J1*VClphOrig*(G*J1)')*inv((-G*J2)');
% else
%    disp('error: incorrect number of inputs for "eas2lph"');
% end;