diff --git a/simpegEM/Analytics/FDEMcasing.py b/simpegEM/Analytics/FDEMcasing.py new file mode 100644 index 00000000..447f4986 --- /dev/null +++ b/simpegEM/Analytics/FDEMcasing.py @@ -0,0 +1,91 @@ +from SimPEG import Utils, np +from scipy.constants import mu_0, epsilon_0 +from simpegEM.Utils.EMUtils import k + +def getKc(freq,sigma,a,b,mu=mu_0,eps=epsilon_0): + a = float(a) + b = float(b) + return 2. * np.sqrt(b / a) * np.exp(-1j*k(freq,sigma,mu,eps)*(b-a)) + +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) + + 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) + + 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.): + HertzZ = _getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment) + + 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) + + 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.): + HertzZ = _getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment) + + nobs = obsloc.shape[0] + dxyz = obsloc - np.c_[np.ones(nobs)]*np.r_[srcloc] + + r2z2 = _r2(dxyz) + sqrtr2z2 = np.sqrt(r2z2) + k2 = k(freq,sigma[2],mu,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.): + HertzZ = _getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment) + dHertzZdr = _getCasingHertzMagDipoleDeriv_r(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment) + + nobs = obsloc.shape[0] + dxyz = obsloc - np.c_[np.ones(nobs)]*np.r_[srcloc] + + r2 = _r2(dxyz[:,:2]) + r = np.sqrt(r2) + z = dxyz[:,2] + sqrtr2z2 = np.sqrt(r2 + z**2) + k2 = k(freq,sigma[2],mu,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.): + HertzZ = _getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment) + dHertzZdz = _getCasingHertzMagDipoleDeriv_z(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment) + + nobs = obsloc.shape[0] + dxyz = obsloc - np.c_[np.ones(nobs)]*np.r_[srcloc] + + r2 = _r2(dxyz[:,:2]) + r = np.sqrt(r2) + z = dxyz[:,2] + sqrtr2z2 = np.sqrt(r2 + z**2) + k2 = k(freq,sigma[2],mu,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): + 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): + 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): + d2HertzZdz2 = _getCasingHertzMagDipole2Deriv_z_z(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment) + k2 = k(freq,sigma[2],mu,eps) + HertzZ(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment) + return d2HertzZdz2 + k2**2 * HertzZ \ No newline at end of file diff --git a/simpegEM/Analytics/__init__.py b/simpegEM/Analytics/__init__.py index 5828ba5c..5b7a8851 100644 --- a/simpegEM/Analytics/__init__.py +++ b/simpegEM/Analytics/__init__.py @@ -1,2 +1,3 @@ from TDEM import hzAnalyticDipoleT from FDEM import hzAnalyticDipoleF +from FDEMcasing import * diff --git a/simpegEM/Tests/test_FDEMCasing.py b/simpegEM/Tests/test_FDEMCasing.py new file mode 100644 index 00000000..21135a09 --- /dev/null +++ b/simpegEM/Tests/test_FDEMCasing.py @@ -0,0 +1,59 @@ +from SimPEG import Tests, Utils, np +import simpegEM.Analytics.FDEMcasing as Casing +import unittest + + +n = 50. +freq = 1. +a = 5e-2 +b = a + 1e-2 +sigma = np.r_[1., 5.5e6, 1e-1] +srcloc = np.r_[0., 0., 0.] +xobs = np.random.rand(n)+10. +yobs = np.zeros(n) +zobs = np.random.randn(n) +plotit = False + +def CasingMagDipoleDeriv_r(x): + obsloc = np.vstack([x, yobs, zobs]).T + + f = Casing._getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b) + g = Utils.sdiag(Casing._getCasingHertzMagDipoleDeriv_r(srcloc,obsloc,freq,sigma,a,b)) + + return f,g + +def CasingMagDipoleDeriv_z(z): + obsloc = np.vstack([xobs, yobs, z]).T + + f = Casing._getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b) + g = Utils.sdiag(Casing._getCasingHertzMagDipoleDeriv_z(srcloc,obsloc,freq,sigma,a,b)) + + return f,g + +def CasingMagDipole2Deriv_z_r(x): + obsloc = np.vstack([x, yobs, zobs]).T + + f = Casing._getCasingHertzMagDipoleDeriv_z(srcloc,obsloc,freq,sigma,a,b) + g = Utils.sdiag(Casing._getCasingHertzMagDipole2Deriv_z_r(srcloc,obsloc,freq,sigma,a,b)) + + return f,g + +def CasingMagDipole2Deriv_z_z(z): + obsloc = np.vstack([xobs, yobs, z]).T + + f = Casing._getCasingHertzMagDipoleDeriv_z(srcloc,obsloc,freq,sigma,a,b) + g = Utils.sdiag(Casing._getCasingHertzMagDipole2Deriv_z_z(srcloc,obsloc,freq,sigma,a,b)) + + return f,g + + + +class Casing_DerivTest(unittest.TestCase): + Tests.checkDerivative(CasingMagDipoleDeriv_r,np.ones(n)*10+np.random.randn(n),plotIt=False) + Tests.checkDerivative(CasingMagDipoleDeriv_z,np.random.randn(n),plotIt=False) + Tests.checkDerivative(CasingMagDipole2Deriv_z_r,np.ones(n)*10+np.random.randn(n),plotIt=False) + Tests.checkDerivative(CasingMagDipole2Deriv_z_z,np.random.randn(n),plotIt=False) + + +if __name__ == '__main__': + unittest.main() \ No newline at end of file diff --git a/simpegEM/Tests/test_FieldsObject.py b/simpegEM/Tests/test_FieldsObject.py index 251924b4..4d9cb9ef 100644 --- a/simpegEM/Tests/test_FieldsObject.py +++ b/simpegEM/Tests/test_FieldsObject.py @@ -11,7 +11,7 @@ class FieldsTest(unittest.TestCase): srcLoc = np.r_[0,0,0.] rxList0 = EM.FDEM.RxFDEM(XYZ, 'exi') Src0 = EM.FDEM.SrcFDEM_MagDipole(srcLoc, 3., [rxList0]) - rxList1 = EM.FDEM.RxFDEM(XYZ, 'bxi') + rxList1 = EM.FDEM.RxFDEM(XYZ, 'bxi') Src1 = EM.FDEM.SrcFDEM_MagDipole(srcLoc, 3., [rxList1]) rxList2 = EM.FDEM.RxFDEM(XYZ, 'bxi') Src2 = EM.FDEM.SrcFDEM_MagDipole(srcLoc, 2., [rxList2])