Python SimPEG类代码示例

def analytic1DModelSource(mesh,freq,sigma_1d):
    '''
        Function that calculates and return background fields

        :param Simpeg mesh object mesh: Holds information on the discretization
        :param float freq: The frequency to solve at
        :param np.array sigma_1d: Background model of conductivity to base the calculations on, 1d model.
        :rtype: numpy.ndarray (mesh.nE,2)
        :return: eBG_bp, E fields for the background model at both polarizations.

    '''
    # import
    from SimPEG.MT.Utils import getEHfields
    # Get a 1d solution for a halfspace background
    if mesh.dim == 1:
        mesh1d = mesh
    elif mesh.dim == 2:
        mesh1d = simpeg.Mesh.TensorMesh([mesh.hy],np.array([mesh.x0[1]]))
    elif mesh.dim == 3:
        mesh1d = simpeg.Mesh.TensorMesh([mesh.hz],np.array([mesh.x0[2]]))

    # # Note: Everything is using e^iwt
    Eu, Ed, _, _ = getEHfields(mesh1d,sigma_1d,freq,mesh.vectorNz)
    # Make the fields into a dictionary of location and the fields
    e0_1d = Eu+Ed
    E1dFieldDict = dict(zip(mesh.vectorNz,e0_1d))
    if mesh.dim == 1:
        eBG_px = simpeg.mkvc(e0_1d,2)
        eBG_py = -simpeg.mkvc(e0_1d,2) # added a minus to make the results in the correct quadrents.
    elif mesh.dim == 2:
        ex_px = np.zeros(mesh.vnEx,dtype=complex)
        ey_px = np.zeros((mesh.nEy,1),dtype=complex)
        for i in np.arange(mesh.vnEx[0]):
            ex_px[i,:] = -e0_1d
        eBG_px = np.vstack((simpeg.Utils.mkvc(ex_px,2),ey_px))
        # Setup y (north) polarization (_py)
        ex_py = np.zeros((mesh.nEx,1), dtype='complex128')
        ey_py = np.zeros(mesh.vnEy, dtype='complex128')
        # Assign the source to ey_py
        for i in np.arange(mesh.vnEy[0]):
            ey_py[i,:] = e0_1d
        # ey_py[1:-1,1:-1,1:-1] = 0
        eBG_py = np.vstack((ex_py,simpeg.Utils.mkvc(ey_py,2),ez_py))
    elif mesh.dim == 3:
        # Setup x (east) polarization (_x)
        ex_px = -np.array([E1dFieldDict[i] for i in mesh.gridEx[:,2]]).reshape(-1,1)
        ey_px = np.zeros((mesh.nEy,1),dtype=complex)
        ez_px = np.zeros((mesh.nEz,1),dtype=complex)
        # Construct the full fields
        eBG_px = np.vstack((ex_px,ey_px,ez_px))
        # Setup y (north) polarization (_py)
        ex_py = np.zeros((mesh.nEx,1), dtype='complex128')
        ey_py = np.array([E1dFieldDict[i] for i in mesh.gridEy[:,2]]).reshape(-1,1)
        ez_py = np.zeros((mesh.nEz,1), dtype='complex128')
        # Construct the full fields
        eBG_py = np.vstack((ex_py,simpeg.Utils.mkvc(ey_py,2),ez_py))

    # Return the electric fields
    eBG_bp = np.hstack((eBG_px,eBG_py))
    return eBG_bp

def resampleNSEMdataAtFreq(NSEMdata, freqs):
    """
    Function to resample NSEMdata at set of frequencies

    """

    # Make a rec array
    NSEMrec = NSEMdata.toRecArray().data

    # Find unique locations
    uniLoc = np.unique(NSEMrec[['x','y','z']])
    uniFreq = NSEMdata.survey.freqs
    # Get the comps
    dNames = NSEMrec.dtype

    # Loop over all the locations and interpolate
    for loc in uniLoc:
        # Find the index of the station
        ind = np.sqrt(np.sum((rec_to_ndarr(NSEMrec[['x','y','z']]) - rec_to_ndarr(loc))**2,axis=1)) < 1. # Find dist of 1 m accuracy
        # Make a temporary recArray and interpolate all the components
        tArrRec = np.concatenate((simpeg.mkvc(freqs,2),np.ones((len(freqs),1))*rec_to_ndarr(loc),np.nan*np.ones((len(freqs),12))),axis=1).view(dNames)
        for comp in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi','tzxr','tzxi','tzyr','tzyi']:
            int1d = sciint.interp1d(NSEMrec[ind]['freq'],NSEMrec[ind][comp],bounds_error=False)
            tArrRec[comp] = simpeg.mkvc(int1d(freqs),2)

        # Join together
        try:
            outRecArr = recFunc.stack_arrays((outRecArr,tArrRec))
        except NameError:
            outRecArr = tArrRec

    # Make the NSEMdata and return
    return Data.fromRecArray(outRecArr)

def run(plotIt=True, nFreq=1):

    # Make a mesh
    M = simpeg.Mesh.TensorMesh(
        [
            [(100, 5, -1.5), (100., 10), (100, 5, 1.5)],
            [(100, 5, -1.5), (100., 10), (100, 5, 1.5)],
            [(100, 5, +1.6), (100., 10), (100, 3, 2)]
        ], x0=['C', 'C', -3529.5360]
    )
    # Setup the model
    conds = [1e-2, 1]
    sig = simpeg.Utils.ModelBuilder.defineBlock(
        M.gridCC, [-1000, -1000, -400], [1000, 1000, -200], conds
    )
    sig[M.gridCC[:, 2] > 0] = 1e-8
    sig[M.gridCC[:, 2] < -600] = 1e-1
    sigBG = np.zeros(M.nC) + conds[0]
    sigBG[M.gridCC[:, 2] > 0] = 1e-8

    # Setup the the survey object
    # Receiver locations
    rx_x, rx_y = np.meshgrid(np.arange(-500, 501, 50), np.arange(-500, 501, 50))
    rx_loc = np.hstack((simpeg.Utils.mkvc(rx_x, 2), simpeg.Utils.mkvc(rx_y, 2), np.zeros((np.prod(rx_x.shape), 1))))
    # Make a receiver list
    rxList = []
    for loc in rx_loc:
        # NOTE: loc has to be a (1, 3) np.ndarray otherwise errors accure
        for rx_orientation in ['xx', 'xy', 'yx', 'yy']:
            rxList.append(NSEM.Rx.Point_impedance3D(simpeg.mkvc(loc, 2).T, rx_orientation, 'real'))
            rxList.append(NSEM.Rx.Point_impedance3D(simpeg.mkvc(loc, 2).T, rx_orientation, 'imag'))
        for rx_orientation in ['zx', 'zy']:
            rxList.append(NSEM.Rx.Point_tipper3D(simpeg.mkvc(loc, 2).T, rx_orientation, 'real'))
            rxList.append(NSEM.Rx.Point_tipper3D(simpeg.mkvc(loc, 2).T, rx_orientation, 'imag'))
    # Source list
    srcList = [
        NSEM.Src.Planewave_xy_1Dprimary(rxList, freq)
        for freq in np.logspace(3, -3, nFreq)
    ]
    # Survey MT
    survey = NSEM.Survey(srcList)

    # Setup the problem object
    problem = NSEM.Problem3D_ePrimSec(M, sigma=sig, sigmaPrimary=sigBG)

    problem.pair(survey)
    problem.Solver = Solver

    # Calculate the data
    fields = problem.fields()
    dataVec = survey.eval(fields)

    # Make the data
    mtData = NSEM.Data(survey, dataVec)

    # Add plots
    if plotIt:
        pass

def convert3Dto1Dobject(NSEMdata, rxType3D='yx'):

    # Find the unique locations
    # Need to find the locations
    recDataTemp = NSEMdata.toRecArray().data.flatten()
    # Check if survey.std has been assigned.
    ## NEED TO: write this...
    # Calculte and add the DET of the tensor to the recArray
    if 'det' in rxType3D:
        Zon = (recDataTemp['zxxr']+1j*recDataTemp['zxxi'])*(recDataTemp['zyyr']+1j*recDataTemp['zyyi'])
        Zoff = (recDataTemp['zxyr']+1j*recDataTemp['zxyi'])*(recDataTemp['zyxr']+1j*recDataTemp['zyxi'])
        det = np.sqrt(Zon - Zoff)
        recData = recFunc.append_fields(recDataTemp,['zdetr','zdeti'],[det.real,det.imag] )
    else:
        recData = recDataTemp



    uniLocs = rec_to_ndarr(np.unique(recData[['x','y','z']]))
    mtData1DList = []
    if 'xy' in rxType3D:
        corr = -1 # Shift the data to comply with the quadtrature of the 1d problem
    else:
        corr = 1
    for loc in uniLocs:
        # Make the receiver list
        rx1DList = []
        rx1DList.append(Point_impedance1D(simpeg.mkvc(loc,2).T,'real'))
        rx1DList.append(Point_impedance1D(simpeg.mkvc(loc,2).T,'imag'))
        # Source list
        locrecData = recData[np.sqrt(np.sum( (rec_to_ndarr(recData[['x','y','z']]) - loc )**2,axis=1)) < 1e-5]
        dat1DList = []
        src1DList = []
        for freq in locrecData['freq']:
            src1DList.append(Planewave_xy_1Dprimary(rx1DList,freq))
            for comp  in ['r','i']:
                dat1DList.append( corr * locrecData[rxType3D+comp][locrecData['freq']== freq] )

        # Make the survey
        sur1D = Survey(src1DList)

        # Make the data
        dataVec = np.hstack(dat1DList)
        dat1D = Data(sur1D,dataVec)
        sur1D.dobs = dataVec
        # Need to take NSEMdata.survey.std and split it as well.
        std=0.05
        sur1D.std =  np.abs(sur1D.dobs*std) #+ 0.01*np.linalg.norm(sur1D.dobs)
        mtData1DList.append(dat1D)

    # Return the the list of data.
    return mtData1DList

def setupSurvey(sigmaHalf,tD=True):

    # Frequency
    nFreq = 33
    freqs = np.logspace(3,-3,nFreq)
    # Make the mesh
    ct = 5
    air = meshTensor([(ct,25,1.3)])
    # coreT0 = meshTensor([(ct,15,1.2)])
    # coreT1 = np.kron(meshTensor([(coreT0[-1],15,1.3)]),np.ones((7,)))
    core = np.concatenate( (  np.kron(meshTensor([(ct,15,-1.2)]),np.ones((10,))) , meshTensor([(ct,20)]) ) )
    bot = meshTensor([(core[0],10,-1.3)])
    x0 = -np.array([np.sum(np.concatenate((core,bot)))])
    m1d = simpeg.Mesh.TensorMesh([np.concatenate((bot,core,air))], x0=x0)
    # Make the model
    sigma = np.zeros(m1d.nC) + sigmaHalf
    sigma[m1d.gridCC > 0 ] = 1e-8

    rxList = []
    for rxType in ['z1dr','z1di']:
        rxList.append(simpegmt.SurveyMT.RxMT(simpeg.mkvc(np.array([0.0]),2).T,rxType))
    # Source list
    srcList =[]
    if tD:
        for freq in freqs:
            srcList.append(simpegmt.SurveyMT.srcMT_polxy_1DhomotD(rxList,freq))
    else:
        for freq in freqs:
            srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,freq))

    survey = simpegmt.SurveyMT.SurveyMT(srcList)
    return survey, sigma, m1d

def setup1DSurvey(sigmaHalf, tD=False, structure=False):

    # Frequency
    nFreq = 33
    freqs = np.logspace(3, -3, nFreq)
    # Make the mesh
    ct = 5
    air = meshTensor([(ct, 25, 1.3)])
    # coreT0 = meshTensor([(ct,15,1.2)])
    # coreT1 = np.kron(meshTensor([(coreT0[-1],15,1.3)]),np.ones((7,)))
    core = np.concatenate((np.kron(meshTensor([(ct, 15, -1.2)]), np.ones((10,))), meshTensor([(ct, 20)])))
    bot = meshTensor([(core[0], 20, -1.3)])
    x0 = -np.array([np.sum(np.concatenate((core, bot)))])
    m1d = simpeg.Mesh.TensorMesh([np.concatenate((bot, core, air))], x0=x0)
    # Make the model
    sigma = np.zeros(m1d.nC) + sigmaHalf
    sigma[m1d.gridCC > 0] = 1e-8
    sigmaBack = sigma.copy()
    # Add structure
    if structure:
        shallow = (m1d.gridCC < -200) * (m1d.gridCC > -600)
        deep = (m1d.gridCC < -3000) * (m1d.gridCC > -5000)
        sigma[shallow] = 1
        sigma[deep] = 0.1

    rxList = []
    for rxType in ["z1d", "z1d"]:
        rxList.append(Point_impedance1D(simpeg.mkvc(np.array([0.0]), 2).T, "real"))
        rxList.append(Point_impedance1D(simpeg.mkvc(np.array([0.0]), 2).T, "imag"))
    # Source list
    srcList = []
    if tD:
        for freq in freqs:
            srcList.append(Planewave_xy_1DhomotD(rxList, freq))
    else:
        for freq in freqs:
            srcList.append(Planewave_xy_1Dprimary(rxList, freq))

    survey = Survey(srcList)
    return (survey, sigma, sigmaBack, m1d)

def DerivProjfieldsTest(inputSetup, comp="All", freq=False):

    survey, problem = NSEM.Utils.testUtils.setupSimpegNSEM_ePrimSec(inputSetup, comp, freq)
    print("Derivative test of data projection for eFormulation primary/secondary\n")
    # problem.mapping = simpeg.Maps.ExpMap(problem.mesh)
    # Initate things for the derivs Test
    src = survey.srcList[0]
    np.random.seed(1983)
    u0x = np.random.randn(survey.mesh.nE) + np.random.randn(survey.mesh.nE) * 1j
    u0y = np.random.randn(survey.mesh.nE) + np.random.randn(survey.mesh.nE) * 1j
    u0 = np.vstack((simpeg.mkvc(u0x, 2), simpeg.mkvc(u0y, 2)))
    f0 = problem.fieldsPair(survey.mesh, survey)
    # u0 = np.hstack((simpeg.mkvc(u0_px,2),simpeg.mkvc(u0_py,2)))
    f0[src, "e_pxSolution"] = u0[: len(u0) / 2]  # u0x
    f0[src, "e_pySolution"] = u0[len(u0) / 2 : :]  # u0y

    def fun(u):
        f = problem.fieldsPair(survey.mesh, survey)
        f[src, "e_pxSolution"] = u[: len(u) / 2]
        f[src, "e_pySolution"] = u[len(u) / 2 : :]
        return rx.eval(src, survey.mesh, f), lambda t: rx.evalDeriv(src, survey.mesh, f0, simpeg.mkvc(t, 2))

    return simpeg.Tests.checkDerivative(fun, u0, num=4, plotIt=False, eps=FLR)

def dataMis_AnalyticPrimarySecondary(sigmaHalf):

    # Make the survey
    # Primary secondary
    surveyPS, sigmaPS, mesh = setupSurvey(sigmaHalf,tD=False)
    problemPS = MT.Problem1D.eForm_psField(mesh)
    problemPS.sigmaPrimary = sigmaPS
    problemPS.pair(surveyPS)
    # Analytic data
    dataAnaObj = calculateAnalyticSolution(surveyPS.srcList,mesh,sigmaPS)

    dataPS = surveyPS.dpred(sigmaPS)
    dataAna = simpeg.mkvc(dataAnaObj)
    return np.all((dataPS - dataAna)/dataAna < 2.)

def run(plotIt=True, nFreq=1):
    """
        MT: 3D: Forward
        =======================

        Forward model 3D MT data.

    """

    # Make a mesh
    M = simpeg.Mesh.TensorMesh([[(100,5,-1.5),(100.,10),(100,5,1.5)],[(100,5,-1.5),(100.,10),(100,5,1.5)],[(100,5,1.6),(100.,10),(100,3,2)]], x0=['C','C',-3529.5360])
    # Setup the model
    conds = [1e-2,1]
    sig = simpeg.Utils.ModelBuilder.defineBlock(M.gridCC,[-1000,-1000,-400],[1000,1000,-200],conds)
    sig[M.gridCC[:,2]>0] = 1e-8
    sig[M.gridCC[:,2]<-600] = 1e-1
    sigBG = np.zeros(M.nC) + conds[0]
    sigBG[M.gridCC[:,2]>0] = 1e-8

    ## Setup the the survey object
    # Receiver locations
    rx_x, rx_y = np.meshgrid(np.arange(-500,501,50),np.arange(-500,501,50))
    rx_loc = np.hstack((simpeg.Utils.mkvc(rx_x,2),simpeg.Utils.mkvc(rx_y,2),np.zeros((np.prod(rx_x.shape),1))))
    # Make a receiver list
    rxList = []
    for loc in rx_loc:
        # NOTE: loc has to be a (1,3) np.ndarray otherwise errors accure
        for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi','tzxr','tzxi','tzyr','tzyi']:
            rxList.append(MT.Rx(simpeg.mkvc(loc,2).T,rxType))
    # Source list
    srcList =[]
    for freq in np.logspace(3,-3,nFreq):
        srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq))
    # Survey MT
    survey = MT.Survey(srcList)

    ## Setup the problem object
    problem = MT.Problem3D.eForm_ps(M, sigmaPrimary=sigBG)
    problem.pair(survey)
    problem.Solver = Solver

    # Calculate the data
    fields = problem.fields(sig)
    dataVec = survey.eval(fields)

    # Make the data
    mtData = MT.Data(survey,dataVec)
    # Add plots
    if plotIt:
        pass

def dataMis_AnalyticTotalDomain(sigmaHalf):

    # Make the survey

    # Total domain solution
    surveyTD, sigma, mesh = setupSurvey(sigmaHalf)
    problemTD = MT.Problem1D.eForm_TotalField(mesh)
    problemTD.pair(surveyTD)
    # Analytic data
    dataAnaObj = calculateAnalyticSolution(surveyTD.srcList,mesh,sigma)
    # dataTDObj = MT.DataMT.DataMT(surveyTD, surveyTD.dpred(sigma))
    dataTD = surveyTD.dpred(sigma)
    dataAna = simpeg.mkvc(dataAnaObj)
    return np.all((dataTD - dataAna)/dataAna < 2.)

def DerivProjfieldsTest(inputSetup,comp='All',freq=False):

    survey, problem = setupSimpegMTfwd_eForm_ps(inputSetup,comp,freq)
    print 'Derivative test of data projection for eFormulation primary/secondary\n\n'
    # problem.mapping = simpeg.Maps.ExpMap(problem.mesh)
    # Initate things for the derivs Test
    src = survey.srcList[0]
    rx = src.rxList[0]

    u0x = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j
    u0y = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j
    u0 = np.vstack((simpeg.mkvc(u0x,2),simpeg.mkvc(u0y,2)))
    f0 = problem.fieldsPair(survey.mesh,survey)
    # u0 = np.hstack((simpeg.mkvc(u0_px,2),simpeg.mkvc(u0_py,2)))
    f0[src,'e_pxSolution'] =  u0[:len(u0)/2]#u0x
    f0[src,'e_pySolution'] = u0[len(u0)/2::]#u0y

    def fun(u):
        f = problem.fieldsPair(survey.mesh,survey)
        f[src,'e_pxSolution'] = u[:len(u)/2]
        f[src,'e_pySolution'] = u[len(u)/2::]
        return rx.projectFields(src,survey.mesh,f), lambda t: rx.projectFieldsDeriv(src,survey.mesh,f0,simpeg.mkvc(t,2))

    return simpeg.Tests.checkDerivative(fun, u0, num=3, plotIt=False, eps=FLR)

 def fun(u):
     f = problem.fieldsPair(survey.mesh, survey)
     f[src, "e_pxSolution"] = u[: len(u) / 2]
     f[src, "e_pySolution"] = u[len(u) / 2 : :]
     return rx.eval(src, survey.mesh, f), lambda t: rx.evalDeriv(src, survey.mesh, f0, simpeg.mkvc(t, 2))

def run(plotIt=True):

    # Setup the forward modeling
    # Setting up 1D mesh and conductivity models to forward model data.
    # Frequency
    nFreq = 26
    freqs = np.logspace(2, -3, nFreq)
    # Set mesh parameters
    ct = 10
    air = simpeg.Utils.meshTensor([(ct, 25, 1.4)])
    core = np.concatenate((
        np.kron(simpeg.Utils.meshTensor([(ct, 15, -1.2)]), np.ones((8, ))),
        simpeg.Utils.meshTensor([(ct, 5)])))
    bot = simpeg.Utils.meshTensor([(core[0], 20, -1.4)])
    x0 = -np.array([np.sum(np.concatenate((core, bot)))])
    # Make the model
    m1d = simpeg.Mesh.TensorMesh([np.concatenate((bot, core, air))], x0=x0)

    # Setup model variables
    active = m1d.vectorCCx < 0.
    layer1 = (m1d.vectorCCx < -500.) & (m1d.vectorCCx >= -800.)
    layer2 = (m1d.vectorCCx < -3500.) & (m1d.vectorCCx >= -5000.)
    # Set the conductivity values
    sig_half = 1e-2
    sig_air = 1e-8
    sig_layer1 = .2
    sig_layer2 = .2
    # Make the true model
    sigma_true = np.ones(m1d.nCx) * sig_air
    sigma_true[active] = sig_half
    sigma_true[layer1] = sig_layer1
    sigma_true[layer2] = sig_layer2
    # Extract the model
    m_true = np.log(sigma_true[active])
    # Make the background model
    sigma_0 = np.ones(m1d.nCx) * sig_air
    sigma_0[active] = sig_half
    m_0 = np.log(sigma_0[active])

    # Set the mapping
    actMap = simpeg.Maps.InjectActiveCells(
        m1d, active, np.log(1e-8), nC=m1d.nCx)
    mappingExpAct = simpeg.Maps.ExpMap(m1d) * actMap

    # Setup the layout of the survey, set the sources and the connected receivers
    # Receivers
    rxList = []
    rxList.append(
        NSEM.Rx.Point_impedance1D(simpeg.mkvc(np.array([-0.5]), 2).T, 'real'))
    rxList.append(
        NSEM.Rx.Point_impedance1D(simpeg.mkvc(np.array([-0.5]), 2).T, 'imag'))
    # Source list
    srcList = []
    for freq in freqs:
            srcList.append(NSEM.Src.Planewave_xy_1Dprimary(rxList, freq))
    # Make the survey
    survey = NSEM.Survey(srcList)
    survey.mtrue = m_true

    # Set the problem
    problem = NSEM.Problem1D_ePrimSec(
        m1d, sigmaPrimary=sigma_0, sigmaMap=mappingExpAct)
    problem.pair(survey)

    # Forward model data
    # Project the data
    survey.dtrue = survey.dpred(m_true)
    survey.dobs = (
        survey.dtrue + 0.01 *
        abs(survey.dtrue) *
        np.random.randn(*survey.dtrue.shape))

    # Assign uncertainties
    std = 0.025  # 5% std
    survey.std = np.abs(survey.dobs * std)
    # Assign the data weight
    Wd = 1. / survey.std

    # Setup the inversion proceedure
    # Define a counter
    C = simpeg.Utils.Counter()
    # Optimization
    opt = simpeg.Optimization.InexactGaussNewton(maxIter=25)
    opt.counter = C
    opt.LSshorten = 0.1
    opt.remember('xc')
    # Data misfit
    dmis = simpeg.DataMisfit.l2_DataMisfit(survey)
    dmis.W = Wd
    # Regularization - with a regularization mesh
    regMesh = simpeg.Mesh.TensorMesh([m1d.hx[active]], m1d.x0)
    reg = simpeg.Regularization.Tikhonov(regMesh)
    reg.mrefInSmooth = True
    reg.alpha_s = 1e-2
    reg.alpha_x = 1.
    reg.mrefInSmooth = True
    # Inversion problem
    invProb = simpeg.InvProblem.BaseInvProblem(dmis, reg, opt)
    invProb.counter = C
    # Beta schedule
#.........这里部分代码省略.........

# Load the model to the uniform cell mesh
modelUniCell = simpeg.Utils.meshutils.readUBCTensorModel(modelname,mesh3dCons)

# Load the model to the mesh with padding cells
modelT = simpeg.Utils.meshutils.readUBCTensorModel(modelname,mesh3d)
# Adjust the model to reflect changes in the mesh (fewer aircells)
modMat = mesh3d.r(modelT,'CC','CC','M')
modNewMat = np.ones((50,50,48))*modMat[0,0,0]
modNewMat[:,:,9::] = modMat[:,:,:-9]
modelTD = mesh3d.r(modNewMat,'CC','CC','V')

# Define the data locations
xG,yG = np.meshgrid(np.linspace(-700,700,8),np.linspace(-700,700,8))
zG = np.zeros_like(xG)
locs = np.hstack((simpeg.mkvc(xG.ravel(),2),simpeg.mkvc(yG.ravel(),2),simpeg.mkvc(zG.ravel(),2)))

# Make the receiver list
rxList = []
for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi']:
    rxList.append(simpegmt.SurveyMT.RxMT(locs,rxType))
# Source list
srcList =[]
freqs = np.logspace(4,0,17)
for freq in freqs:
    srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,freq))
# Survey MT
survey = simpegmt.SurveyMT.SurveyMT(srcList)

# Setup the problem object
sigma1d = mesh3d.r(modelTD,'CC','CC','M')[0,0,:] # Use the edge column as a background model

def homo1DModelSource(mesh, freq, sigma_1d):
    """
        Function that calculates and return background fields

        :param Simpeg mesh object mesh: Holds information on the discretization
        :param float freq: The frequency to solve at
        :param np.array sigma_1d: Background model of conductivity to base the calculations on, 1d model.
        :rtype: numpy.ndarray (mesh.nE, 2)
        :return: eBG_bp, E fields for the background model at both polarizations.

    """
    from . import get1DEfields
    # Get a 1d solution for a halfspace background
    if mesh.dim == 1:
        mesh1d = mesh
    elif mesh.dim == 2:
        mesh1d = simpeg.Mesh.TensorMesh([mesh.hy], np.array([mesh.x0[1]]))
    elif mesh.dim == 3:
        mesh1d = simpeg.Mesh.TensorMesh([mesh.hz], np.array([mesh.x0[2]]))

    # # Note: Everything is using e^iwt
    e0_1d = get1DEfields(mesh1d, sigma_1d, freq)
    if mesh.dim == 1:
        eBG_px = simpeg.mkvc(e0_1d, 2)
        eBG_py = -simpeg.mkvc(e0_1d, 2) # added a minus to make the results in the correct quadrents.
    elif mesh.dim == 2:
        ex_px = np.zeros(mesh.vnEx, dtype=complex)
        ey_px = np.zeros((mesh.nEy, 1), dtype=complex)
        for i in np.arange(mesh.vnEx[0]):
            ex_px[i, :] = -e0_1d
        eBG_px = np.vstack((simpeg.Utils.mkvc(ex_px, 2), ey_px))
        # Setup y (north) polarization (_py)
        ex_py = np.zeros((mesh.nEx, 1), dtype='complex128')
        ey_py = np.zeros(mesh.vnEy, dtype='complex128')
        # Assign the source to ey_py
        for i in np.arange(mesh.vnEy[0]):
            ey_py[i, :] = e0_1d
        # ey_py[1:-1, 1:-1, 1:-1] = 0
        eBG_py = np.vstack((ex_py, simpeg.Utils.mkvc(ey_py, 2), ez_py))
    elif mesh.dim == 3:
        # Setup x (east) polarization (_x)
        ex_px = np.zeros(mesh.vnEx, dtype=complex)
        ey_px = np.zeros((mesh.nEy, 1), dtype=complex)
        ez_px = np.zeros((mesh.nEz, 1), dtype=complex)
        # Assign the source to ex_x
        for i in np.arange(mesh.vnEx[0]):
            for j in np.arange(mesh.vnEx[1]):
                ex_px[i, j, :] = -e0_1d
        eBG_px = np.vstack((simpeg.Utils.mkvc(ex_px, 2), ey_px, ez_px))
        # Setup y (north) polarization (_py)
        ex_py = np.zeros((mesh.nEx, 1), dtype='complex128')
        ey_py = np.zeros(mesh.vnEy, dtype='complex128')
        ez_py = np.zeros((mesh.nEz, 1), dtype='complex128')
        # Assign the source to ey_py
        for i in np.arange(mesh.vnEy[0]):
            for j in np.arange(mesh.vnEy[1]):
                ey_py[i, j, :] = e0_1d
        # ey_py[1:-1, 1:-1, 1:-1] = 0
        eBG_py = np.vstack((ex_py, simpeg.Utils.mkvc(ey_py, 2), ez_py))

    # Return the electric fields
    eBG_bp = np.hstack((eBG_px, eBG_py))
    return eBG_bp