

In [1]:

    
%pylab inline
import sys
sys.path.append('../')

import numpy as np
from scipy import sparse as sp
import matplotlib.pyplot as plt
from SimPEG import TensorMesh









    



Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.kernel.zmq.pylab.backend_inline].
For more information, type 'help(pylab)'.

Set up mesh



In [2]:

    
x0 = np.zeros(3)
h1 = np.ones(100)*0.5
h2 = np.ones(100)*0.5
h3 = np.ones(100)*0.5
mesh = TensorMesh([h1,h2,h3],x0)

Define magnetic dipole (MD) source



In [115]:

    
def sourceMD(mesh,dipoleLoc,dimDip=3):
    """
    Calculates the response from magnetic dipole source(s). 
      
    Inputs:
        mesh      - tensorMesh object
        dipoleLoc - an array of [n x 3] dipole locations (xyz)
        dimDip    - dipole dimension (e.g., x = 1,y = 2,z = 3); Default is 3.

    Outputs:
        b - the magnetic field on the face grid
    """
    
    def getMagneticDipolePotential(loc,m,grid,dimAxis):
        # Get the potential of the dipole
        # Returns A(# faces x # transmitters)
        
        nL = grid.shape[0]
        nT = loc.shape[0]
        fullLoc = ones((nL,3))
        A = zeros((nL,nT))
        for ii in range(nT):
            fullM = ones((nL,3))*m[ii,:]
            br = grid - fullLoc*loc[ii,:]
            cp = np.cross(fullM,br)
            r = sqrt(br[:,0]**2 + br[:,1]**2 + br[:,2]**2)
            A[:,ii] = ((1e-6)*cp[:,dimAxis-1])/r**3
        return A
              
    dipoleLoc = np.atleast_2d(dipoleLoc)
    dimDip = np.atleast_2d(dimDip)
    m = np.zeros((dipoleLoc.shape[0],3))
    if (dimDip.shape[0] == 1):
        try:
            m[:,dimDip-1] = 1
        except IndexError:
            print "magneticSource:Error: Dipole dimension should be 1 (x), 2 (y), or 3 (z)."
            raise
    elif (dimDip.shape[0] == dipoleLoc.shape[0]):
        try:
            for jj in range(dipoleLoc.shape[0]):
                m[jj,dimDip[jj] - 1] = 1

        except IndexError:
           print "magneticSource:Error: Dipole dimension should be 1 (x), 2 (y), or 3 (z)."
           raise
    else:
        print "magneticSource:Error: Dipole direction should also be vector of same length as dipole locations."
        raise

    # Get magnetic potential at each set of orthogonal faces:
    Ax = getMagneticDipolePotential(dipoleLoc,m,mesh.gridEx,1)
    Ay = getMagneticDipolePotential(dipoleLoc,m,mesh.gridEy,2)
    Az = getMagneticDipolePotential(dipoleLoc,m,mesh.gridEz,3)

    # Combine potential
    A = concatenate((Ax, Ay, Az),axis=0)
    
    # B = curl A
    CURL = mesh.edgeCurl
    b0 = CURL*A

    return b0

Dipole location and oriention



In [112]:

    
# make 2 dipole locations
loc = np.array([[3.2,27.1,40.1],[42.0,5,10]])
m = np.array((3,1),'i') # vertical dipole

Test code



In [117]:

    
A = sourceMD(mesh,loc,m)
print A.shape
#newL = np.atleast_2d(loc)
#print newM.shape,newL.shape
#mVec = np.zeros((newL.shape[0],3))
#print mVec
#for jj in range(mVec.shape[0]):
#    mVec[jj, m[jj] - 1] = 1
#print "m: ",m,m.shape[0]









    



(3030000L, 2L)

Plot faces



In [98]:

    
plt.colorbar(mesh.plotImage(np.log(mesh.r(A[:,0],'F','Fz','V')),'Fz'))









    



-c:1: RuntimeWarning: invalid value encountered in log
C:\Users\kdavis\AppData\Local\Enthought\Canopy\App\appdata\canopy-1.0.3.1262.win-x86_64\lib\site-packages\matplotlib\figure.py:362: UserWarning: matplotlib is currently using a non-GUI backend, so cannot show the figure
  "matplotlib is currently using a non-GUI backend, "






    Out[98]:





<matplotlib.colorbar.Colorbar instance at 0x0000000008A1AC48>



In [100]:

    
plt.colorbar(mesh.plotImage(np.log(mesh.r(A[:,1],'F','Fz','V')),'Fz'))









    



-c:1: RuntimeWarning: divide by zero encountered in log
C:\Users\kdavis\AppData\Local\Enthought\Canopy\App\appdata\canopy-1.0.3.1262.win-x86_64\lib\site-packages\matplotlib\colors.py:897: RuntimeWarning: invalid value encountered in true_divide
  resdat /= (vmax - vmin)
C:\Users\kdavis\AppData\Local\Enthought\Canopy\App\appdata\canopy-1.0.3.1262.win-x86_64\lib\site-packages\matplotlib\colorbar.py:561: RuntimeWarning: invalid value encountered in greater
  inrange = (ticks > -0.001) & (ticks < 1.001)
C:\Users\kdavis\AppData\Local\Enthought\Canopy\App\appdata\canopy-1.0.3.1262.win-x86_64\lib\site-packages\matplotlib\colorbar.py:561: RuntimeWarning: invalid value encountered in less
  inrange = (ticks > -0.001) & (ticks < 1.001)






    Out[100]:





<matplotlib.colorbar.Colorbar instance at 0x0000000018D2BAC8>






    



C:\Users\kdavis\AppData\Local\Enthought\Canopy\App\appdata\canopy-1.0.3.1262.win-x86_64\lib\site-packages\matplotlib\colors.py:896: RuntimeWarning: invalid value encountered in subtract
  resdat -= vmin
C:\Users\kdavis\AppData\Local\Enthought\Canopy\App\appdata\canopy-1.0.3.1262.win-x86_64\lib\site-packages\matplotlib\colors.py:557: RuntimeWarning: invalid value encountered in less
  cbook._putmask(xa, xa < 0.0, -1)



In [20]:

    
plt.colorbar(mesh.plotImage(np.log(mesh.r(A[:,0],'F','Fy','V')),'Fy'))









    Out[20]:





<matplotlib.colorbar.Colorbar instance at 0x00000000342B2208>



In [ ]: