function [res, phz] = MT1d(f,thick,rho)
% computes the 1D MT apparent resistivity and phase
% usage:
%    [res, phz] = MT1d(f,thick,rho)
% f = vector of frequencies, Hz
% thick = vector of L-1 layer thicknesses
% rho = vector of L layer resistivities
% (L can be anything, - the code will find it from the length of rho)
%
% res = output vector of apparent resistivities
% phz = outupt vector of phases.

    mu = 1.256637062E-06;    
    w = 2*pi*f;
    i = sqrt(-1);
    ki = sqrt(-i*w.*mu./rho(end));    
    Z = w.*mu./ki;
    for j = length(rho)-1:-1:1
        ki = sqrt(-i*w*mu./rho(j));
        Zi = w.*mu./ki;    
        Z = Zi .*(Z + Zi.*tanh(i*ki*thick(j)) ) ./ ( Zi +  Z.*tanh(i*ki*thick(j)) );
    end
    res = Z.*conj(Z) ./ (mu*2*pi*f);
    phz = 180/pi*angle(Z);

