diff --git a/MT1D/mt1DAniAnalyticFieldEspMu.m b/MT1D/mt1DAniAnalyticFieldEspMu.m
new file mode 100644
index 0000000..cca3ded
--- /dev/null
+++ b/MT1D/mt1DAniAnalyticFieldEspMu.m
@@ -0,0 +1,202 @@
+function [Ex, Ey, Ez, Hx, Hy]=mt1DAniAnalyticFieldEspMu(freqs,sigma,epsr,mur,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);
+MU0 = 4*pi*1e-7;
+MU  = MU0 * mur;
+EPS0  =  8.854187817620389 * 1e-12 ;
+EPS  =  epsr*EPS0;
+    % assume halfspace at the bottom of model
+
+sigma    = [sigma; sigma(end:end, :)];
+EPS     = [ EPS;  EPS(end:end, :)];
+MU       = [ MU;  MU(end:end, :)];
+    % effective azimuthal anisotropic conductivity
+A = getEffSigma(sigma,EPS,freqs);% get Axx, Axy, Ayx
+nFreq  = length(freqs);
+nLayer = size(sigma, 1);
+
+Ex = zeros(nLayer, nFreq);
+Ey = zeros(nLayer, nFreq);
+% Ez = zeros(nLayer-1, 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()
+
+
+
+for iFreq=1:nFreq
+    omega = 2 * pi * freqs(iFreq);
+    iom = 1i*omega*MU;
+    
+    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(j)) * sqrt(A1);
+        km = sqrt(-iom(j)) * sqrt(A2);
+    
+        % xip, xim, Qp, and Qm for each layer are saved for later use. 
+        xip(j) = -kp/iom(j);
+        xim(j) = -km/iom(j);
+    
+        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(j-1) * xip(j-1);
+            km = -iom(j-1) * 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
+     Ez(j, iFreq) =   Ez(j-1);
+
+end
+
+
+
+