TDEMDipolarfields.py 文件源码

def E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, t, current=1., length=1., orientation='X', kappa=0., epsr=1.):

    """
        Computing the analytic electric fields (E) from an electrical dipole in a wholespace
        - You have the option of computing E for multiple times at a single reciever location
          or a single time at multiple locations

        :param numpy.array XYZ: reciever locations at which to evaluate E
        :param numpy.array srcLoc: [x,y,z] triplet defining the location of the electric dipole source
        :param float sig: value specifying the conductivity (S/m) of the wholespace
        :param numpy.array t: array of times at which to measure
        :param float current: size of the injected current (A), default is 1.0 A
        :param float length: length of the dipole (m), default is 1.0 m
        :param str orientation: orientation of dipole: 'X', 'Y', or 'Z'
        :param float kappa: magnetic susceptiblity value (unitless), default is 0.
        :param float epsr: relative permitivitty value (unitless),  default is 1.0
        :rtype: numpy.array
        :return: Ex, Ey, Ez: arrays containing all 3 components of E evaluated at the specified locations and times.
    """

    mu = mu_0*(1+kappa)
    epsilon = epsilon_0*epsr

    XYZ = Utils.asArray_N_x_Dim(XYZ, 3)
    # Check
    if XYZ.shape[0] > 1 & t.shape[0] > 1:
        raise Exception("I/O type error: For multiple field locations only a single time can be specified.")

    dx = XYZ[:,0]-srcLoc[0]
    dy = XYZ[:,1]-srcLoc[1]
    dz = XYZ[:,2]-srcLoc[2]

    r  = np.sqrt( dx**2. + dy**2. + dz**2.)
    theta = np.sqrt((mu*sig)/(4*t))

    front = current * length / (4.* pi * sig * r**3)
    mid   = 3 * erf(theta*r) - (4/np.sqrt(pi) * (theta)**3 * r**3 + 6/np.sqrt(pi) * theta * r) * np.exp(-(theta)**2 * (r)**2)
    extra = (erf(theta*r) - (4/np.sqrt(pi) * (theta)**3 * r**3 + 2/np.sqrt(pi) * theta * r) * np.exp(-(theta)**2 * (r)**2))

    if orientation.upper() == 'X':
        Ex = front*(dx**2 / r**2)*mid - front*extra
        Ey = front*(dx*dy  / r**2)*mid
        Ez = front*(dx*dz  / r**2)*mid
        return Ex, Ey, Ez

    elif orientation.upper() == 'Y':
        #  x--> y, y--> z, z-->x
        Ey = front*(dy**2 / r**2)*mid - front*extra
        Ez = front*(dy*dz  / r**2)*mid
        Ex = front*(dy*dx  / r**2)*mid
        return Ex, Ey, Ez

    elif orientation.upper() == 'Z':
        # x --> z, y --> x, z --> y
        Ez = front*(dz**2 / r**2)*mid - front*extra
        Ex = front*(dz*dx  / r**2)*mid
        Ey = front*(dz*dy  / r**2)*mid
        return Ex, Ey, Ez