function [Ex, Ey, Ez, Hx, Hy]=mt1DAniAnalyticField(freqs,sigma,zNode,eTop)
%input:
% fre:are the computational frequencies
% sigma includes the conductivity and Euler angle of the calculation element
% znode  is the node coordinate
% eTop is Calculation mode: XY,eTop=[1;0]
%                           YX,eTop=[0;1]
%output:
%Ex is the boundary electric field along the x direction
%Ey is the boundary electric field along the y direction
%Ez is the boundary electric field along the z direction
%Hx is the boundary magnetic field along the x direction
%Hy is the boundary magnetic field along the y direction
if size(sigma,1)~=length(zNode)-1
    fprintf('error')
end
h   = diff(zNode);
    % assume halfspace at the bottom of model
sigma    = [sigma; sigma(end:end, :)];
    % effective azimuthal anisotropic conductivity
A= getEffSigma(sigma);% get Axx, Axy, Ayx
nFreq  = length(freqs);
nLayer = size(sigma, 1);

Ex = zeros(nLayer, nFreq);
Ey = zeros(nLayer, nFreq);
Ez = zeros(nLayer-1, nFreq);
Hx = zeros(nLayer, nFreq);
Hy = zeros(nLayer, nFreq);

sh = @(x) 0.5*(exp(x)-exp(-x));   % sinh()
ch = @(x) 0.5*(exp(x)+exp(-x));   % cosh()

MU0 = 4*pi*1e-7;

for iFreq=1:nFreq
    omega = 2 * pi * freqs(iFreq);
    iom = 1i*omega*MU0;
    
    Sprod = zeros(4, 4, nLayer);
    Sprod(:,:,end) = eye(4);
    
    % loop over layers to get the accumulative product of S (i.e. Sprod). 
    for j = nLayer:-1:1
        
        ada = A.xx(j) + A.yy(j);
        add = A.xx(j) - A.yy(j);

        if add>=0
            dA12 =  sqrt( add^2 + 4*A.xy(j)^2 );
        elseif add<0
            dA12 = -sqrt( add^2 + 4*A.xy(j)^2 );
        end

        A1 = 0.5 * (ada + dA12);
        A2 = 0.5 * (ada - dA12);

        kp = sqrt(-iom) * sqrt(A1);
        km = sqrt(-iom) * sqrt(A2);
    
        % xip, xim, Qp, and Qm for each layer are saved for later use. 
        xip(j) = -kp/iom;
        xim(j) = -km/iom;
    
        if A.xy(j) == 0
            QpNeg(j) = 0;
            QmNegRec(j) = 0;
        else
            QpNeg(j) = 2 * A.xy(j) / (add+dA12);          % -Qp
            QmNegRec(j) = 0.5 * (add-dA12) / A.xy(j);     % -1/Qm
        end
    
        QpDerQm = QpNeg(j) * QmNegRec(j);                 % Q_p/Q_m
        dq = 1-QpDerQm;                                   % 1-Q_p/Q_m
    
        % M of the basement
        if j==nLayer
            
            % 1st and 3rd columns won't be used, so set them to zeros
            M = [ 0                  1     0          QmNegRec(j); ...
                  0           QpNeg(j)     0                    1; ...
                  0   -xip(j)*QpNeg(j)     0              -xim(j); ...
                  0             xip(j)     0   xim(j)*QmNegRec(j) ];
            
            continue;
        end
        
        sp = sh(kp*h(j));
        sm = sh(km*h(j));
        cp = ch(kp*h(j));
        cm = ch(km*h(j));
        
        S(1,1) =  cp-QpDerQm*cm;
        S(1,2) = -QmNegRec(j)*(cp-cm);
        S(1,3) =  QmNegRec(j)*(sp/xip(j)-sm/xim(j));
        S(1,4) =  sp/xip(j)-QpDerQm*sm/xim(j);
        S(2,1) =  QpNeg(j)*(cp-cm);
        S(2,2) = -QpDerQm*cp+cm;
        S(2,3) =  QpDerQm*sp/xip(j)-sm/xim(j);
        S(2,4) =  QpNeg(j)*(sp/xip(j)-sm/xim(j));
        S(3,1) = -QpNeg(j)*(xip(j)*sp-xim(j)*sm);
        S(3,2) =  QpDerQm*xip(j)*sp-xim(j)*sm;
        S(3,3) = -QpDerQm*cp+cm;
        S(3,4) = -QpNeg(j)*(cp-cm);
        S(4,1) =  xip(j)*sp-QpDerQm*xim(j)*sm;
        S(4,2) = -QmNegRec(j)*(xip(j)*sp-xim(j)*sm);
        S(4,3) =  QmNegRec(j)*(cp-cm);
        S(4,4) =  cp-QpDerQm*cm;
        
        S = S/dq;
        
%         disp(['Layer#: ', num2str(j)]);
%         rcond(S)
        
        Sprod(:,:,j) = S * Sprod(:,:,j+1);
        
    end  % nLayer
    
    SM = Sprod(:,:,1) * M;
    
    % assume Ex0 & Ey0 are known, compute Cm and Dm of the bottom layer
    Ex0 = eTop(1);
    Ey0 = eTop(2);
    G = [SM(1,2)  SM(1,4); SM(2,2)  SM(2,4)];
    
    % CD = inv(G) * [Ex0; Ey0];
    CD = zeros(2,1);
    CD(1) = ( G(2,2)*Ex0 - G(1,2)*Ey0)/det(G);
    CD(2) = (-G(2,1)*Ex0 + G(1,1)*Ey0)/det(G);
    
%     % Alternatively, one can assume Hx0 & Hy0 are known
%     Hx0 = 1.;
%     Hy0 = 1.;
%     G = [SM(3,2)  SM(3,4); SM(4,2)  SM(4,4)];
%     CD = inv(G) * [Hx0; Hy0];
    
    % the full C of the bottom layer
    Cbot = [0; CD(1); 0; CD(2)];
    
    MC = M * Cbot;
    
    for j=1:nLayer
        F = Sprod(:,:,j) * MC;
        Ex(j, iFreq) = F(1);
        Ey(j, iFreq) = F(2);
        Hx(j, iFreq) = F(3);
        Hy(j, iFreq) = F(4);
        
        if j>1
            kp = -iom * xip(j-1);
            km = -iom * xim(j-1);
            QpDerQm = QpNeg(j-1)*QmNegRec(j-1);
            dq = 1-QpDerQm;
            dz = 0.5 * h(j-1);
        
            sp = sh(kp*dz);
            sm = sh(km*dz);
            cp = ch(kp*dz);
            cm = ch(km*dz);
        
			S(1,1) =  cp-QpDerQm*cm;
			S(1,2) = -QmNegRec(j-1)*(cp-cm);
			S(1,3) =  QmNegRec(j-1)*(sp/xip(j-1)-sm/xim(j-1));
			S(1,4) =  sp/xip(j-1)-QpDerQm*sm/xim(j-1);
			S(2,1) =  QpNeg(j-1)*(cp-cm);
			S(2,2) = -QpDerQm*cp+cm;
			S(2,3) =  QpDerQm*sp/xip(j-1)-sm/xim(j-1);
			S(2,4) =  QpNeg(j-1)*(sp/xip(j-1)-sm/xim(j-1));
			S(3,1) = -QpNeg(j-1)*(xip(j-1)*sp-xim(j-1)*sm);
			S(3,2) =  QpDerQm*xip(j-1)*sp-xim(j-1)*sm;
			S(3,3) = -QpDerQm*cp+cm;
			S(3,4) = -QpNeg(j-1)*(cp-cm);
			S(4,1) =  xip(j-1)*sp-QpDerQm*xim(j-1)*sm;
			S(4,2) = -QmNegRec(j-1)*(xip(j-1)*sp-xim(j-1)*sm);
			S(4,3) =  QmNegRec(j-1)*(cp-cm);
			S(4,4) =  cp-QpDerQm*cm;
        
            S = S/dq;

            % computes E-fields at layer center
            SF = S * F;
            Exc = SF(1);
            Eyc = SF(2);
            Ez(j-1, iFreq) = -(sigma(j-1, 5)*Exc + sigma(j-1, 6)*Eyc)/sigma(j-1, 3);
                 
        end
    end

end