From 28d67e311205ed1ba12b24f06a891df8180eb514 Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Fri, 27 May 2016 11:26:47 -0700 Subject: [PATCH 01/11] Start of ED !! --- SimPEG/EM/Analytics/FDEM_fields.py | 51 ++++++++++++++++++++++++++++++ 1 file changed, 51 insertions(+) create mode 100644 SimPEG/EM/Analytics/FDEM_fields.py diff --git a/SimPEG/EM/Analytics/FDEM_fields.py b/SimPEG/EM/Analytics/FDEM_fields.py new file mode 100644 index 00000000..adf10e4e --- /dev/null +++ b/SimPEG/EM/Analytics/FDEM_fields.py @@ -0,0 +1,51 @@ +from __future__ import division +import numpy as np +from scipy.constants import mu_0, pi +from scipy.special import erf +from SimPEG import Utils + + +def ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): + 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*sig ) + kr = k*r + + front = current * length / (4. * np.pi * sig * r**3) * np.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': + Ex = front*((dx**2 / r**2)*mid + (k**2 * r**2 -1j*k*r-1.)) + 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 + (k**2 * r**2 -1j*k*r-1.)) + 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 + (k**2 * r**2 -1j*k*r-1.)) + Ex = front*(dz*dx / r**2)*mid + Ey = front*(dz*dy / r**2)*mid + return Ex, Ey, Ez + # return Ey, Ez, Ex + +def A_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): +def J_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): +def H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): +def B_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): + mu*H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0) From 1960b52dfd58a9fc81c225b5ad005d98ee9ca74a Mon Sep 17 00:00:00 2001 From: micmitch Date: Fri, 27 May 2016 14:56:48 -0700 Subject: [PATCH 02/11] First stab at analytic functions for the fields from a harmonic electric dipole source. Not sure about the exception that I try to throw if multiple frequencies and multiple evaluation locations are both specified. --- SimPEG/EM/Analytics/FDEM_fields.py | 119 +++++++++++++++++++++++++---- 1 file changed, 106 insertions(+), 13 deletions(-) diff --git a/SimPEG/EM/Analytics/FDEM_fields.py b/SimPEG/EM/Analytics/FDEM_fields.py index adf10e4e..abf956e8 100644 --- a/SimPEG/EM/Analytics/FDEM_fields.py +++ b/SimPEG/EM/Analytics/FDEM_fields.py @@ -5,24 +5,27 @@ from scipy.special import erf from SimPEG import Utils -def ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): +def E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): + epsilon = 8.854187817*(10.**-12) + omega = 2.*np.pi*f + XYZ = Utils.asArray_N_x_Dim(XYZ, 3) + # Check + if XYZ.shape[0] > 1 & f.length > 1 + except: + print "I/O type error: For multiple field locations only a single frequency 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.) - k = np.sqrt( -1j*2.*np.pi*f*mu*sig ) - kr = k*r + # k = np.sqrt( -1j*2.*np.pi*f*mu*sig ) + k = np.sqrt( omega**2. *mu*epsilon -1j*omega*mu*sig ) - front = current * length / (4. * np.pi * sig * r**3) * np.exp(-1j*k*r) + front = current * length / (4.*np.pi*sig* r**3) * np.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': Ex = front*((dx**2 / r**2)*mid + (k**2 * r**2 -1j*k*r-1.)) Ey = front*(dx*dy / r**2)*mid @@ -42,10 +45,100 @@ def ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orienta Ex = front*(dz*dx / r**2)*mid Ey = front*(dz*dy / r**2)*mid return Ex, Ey, Ez - # return Ey, Ez, Ex + +def J_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): + Ex, Ey, Ez = E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0) + Jx = sig*Ex + Jy = sig*Ey + Jz = sig*Ez + return Jx, Jy, Jz + + +def H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): + epsilon = 8.854187817*(10.**-12) + omega = 2.*np.pi*f + + XYZ = Utils.asArray_N_x_Dim(XYZ, 3) + # Check + if XYZ.shape[0] > 1 & f.length > 1 + except: + print "I/O type error: For multiple field locations only a single frequency 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.) + # k = np.sqrt( -1j*2.*np.pi*f*mu*sig ) + k = np.sqrt( omega**2. *mu*epsilon -1j*omega*mu*sig ) + + front = current * length / (4.*np.pi* r**2) * (-1j*k*r + 1) * np.exp(-1j*k*r) + + if orientation.upper() == 'X': + Hy = front*(-dz / r) + Hz = front*(dy / r) + Hx = np.zeros_like(Hy) + return Hx, Hy, Hz + + elif orientation.upper() == 'Y': + Hx = front*(dz / r) + Hz = front*(-dx / r) + Hy = np.zeros_like(Hx) + return Hx, Hy, Hz + + elif orientation.upper() == 'Z': + Hx = front*(-dy / r) + Hy = front*(dx / r) + Hz = np.zeros_like(Hx) + return Hx, Hy, Hz + + +def B_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): + Hx, Hy, Hz = H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0) + Bx = mu*Hx + By = mu*Hy + Bz = mu*Hz + return Bx, By, Bz + def A_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): -def J_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): -def H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): -def B_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): - mu*H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0) + epsilon = 8.854187817*(10.**-12) + omega = 2.*np.pi*f + + XYZ = Utils.asArray_N_x_Dim(XYZ, 3) + # Check + if XYZ.shape[0] > 1 & f.length > 1 + except: + print "I/O type error: For multiple field locations only a single frequency 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.) + k = np.sqrt( omega**2. *mu*epsilon -1j*omega*mu*sig ) + + front = current * length / (4.*np.pi*r) + + if orientation.upper() == 'X': + Ax = front*np.exp(-1j*k*r) + Ay = np.zeros_like(Ax) + Az = np.zeros_like(Ax) + return Ax, Ay, Az + + elif orientation.upper() == 'Y': + Ay = front*np.exp(-1j*k*r) + Ax = np.zeros_like(Ay) + Az = np.zeros_like(Ay) + return Ax, Ay, Az + + elif orientation.upper() == 'Z': + Az = front*np.exp(-1j*k*r) + Ax = np.zeros_like(Ay) + Ay = np.zeros_like(Ay) + return Ax, Ay, Az + + + + + From a3a5c860088bfb58d388c47b6d2d26c933573e49 Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Sat, 11 Jun 2016 21:58:03 +0200 Subject: [PATCH 03/11] Fix couple bugs in FDEM analytics --- SimPEG/EM/Analytics/FDEM_fields.py | 21 +++++++++------------ SimPEG/EM/Analytics/__init__.py | 1 + 2 files changed, 10 insertions(+), 12 deletions(-) diff --git a/SimPEG/EM/Analytics/FDEM_fields.py b/SimPEG/EM/Analytics/FDEM_fields.py index abf956e8..e0ac1b34 100644 --- a/SimPEG/EM/Analytics/FDEM_fields.py +++ b/SimPEG/EM/Analytics/FDEM_fields.py @@ -11,9 +11,8 @@ def E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., XYZ = Utils.asArray_N_x_Dim(XYZ, 3) # Check - if XYZ.shape[0] > 1 & f.length > 1 - except: - print "I/O type error: For multiple field locations only a single frequency can be specified." + if XYZ.shape[0] > 1 & f.shape[0] > 1: + raise Exception("I/O type error: For multiple field locations only a single frequency can be specified.") dx = XYZ[:,0]-srcLoc[0] dy = XYZ[:,1]-srcLoc[1] @@ -51,7 +50,7 @@ def J_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., Jx = sig*Ex Jy = sig*Ey Jz = sig*Ez - return Jx, Jy, Jz + return Jx, Jy, Jz def H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): @@ -60,10 +59,9 @@ def H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., XYZ = Utils.asArray_N_x_Dim(XYZ, 3) # Check - if XYZ.shape[0] > 1 & f.length > 1 - except: - print "I/O type error: For multiple field locations only a single frequency can be specified." - + if XYZ.shape[0] > 1 & f.shape[0] > 1: + raise Exception("I/O type error: For multiple field locations only a single frequency can be specified.") + dx = XYZ[:,0]-srcLoc[0] dy = XYZ[:,1]-srcLoc[1] dz = XYZ[:,2]-srcLoc[2] @@ -99,7 +97,7 @@ def B_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., By = mu*Hy Bz = mu*Hz return Bx, By, Bz - + def A_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): epsilon = 8.854187817*(10.**-12) @@ -107,9 +105,8 @@ def A_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., XYZ = Utils.asArray_N_x_Dim(XYZ, 3) # Check - if XYZ.shape[0] > 1 & f.length > 1 - except: - print "I/O type error: For multiple field locations only a single frequency can be specified." + if XYZ.shape[0] > 1 & f.shape[0] > 1: + raise Exception("I/O type error: For multiple field locations only a single frequency can be specified.") dx = XYZ[:,0]-srcLoc[0] dy = XYZ[:,1]-srcLoc[1] diff --git a/SimPEG/EM/Analytics/__init__.py b/SimPEG/EM/Analytics/__init__.py index 9df2aef7..9c5e5e28 100644 --- a/SimPEG/EM/Analytics/__init__.py +++ b/SimPEG/EM/Analytics/__init__.py @@ -2,3 +2,4 @@ from TDEM import hzAnalyticDipoleT from FDEM import hzAnalyticDipoleF from FDEMcasing import * from DC import DCAnalyticHalf, DCAnalyticSphere +from FDEM_fields import * From 5e3c1da8e2626db11b14578b489faa17e0bad90a Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Mon, 13 Jun 2016 16:41:59 -0700 Subject: [PATCH 04/11] add place holder for galvanic and inductive electric fields... --- SimPEG/EM/Analytics/FDEM_fields.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/SimPEG/EM/Analytics/FDEM_fields.py b/SimPEG/EM/Analytics/FDEM_fields.py index e0ac1b34..4a521ba1 100644 --- a/SimPEG/EM/Analytics/FDEM_fields.py +++ b/SimPEG/EM/Analytics/FDEM_fields.py @@ -5,6 +5,9 @@ from scipy.special import erf from SimPEG import Utils +# def E_galvanic_from_ElectricDipoleWholeSpaced +# def E_inductive_from_ElectricDipoleWholeSpaced + def E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): epsilon = 8.854187817*(10.**-12) omega = 2.*np.pi*f @@ -45,6 +48,10 @@ def E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., Ey = front*(dz*dy / r**2)*mid return Ex, Ey, Ez + +# def J_galvanic_from_ElectricDipoleWholeSpaced +# def J_inductive_from_ElectricDipoleWholeSpaced + def J_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): Ex, Ey, Ez = E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0) Jx = sig*Ex From 093f441331c034e8dfb27ce91d4495c13e690efb Mon Sep 17 00:00:00 2001 From: micmitch Date: Mon, 13 Jun 2016 17:35:52 -0700 Subject: [PATCH 05/11] Added functions to split electric field into "galvanic" and "inductive" portions. --- SimPEG/EM/Analytics/FDEM_fields.py | 102 +++++++++++++++++++++++++++-- 1 file changed, 96 insertions(+), 6 deletions(-) diff --git a/SimPEG/EM/Analytics/FDEM_fields.py b/SimPEG/EM/Analytics/FDEM_fields.py index 4a521ba1..ac18da7f 100644 --- a/SimPEG/EM/Analytics/FDEM_fields.py +++ b/SimPEG/EM/Analytics/FDEM_fields.py @@ -4,10 +4,6 @@ from scipy.constants import mu_0, pi from scipy.special import erf from SimPEG import Utils - -# def E_galvanic_from_ElectricDipoleWholeSpaced -# def E_inductive_from_ElectricDipoleWholeSpaced - def E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): epsilon = 8.854187817*(10.**-12) omega = 2.*np.pi*f @@ -49,8 +45,86 @@ def E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., return Ex, Ey, Ez -# def J_galvanic_from_ElectricDipoleWholeSpaced -# def J_inductive_from_ElectricDipoleWholeSpaced +def E_galvanic_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): + epsilon = 8.854187817*(10.**-12) + omega = 2.*np.pi*f + + XYZ = Utils.asArray_N_x_Dim(XYZ, 3) + # Check + if XYZ.shape[0] > 1 & f.shape[0] > 1: + raise Exception("I/O type error: For multiple field locations only a single frequency 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.) + # k = np.sqrt( -1j*2.*np.pi*f*mu*sig ) + k = np.sqrt( omega**2. *mu*epsilon -1j*omega*mu*sig ) + + front = current * length / (4.*np.pi*sig* r**3) * np.exp(-1j*k*r) + mid = -k**2 * r**2 + 3*1j*k*r + 3 + + if orientation.upper() == 'X': + Ex_galvanic = front*((dx**2 / r**2)*mid + (-1j*k*r-1.)) + Ey_galvanic = front*(dx*dy / r**2)*mid + Ez_galvanic = front*(dx*dz / r**2)*mid + return Ex_galvanic, Ey_galvanic, Ez_galvanic + + elif orientation.upper() == 'Y': + # x--> y, y--> z, z-->x + Ey_galvanic = front*((dy**2 / r**2)*mid + (-1j*k*r-1.)) + Ez_galvanic = front*(dy*dz / r**2)*mid + Ex_galvanic = front*(dy*dx / r**2)*mid + return Ex_galvanic, Ey_galvanic, Ez_galvanic + + elif orientation.upper() == 'Z': + # x --> z, y --> x, z --> y + Ez_galvanic = front*((dz**2 / r**2)*mid + (-1j*k*r-1.)) + Ex_galvanic = front*(dz*dx / r**2)*mid + Ey_galvanic = front*(dz*dy / r**2)*mid + return Ex_galvanic, Ey_galvanic, Ez_galvanic + + +def E_inductive_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): + epsilon = 8.854187817*(10.**-12) + omega = 2.*np.pi*f + + XYZ = Utils.asArray_N_x_Dim(XYZ, 3) + # Check + if XYZ.shape[0] > 1 & f.shape[0] > 1: + raise Exception("I/O type error: For multiple field locations only a single frequency 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.) + # k = np.sqrt( -1j*2.*np.pi*f*mu*sig ) + k = np.sqrt( omega**2. *mu*epsilon -1j*omega*mu*sig ) + + front = current * length / (4.*np.pi*sig* r**3) * np.exp(-1j*k*r) + + if orientation.upper() == 'X': + Ex_inductive = front*(k**2 * r**2) + Ey_inductive = 0 + Ez_inductive = 0 + return Ex_inductive, Ey_inductive, Ez_inductive + + elif orientation.upper() == 'Y': + # x--> y, y--> z, z-->x + Ey_inductive = front*(k**2 * r**2) + Ez_inductive = 0 + Ex_inductive = 0 + return Ex_inductive, Ey_inductive, Ez_inductive + + elif orientation.upper() == 'Z': + # x --> z, y --> x, z --> y + Ez_inductive = front*(k**2 * r**2) + Ex_inductive = 0 + Ey_inductive = 0 + return Ex_inductive, Ey_inductive, Ez_inductive + def J_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): Ex, Ey, Ez = E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0) @@ -60,6 +134,22 @@ def J_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., return Jx, Jy, Jz +def J_galvanic_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): + Ex_galvanic, Ey_galvanic, Ez_galvanic = E_galvanic_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0) + Jx_galvanic = sig*Ex_galvanic + Jy_galvanic = sig*Ey_galvanic + Jz_galvanic = sig*Ez_galvanic + return Jx_galvanic, Jy_galvanic, Jz_galvanic + + +def J_inductive_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): + Ex_inductive, Ey_inductive, Ez_inductive = E_inductive_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0) + Jx_inductive = sig*Ex_inductive + Jy_inductive = sig*Ey_inductive + Jz_inductive = sig*Ez_inductive + return Jx_inductive, Jy_inductive, Jz_inductive + + def H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): epsilon = 8.854187817*(10.**-12) omega = 2.*np.pi*f From b9d30af4a88903da8941e4afa6355cece6c939eb Mon Sep 17 00:00:00 2001 From: micmitch Date: Mon, 13 Jun 2016 17:47:46 -0700 Subject: [PATCH 06/11] Changed \sigma to \hat{\sigma} = \sigma + i \omega \epsilon in the E field calculations. --- SimPEG/EM/Analytics/FDEM_fields.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/SimPEG/EM/Analytics/FDEM_fields.py b/SimPEG/EM/Analytics/FDEM_fields.py index ac18da7f..467407af 100644 --- a/SimPEG/EM/Analytics/FDEM_fields.py +++ b/SimPEG/EM/Analytics/FDEM_fields.py @@ -7,6 +7,7 @@ from SimPEG import Utils def E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): epsilon = 8.854187817*(10.**-12) omega = 2.*np.pi*f + sig_hat = sig + 1j*omega*epsilon XYZ = Utils.asArray_N_x_Dim(XYZ, 3) # Check @@ -21,7 +22,7 @@ def E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., # k = np.sqrt( -1j*2.*np.pi*f*mu*sig ) k = np.sqrt( omega**2. *mu*epsilon -1j*omega*mu*sig ) - front = current * length / (4.*np.pi*sig* r**3) * np.exp(-1j*k*r) + front = current * length / (4.*np.pi*sig_hat* r**3) * np.exp(-1j*k*r) mid = -k**2 * r**2 + 3*1j*k*r + 3 if orientation.upper() == 'X': @@ -48,6 +49,7 @@ def E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., def E_galvanic_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): epsilon = 8.854187817*(10.**-12) omega = 2.*np.pi*f + sig_hat = sig + 1j*omega*epsilon XYZ = Utils.asArray_N_x_Dim(XYZ, 3) # Check @@ -62,7 +64,7 @@ def E_galvanic_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., l # k = np.sqrt( -1j*2.*np.pi*f*mu*sig ) k = np.sqrt( omega**2. *mu*epsilon -1j*omega*mu*sig ) - front = current * length / (4.*np.pi*sig* r**3) * np.exp(-1j*k*r) + front = current * length / (4.*np.pi*sig_hat* r**3) * np.exp(-1j*k*r) mid = -k**2 * r**2 + 3*1j*k*r + 3 if orientation.upper() == 'X': @@ -89,6 +91,7 @@ def E_galvanic_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., l def E_inductive_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): epsilon = 8.854187817*(10.**-12) omega = 2.*np.pi*f + sig_hat = sig + 1j*omega*epsilon XYZ = Utils.asArray_N_x_Dim(XYZ, 3) # Check @@ -103,7 +106,7 @@ def E_inductive_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., # k = np.sqrt( -1j*2.*np.pi*f*mu*sig ) k = np.sqrt( omega**2. *mu*epsilon -1j*omega*mu*sig ) - front = current * length / (4.*np.pi*sig* r**3) * np.exp(-1j*k*r) + front = current * length / (4.*np.pi*sig_hat* r**3) * np.exp(-1j*k*r) if orientation.upper() == 'X': Ex_inductive = front*(k**2 * r**2) From e815ddaec79d59e142a9a3899d0fc1bac7434340 Mon Sep 17 00:00:00 2001 From: micmitch Date: Tue, 14 Jun 2016 15:15:53 -0700 Subject: [PATCH 07/11] Removed d typos from the end of function names. --- SimPEG/EM/Analytics/FDEM_fields.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/SimPEG/EM/Analytics/FDEM_fields.py b/SimPEG/EM/Analytics/FDEM_fields.py index 467407af..bd1ee368 100644 --- a/SimPEG/EM/Analytics/FDEM_fields.py +++ b/SimPEG/EM/Analytics/FDEM_fields.py @@ -46,7 +46,7 @@ def E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., return Ex, Ey, Ez -def E_galvanic_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): +def E_galvanic_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): epsilon = 8.854187817*(10.**-12) omega = 2.*np.pi*f sig_hat = sig + 1j*omega*epsilon @@ -88,7 +88,7 @@ def E_galvanic_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., l return Ex_galvanic, Ey_galvanic, Ez_galvanic -def E_inductive_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): +def E_inductive_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): epsilon = 8.854187817*(10.**-12) omega = 2.*np.pi*f sig_hat = sig + 1j*omega*epsilon @@ -137,7 +137,7 @@ def J_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., return Jx, Jy, Jz -def J_galvanic_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): +def J_galvanic_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): Ex_galvanic, Ey_galvanic, Ez_galvanic = E_galvanic_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0) Jx_galvanic = sig*Ex_galvanic Jy_galvanic = sig*Ey_galvanic @@ -145,7 +145,7 @@ def J_galvanic_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., l return Jx_galvanic, Jy_galvanic, Jz_galvanic -def J_inductive_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): +def J_inductive_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): Ex_inductive, Ey_inductive, Ez_inductive = E_inductive_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0) Jx_inductive = sig*Ex_inductive Jy_inductive = sig*Ey_inductive From f0944362c8afded6dc263629681775705fb73256 Mon Sep 17 00:00:00 2001 From: micmitch Date: Tue, 14 Jun 2016 15:22:36 -0700 Subject: [PATCH 08/11] Silly mistake... needed zero arrays instead of scalars. --- SimPEG/EM/Analytics/FDEM_fields.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/SimPEG/EM/Analytics/FDEM_fields.py b/SimPEG/EM/Analytics/FDEM_fields.py index bd1ee368..ca250ad6 100644 --- a/SimPEG/EM/Analytics/FDEM_fields.py +++ b/SimPEG/EM/Analytics/FDEM_fields.py @@ -110,22 +110,22 @@ def E_inductive_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., l if orientation.upper() == 'X': Ex_inductive = front*(k**2 * r**2) - Ey_inductive = 0 - Ez_inductive = 0 + Ey_inductive = np.zeros_like(Ex_inductive) + Ez_inductive = np.zeros_like(Ex_inductive) return Ex_inductive, Ey_inductive, Ez_inductive elif orientation.upper() == 'Y': # x--> y, y--> z, z-->x Ey_inductive = front*(k**2 * r**2) - Ez_inductive = 0 - Ex_inductive = 0 + Ez_inductive = np.zeros_like(Ey_inductive) + Ex_inductive = np.zeros_like(Ey_inductive) return Ex_inductive, Ey_inductive, Ez_inductive elif orientation.upper() == 'Z': # x --> z, y --> x, z --> y Ez_inductive = front*(k**2 * r**2) - Ex_inductive = 0 - Ey_inductive = 0 + Ex_inductive = np.zeros_like(Ez_inductive) + Ey_inductive = np.zeros_like(Ez_inductive) return Ex_inductive, Ey_inductive, Ez_inductive From e1ba80883d1dcaaac8f7804812ec216b8b974fea Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Thu, 23 Jun 2016 09:10:50 -0700 Subject: [PATCH 09/11] Incorporate Lindsey's suggestoins --- .../{FDEM_fields.py => FDEMDipolarfields.py} | 125 +++++++++++++----- 1 file changed, 93 insertions(+), 32 deletions(-) rename SimPEG/EM/Analytics/{FDEM_fields.py => FDEMDipolarfields.py} (69%) diff --git a/SimPEG/EM/Analytics/FDEM_fields.py b/SimPEG/EM/Analytics/FDEMDipolarfields.py similarity index 69% rename from SimPEG/EM/Analytics/FDEM_fields.py rename to SimPEG/EM/Analytics/FDEMDipolarfields.py index ca250ad6..ef724542 100644 --- a/SimPEG/EM/Analytics/FDEM_fields.py +++ b/SimPEG/EM/Analytics/FDEMDipolarfields.py @@ -1,13 +1,24 @@ from __future__ import division import numpy as np -from scipy.constants import mu_0, pi +from scipy.constants import mu_0, pi, epsilon_0 from scipy.special import erf from SimPEG import Utils -def E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): - epsilon = 8.854187817*(10.**-12) - omega = 2.*np.pi*f - sig_hat = sig + 1j*omega*epsilon +omega = lambda f: 2.*np.pi*f +# TODO: +# r = lambda dx, dy, dz: np.sqrt( dx**2. + dy**2. + dz**2.) +# k = lambda f, mu, epsilon, sig: np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig ) + +def E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=0., epsr=1.): + + """ + Computing Analytic Electric fields from Electrical Dipole in a Wholespace + TODO: + Add description of parameters + """ + mu = mu_0*(1+kappa) + epsilon = epsilon_0*epsr + sig_hat = sig + 1j*omeg*epsilon XYZ = Utils.asArray_N_x_Dim(XYZ, 3) # Check @@ -20,7 +31,7 @@ def E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., r = np.sqrt( dx**2. + dy**2. + dz**2.) # k = np.sqrt( -1j*2.*np.pi*f*mu*sig ) - k = np.sqrt( omega**2. *mu*epsilon -1j*omega*mu*sig ) + k = np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig ) front = current * length / (4.*np.pi*sig_hat* r**3) * np.exp(-1j*k*r) mid = -k**2 * r**2 + 3*1j*k*r + 3 @@ -46,10 +57,16 @@ def E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., return Ex, Ey, Ez -def E_galvanic_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): - epsilon = 8.854187817*(10.**-12) - omega = 2.*np.pi*f - sig_hat = sig + 1j*omega*epsilon +def E_galvanic_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.): + + """ + Computing Galvanic portion of Electric fields from Electrical Dipole in a Wholespace + TODO: + Add description of parameters + """ + mu = mu_0*(1+kappa) + epsilon = epsilon_0*epsr + sig_hat = sig + 1j*omeg*epsilon XYZ = Utils.asArray_N_x_Dim(XYZ, 3) # Check @@ -62,7 +79,7 @@ def E_galvanic_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., le r = np.sqrt( dx**2. + dy**2. + dz**2.) # k = np.sqrt( -1j*2.*np.pi*f*mu*sig ) - k = np.sqrt( omega**2. *mu*epsilon -1j*omega*mu*sig ) + k = np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig ) front = current * length / (4.*np.pi*sig_hat* r**3) * np.exp(-1j*k*r) mid = -k**2 * r**2 + 3*1j*k*r + 3 @@ -88,10 +105,16 @@ def E_galvanic_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., le return Ex_galvanic, Ey_galvanic, Ez_galvanic -def E_inductive_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): - epsilon = 8.854187817*(10.**-12) - omega = 2.*np.pi*f - sig_hat = sig + 1j*omega*epsilon +def E_inductive_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.): + + """ + Computing Inductive portion of Electric fields from Electrical Dipole in a Wholespace + TODO: + Add description of parameters + """ + mu = mu_0*(1+kappa) + epsilon = epsilon_0*epsr + sig_hat = sig + 1j*omeg*epsilon XYZ = Utils.asArray_N_x_Dim(XYZ, 3) # Check @@ -104,7 +127,7 @@ def E_inductive_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., l r = np.sqrt( dx**2. + dy**2. + dz**2.) # k = np.sqrt( -1j*2.*np.pi*f*mu*sig ) - k = np.sqrt( omega**2. *mu*epsilon -1j*omega*mu*sig ) + k = np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig ) front = current * length / (4.*np.pi*sig_hat* r**3) * np.exp(-1j*k*r) @@ -129,34 +152,60 @@ def E_inductive_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., l return Ex_inductive, Ey_inductive, Ez_inductive -def J_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): - Ex, Ey, Ez = E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0) +def J_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.): + + """ + Computing Current densities from Electrical Dipole in a Wholespace + TODO: + Add description of parameters + """ + + Ex, Ey, Ez = E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.) Jx = sig*Ex Jy = sig*Ey Jz = sig*Ez return Jx, Jy, Jz -def J_galvanic_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): - Ex_galvanic, Ey_galvanic, Ez_galvanic = E_galvanic_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0) +def J_galvanic_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.): + + """ + Computing Galvanic portion of Current densities from Electrical Dipole in a Wholespace + TODO: + Add description of parameters + """ + + Ex_galvanic, Ey_galvanic, Ez_galvanic = E_galvanic_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.) Jx_galvanic = sig*Ex_galvanic Jy_galvanic = sig*Ey_galvanic Jz_galvanic = sig*Ez_galvanic return Jx_galvanic, Jy_galvanic, Jz_galvanic -def J_inductive_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): - Ex_inductive, Ey_inductive, Ez_inductive = E_inductive_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0) +def J_inductive_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.): + + """ + Computing Inductive portion of Current densities from Electrical Dipole in a Wholespace + TODO: + Add description of parameters + """ + + Ex_inductive, Ey_inductive, Ez_inductive = E_inductive_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.) Jx_inductive = sig*Ex_inductive Jy_inductive = sig*Ey_inductive Jz_inductive = sig*Ez_inductive return Jx_inductive, Jy_inductive, Jz_inductive -def H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): - epsilon = 8.854187817*(10.**-12) - omega = 2.*np.pi*f +def H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.): + """ + Computing Magnetic fields from Electrical Dipole in a Wholespace + TODO: + Add description of parameters + """ + mu = mu_0*(1+kappa) + epsilon = epsilon_0*epsr XYZ = Utils.asArray_N_x_Dim(XYZ, 3) # Check if XYZ.shape[0] > 1 & f.shape[0] > 1: @@ -168,7 +217,7 @@ def H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., r = np.sqrt( dx**2. + dy**2. + dz**2.) # k = np.sqrt( -1j*2.*np.pi*f*mu*sig ) - k = np.sqrt( omega**2. *mu*epsilon -1j*omega*mu*sig ) + k = np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig ) front = current * length / (4.*np.pi* r**2) * (-1j*k*r + 1) * np.exp(-1j*k*r) @@ -191,18 +240,30 @@ def H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., return Hx, Hy, Hz -def B_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): - Hx, Hy, Hz = H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0) +def B_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.): + + """ + Computing Magnetic flux densites from Electrical Dipole in a Wholespace + TODO: + Add description of parameters + """ + + Hx, Hy, Hz = H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.) Bx = mu*Hx By = mu*Hy Bz = mu*Hz return Bx, By, Bz -def A_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): - epsilon = 8.854187817*(10.**-12) - omega = 2.*np.pi*f +def A_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.): + """ + Computing Electric vector potentials from Electrical Dipole in a Wholespace + TODO: + Add description of parameters + """ + mu = mu_0*(1+kappa) + epsilon = epsilon_0*epsr XYZ = Utils.asArray_N_x_Dim(XYZ, 3) # Check if XYZ.shape[0] > 1 & f.shape[0] > 1: @@ -213,7 +274,7 @@ def A_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., dz = XYZ[:,2]-srcLoc[2] r = np.sqrt( dx**2. + dy**2. + dz**2.) - k = np.sqrt( omega**2. *mu*epsilon -1j*omega*mu*sig ) + k = np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig ) front = current * length / (4.*np.pi*r) From 8b44f8d96b91dec91a1f0d5dc1cad1c96e59cf43 Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Thu, 23 Jun 2016 10:14:13 -0700 Subject: [PATCH 10/11] change FDEM_fields.py to FDEMDipolarfield.py --- SimPEG/EM/Analytics/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/EM/Analytics/__init__.py b/SimPEG/EM/Analytics/__init__.py index 9c5e5e28..8331501d 100644 --- a/SimPEG/EM/Analytics/__init__.py +++ b/SimPEG/EM/Analytics/__init__.py @@ -2,4 +2,4 @@ from TDEM import hzAnalyticDipoleT from FDEM import hzAnalyticDipoleF from FDEMcasing import * from DC import DCAnalyticHalf, DCAnalyticSphere -from FDEM_fields import * +from FDEMDipolarfields import * From 0763925743e4ff26ce2e4831824870b78ab5e0c7 Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Thu, 23 Jun 2016 14:12:49 -0700 Subject: [PATCH 11/11] fix minor bugs in analytics --- SimPEG/EM/Analytics/FDEMDipolarfields.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/SimPEG/EM/Analytics/FDEMDipolarfields.py b/SimPEG/EM/Analytics/FDEMDipolarfields.py index ef724542..7b79a680 100644 --- a/SimPEG/EM/Analytics/FDEMDipolarfields.py +++ b/SimPEG/EM/Analytics/FDEMDipolarfields.py @@ -18,7 +18,7 @@ def E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., """ mu = mu_0*(1+kappa) epsilon = epsilon_0*epsr - sig_hat = sig + 1j*omeg*epsilon + sig_hat = sig + 1j*omega(f)*epsilon XYZ = Utils.asArray_N_x_Dim(XYZ, 3) # Check @@ -66,7 +66,7 @@ def E_galvanic_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., le """ mu = mu_0*(1+kappa) epsilon = epsilon_0*epsr - sig_hat = sig + 1j*omeg*epsilon + sig_hat = sig + 1j*omega(f)*epsilon XYZ = Utils.asArray_N_x_Dim(XYZ, 3) # Check