diff --git a/simpegEM/Analytics/FDEM.py b/simpegEM/Analytics/FDEM.py index 3d02b589..36a2de63 100644 --- a/simpegEM/Analytics/FDEM.py +++ b/simpegEM/Analytics/FDEM.py @@ -111,3 +111,30 @@ def AnalyticMagDipoleWholeSpace(XYZ, srcLoc, sig, f, m=1., orientation='X'): return Bx, By, Bz + +def ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, m=1., orientation='X'): + XYZ = Utils.asArray_N_x_Dim(XYZ, 3) + + dx = XYZ[:,0]-srcLoc[0] + dy = XYZ[:,1]-srcLoc[1] + dz = XYZ[:,2]-srcLoc[2] + + r = np.sqrt( dx**2. + dy**2. + dz**2.) + k = np.sqrt( -1j*2.*np.pi*f*mu_0*sig ) + kr = k*r + + front = moment / (4. * np.pi * sig * r**3) * exp(-1j*k*r) + mid = -k**2 * r**2 + 3*1j*k*r + 3 + + Ex = front*((dx**2 / r**2)*mid + (k**2 * r**2 -1j*k*r)) + Ey = front*(dx*dy / r**2)*mid + Ez = front*(dx*dz / r**2)*mid + + if orientation.upper() == 'X': + return Ex, Ey, Ez + + elif orientation.upper() == 'Y': + return Ez, Ex, Ey + + elif orientation.upper() == 'Z': + return Ey, Ez, Ex \ No newline at end of file diff --git a/simpegEM/Analytics/FDEMcasing.py b/simpegEM/Analytics/FDEMcasing.py index 547d40ac..21f42cc3 100644 --- a/simpegEM/Analytics/FDEMcasing.py +++ b/simpegEM/Analytics/FDEMcasing.py @@ -10,20 +10,20 @@ def getKc(freq,sigma,a,b,mu=mu_0,eps=epsilon_0): def _r2(xyz): return np.sum(xyz**2,1) -def _getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu=mu_0,eps=epsilon_0,moment=1.): - Kc1 = getKc(freq,sigma[1],a,b,mu,eps) +def _getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1.): + Kc1 = getKc(freq,sigma[1],a,b,mu[1],eps) nobs = obsloc.shape[0] dxyz = obsloc - np.c_[np.ones(nobs)]*np.r_[srcloc] r2 = _r2(dxyz[:,:2]) sqrtr2z2 = np.sqrt(r2 + dxyz[:,2]**2) - k2 = k(freq,sigma[2],mu,eps) + k2 = k(freq,sigma[2],mu[2],eps) return Kc1 * moment / (4.*np.pi) *np.exp(-1j*k2*sqrtr2z2) / sqrtr2z2 -def _getCasingHertzMagDipoleDeriv_r(srcloc,obsloc,freq,sigma,a,b,mu=mu_0,eps=epsilon_0,moment=1.): +def _getCasingHertzMagDipoleDeriv_r(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1.): HertzZ = _getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment) nobs = obsloc.shape[0] @@ -31,12 +31,12 @@ def _getCasingHertzMagDipoleDeriv_r(srcloc,obsloc,freq,sigma,a,b,mu=mu_0,eps=eps r2 = _r2(dxyz[:,:2]) sqrtr2z2 = np.sqrt(r2 + dxyz[:,2]**2) - k2 = k(freq,sigma[2],mu,eps) + k2 = k(freq,sigma[2],mu[2],eps) return -HertzZ * np.sqrt(r2) / sqrtr2z2 * (1j*k2 + 1./ sqrtr2z2) -def _getCasingHertzMagDipoleDeriv_z(srcloc,obsloc,freq,sigma,a,b,mu=mu_0,eps=epsilon_0,moment=1.): +def _getCasingHertzMagDipoleDeriv_z(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1.): HertzZ = _getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment) nobs = obsloc.shape[0] @@ -44,11 +44,11 @@ def _getCasingHertzMagDipoleDeriv_z(srcloc,obsloc,freq,sigma,a,b,mu=mu_0,eps=eps r2z2 = _r2(dxyz) sqrtr2z2 = np.sqrt(r2z2) - k2 = k(freq,sigma[2],mu,eps) + k2 = k(freq,sigma[2],mu[2],eps) return -HertzZ*dxyz[:,2] /sqrtr2z2 * (1j*k2 + 1./sqrtr2z2) -def _getCasingHertzMagDipole2Deriv_z_r(srcloc,obsloc,freq,sigma,a,b,mu=mu_0,eps=epsilon_0,moment=1.): +def _getCasingHertzMagDipole2Deriv_z_r(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1.): HertzZ = _getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment) dHertzZdr = _getCasingHertzMagDipoleDeriv_r(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment) @@ -59,11 +59,11 @@ def _getCasingHertzMagDipole2Deriv_z_r(srcloc,obsloc,freq,sigma,a,b,mu=mu_0,eps= r = np.sqrt(r2) z = dxyz[:,2] sqrtr2z2 = np.sqrt(r2 + z**2) - k2 = k(freq,sigma[2],mu,eps) + k2 = k(freq,sigma[2],mu[2],eps) return dHertzZdr*(-z/sqrtr2z2)*(1j*k2+1./sqrtr2z2) + HertzZ*(z*r/sqrtr2z2**3)*(1j*k2 + 2./sqrtr2z2) -def _getCasingHertzMagDipole2Deriv_z_z(srcloc,obsloc,freq,sigma,a,b,mu=mu_0,eps=epsilon_0,moment=1.): +def _getCasingHertzMagDipole2Deriv_z_z(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1.): HertzZ = _getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment) dHertzZdz = _getCasingHertzMagDipoleDeriv_z(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment) @@ -74,24 +74,24 @@ def _getCasingHertzMagDipole2Deriv_z_z(srcloc,obsloc,freq,sigma,a,b,mu=mu_0,eps= r = np.sqrt(r2) z = dxyz[:,2] sqrtr2z2 = np.sqrt(r2 + z**2) - k2 = k(freq,sigma[2],mu,eps) + k2 = k(freq,sigma[2],mu[2],eps) return (dHertzZdz*z + HertzZ)/sqrtr2z2*(-1j*k2 - 1./sqrtr2z2) + HertzZ*z/sqrtr2z2**3*(1j*k2*z + 2.*z/sqrtr2z2) -def getCasingEphiMagDipole(srcloc,obsloc,freq,sigma,a,b,mu=mu_0,eps=epsilon_0,moment=1): +def getCasingEphiMagDipole(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1): return 1j * omega(freq) * mu * _getCasingHertzMagDipoleDeriv_r(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment) -def getCasingHrMagDipole(srcloc,obsloc,freq,sigma,a,b,mu=mu_0,eps=epsilon_0,moment=1): +def getCasingHrMagDipole(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1): return _getCasingHertzMagDipole2Deriv_z_r(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment) -def getCasingHzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu=mu_0,eps=epsilon_0,moment=1): +def getCasingHzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1): d2HertzZdz2 = _getCasingHertzMagDipole2Deriv_z_z(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment) - k2 = k(freq,sigma[2],mu,eps) + k2 = k(freq,sigma[2],mu[2],eps) HertzZ = _getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment) return d2HertzZdz2 + k2**2 * HertzZ -def getCasingBrMagDipole(srcloc,obsloc,freq,sigma,a,b,mu=mu_0,eps=epsilon_0,moment=1): +def getCasingBrMagDipole(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1): return mu_0 * getCasingHrMagDipole(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment) -def getCasingBzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu=mu_0,eps=epsilon_0,moment=1): +def getCasingBzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1): return mu_0 * getCasingHzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment) \ No newline at end of file