%% mt1DAnalyticFieldsEpsMu
% modified form EM3DANI/mt1DFwd.jl at master · CUG-EMI/EM3DANI · GitHub
% computes analytic fields for 1D layered model.
%     1) time dependence: e^{iwt};
%    eField = complex(sum(eLayer,1)');
%    hField = complex(sum(hLayer)') ;
%    time dependence:e^{iwt};
%    eField = complex(sum(eLayer,1).');
%    hField = complex(sum(hLayer).') ;
%     2) wavenumber: k^2 = -iwu/p + uew^2;
%     3) for Ex-Hy.
%  Assuming a halfspace with the same conductivity as the bottom layer.
%  Note: length(sigma) = length(zNode)-1
% Input:
%     freq    :: Real     - a single frequency value.
%     sigma   :: Vector   - conductivity of layered model.
%     zNode   :: Vector   - depth of top of each layer.
%     srcType :: String   - "E" or "H", determines the top E- or H-field is source.
%     fTop    :: Complex  - the given top boundary E-field or H-field value.
%     compH   :: Bool     - indicates whether compute the H-field or not.
% Output:
%      eField :: Vector   - electric fields at the top of each layer.
%      hField :: Vector   - magnetic fields at the top of each layer.
%


function  [eField, hField] = mt1DAnalyticFieldsEpsMu(freq,sigma,epsr,mur,zNode,fTop,srcType, compH)

% check if layer conductivity and depth array have the correct size
if size(sigma,1)~=length(zNode)-1
    error("The dimension of layer conductivity and depth mismatch!")
end

% physical constant
MU0   = 4 * pi * 1e-7 ;
EPS0  =  8.854187817620389 * 1e-12 ;
omega = 2 * pi * freq;
EPS   =  epsr*EPS0 ;
MU   = MU0 * mur;
omu   = omega .* [MU;MU(end)] ;

nLayer = length(zNode);
zLen   = diff(zNode);

% 1st, compute the impedance at the top layer using the well-known recurrence formula for impedances.

% wave number of all layers
% assume halfspace at the bottom of model: [sigma; sigma[end]]
% ka = sqrt.(MU0*EPS0*omega^2 - 1im*omu*sigma)

% ka = sqrt(-1i*omu.*[sigma; sigma(end)]) ;

ka = sqrt([MU;MU(end)].*[EPS;EPS(end)]*omega^2 -1i*omu.*[sigma; sigma(end)]) ;

% intrinsic impedance of the bottom layer
% ztmp = sqr[MU;MU(end)].*[EPS;EPS(end)]*omega^2 t(-1i*omu*(1/sigma(end)));
ztmp = omega * MU(end) / ka(end) ;

for j = nLayer-1:-1:1
    k = ka(j);
    zp = omega * MU(j) / k ;
    ztmp = zp * (ztmp + zp * tanh(k*zLen(j)*1i)) / (zp + ztmp * tanh(k*zLen(j)*1i)) ;
end
z1 = ztmp ;

% Up- and down-going electric fileds of each layer boundary
eLayer = zeros(2, nLayer);
hLayer = zeros(2, nLayer);


% Up- and down-going components of the top layer, respectively
k1 = ka(1)   ;
omu1 =omu(1) ;
if srcType == "E"          %E1 = c
    eLayer(1, 1) = 0.5 * fTop * (1 - omu1 / (z1*k1)) ;
    eLayer(2, 1) = 0.5 * fTop * (1 + omu1 / (z1*k1)) ;
elseif srcType == "H"     % H1 = c
    hLayer(1, 1) = 0.5 * fTop * (1 - z1*k1 / omu1) ;
    hLayer(2, 1) = 0.5 * fTop * (1 + z1*k1 / omu1) ;

    eLayer(1, 1) = -omu1/k1 * hLayer(1, 1) ;
    eLayer(2, 1) =  omu1/k1 * hLayer(2, 1) ;
end

% 2nd, propagate the EM fields from top to bottom
for j = 1:nLayer-1
    kr = ka(j) / ka(j+1) ;
    pInv = 0.5 * [1+kr 1-kr; 1-kr 1+kr] ;
    eUD = [exp(ka(j)*zLen(j)*1i) 0; 0 exp(-ka(j)*zLen(j)*1i)] ;
    eLayer(:, j+1) = pInv * eUD * eLayer(:, j) ;
    %check if overflow happens
    e2 = abs(eLayer(1, j+1) + eLayer(2, j+1)) ;
    e1 = abs(eLayer(1, j)   + eLayer(2, j) ) ;
    if e2-e1>0 || isnan(e2)
        eLayer(:, j+1:end) = 0.0 ;
        warning 'check if overflow happens'
        break;
    end
end

eField = complex(sum(eLayer,1)');

if compH
    hLayer = [ eLayer(1:1,:) * sparse(diag(-ka)) ./omu';...
        eLayer(2:2,:) * sparse(diag(ka))  ./omu' ] ;

    hField = complex(sum(hLayer)') ;
end

end


%