def getImpedance(m1d,sigma,freq):

    # Initiate the impedances
    Z1d = np.empty(len(freq) , dtype='complex')
    h = m1d.hx   #vectorNx[:-1]
    # Start the process
    for nrFr, fr in enumerate(freq):
        om = 2*np.pi*fr
        Zall = np.empty(len(h)+1,dtype='complex')
        # Calculate the impedance for the bottom layer
        Zall[0] = (mu_0*om)/np.sqrt(mu_0*eps_0*(om)**2 - 1j*mu_0*sigma[0]*om)

        for nr,hi in enumerate(h):
            # Calculate the wave number
            # print nr,sigma[nr]
            k = np.sqrt(mu_0*eps_0*om**2 - 1j*mu_0*sigma[nr]*om)
            Z = (mu_0*om)/k

            Zall[nr+1] = Z *((Zall[nr] + Z*np.tanh(1j*k*hi))/(Z + Zall[nr]*np.tanh(1j*k*hi)))

        #pdb.set_trace()
        Z1d[nrFr] = Zall[-1]

    return Z1d
Analytic solution for MT 1D layered earth. Returns the impedance at the surface.
dataset/ETHPy150Open simpeg/simpeg/SimPEG/MT/Utils/MT1Danalytic.py/getImpedance