From 18357a11da01a85b0b403074f9b8962662f8b51e Mon Sep 17 00:00:00 2001
From: GudniRos <gkrosen@gmail.com>
Date: Fri, 15 Apr 2016 14:48:58 -0700
Subject: [PATCH 01/14] Fixing analytic function

---
 SimPEG/MT/Utils/MT1Danalytic.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/SimPEG/MT/Utils/MT1Danalytic.py b/SimPEG/MT/Utils/MT1Danalytic.py
index 1fc177f4..d2ed7a55 100644
--- a/SimPEG/MT/Utils/MT1Danalytic.py
+++ b/SimPEG/MT/Utils/MT1Danalytic.py
@@ -3,7 +3,7 @@
 import numpy as np, SimPEG as simpeg
 from scipy.constants import mu_0, epsilon_0 as eps_0
 
-def getEHfields(m1d,sigma,freq,zd,scaleUD=True):
+def getEHfields(m1d,sigma,freq,zd,scaleUD=True,scaleValue=1):
     '''Analytic solution for MT 1D layered earth. Returns E and H fields.
 
     :param SimPEG.mesh, object m1d: Mesh object with the 1D spatial information.
@@ -12,7 +12,7 @@ def getEHfields(m1d,sigma,freq,zd,scaleUD=True):
     :param numpy array, vector zd: location to calculate EH fields at
     :param bollean, scaleUD: scales the output to be 1 at the top, increases numeracal stability.
 
-    Assumes a halfspace with the same conductive as the last cell below.
+    Assumes a halfspace with the same conductive as the deepest cell.
 
     '''
     # Note add an error check for the mesh and sigma are the same size.
@@ -29,7 +29,7 @@ def getEHfields(m1d,sigma,freq,zd,scaleUD=True):
 
     # Initiate the propagation matrix, in the order down up.
     UDp = np.zeros((2,m1d.nC+1),dtype=complex)
-    UDp[1,0] = 1. # Set the wave amplitude as 1 into the half-space at the bottom of the mesh
+    UDp[1,0] = scaleValue # Set the wave amplitude as 1 into the half-space at the bottom of the mesh
     # Loop over all the layers, starting at the bottom layer
     for lnr, h in enumerate(m1d.hx): # lnr-number of layer, h-thickness of the layer
         # Calculate
@@ -53,7 +53,7 @@ def getEHfields(m1d,sigma,freq,zd,scaleUD=True):
             if np.any(np.isnan(scaleVal)):
                 # If there is a nan (thickness very great), rebuild the move up cell
                 scaleVal = np.zeros_like(UDp[:,lnr+1::-1],dtype=complex)
-                scaleVal[1,0] = 1.
+                scaleVal[1,0] = scaleValue
 
             UDp[:,lnr+1::-1] = scaleVal
 

From 22f0a742e732277fcd518da0e1b2812ad921eb93 Mon Sep 17 00:00:00 2001
From: GudniRos <gkrosen@gmail.com>
Date: Mon, 2 May 2016 05:16:56 -0700
Subject: [PATCH 02/14] Fixed UBC mesh read in function

---
 SimPEG/Mesh/MeshIO.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/SimPEG/Mesh/MeshIO.py b/SimPEG/Mesh/MeshIO.py
index 3ae4d88f..2b0de9e7 100644
--- a/SimPEG/Mesh/MeshIO.py
+++ b/SimPEG/Mesh/MeshIO.py
@@ -21,7 +21,7 @@ class TensorMeshIO(object):
                 if '*' in seg:
                     st = seg
                     sp = seg.split('*')
-                    re = np.array(sp[0],dtype=int)*(' ' + sp[1])
+                    re = int(sp[0])*(' ' + sp[1])
                     line = line.replace(st,re.strip())
             return np.array(line.split(),dtype=float)
 

From b1569e57341bac56b5c770fda038d9a53ba98feb Mon Sep 17 00:00:00 2001
From: GudniRos <gkrosen@gmail.com>
Date: Tue, 10 May 2016 14:00:48 -0700
Subject: [PATCH 03/14] Changes to MT1D analytic example

---
 .../Examples/MT_1D_analytic_nlayer_Earth.py   | 428 ++++++++++++++++++
 1 file changed, 428 insertions(+)
 create mode 100644 SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py

diff --git a/SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py b/SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py
new file mode 100644
index 00000000..62ab2a68
--- /dev/null
+++ b/SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py
@@ -0,0 +1,428 @@
+from scipy.constants import epsilon_0, mu_0
+import matplotlib.pyplot as plt
+import numpy as np
+from ipywidgets import *
+from SimPEG.EM.Utils import k, omega
+
+"""
+MT1D: n layered earth problem
+*****************************
+
+Author: Thibaut Astic
+Contact: thast@eos.ubc.ca
+Date: January 2016
+
+This code compute the analytic response of a n-layered Earth to a plane wave (Magneto-Tellurics).
+
+We start by looking at Maxwell's equations in the electric
+field \\\(\\\mathbf{E}\\) and the magnetic flux
+\\\(\\\mathbf{H}\\) to write the wave equations
+\\(\\ \nabla ^2  \mathbf{E_x} + k^2 \mathbf{E_x} = 0 \\) &
+\\(\\ \nabla ^2  \mathbf{H_y} + k^2 \mathbf{H_y} = 0 \\)
+
+Then solving the equations in each layer "j" between z_{j-1} and z_j in the form of
+\\(\\ E_{x,j} (z) = U_j e^{i k (z-z_{j-1})} + D_j e^{-i k (z-z_{j-1})} \\)
+\\(\\ H_{y,j} (z) = \frac{1}{Z_j} (D_j e^{-i k (z-z_{j-1})} - U_j e^{i k (z-z_{j-1})}) \\)
+
+With U and D the Up and Down components of the E-field.
+
+The iteration from one layer to another is ensure by:
+
+\\(\\  \left(\begin{matrix} E_{x,j} \\ H_{y,j} \end{matrix} \right) =
+    P_j T_j P^{-1}_J \left(\begin{matrix} E_{x,j+1} \\ H_{y,j+1} \end{matrix} \right) \\)
+
+And the Boundary Condition is set for the E-field in the last layer, with no Up component (=0)
+and only a down component (=1 then normalized by the highest amplitude to ensure numeric stability)
+
+The layer 0 is assumed to be the air layer.
+
+"""
+
+#Define a frquency range for a survey
+frange = lambda minfreq, maxfreq, step: np.logspace(minfreq,maxfreq,num = step, base = 10.)
+
+#Functions to create random physical Perties for a n-layered earth
+thick = lambda minthick, maxthick, nlayer: np.append(np.array([1.2*10.**5]),
+                                                     np.ndarray.round(minthick + (maxthick-minthick)* np.random.rand(nlayer-1,1)
+                                                            ,decimals =1))
+
+sig = lambda minsig, maxsig, nlayer: np.append(np.array([0.]),
+                                               np.ndarray.round(10.**minsig + (10.**maxsig-10.**minsig)* np.random.rand(nlayer,1)
+                                                      ,decimals=3))
+
+mu  = lambda minmu, maxmu, nlayer: np.append(np.array([1.]),
+                                             np.ndarray.round(minmu + (maxmu-minmu)* np.random.rand(nlayer,1)
+                                                    ,decimals=1))
+
+eps = lambda mineps, maxeps, nlayer: np.append(np.array([1.]),
+                                               np.ndarray.round(mineps + (maxeps-mineps)* np.random.rand(nlayer,1)
+                                                                ,decimals=1))
+
+#Evaluate Impedance Z of a layer
+ImpZ = lambda f, mu, k: omega(f)*mu*mu_0/k
+
+#Complex Cole-Cole Conductivity - EM utils
+PCC= lambda siginf,m,t,c,f: siginf*(1.-(m/(1.+(1j*omega(f)*t)**c)))
+
+#Converted thickness array into top of layer array
+top = lambda thick: np.cumsum(thick)
+
+#Propagation Matrix and theirs inverses
+
+#matrix T for transition of Up and Down components accross a layer
+T = lambda h,k: np.matrix([[np.exp(1j*k*h),0.],[0.,np.exp(-1j*k*h)]],dtype='complex_')
+
+Tinv = lambda h,k: np.matrix([[np.exp(-1j*k*h),0.],[0.,np.exp(1j*k*h)]],dtype='complex_')
+
+#transition of Up and Down components accross a layer
+UD_Z = lambda UD,z,zj,k : T((z-zj),k)*UD
+
+
+#matrix P relating Up and Down components with E and H fields
+P = lambda z: np.matrix([[1.,1,],[-1./z,1./z]],dtype='complex_')
+
+Pinv = lambda z: np.matrix([[1.,-z],[1.,z]],dtype='complex_')/2.
+
+
+#Time Variation of E and H
+E_ZT = lambda U,D,f,t : np.exp(1j*omega(f)*t)*(U+D)
+H_ZT = lambda U,D,Z,f,t : (1./Z)*np.exp(1j*omega(f)*t)*(D-U)
+
+#Plot the configuration of the problem
+def PlotConfiguration(thick,sig,eps,mu,ax,widthg,z):
+
+    topn = top(thick)
+    widthn = np.arange(-widthg,widthg+widthg/10.,widthg/10.)
+
+    ax.set_ylim([z.min(),z.max()])
+    ax.set_xlim([-widthg,widthg])
+
+    ax.set_ylabel("Depth (m)", fontsize=16.)
+    ax.yaxis.tick_right()
+    ax.yaxis.set_label_position("right")
+
+    #define filling for the different layers
+    hatches=['/'  , '+', 'x', '|' , '\\', '-'  , 'o' , 'O' , '.' , '*' ]
+
+    #Write the physical properties of air
+    ax.annotate(("Air, $\sigma$ =%1.0f mS/m")%(sig[0]*10**(3)),
+            xy=(-widthg/2., -np.abs(z.max())/2.), xycoords='data',
+            xytext=(-widthg/2., -np.abs(z.max())/2.), textcoords='data',
+            fontsize=14.)
+
+    ax.annotate(("$\epsilon_r$= %1i")%(eps[0]),
+            xy=(-widthg/2., -np.abs(z.max())/3.), xycoords='data',
+            xytext=(-widthg/2., -np.abs(z.max())/3.), textcoords='data',
+            fontsize=14.)
+
+    ax.annotate(("$\mu_r$= %1i")%(mu[0]),
+            xy=(-widthg/2., -np.abs(z.max())/3.), xycoords='data',
+            xytext=(0, -np.abs(z.max())/3.), textcoords='data',
+            fontsize=14.)
+
+    #Write the physical properties of the differents layers up to the (n-1)-th and fill it with pattern
+    for i in range(1,len(topn)-1,1):
+        if topn[i] == topn[i+1]:
+            pass
+        else:
+            ax.annotate(("$\sigma$ =%3.3f mS/m")%(sig[i]*10**(3)),
+                xy=(0., (2.*topn[i]+topn[i+1])/3), xycoords='data',
+                xytext=(0., (2.*topn[i]+topn[i+1])/3), textcoords='data',
+                fontsize=14.)
+
+            ax.annotate(("$\epsilon_r$= %1i")%(eps[i]),
+                xy=(-widthg/1.1, (2.*topn[i]+topn[i+1])/3), xycoords='data',
+                xytext=(-widthg/1.1, (2.*topn[i]+topn[i+1])/3), textcoords='data',
+                fontsize=14.)
+
+            ax.annotate(("$\mu_r$= %1.2f")%(mu[i]),
+                xy=(-widthg/2., (2.*topn[i]+topn[i+1])/3), xycoords='data',
+                xytext=(-widthg/2., (2.*topn[i]+topn[i+1])/3), textcoords='data',
+                fontsize=14.)
+
+            ax.plot(widthn,topn[i]*np.ones_like(widthn),color='black')
+            ax.fill_between(widthn,topn[i],topn[i+1],alpha=0.3,color="none",edgecolor='black', hatch=hatches[(i-1)%10])
+
+    #Write the physical properties of the n-th layer and fill it with pattern
+    ax.plot(widthn,topn[-1]*np.ones_like(widthn),color='black')
+    ax.fill_between(widthn,topn[-1],z.max(),alpha=0.3,color="none",edgecolor='black', hatch=hatches[(len(topn)-2)%10])
+
+    ax.annotate(("$\sigma$ =%3.3f mS/m")%(sig[-1]*10**(3)),
+            xy=(0., (2.*topn[-1]+z.max())/3), xycoords='data',
+            xytext=(0., (2.*topn[-1]+z.max())/3), textcoords='data',
+            fontsize=14.)
+
+    ax.annotate(("$\epsilon_r$= %1i")%(eps[-1]),
+            xy=(-widthg/1.1, (2.*topn[-1]+z.max())/3), xycoords='data',
+            xytext=(-widthg/1.1, (2.*topn[-1]+z.max())/3), textcoords='data',
+            fontsize=14.)
+
+    ax.annotate(("$\mu_r$= %1.2f")%(mu[-1]),
+            xy=(-widthg/2., (2.*topn[-1]+z.max())/3), xycoords='data',
+            xytext=(-widthg/2., (2.*topn[-1]+z.max())/3), textcoords='data',
+            fontsize=14.)
+
+    #plot Trees!
+    ax.annotate("",
+            xy=(widthg/2., -1.*z.max()/5.), xycoords='data',
+            xytext=(widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.2,head_length=1.2',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(widthg/2., -3./4.*z.max()/5.), xycoords='data',
+            xytext=(widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.4,head_length=1.4',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(widthg/2., -1./2.*z.max()/5.), xycoords='data',
+            xytext=(widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.6,head_length=1.6',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(1.2*widthg/2., -1.*z.max()/5.), xycoords='data',
+            xytext=(1.2*widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.2,head_length=1.2',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(1.2*widthg/2., -3./4.*z.max()/5.), xycoords='data',
+            xytext=(1.2*widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.4,head_length=1.4',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(1.2*widthg/2., -1./2.*z.max()/5.), xycoords='data',
+            xytext=(1.2*widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.6,head_length=1.6',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(1.5*widthg/2., -1.*z.max()/5.), xycoords='data',
+            xytext=(1.5*widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.2,head_length=1.2',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(1.5*widthg/2., -3./4.*z.max()/5.), xycoords='data',
+            xytext=(1.5*widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.4,head_length=1.4',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(1.5*widthg/2., -1./2.*z.max()/5.), xycoords='data',
+            xytext=(1.5*widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.6,head_length=1.6',color='green',linewidth=2.)
+            )
+
+
+    ax.invert_yaxis()
+
+    return ax
+
+#Propagate Up and Down component for a certain frequency & evaluate E and H field
+
+def Propagate(f,H,sig,chg,taux,c,mu,eps,n):
+
+    sigcm = np.zeros_like(sig,dtype='complex_')
+
+    for j in range(1,len(sig)):
+        sigcm[j]=PCC(sig[j],chg[j],taux[j],c[j],f)
+
+    K = k(f, sigcm, mu, eps)
+    Z = ImpZ(f,mu,K)
+
+    EH = np.matrix(np.zeros((2,n+1),dtype = 'complex_'),dtype = 'complex_')
+    UD = np.matrix(np.zeros((2,n+1),dtype = 'complex_'),dtype = 'complex_')
+
+    UD[1,-1] = 1.
+
+    for i in range(-2,-(n+2),-1):
+
+        UD[:,i] = Tinv(H[i+1],K[i])*Pinv(Z[i])*P(Z[i+1])*UD[:,i+1]
+        UD = UD/((np.abs(UD[0,:]+UD[1,:])).max())
+
+    for j in range(0,n+1):
+        EH[:,j] = np.matrix([[1.,1,],[-1./Z[j],1./Z[j]]])*UD[:,j]
+
+    return UD, EH, Z ,K
+
+
+#Evaluate the apparent resistivity and phase for a frequency range
+def appres(F,H,sig,chg,taux,c,mu,eps,n):
+
+    Res = np.zeros_like(F)
+    Phase = np.zeros_like(F)
+    App_ImpZ= np.zeros_like(F,dtype='complex_')
+
+    for i in range(0,len(F)):
+
+        UD,EH,Z ,K = Propagate(F[i],H,sig,chg,taux,c,mu,eps,n)
+
+        App_ImpZ[i] = EH[0,1]/EH[1,1]
+
+        Res[i] = np.abs(App_ImpZ[i])**2./(mu_0*omega(F[i]))
+        Phase[i] = np.angle(App_ImpZ[i], deg = True)
+
+    return Res,Phase
+
+#Evaluate Up, Down components, E and H field, for a frequency range,
+#a discretized depth range and a time range (use to calculate envelope)
+def calculateEHzt(F,H,sig,chg,taux,c,mu,eps,n,zsample,tsample):
+
+    topc = top(H)
+
+    layer = np.zeros(len(zsample),dtype=np.int)-1
+
+    Exzt = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
+    Hyzt = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
+    Uz = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
+    Dz = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
+    UDaux = np.matrix(np.zeros((2,len(zsample)),dtype = 'complex_'),dtype = 'complex_')
+
+    for i in range(0,n+1,1):
+        layer = layer+(zsample>=topc[i])*1
+
+    for j in range(0,len(F)):
+
+        UD,EH,Z ,K = Propagate(F[j],H,sig,chg,taux,c,mu,eps,n)
+
+        for p in range(0,len(zsample)):
+
+            UDaux[:,p] = UD_Z(UD[:,layer[p]],zsample[p],topc[layer[p]],K[layer[p]])
+
+            for q in range(0,len(tsample)):
+
+                Exzt[p,q]  = Exzt[p,q] + E_ZT(UDaux[0,p],UDaux[1,p],F[j],tsample[q])/len(F)
+                Hyzt[p,q] = Hyzt[p,q] + H_ZT(UDaux[0,p],UDaux[1,p],Z[layer[p]],F[j],tsample[q])/len(F)
+                Uz[p,q] = Uz[p,q] + UDaux[0,p]*np.exp(1j*omega(F[j])*tsample[q])/len(F)
+                Dz[p,q] = Dz[p,q] + UDaux[1,p]*np.exp(1j*omega(F[j])*tsample[q])/len(F)
+
+    return  Exzt,Hyzt,Uz,Dz,UDaux,layer
+
+
+#Function to Plot Apparent Resistivity and Phase
+def PlotAppRes(F,H,sig,chg,taux,c,mu,eps,n,fenvelope,PlotEnvelope):
+
+    Res, Phase = appres(F,H,sig,chg,taux,c,mu,eps,n)
+
+    fig,ax = plt.subplots(1,2,figsize=(16,10))
+
+    ax[0].scatter(Res,F,color='black')
+    ax[0].set_xscale('Log')
+    ax[0].set_yscale('Log')
+    ax[0].set_xlim([10.**(np.log10(Res.min())-1.),10.**(np.log10(Res.max())+1.)])
+    ax[0].set_ylim([F.min(),F.max()])
+    ax[0].set_xlabel('Apparent Resistivity (Ohm*m)',fontsize=16.,color="black")
+    ax[0].set_ylabel('Frequency (Hz)',fontsize=16.)
+    ax[0].grid(which='major')
+
+    ax0 = ax[0].twiny()
+
+    ax0.set_xlim([0.,90.])
+    ax0.set_ylim([F.min(),F.max()])
+    ax0.scatter(Phase,F,color='purple')
+    ax0.set_xlabel('Phase (Degrees)',fontsize=16.,color="purple")
+
+    zc=np.arange(-(H[1:].max()+10)*n,(H[1:].max()+10)*n,10.)
+
+    ax[0].tick_params(labelsize=16)
+    ax[1].tick_params(labelsize=16)
+    ax0.tick_params(labelsize=16)
+
+    if PlotEnvelope:
+
+        widthn=np.logspace(np.log10(Res.min())-1., np.log10(Res.max())+1., num=100, endpoint=True, base=10.0)
+        fenvelope1n=np.ones(100)*fenvelope
+        ax[0].plot(widthn,fenvelope1n,linestyle='dashed',color='black')
+
+        tc=np.arange(0.,1./fenvelope,0.01/(fenvelope))
+        Exzt,Hyzt,Uz,Dz,UDaux,layer = calculateEHzt(np.array([fenvelope]),H,sig,chg,taux,c,mu,eps,n,zc,tc)
+
+        ax1=ax[1].twiny()
+
+        ax[1].tick_params(labelsize=16)
+        ax1.tick_params(labelsize=16)
+
+        ax[1].set_xlabel('Amplitude Electric Field E (V/m)',color='blue',fontsize=16)
+
+        ax1.set_xlabel('Amplitude Magnetic Field H (A/m)',color='red',fontsize=16)
+
+        ax[1].fill_betweenx(zc,np.squeeze(np.asarray(np.real(Exzt.min(axis=1)))),
+                      np.squeeze(np.asarray(np.real(Exzt.max(axis=1)))),
+                      color='blue', alpha=0.1)
+
+        ax1.fill_betweenx(zc,np.squeeze(np.asarray(np.real(Hyzt.min(axis=1)))),
+                      np.squeeze(np.asarray(np.real(Hyzt.max(axis=1)))),
+                      color='red', alpha=0.1)
+
+        ax[1] = PlotConfiguration(H,sig,eps,mu,ax[1],(1.5*np.abs(Exzt).max()),zc)
+        ax1.set_xlim([-1.5*np.abs(Hyzt).max(),1.5*np.abs(Hyzt).max()])
+        ax1.set_xlim([-1.5*np.abs(Hyzt).max(),1.5*np.abs(Hyzt).max()])
+    else:
+        print 'No envelop (if True, might be slow)'
+        ax[1] = PlotConfiguration(H,sig,eps,mu,ax[1],1.,zc)
+        ax[1].get_xaxis().set_ticks([])
+
+    plt.show()
+
+#Interactive MT for Notebook
+def PlotAppRes3LayersInteract(h1,h2,sigl1,sigl2,sigl3,mul1,mul2,mul3,epsl1,epsl2,epsl3,PlotEnvelope,F_Envelope):
+
+    frangn=frange(-5,5,100.)
+    sig3= np.array([0.,0.001,0.1, 0.001])
+    thick3 = np.array([120000.,50.,50.])
+    eps3=np.array([1.,1.,1.,1])
+    mu3=np.array([1.,1.,1.,1])
+    chg3=np.array([0.,0.1,0.,0.2])
+    chg3_0=np.array([0.,0.1,0.,0.])
+    taux3=np.array([0.,0.1,0.,0.1])
+    c3=np.array([1.,1.,1.,1.])
+
+    sig3[1]=sigl1
+    sig3[1]=10.**sig3[1]
+    sig3[2]=sigl2
+    sig3[2]=10.**sig3[2]
+    sig3[3]=sigl3
+    sig3[3]=10.**sig3[3]
+    mu3[1]=mul1
+    mu3[2]=mul2
+    mu3[3]=mul3
+    eps3[1]=epsl1
+    eps3[2]=epsl2
+    eps3[3]=epsl3
+    thick3[1]=h1
+    thick3[2]=h2
+
+    PlotAppRes(frangn,thick3,sig3,chg3_0,taux3,c3,mu3,eps3,3,F_Envelope,PlotEnvelope)
+
+
+def run(n,plotIt=True):
+    # something to make a plot
+
+    F = frange(-5.,5.,20)
+    H = thick(50.,100.,n)
+    sign = sig(-5.,0.,n)
+    mun = mu(1.,2.,n)
+    epsn = eps(1.,9.,n)
+    chg = np.zeros_like(sign)
+    taux = np.zeros_like(sign)
+    c = np.zeros_like(sign)
+
+    Res, Phase = appres(F,H,sign,chg,taux,c,mun,epsn,n)
+
+    if plotIt:
+
+        PlotAppRes(F, H, sign, chg, taux, c, mun, epsn, n, fenvelope=1000., PlotEnvelope=True)
+
+    return Res, Phase
+
+if __name__ == '__main__':
+    run(3)
+
+
+
+
+

From 1ab67b5790905052740273ef0062be9e184090d6 Mon Sep 17 00:00:00 2001
From: GudniRos <gkrosen@gmail.com>
Date: Tue, 10 May 2016 14:41:20 -0700
Subject: [PATCH 04/14] Updated MT-->NSEM for all the classes. Changed
 Vertical1DMap --> SurjectVertical1D

---
 SimPEG/MT/__init__.py                         |  5 --
 SimPEG/{MT/FieldsMT.py => NSEM/FieldsNSEM.py} | 26 ++++-----
 SimPEG/{MT/BaseMT.py => NSEM/NSEM.py}         | 16 ++---
 SimPEG/{MT => NSEM}/Problem1D/Probs.py        | 24 ++++----
 SimPEG/{MT => NSEM}/Problem1D/__init__.py     |  0
 SimPEG/{MT => NSEM}/Problem2D/Probs.py        |  0
 SimPEG/{MT => NSEM}/Problem2D/__init__.py     |  0
 SimPEG/{MT => NSEM}/Problem3D/Probs.py        | 14 ++---
 SimPEG/{MT => NSEM}/Problem3D/__init__.py     |  0
 SimPEG/{MT/SrcMT.py => NSEM/SrcNSEM.py}       | 26 ++++-----
 SimPEG/{MT/SurveyMT.py => NSEM/SurveyNSEM.py} | 28 ++++-----
 SimPEG/{MT => NSEM}/Utils/MT1Danalytic.py     |  0
 SimPEG/{MT => NSEM}/Utils/MT1Dsolutions.py    |  0
 SimPEG/{MT => NSEM}/Utils/__init__.py         |  0
 SimPEG/{MT => NSEM}/Utils/dataUtils.py        | 58 +++++++++----------
 SimPEG/{MT => NSEM}/Utils/ediFilesUtils.py    |  2 +-
 SimPEG/{MT => NSEM}/Utils/plotDataTypes.py    |  0
 SimPEG/{MT => NSEM}/Utils/sourceUtils.py      |  4 +-
 SimPEG/{MT => NSEM}/Utils/srcUtils.py         |  2 +-
 SimPEG/NSEM/__init__.py                       |  5 ++
 tests/mt/test_ApparentResistivityAnalytic.py  |  4 +-
 ...test_Problem1D_againstAnalyticHalfspace.py | 24 ++++----
 .../mt/test_Problem1D_totalDvsPSvsAnalytic.py | 28 ++++-----
 tests/mt/test_Problem3D_againstAnalytic.py    | 36 ++++++------
 24 files changed, 151 insertions(+), 151 deletions(-)
 delete mode 100644 SimPEG/MT/__init__.py
 rename SimPEG/{MT/FieldsMT.py => NSEM/FieldsNSEM.py} (95%)
 rename SimPEG/{MT/BaseMT.py => NSEM/NSEM.py} (91%)
 rename SimPEG/{MT => NSEM}/Problem1D/Probs.py (93%)
 rename SimPEG/{MT => NSEM}/Problem1D/__init__.py (100%)
 rename SimPEG/{MT => NSEM}/Problem2D/Probs.py (100%)
 rename SimPEG/{MT => NSEM}/Problem2D/__init__.py (100%)
 rename SimPEG/{MT => NSEM}/Problem3D/Probs.py (92%)
 rename SimPEG/{MT => NSEM}/Problem3D/__init__.py (100%)
 rename SimPEG/{MT/SrcMT.py => NSEM/SrcNSEM.py} (91%)
 rename SimPEG/{MT/SurveyMT.py => NSEM/SurveyNSEM.py} (97%)
 rename SimPEG/{MT => NSEM}/Utils/MT1Danalytic.py (100%)
 rename SimPEG/{MT => NSEM}/Utils/MT1Dsolutions.py (100%)
 rename SimPEG/{MT => NSEM}/Utils/__init__.py (100%)
 rename SimPEG/{MT => NSEM}/Utils/dataUtils.py (88%)
 rename SimPEG/{MT => NSEM}/Utils/ediFilesUtils.py (99%)
 rename SimPEG/{MT => NSEM}/Utils/plotDataTypes.py (100%)
 rename SimPEG/{MT => NSEM}/Utils/sourceUtils.py (98%)
 rename SimPEG/{MT => NSEM}/Utils/srcUtils.py (97%)
 create mode 100644 SimPEG/NSEM/__init__.py

diff --git a/SimPEG/MT/__init__.py b/SimPEG/MT/__init__.py
deleted file mode 100644
index b9f5536b..00000000
--- a/SimPEG/MT/__init__.py
+++ /dev/null
@@ -1,5 +0,0 @@
-import Utils
-from SurveyMT import Rx, Survey, Data
-from FieldsMT import Fields1D_e, Fields3D_e
-import Problem1D, Problem2D, Problem3D
-import SrcMT
\ No newline at end of file
diff --git a/SimPEG/MT/FieldsMT.py b/SimPEG/NSEM/FieldsNSEM.py
similarity index 95%
rename from SimPEG/MT/FieldsMT.py
rename to SimPEG/NSEM/FieldsNSEM.py
index 86ac7343..72bf3f6b 100644
--- a/SimPEG/MT/FieldsMT.py
+++ b/SimPEG/NSEM/FieldsNSEM.py
@@ -7,15 +7,15 @@ from SimPEG.EM.Utils import omega
 ##############
 ### Fields ###
 ##############
-class BaseMTFields(Problem.Fields):
-    """Field Storage for a MT survey."""
+class BaseNSEMFields(Problem.Fields):
+    """Field Storage for a NSEM survey."""
     knownFields = {}
     dtype = complex
 
 
-class Fields1D_e(BaseMTFields):
+class Fields1D_e(BaseNSEMFields):
     """
-    Fields storage for the 1D MT solution.
+    Fields storage for the 1D NSEM solution.
     """
     knownFields = {'e_1dSolution':'F'}
     aliasFields = {
@@ -28,7 +28,7 @@ class Fields1D_e(BaseMTFields):
                   }
 
     def __init__(self,mesh,survey,**kwargs):
-        BaseMTFields.__init__(self,mesh,survey,**kwargs)
+        BaseNSEMFields.__init__(self,mesh,survey,**kwargs)
 
     def _ePrimary(self, eSolution, srcList):
         ePrimary = np.zeros_like(eSolution)
@@ -99,7 +99,7 @@ class Fields1D_e(BaseMTFields):
         """
         Derivative of the fields object wrt u.
 
-        :param MTsrc src: MT source
+        :param NSEMsrc src: NSEM source
         :param numpy.ndarray v: random vector of f_sol.size
         This function stacks the fields derivatives appropriately
 
@@ -120,9 +120,9 @@ class Fields1D_e(BaseMTFields):
         """
         return None
 
-class Fields3D_e(BaseMTFields):
+class Fields3D_e(BaseNSEMFields):
     """
-    Fields storage for the 3D MT solution. Labels polarizations by px and py.
+    Fields storage for the 3D NSEM solution. Labels polarizations by px and py.
 
         :param SimPEG object mesh: The solution mesh
         :param SimPEG object survey: A survey object
@@ -147,7 +147,7 @@ class Fields3D_e(BaseMTFields):
                   }
 
     def __init__(self,mesh,survey,**kwargs):
-        BaseMTFields.__init__(self,mesh,survey,**kwargs)
+        BaseNSEMFields.__init__(self,mesh,survey,**kwargs)
 
     def _e_pxPrimary(self, e_pxSolution, srcList):
         e_pxPrimary = np.zeros_like(e_pxSolution)
@@ -228,7 +228,7 @@ class Fields3D_e(BaseMTFields):
         b = (C * e_pxSolution)
         for i, src in enumerate(srcList):
             b[:,i] *= - 1./(1j*omega(src.freq))
-            # There is no magnetic source in the MT problem
+            # There is no magnetic source in the NSEM problem
             # S_m, _ = src.eval(self.survey.prob)
             # if S_m is not None:
             #     b[:,i] += 1./(1j*omega(src.freq)) * S_m
@@ -239,7 +239,7 @@ class Fields3D_e(BaseMTFields):
         b = (C * e_pySolution)
         for i, src in enumerate(srcList):
             b[:,i] *= - 1./(1j*omega(src.freq))
-            # There is no magnetic source in the MT problem
+            # There is no magnetic source in the NSEM problem
             # S_m, _ = src.eval(self.survey.prob)
             # if S_m is not None:
             #     b[:,i] += 1./(1j*omega(src.freq)) * S_m
@@ -302,7 +302,7 @@ class Fields3D_e(BaseMTFields):
         """
         Derivative of the fields object wrt u.
 
-        :param MTsrc src: MT source
+        :param NSEMsrc src: NSEM source
         :param numpy.ndarray v: random vector of f_sol.size
         This function stacks the fields derivatives appropriately
 
@@ -319,7 +319,7 @@ class Fields3D_e(BaseMTFields):
         """
         Derivative of the fields object wrt u.
 
-        :param MTsrc src: MT source
+        :param NSEMsrc src: NSEM source
         :param numpy.ndarray v: random vector of f_sol.size
         This function stacks the fields derivatives appropriately
 
diff --git a/SimPEG/MT/BaseMT.py b/SimPEG/NSEM/NSEM.py
similarity index 91%
rename from SimPEG/MT/BaseMT.py
rename to SimPEG/NSEM/NSEM.py
index c201dfb0..e9915f5d 100644
--- a/SimPEG/MT/BaseMT.py
+++ b/SimPEG/NSEM/NSEM.py
@@ -1,10 +1,10 @@
 from SimPEG import SolverLU as SimpegSolver, PropMaps, Utils, mkvc, sp, np
 from SimPEG.EM.FDEM.FDEM import BaseFDEMProblem
-from SurveyMT import Survey, Data
-from FieldsMT import BaseMTFields
+from SurveyNSEM import Survey, Data
+from FieldsNSEM import BaseNSEMFields
 
 
-class BaseMTProblem(BaseFDEMProblem):
+class BaseNSEMProblem(BaseFDEMProblem):
     """
         Base class for all Natural source problems.
     """
@@ -15,7 +15,7 @@ class BaseMTProblem(BaseFDEMProblem):
     # Set the default pairs of the problem
     surveyPair = Survey
     dataPair = Data
-    fieldsPair = BaseMTFields
+    fieldsPair = BaseNSEMFields
 
     # Set the solver
     Solver = SimpegSolver
@@ -33,8 +33,8 @@ class BaseMTProblem(BaseFDEMProblem):
 
             :param numpy.ndarray m (nC, 1) - conductive model
             :param numpy.ndarray v (nC, 1) - random vector
-            :param MTfields object (optional) - MT fields object, if not given it is calculated
-            :rtype: MTdata object
+            :param NSEMfields object (optional) - NSEM fields object, if not given it is calculated
+            :rtype: NSEMdata object
             :return: Data sensitivities wrt m
         """
 
@@ -80,8 +80,8 @@ class BaseMTProblem(BaseFDEMProblem):
 
             :param numpy.ndarray m (nC, 1) - conductive model
             :param numpy.ndarray v (nD, 1) - vector
-            :param MTfields object u (optional) - MT fields object, if not given it is calculated
-            :rtype: MTdata object
+            :param NSEMfields object u (optional) - NSEM fields object, if not given it is calculated
+            :rtype: NSEMdata object
             :return: Data sensitivities wrt m
         """
 
diff --git a/SimPEG/MT/Problem1D/Probs.py b/SimPEG/NSEM/Problem1D/Probs.py
similarity index 93%
rename from SimPEG/MT/Problem1D/Probs.py
rename to SimPEG/NSEM/Problem1D/Probs.py
index eb7aa331..1ac6771d 100644
--- a/SimPEG/MT/Problem1D/Probs.py
+++ b/SimPEG/NSEM/Problem1D/Probs.py
@@ -1,17 +1,17 @@
 from SimPEG.EM.Utils import omega
 from SimPEG import mkvc
 from scipy.constants import mu_0
-from SimPEG.MT.BaseMT import BaseMTProblem
-from SimPEG.MT.SurveyMT import Survey, Data
-from SimPEG.MT.FieldsMT import Fields1D_e
-from SimPEG.MT.Utils.MT1Danalytic import getEHfields
+from SimPEG.NSEM.NSEM import BaseNSEMProblem
+from SimPEG.NSEM.SurveyNSEM import Survey, Data
+from SimPEG.NSEM.FieldsNSEM import Fields1D_e
+from SimPEG.NSEM.Utils.MT1Danalytic import getEHfields
 import numpy as np
 import multiprocessing, sys, time
 
 
-class eForm_psField(BaseMTProblem):
+class eForm_psField(BaseNSEMProblem):
     """
-    A MT problem soving a e formulation and primary/secondary fields decomposion.
+    A NSEM problem soving a e formulation and primary/secondary fields decomposion.
 
     By eliminating the magnetic flux density using
 
@@ -30,14 +30,14 @@ class eForm_psField(BaseMTProblem):
 
 
     """
-    # From FDEMproblem: Used to project the fields. Currently not used for MTproblem.
+    # From FDEMproblem: Used to project the fields. Currently not used for NSEMproblem.
     _fieldType = 'e_1d'
     _eqLocs    = 'EF'
     _sigmaPrimary = None
 
 
     def __init__(self, mesh, **kwargs):
-        BaseMTProblem.__init__(self, mesh, **kwargs)
+        BaseNSEMProblem.__init__(self, mesh, **kwargs)
         self.fieldsPair = Fields1D_e
         # self._sigmaPrimary = sigmaPrimary
     @property
@@ -167,9 +167,9 @@ class eForm_psField(BaseMTProblem):
 # Fields class corresponding to the fields
 # Update Jvec and Jtvec to include all the derivatives components
 # Other things ...
-class eForm_TotalField(BaseMTProblem):
+class eForm_TotalField(BaseNSEMProblem):
     """
-    A MT problem solving a e formulation and a Total bondary domain decompostion.
+    A NSEM problem solving a e formulation and a Total bondary domain decompostion.
 
     Solves the equation:
 
@@ -178,13 +178,13 @@ class eForm_TotalField(BaseMTProblem):
 
     """
 
-    # From FDEMproblem: Used to project the fields. Currently not used for MTproblem.
+    # From FDEMproblem: Used to project the fields. Currently not used for NSEMproblem.
     _fieldType = 'e'
     _eqLocs    = 'EF'
 
 
     def __init__(self, mesh, **kwargs):
-        BaseMTProblem.__init__(self, mesh, **kwargs)
+        BaseNSEMProblem.__init__(self, mesh, **kwargs)
     @property
     def MeMui(self):
         """
diff --git a/SimPEG/MT/Problem1D/__init__.py b/SimPEG/NSEM/Problem1D/__init__.py
similarity index 100%
rename from SimPEG/MT/Problem1D/__init__.py
rename to SimPEG/NSEM/Problem1D/__init__.py
diff --git a/SimPEG/MT/Problem2D/Probs.py b/SimPEG/NSEM/Problem2D/Probs.py
similarity index 100%
rename from SimPEG/MT/Problem2D/Probs.py
rename to SimPEG/NSEM/Problem2D/Probs.py
diff --git a/SimPEG/MT/Problem2D/__init__.py b/SimPEG/NSEM/Problem2D/__init__.py
similarity index 100%
rename from SimPEG/MT/Problem2D/__init__.py
rename to SimPEG/NSEM/Problem2D/__init__.py
diff --git a/SimPEG/MT/Problem3D/Probs.py b/SimPEG/NSEM/Problem3D/Probs.py
similarity index 92%
rename from SimPEG/MT/Problem3D/Probs.py
rename to SimPEG/NSEM/Problem3D/Probs.py
index 512362b8..cb217abe 100644
--- a/SimPEG/MT/Problem3D/Probs.py
+++ b/SimPEG/NSEM/Problem3D/Probs.py
@@ -1,16 +1,16 @@
 from SimPEG import Survey, Problem, Utils, Models, np, sp, mkvc, SolverLU as SimpegSolver
 from SimPEG.EM.Utils import omega
 from scipy.constants import mu_0
-from SimPEG.MT.BaseMT import BaseMTProblem
-from SimPEG.MT.SurveyMT import Survey, Data
-from SimPEG.MT.FieldsMT import Fields3D_e
+from SimPEG.NSEM.NSEM import BaseNSEMProblem
+from SimPEG.NSEM.SurveyNSEM import Survey, Data
+from SimPEG.NSEM.FieldsNSEM import Fields3D_e
 import multiprocessing, sys, time
 
 
 
-class eForm_ps(BaseMTProblem):
+class eForm_ps(BaseNSEMProblem):
     """
-    A MT problem solving a e formulation and a primary/secondary fields decompostion.
+    A NSEM problem solving a e formulation and a primary/secondary fields decompostion.
 
     By eliminating the magnetic flux density using
 
@@ -29,14 +29,14 @@ class eForm_ps(BaseMTProblem):
 
     """
 
-    # From FDEMproblem: Used to project the fields. Currently not used for MTproblem.
+    # From FDEMproblem: Used to project the fields. Currently not used for NSEMproblem.
     _fieldType = 'e'
     _eqLocs    = 'FE'
     fieldsPair = Fields3D_e
     _sigmaPrimary = None
 
     def __init__(self, mesh, **kwargs):
-        BaseMTProblem.__init__(self, mesh, **kwargs)
+        BaseNSEMProblem.__init__(self, mesh, **kwargs)
 
     @property
     def sigmaPrimary(self):
diff --git a/SimPEG/MT/Problem3D/__init__.py b/SimPEG/NSEM/Problem3D/__init__.py
similarity index 100%
rename from SimPEG/MT/Problem3D/__init__.py
rename to SimPEG/NSEM/Problem3D/__init__.py
diff --git a/SimPEG/MT/SrcMT.py b/SimPEG/NSEM/SrcNSEM.py
similarity index 91%
rename from SimPEG/MT/SrcMT.py
rename to SimPEG/NSEM/SrcNSEM.py
index 70698c88..3365fc6c 100644
--- a/SimPEG/MT/SrcMT.py
+++ b/SimPEG/NSEM/SrcNSEM.py
@@ -11,9 +11,9 @@ import sys
 ###   Sources ###
 #################
 
-class BaseMTSrc(FDEMBaseSrc):
+class BaseNSEMSrc(FDEMBaseSrc):
     '''
-    Sources for the MT problem.
+    Sources for the NSEM problem.
     Use the SimPEG BaseSrc, since the source fields share properties with the transmitters.
 
     :param float freq: The frequency of the source
@@ -29,28 +29,28 @@ class BaseMTSrc(FDEMBaseSrc):
         FDEMBaseSrc.__init__(self, rxList)
 
 # 1D sources
-class polxy_1DhomotD(BaseMTSrc):
+class polxy_1DhomotD(BaseNSEMSrc):
     """
-    MT source for both polarizations (x and y) for the total Domain.
+    NSEM source for both polarizations (x and y) for the total Domain.
 
     It calculates fields calculated based on conditions on the boundary of the domain.
     """
     def __init__(self, rxList, freq):
-        BaseMTSrc.__init__(self, rxList, freq)
+        BaseNSEMSrc.__init__(self, rxList, freq)
 
 
     # TODO: need to add the  primary fields calc and source terms into the problem.
 
 # Need to implement such that it works for all dims.
-class polxy_1Dprimary(BaseMTSrc):
+class polxy_1Dprimary(BaseNSEMSrc):
     """
-    MT source for both polarizations (x and y) given a 1D primary models.
+    NSEM source for both polarizations (x and y) given a 1D primary models.
     It assigns fields calculated from the 1D model as fields in the full space of the problem.
     """
     def __init__(self, rxList, freq):
         # assert mkvc(self.mesh.hz.shape,1) == mkvc(sigma1d.shape,1),'The number of values in the 1D background model does not match the number of vertical cells (hz).'
         self.sigma1d = None
-        BaseMTSrc.__init__(self, rxList, freq)
+        BaseNSEMSrc.__init__(self, rxList, freq)
         # Hidden property of the ePrimary
         self._ePrimary = None
 
@@ -86,7 +86,7 @@ class polxy_1Dprimary(BaseMTSrc):
         Get the electrical field source
         """
         e_p = self.ePrimary(problem)
-        Map_sigma_p = Maps.Vertical1DMap(problem.mesh)
+        Map_sigma_p = Maps.SurjectVertical1D(problem.mesh)
         sigma_p = Map_sigma_p._transform(self.sigma1d)
         # Make mass matrix
         # Note: M(sig) - M(sig_p) = M(sig - sig_p)
@@ -128,15 +128,15 @@ class polxy_1Dprimary(BaseMTSrc):
             # v should be nC size
             return MsigmaDeriv * v
 
-class polxy_3Dprimary(BaseMTSrc):
+class polxy_3Dprimary(BaseNSEMSrc):
     """
-    MT source for both polarizations (x and y) given a 3D primary model. It assigns fields calculated from the 1D model
+    NSEM source for both polarizations (x and y) given a 3D primary model. It assigns fields calculated from the 1D model
     as fields in the full space of the problem.
     """
     def __init__(self, rxList, freq):
         # assert mkvc(self.mesh.hz.shape,1) == mkvc(sigma1d.shape,1),'The number of values in the 1D background model does not match the number of vertical cells (hz).'
         self.sigmaPrimary = None
-        BaseMTSrc.__init__(self, rxList, freq)
+        BaseNSEMSrc.__init__(self, rxList, freq)
         # Hidden property of the ePrimary
         self._ePrimary = None
 
@@ -163,7 +163,7 @@ class polxy_3Dprimary(BaseMTSrc):
         Get the electrical field source
         """
         e_p = self.ePrimary(problem)
-        Map_sigma_p = Maps.Vertical1DMap(problem.mesh)
+        Map_sigma_p = Maps.SurjectVertical1D(problem.mesh)
         sigma_p = Map_sigma_p._transform(self.sigma1d)
         # Make mass matrix
         # Note: M(sig) - M(sig_p) = M(sig - sig_p)
diff --git a/SimPEG/MT/SurveyMT.py b/SimPEG/NSEM/SurveyNSEM.py
similarity index 97%
rename from SimPEG/MT/SurveyMT.py
rename to SimPEG/NSEM/SurveyNSEM.py
index 0ec91a0e..29026eb8 100644
--- a/SimPEG/MT/SurveyMT.py
+++ b/SimPEG/NSEM/SurveyNSEM.py
@@ -4,7 +4,7 @@ from SimPEG.EM.Utils import omega
 from scipy.constants import mu_0
 from numpy.lib import recfunctions as recFunc
 from Utils import rec2ndarr
-import SrcMT
+import SrcNSEM
 import sys
 
 #################
@@ -63,9 +63,9 @@ class Rx(SimPEGsurvey.BaseRx):
         '''
         Project the fields to natural source data.
 
-            :param SrcMT src: The source of the fields to project
+            :param SrcNSEM src: The source of the fields to project
             :param SimPEG.Mesh mesh:
-            :param FieldsMT f: Natural source fields object to project
+            :param FieldsNSEM f: Natural source fields object to project
         '''
 
         ## NOTE: Assumes that e is on t
@@ -143,9 +143,9 @@ class Rx(SimPEGsurvey.BaseRx):
         """
         The derivative of the projection wrt u
 
-        :param MTsrc src: MT source
+        :param NSEMsrc src: NSEM source
         :param TensorMesh mesh: Mesh defining the topology of the problem
-        :param MTfields f: MT fields object of the source
+        :param NSEMfields f: NSEM fields object of the source
         :param numpy.ndarray v: Random vector of size
         """
 
@@ -390,12 +390,12 @@ class Rx(SimPEGsurvey.BaseRx):
 #################
 class Survey(SimPEGsurvey.BaseSurvey):
     """
-        Survey class for MT. Contains all the sources associated with the survey.
+        Survey class for NSEM. Contains all the sources associated with the survey.
 
         :param list srcList: List of sources associated with the survey
 
     """
-    srcPair = SrcMT.BaseMTSrc
+    srcPair = SrcNSEM.BaseNSEMSrc
 
     def __init__(self, srcList, **kwargs):
         # Sort these by frequency
@@ -443,7 +443,7 @@ class Survey(SimPEGsurvey.BaseSurvey):
 #################
 class Data(SimPEGsurvey.Data):
     '''
-    Data class for MTdata. Stores the data vector indexed by the survey.
+    Data class for NSEMdata. Stores the data vector indexed by the survey.
 
     :param SimPEG survey object survey:
     :param v vector of the data in order matching of the survey
@@ -461,7 +461,7 @@ class Data(SimPEGsurvey.Data):
 
     def toRecArray(self,returnType='RealImag'):
         '''
-        Function that returns a numpy.recarray for a SimpegMT impedance data object.
+        Function that returns a numpy.recarray for a SimpegNSEM impedance data object.
 
         :param str returnType: Switches between returning a rec array where the impedance is split to real and imaginary ('RealImag') or is a complex ('Complex')
 
@@ -483,7 +483,7 @@ class Data(SimPEGsurvey.Data):
                 locs = np.hstack((np.array([[0.0]]),locs))
             tArrRec = np.concatenate((src.freq*np.ones((locs.shape[0],1)),locs,np.nan*np.ones((locs.shape[0],12))),axis=1).view(dtRI)
             # np.array([(src.freq,rx.locs[0,0],rx.locs[0,1],rx.locs[0,2],np.nan ,np.nan ,np.nan ,np.nan ,np.nan ,np.nan ,np.nan ,np.nan ) for rx in src.rxList],dtype=dtRI)
-            # Get the type and the value for the DataMT object as a list
+            # Get the type and the value for the DataNSEM object as a list
             typeList = [[rx.rxType.replace('z1d','zyx'),self[src,rx]] for rx in src.rxList]
             # Insert the values to the temp array
             for nr,(key,val) in enumerate(typeList):
@@ -517,17 +517,17 @@ class Data(SimPEGsurvey.Data):
     @classmethod
     def fromRecArray(cls, recArray, srcType='primary'):
         """
-        Class method that reads in a numpy record array to MTdata object.
+        Class method that reads in a numpy record array to NSEMdata object.
 
         Only imports the impedance data.
 
         """
         if srcType=='primary':
-            src = SrcMT.polxy_1Dprimary
+            src = SrcNSEM.polxy_1Dprimary
         elif srcType=='total':
-            src = SrcMT.polxy_1DhomotD
+            src = SrcNSEM.polxy_1DhomotD
         else:
-            raise NotImplementedError('{:s} is not a valid source type for MTdata')
+            raise NotImplementedError('{:s} is not a valid source type for NSEMdata')
 
         # Find all the frequencies in recArray
         uniFreq = np.unique(recArray['freq'])
diff --git a/SimPEG/MT/Utils/MT1Danalytic.py b/SimPEG/NSEM/Utils/MT1Danalytic.py
similarity index 100%
rename from SimPEG/MT/Utils/MT1Danalytic.py
rename to SimPEG/NSEM/Utils/MT1Danalytic.py
diff --git a/SimPEG/MT/Utils/MT1Dsolutions.py b/SimPEG/NSEM/Utils/MT1Dsolutions.py
similarity index 100%
rename from SimPEG/MT/Utils/MT1Dsolutions.py
rename to SimPEG/NSEM/Utils/MT1Dsolutions.py
diff --git a/SimPEG/MT/Utils/__init__.py b/SimPEG/NSEM/Utils/__init__.py
similarity index 100%
rename from SimPEG/MT/Utils/__init__.py
rename to SimPEG/NSEM/Utils/__init__.py
diff --git a/SimPEG/MT/Utils/dataUtils.py b/SimPEG/NSEM/Utils/dataUtils.py
similarity index 88%
rename from SimPEG/MT/Utils/dataUtils.py
rename to SimPEG/NSEM/Utils/dataUtils.py
index 8df9d5d0..fa97aa9a 100644
--- a/SimPEG/MT/Utils/dataUtils.py
+++ b/SimPEG/NSEM/Utils/dataUtils.py
@@ -5,25 +5,25 @@ import numpy.lib.recfunctions as recFunc
 from scipy.constants import mu_0
 from scipy import interpolate as sciint
 
-def getAppRes(MTdata):
+def getAppRes(NSEMdata):
     # Make impedance
     zList = []
-    for src in MTdata.survey.srcList:
+    for src in NSEMdata.survey.srcList:
         zc = [src.freq]
         for rx in src.rxList:
             if 'i' in rx.rxType:
                 m=1j
             else:
                 m = 1
-            zc.append(m*MTdata[src,rx])
+            zc.append(m*NSEMdata[src,rx])
         zList.append(zc)
     return [appResPhs(zList[i][0],np.sum(zList[i][1:3])) for i in np.arange(len(zList))]
 
-def rotateData(MTdata,rotAngle):
+def rotateData(NSEMdata,rotAngle):
     '''
     Function that rotates clockwist by rotAngle (- negative for a counter-clockwise rotation)
     '''
-    recData = MTdata.toRecArray('Complex')
+    recData = NSEMdata.toRecArray('Complex')
     impData = rec2ndarr(recData[['zxx','zxy','zyx','zyy']],complex)
     # Make the rotation matrix
     # c,s,zxx,zxy,zyx,zyy = sympy.symbols('c,s,zxx,zxy,zyx,zyy')
@@ -40,8 +40,8 @@ def rotateData(MTdata,rotAngle):
     for nr,comp in enumerate(['zxx','zxy','zyx','zyy']):
         outRec[comp] = rotData[:,nr]
 
-    from SimPEG import MT
-    return MT.Data.fromRecArray(outRec)
+    from SimPEG import NSEM
+    return NSEM.Data.fromRecArray(outRec)
 
 
 def appResPhs(freq,z):
@@ -57,10 +57,10 @@ def rec2ndarr(x,dt=float):
     return x.view((dt, len(x.dtype.names)))
 
 def makeAnalyticSolution(mesh,model,elev,freqs):
-    from SimPEG import MT
+    from SimPEG import NSEM
     data1D = []
     for freq in freqs:
-        anaEd, anaEu, anaHd, anaHu = MT.Utils.MT1Danalytic.getEHfields(mesh,model,freq,elev)
+        anaEd, anaEu, anaHd, anaHu = NSEM.Utils.MT1Danalytic.getEHfields(mesh,model,freq,elev)
         anaE = anaEd+anaEu
         anaH = anaHd+anaHu
 
@@ -71,7 +71,7 @@ def makeAnalyticSolution(mesh,model,elev,freqs):
     return dataRec
 
 def plotMT1DModelData(problem,models,symList=None):
-    from SimPEG import MT
+    from SimPEG import NSEM
     # Setup the figure
     fontSize = 15
 
@@ -214,11 +214,11 @@ def printTime():
     import time
     print time.strftime("%a, %d %b %Y %H:%M:%S +0000", time.localtime())
 
-def convert3Dto1Dobject(MTdata,rxType3D='zyx'):
-    from SimPEG import MT
+def convert3Dto1Dobject(NSEMdata,rxType3D='zyx'):
+    from SimPEG import NSEM
     # Find the unique locations
     # Need to find the locations
-    recDataTemp = MTdata.toRecArray()
+    recDataTemp = NSEMdata.toRecArray()
     # Check if survey.std has been assigned.
     ## NEED TO: write this...
     # Calculte and add the DET of the tensor to the recArray
@@ -240,24 +240,24 @@ def convert3Dto1Dobject(MTdata,rxType3D='zyx'):
         # Make the receiver list
         rx1DList = []
         for rxType in ['z1dr','z1di']:
-            rx1DList.append(MT.Rx(simpeg.mkvc(loc,2).T,rxType))
+            rx1DList.append(NSEM.Rx(simpeg.mkvc(loc,2).T,rxType))
         # Source list
         locrecData = recData[np.sqrt(np.sum( (rec2ndarr(recData[['x','y','z']]).data - loc )**2,axis=1)) < 1e-5]
         dat1DList = []
         src1DList = []
         for freq in locrecData['freq']:
-            src1DList.append(MT.SrcMT.src_polxy_1Dprimary(rx1DList,freq))
+            src1DList.append(NSEM.SrcNSEM.src_polxy_1Dprimary(rx1DList,freq))
             for comp  in ['r','i']:
                 dat1DList.append( corr * locrecData[rxType3D+comp][locrecData['freq']== freq].data )
 
         # Make the survey
-        sur1D = MT.Survey(src1DList)
+        sur1D = NSEM.Survey(src1DList)
 
         # Make the data
         dataVec = np.hstack(dat1DList)
-        dat1D = MT.Data(sur1D,dataVec)
+        dat1D = NSEM.Data(sur1D,dataVec)
         sur1D.dobs = dataVec
-        # Need to take MTdata.survey.std and split it as well.
+        # 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)
@@ -265,29 +265,29 @@ def convert3Dto1Dobject(MTdata,rxType3D='zyx'):
     # Return the the list of data.
     return mtData1DList
 
-def resampleMTdataAtFreq(MTdata,freqs):
+def resampleNSEMdataAtFreq(NSEMdata,freqs):
     """
-    Function to resample MTdata at set of frequencies
+    Function to resample NSEMdata at set of frequencies
 
     """
-    from SimPEG import MT
+    from SimPEG import NSEM
     # Make a rec array
-    MTrec = MTdata.toRecArray().data
+    NSEMrec = NSEMdata.toRecArray().data
 
     # Find unique locations
-    uniLoc = np.unique(MTrec[['x','y','z']])
-    uniFreq = MTdata.survey.freqs
+    uniLoc = np.unique(NSEMrec[['x','y','z']])
+    uniFreq = NSEMdata.survey.freqs
     # Get the comps
-    dNames = MTrec.dtype
+    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((rec2ndarr(MTrec[['x','y','z']]) - rec2ndarr(loc))**2,axis=1)) < 1. # Find dist of 1 m accuracy
+        ind = np.sqrt(np.sum((rec2ndarr(NSEMrec[['x','y','z']]) - rec2ndarr(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))*rec2ndarr(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(MTrec[ind]['freq'],MTrec[ind][comp],bounds_error=False)
+            int1d = sciint.interp1d(NSEMrec[ind]['freq'],NSEMrec[ind][comp],bounds_error=False)
             tArrRec[comp] = simpeg.mkvc(int1d(freqs),2)
 
         # Join together
@@ -296,5 +296,5 @@ def resampleMTdataAtFreq(MTdata,freqs):
         except NameError as e:
             outRecArr = tArrRec
 
-    # Make the MTdata and return
-    return MT.Data.fromRecArray(outRecArr)
+    # Make the NSEMdata and return
+    return NSEM.Data.fromRecArray(outRecArr)
diff --git a/SimPEG/MT/Utils/ediFilesUtils.py b/SimPEG/NSEM/Utils/ediFilesUtils.py
similarity index 99%
rename from SimPEG/MT/Utils/ediFilesUtils.py
rename to SimPEG/NSEM/Utils/ediFilesUtils.py
index d8f5e0ad..9922469c 100644
--- a/SimPEG/MT/Utils/ediFilesUtils.py
+++ b/SimPEG/NSEM/Utils/ediFilesUtils.py
@@ -2,7 +2,7 @@
 from SimPEG import mkvc
 from scipy.constants import mu_0
 from numpy.lib import recfunctions as recFunc
-from SimPEG.MT.Utils.dataUtils import rec2ndarr
+from SimPEG.NSEM.Utils.dataUtils import rec2ndarr
 
 # Import modules
 import numpy as np
diff --git a/SimPEG/MT/Utils/plotDataTypes.py b/SimPEG/NSEM/Utils/plotDataTypes.py
similarity index 100%
rename from SimPEG/MT/Utils/plotDataTypes.py
rename to SimPEG/NSEM/Utils/plotDataTypes.py
diff --git a/SimPEG/MT/Utils/sourceUtils.py b/SimPEG/NSEM/Utils/sourceUtils.py
similarity index 98%
rename from SimPEG/MT/Utils/sourceUtils.py
rename to SimPEG/NSEM/Utils/sourceUtils.py
index 24b36bfe..fc8e6d81 100644
--- a/SimPEG/MT/Utils/sourceUtils.py
+++ b/SimPEG/NSEM/Utils/sourceUtils.py
@@ -12,7 +12,7 @@ def homo1DModelSource(mesh,freq,sigma_1d):
 
     '''
     # import
-    from SimPEG.MT.Utils import get1DEfields
+    from SimPEG.NSEM.Utils import get1DEfields
     # Get a 1d solution for a halfspace background
     if mesh.dim == 1:
         mesh1d = mesh
@@ -77,7 +77,7 @@ def analytic1DModelSource(mesh,freq,sigma_1d):
 
     '''
     # import
-    from SimPEG.MT.Utils import getEHfields
+    from SimPEG.NSEM.Utils import getEHfields
     # Get a 1d solution for a halfspace background
     if mesh.dim == 1:
         mesh1d = mesh
diff --git a/SimPEG/MT/Utils/srcUtils.py b/SimPEG/NSEM/Utils/srcUtils.py
similarity index 97%
rename from SimPEG/MT/Utils/srcUtils.py
rename to SimPEG/NSEM/Utils/srcUtils.py
index a98cc5b2..9bf141d8 100644
--- a/SimPEG/MT/Utils/srcUtils.py
+++ b/SimPEG/NSEM/Utils/srcUtils.py
@@ -15,7 +15,7 @@ def homo1DModelSource(mesh,freq,m_back):
     '''
 
     # import
-    from SimPEG.MT.Utils import get1DEfields
+    from SimPEG.NSEM.Utils import get1DEfields
     # Get a 1d solution for a halfspace background
     mesh1d = simpeg.Mesh.TensorMesh([mesh.hz],np.array([mesh.x0[2]]))
     # Note: Everything is using e^iwt
diff --git a/SimPEG/NSEM/__init__.py b/SimPEG/NSEM/__init__.py
new file mode 100644
index 00000000..98323827
--- /dev/null
+++ b/SimPEG/NSEM/__init__.py
@@ -0,0 +1,5 @@
+import Utils
+from SurveyNSEM import Rx, Survey, Data
+from FieldsNSEM import Fields1D_e, Fields3D_e
+import Problem1D, Problem2D, Problem3D
+import SrcNSEM
\ No newline at end of file
diff --git a/tests/mt/test_ApparentResistivityAnalytic.py b/tests/mt/test_ApparentResistivityAnalytic.py
index 2a3b1ba9..44bc3c23 100644
--- a/tests/mt/test_ApparentResistivityAnalytic.py
+++ b/tests/mt/test_ApparentResistivityAnalytic.py
@@ -1,6 +1,6 @@
 import unittest
 from SimPEG import *
-from SimPEG import MT
+from SimPEG import NSEM
 
 TOL = 1e-6
 
@@ -20,7 +20,7 @@ def appResNorm(sigmaHalf):
     freqs = np.logspace(4,-4,nFreq)
     Z = []
     for freq in freqs:
-        Ed, Eu, Hd, Hu = MT.Utils.getEHfields(m1d,sigma,freq,np.array([200]))
+        Ed, Eu, Hd, Hu = NSEM.Utils.getEHfields(m1d,sigma,freq,np.array([200]))
         Z.append((Ed + Eu)/(Hd + Hu))
 
     Zarr = np.concatenate(Z)
diff --git a/tests/mt/test_Problem1D_againstAnalyticHalfspace.py b/tests/mt/test_Problem1D_againstAnalyticHalfspace.py
index b023476e..1951271b 100644
--- a/tests/mt/test_Problem1D_againstAnalyticHalfspace.py
+++ b/tests/mt/test_Problem1D_againstAnalyticHalfspace.py
@@ -1,6 +1,6 @@
 import unittest
 import SimPEG as simpeg
-from SimPEG import MT
+from SimPEG import NSEM
 from SimPEG.Utils import meshTensor
 import numpy as np
 # Define the tolerances
@@ -28,34 +28,34 @@ def setupSurvey(sigmaHalf,tD=True):
 
     rxList = []
     for rxType in ['z1dr','z1di']:
-        rxList.append(MT.Rx(simpeg.mkvc(np.array([0.0]),2).T,rxType))
+        rxList.append(NSEM.Rx(simpeg.mkvc(np.array([0.0]),2).T,rxType))
     # Source list
     srcList =[]
     if tD:
         for freq in freqs:
-            srcList.append(MT.SrcMT.polxy_1DhomotD(rxList,freq))
+            srcList.append(NSEM.SrcNSEM.polxy_1DhomotD(rxList,freq))
     else:
         for freq in freqs:
-            srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq))
+            srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,freq))
 
-    survey = MT.Survey(srcList)
+    survey = NSEM.Survey(srcList)
     return survey, sigma, m1d
 
-def getAppResPhs(MTdata):
+def getAppResPhs(NSEMdata):
     # Make impedance
     def appResPhs(freq,z):
         app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2
         app_phs = np.arctan2(z.imag,z.real)*(180/np.pi)
         return app_res, app_phs
     zList = []
-    for src in MTdata.survey.srcList:
+    for src in NSEMdata.survey.srcList:
         zc = [src.freq]
         for rx in src.rxList:
             if 'i' in rx.rxType:
                 m=1j
             else:
                 m = 1
-            zc.append(m*MTdata[src,rx])
+            zc.append(m*NSEMdata[src,rx])
         zList.append(zc)
     return [appResPhs(zList[i][0],np.sum(zList[i][1:3])) for i in np.arange(len(zList))]
 
@@ -63,7 +63,7 @@ def appRes_TotalFieldNorm(sigmaHalf):
 
     # Make the survey
     survey, sigma, mesh = setupSurvey(sigmaHalf)
-    problem = MT.Problem1D.eForm_TotalField(mesh)
+    problem = NSEM.Problem1D.eForm_TotalField(mesh)
     problem.pair(survey)
 
     # Get the fields
@@ -81,7 +81,7 @@ def appPhs_TotalFieldNorm(sigmaHalf):
 
     # Make the survey
     survey, sigma, mesh = setupSurvey(sigmaHalf)
-    problem = MT.Problem1D.eForm_TotalField(mesh)
+    problem = NSEM.Problem1D.eForm_TotalField(mesh)
     problem.pair(survey)
 
     # Get the fields
@@ -99,7 +99,7 @@ def appRes_psFieldNorm(sigmaHalf):
 
     # Make the survey
     survey, sigma, mesh = setupSurvey(sigmaHalf,False)
-    problem = MT.Problem1D.eForm_psField(mesh, sigmaPrimary = sigma)
+    problem = NSEM.Problem1D.eForm_psField(mesh, sigmaPrimary = sigma)
     problem.pair(survey)
 
     # Get the fields
@@ -117,7 +117,7 @@ def appPhs_psFieldNorm(sigmaHalf):
 
     # Make the survey
     survey, sigma, mesh = setupSurvey(sigmaHalf,False)
-    problem = MT.Problem1D.eForm_psField(mesh, sigmaPrimary = sigma)
+    problem = NSEM.Problem1D.eForm_psField(mesh, sigmaPrimary = sigma)
     problem.pair(survey)
 
     # Get the fields
diff --git a/tests/mt/test_Problem1D_totalDvsPSvsAnalytic.py b/tests/mt/test_Problem1D_totalDvsPSvsAnalytic.py
index 2a2e82fd..09c68400 100644
--- a/tests/mt/test_Problem1D_totalDvsPSvsAnalytic.py
+++ b/tests/mt/test_Problem1D_totalDvsPSvsAnalytic.py
@@ -1,6 +1,6 @@
 import unittest
 import SimPEG as simpeg
-from SimPEG import MT
+from SimPEG import NSEM
 from SimPEG.Utils import meshTensor
 import numpy as np
 # Define the tolerances
@@ -34,43 +34,43 @@ def setupSurvey(sigmaHalf,tD=True):
 
     rxList = []
     for rxType in ['z1dr','z1di']:
-        rxList.append(MT.Rx(simpeg.mkvc(np.array([0.0]),2).T,rxType))
+        rxList.append(NSEM.Rx(simpeg.mkvc(np.array([0.0]),2).T,rxType))
     # Source list
     srcList =[]
     if tD:
         for freq in freqs:
-            srcList.append(MT.SrcMT.polxy_1DhomotD(rxList,freq))
+            srcList.append(NSEM.SrcNSEM.polxy_1DhomotD(rxList,freq))
     else:
         for freq in freqs:
-            srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq))
+            srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,freq))
 
-    survey = MT.Survey(srcList)
+    survey = NSEM.Survey(srcList)
     return survey, sigma, m1d
 
-def getAppResPhs(MTdata):
+def getAppResPhs(NSEMdata):
     # Make impedance
     def appResPhs(freq,z):
         app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2
         app_phs = np.arctan2(z.imag,z.real)*(180/np.pi)
         return app_res, app_phs
     zList = []
-    for src in MTdata.survey.srcList:
+    for src in NSEMdata.survey.srcList:
         zc = [src.freq]
         for rx in src.rxList:
             if 'i' in rx.rxType:
                 m=1j
             else:
                 m = 1
-            zc.append(m*MTdata[src,rx])
+            zc.append(m*NSEMdata[src,rx])
         zList.append(zc)
     return [appResPhs(zList[i][0],np.sum(zList[i][1:3])) for i in np.arange(len(zList))]
 
 def calculateAnalyticSolution(srcList,mesh,model):
-    surveyAna = MT.Survey(srcList)
-    data1D = MT.Data(surveyAna)
+    surveyAna = NSEM.Survey(srcList)
+    data1D = NSEM.Data(surveyAna)
     for src in surveyAna.srcList:
         elev = src.rxList[0].locs[0]
-        anaEd, anaEu, anaHd, anaHu = MT.Utils.MT1Danalytic.getEHfields(mesh,model,src.freq,elev)
+        anaEd, anaEu, anaHd, anaHu = NSEM.Utils.NSEM1Danalytic.getEHfields(mesh,model,src.freq,elev)
         anaE = anaEd+anaEu
         anaH = anaHd+anaHu
         # Scale the solution
@@ -87,11 +87,11 @@ def dataMis_AnalyticTotalDomain(sigmaHalf):
 
     # Total domain solution
     surveyTD, sigma, mesh = setupSurvey(sigmaHalf)
-    problemTD = MT.Problem1D.eForm_TotalField(mesh)
+    problemTD = NSEM.Problem1D.eForm_TotalField(mesh)
     problemTD.pair(surveyTD)
     # Analytic data
     dataAnaObj = calculateAnalyticSolution(surveyTD.srcList,mesh,sigma)
-    # dataTDObj = MT.DataMT.DataMT(surveyTD, surveyTD.dpred(sigma))
+    # dataTDObj = NSEM.DataNSEM.DataNSEM(surveyTD, surveyTD.dpred(sigma))
     dataTD = surveyTD.dpred(sigma)
     dataAna = simpeg.mkvc(dataAnaObj)
     return np.all((dataTD - dataAna)/dataAna < 2.)
@@ -109,7 +109,7 @@ def dataMis_AnalyticPrimarySecondary(sigmaHalf):
     # Make the survey
     # Primary secondary
     surveyPS, sigmaPS, mesh = setupSurvey(sigmaHalf,tD=False)
-    problemPS = MT.Problem1D.eForm_psField(mesh)
+    problemPS = NSEM.Problem1D.eForm_psField(mesh)
     problemPS.sigmaPrimary = sigmaPS
     problemPS.pair(surveyPS)
     # Analytic data
diff --git a/tests/mt/test_Problem3D_againstAnalytic.py b/tests/mt/test_Problem3D_againstAnalytic.py
index 602dbaeb..9ad11a33 100644
--- a/tests/mt/test_Problem3D_againstAnalytic.py
+++ b/tests/mt/test_Problem3D_againstAnalytic.py
@@ -3,7 +3,7 @@ from glob import glob
 import numpy as np, sys, os, time, scipy, subprocess
 import SimPEG as simpeg
 import unittest
-from SimPEG import MT
+from SimPEG import NSEM
 from SimPEG.Utils import meshTensor
 from scipy.constants import mu_0
 
@@ -98,40 +98,40 @@ def twoLayer(conds):
 
 
 
-def setupSimpegMTfwd_eForm_ps(inputSetup,comp='Imp',singleFreq=False,expMap=True):
+def setupSimpegNSEMfwd_eForm_ps(inputSetup,comp='Imp',singleFreq=False,expMap=True):
     M,freqs,sig,sigBG,rx_loc = inputSetup
     # Make a receiver list
     rxList = []
     if comp == 'All':
         for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi','tzxr','tzxi','tzyr','tzyi']:
-            rxList.append(MT.Rx(rx_loc,rxType))
+            rxList.append(NSEM.Rx(rx_loc,rxType))
     elif comp == 'Imp':
         for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi']:
-            rxList.append(MT.Rx(rx_loc,rxType))
+            rxList.append(NSEM.Rx(rx_loc,rxType))
     elif comp == 'Tip':
         for rxType in ['tzxr','tzxi','tzyr','tzyi']:
-            rxList.append(MT.Rx(rx_loc,rxType))
+            rxList.append(NSEM.Rx(rx_loc,rxType))
     else:
-        rxList.append(MT.Rx(rx_loc,comp))
+        rxList.append(NSEM.Rx(rx_loc,comp))
     # Source list
     srcList =[]
 
     if singleFreq:
-        srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,singleFreq))
+        srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,singleFreq))
     else:
         for freq in freqs:
-            srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq))
-    # Survey MT
-    survey = MT.Survey(srcList)
+            srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,freq))
+    # Survey NSEM
+    survey = NSEM.Survey(srcList)
 
     ## Setup the problem object
     sigma1d = M.r(sigBG,'CC','CC','M')[0,0,:]
     if expMap:
-        problem = MT.Problem3D.eForm_ps(M,sigmaPrimary= np.log(sigma1d) )
+        problem = NSEM.Problem3D.eForm_ps(M,sigmaPrimary= np.log(sigma1d) )
         problem.mapping = simpeg.Maps.ExpMap(problem.mesh)
         problem.curModel = np.log(sig)
     else:
-        problem = MT.Problem3D.eForm_ps(M,sigmaPrimary= sigma1d)
+        problem = NSEM.Problem3D.eForm_ps(M,sigmaPrimary= sigma1d)
         problem.curModel = sig
     problem.pair(survey)
     problem.verbose = False
@@ -143,18 +143,18 @@ def setupSimpegMTfwd_eForm_ps(inputSetup,comp='Imp',singleFreq=False,expMap=True
 
     return (survey, problem)
 
-def getAppResPhs(MTdata):
+def getAppResPhs(NSEMdata):
     # Make impedance
     def appResPhs(freq,z):
         app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2
         app_phs = np.arctan2(z.imag,z.real)*(180/np.pi)
         return app_res, app_phs
-    recData = MTdata.toRecArray('Complex')
+    recData = NSEMdata.toRecArray('Complex')
     return appResPhs(recData['freq'],recData['zxy']), appResPhs(recData['freq'],recData['zyx'])
 
 def JvecAdjointTest(inputSetup,comp='All',freq=False):
     (M, freqs, sig, sigBG, rx_loc) = inputSetup
-    survey, problem = setupSimpegMTfwd_eForm_ps(inputSetup,comp='All',singleFreq=freq)
+    survey, problem = setupSimpegNSEMfwd_eForm_ps(inputSetup,comp='All',singleFreq=freq)
     print 'Adjoint test of eForm primary/secondary for {:s} comp at {:s}\n'.format(comp,str(survey.freqs))
 
     m  = sig
@@ -174,7 +174,7 @@ def JvecAdjointTest(inputSetup,comp='All',freq=False):
 # Test the Jvec derivative
 def DerivJvecTest(inputSetup,comp='All',freq=False,expMap=True):
     (M, freqs, sig, sigBG, rx_loc) = inputSetup
-    survey, problem = setupSimpegMTfwd_eForm_ps(inputSetup,comp=comp,singleFreq=freq,expMap=expMap)
+    survey, problem = setupSimpegNSEMfwd_eForm_ps(inputSetup,comp=comp,singleFreq=freq,expMap=expMap)
     print 'Derivative test of Jvec for eForm primary/secondary for {:s} comp at {:s}\n'.format(comp,survey.freqs)
     # problem.mapping = simpeg.Maps.ExpMap(problem.mesh)
     # problem.sigmaPrimary = np.log(sigBG)
@@ -191,7 +191,7 @@ def DerivJvecTest(inputSetup,comp='All',freq=False,expMap=True):
 
 def DerivProjfieldsTest(inputSetup,comp='All',freq=False):
 
-    survey, problem = setupSimpegMTfwd_eForm_ps(inputSetup,comp,freq)
+    survey, problem = setupSimpegNSEMfwd_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
@@ -220,7 +220,7 @@ def appResPhsHalfspace_eFrom_ps_Norm(sigmaHalf,appR=True,expMap=False):
     else:
         label = 'phase'
     # Make the survey and the problem
-    survey, problem = setupSimpegMTfwd_eForm_ps(halfSpace(sigmaHalf),expMap=expMap)
+    survey, problem = setupSimpegNSEMfwd_eForm_ps(halfSpace(sigmaHalf),expMap=expMap)
     print 'Apperent {:s} test of eFormulation primary/secondary at {:g}\n\n'.format(label,sigmaHalf)
 
     data = problem.dataPair(survey,survey.dpred(problem.curModel))

From 6165e619edc8cfdefcfc3048603be7887dcf81f0 Mon Sep 17 00:00:00 2001
From: GudniRos <gkrosen@gmail.com>
Date: Tue, 10 May 2016 16:07:45 -0700
Subject: [PATCH 05/14] Moving Problem1D/3D folders to NSEM file Looking at
 making Jvec more general.

---
 SimPEG/NSEM/FieldsNSEM.py         | 130 ++++++++-
 SimPEG/NSEM/NSEM.py               | 438 +++++++++++++++++++++++++++++-
 SimPEG/NSEM/Problem1D/Probs.py    | 291 --------------------
 SimPEG/NSEM/Problem1D/__init__.py |   1 -
 SimPEG/NSEM/Problem2D/Probs.py    |   0
 SimPEG/NSEM/Problem2D/__init__.py |   1 -
 SimPEG/NSEM/Problem3D/Probs.py    | 138 ----------
 SimPEG/NSEM/Problem3D/__init__.py |   1 -
 8 files changed, 559 insertions(+), 441 deletions(-)
 delete mode 100644 SimPEG/NSEM/Problem1D/Probs.py
 delete mode 100644 SimPEG/NSEM/Problem1D/__init__.py
 delete mode 100644 SimPEG/NSEM/Problem2D/Probs.py
 delete mode 100644 SimPEG/NSEM/Problem2D/__init__.py
 delete mode 100644 SimPEG/NSEM/Problem3D/Probs.py
 delete mode 100644 SimPEG/NSEM/Problem3D/__init__.py

diff --git a/SimPEG/NSEM/FieldsNSEM.py b/SimPEG/NSEM/FieldsNSEM.py
index 72bf3f6b..a13a0e0f 100644
--- a/SimPEG/NSEM/FieldsNSEM.py
+++ b/SimPEG/NSEM/FieldsNSEM.py
@@ -4,6 +4,7 @@ import sys
 from numpy.lib import recfunctions as recFunc
 from SimPEG.EM.Utils import omega
 
+
 ##############
 ### Fields ###
 ##############
@@ -12,8 +13,10 @@ class BaseNSEMFields(Problem.Fields):
     knownFields = {}
     dtype = complex
 
-
-class Fields1D_e(BaseNSEMFields):
+###########
+# 1D Fields
+###########
+class Fields1D_ePrimSec(BaseNSEMFields):
     """
     Fields storage for the 1D NSEM solution.
     """
@@ -44,6 +47,118 @@ class Fields1D_e(BaseNSEMFields):
     def _e(self, eSolution, srcList):
         return self._ePrimary(eSolution,srcList) + self._eSecondary(eSolution,srcList)
 
+    def _eDeriv_u(self, src, du_dm_v, adjoint = False):
+
+
+        return Utils.Identity()*du_dm_v
+
+    def _eDeriv_m(self, src, v, adjoint = False):
+        # assuming primary does not depend on the model
+        return Utils.Zero()
+
+    def _bPrimary(self, eSolution, srcList):
+        bPrimary = np.zeros([self.survey.mesh.nE,eSolution.shape[1]], dtype = complex)
+        for i, src in enumerate(srcList):
+            bp = src.bPrimary(self.survey.prob)
+            if bp is not None:
+                bPrimary[:,i] += bp[:,-1]
+        return bPrimary
+
+    def _bSecondary(self, eSolution, srcList):
+        C = self.mesh.nodalGrad
+        b = (C * eSolution)
+        for i, src in enumerate(srcList):
+            b[:,i] *= - 1./(1j*omega(src.freq))
+            # There is no magnetic source in the MT problem
+            # S_m, _ = src.eval(self.survey.prob)
+            # if S_m is not None:
+            #     b[:,i] += 1./(1j*omega(src.freq)) * S_m
+        return b
+
+    def _b(self, eSolution, srcList):
+        return self._bPrimary(eSolution, srcList) + self._bSecondary(eSolution, srcList)
+
+    def _bSecondaryDeriv_u(self, src, v, adjoint = False):
+        C = self.mesh.nodalGrad
+        if adjoint:
+            return - 1./(1j*omega(src.freq)) * (C.T * v)
+        return - 1./(1j*omega(src.freq)) * (C * v)
+
+    def _bSecondaryDeriv_m(self, src, v, adjoint = False):
+        # Doesn't depend on m
+        # _, S_eDeriv = src.evalDeriv(self.survey.prob, adjoint)
+        # S_eDeriv = S_eDeriv(v)
+        # if S_eDeriv is not None:
+        #     return 1./(1j * omega(src.freq)) * S_eDeriv
+        return None
+
+    def _bDeriv_u(self, src, v, adjoint=False):
+        # Primary does not depend on u
+        return self._bSecondaryDeriv_u(src, v, adjoint)
+
+    def _bDeriv_m(self, src, v, adjoint=False):
+        # Assuming the primary does not depend on the model
+        return self._bSecondaryDeriv_m(src, v, adjoint)
+
+    def _fDeriv_u(self, src, v, adjoint=False):
+        """
+        Derivative of the fields object wrt u.
+
+        :param NSEMsrc src: NSEM source
+        :param numpy.ndarray v: random vector of f_sol.size
+        This function stacks the fields derivatives appropriately
+
+        return a vector of size (nreEle+nrbEle)
+        """
+
+        de_du = v #Utils.spdiag(np.ones((self.nF,)))
+        db_du = self._bDeriv_u(src, v, adjoint)
+        # Return the stack
+        # This doesn't work...
+        return np.vstack((de_du,db_du))
+
+    def _fDeriv_m(self, src, v, adjoint=False):
+        """
+        Derivative of the fields object wrt m.
+
+        This function stacks the fields derivatives appropriately
+        """
+        return None
+
+
+class Fields1D_eTotal(BaseNSEMFields):
+    """
+    Fields storage for the 1D NSEM solution solved with for a total domain formulation.
+
+    Used in conjuction with Problem1D_eTotal.
+    """
+    knownFields = {'e_1dSolution':'F'}
+    aliasFields = {
+                    'e_1d' : ['e_1dSolution','F','_e'],
+                    'e_1dPrimary' : ['e_1dSolution','F','_ePrimary'],
+                    'e_1dSecondary' : ['e_1dSolution','F','_eSecondary'],
+                    'b_1d' : ['e_1dSolution','E','_b'],
+                    'b_1dPrimary' : ['e_1dSolution','E','_bPrimary'],
+                    'b_1dSecondary' : ['e_1dSolution','E','_bSecondary']
+                  }
+
+    def __init__(self,mesh,survey,**kwargs):
+        BaseNSEMFields.__init__(self,mesh,survey,**kwargs)
+
+    def _ePrimary(self, eSolution, srcList):
+        ePrimary = np.zeros_like(eSolution)
+        for i, src in enumerate(srcList):
+            ep = src.ePrimary(self.survey.prob)
+            if ep is not None:
+                ePrimary[:,i] = ep[:,-1]
+        return ePrimary
+
+    def _eSecondary(self, eSolution, srcList):
+        return eSolution
+
+    def _e(self, eSolution, srcList):
+        return self._ePrimary(eSolution,srcList) + self._eSecondary(eSolution,srcList)
+
     def _eDeriv_u(self, src, v, adjoint = False):
         return v
 
@@ -120,7 +235,16 @@ class Fields1D_e(BaseNSEMFields):
         """
         return None
 
-class Fields3D_e(BaseNSEMFields):
+
+###########
+# 2D Fields
+###########
+
+
+###########
+# 3D Fields
+###########
+class Fields3D_ePrimSec(BaseNSEMFields):
     """
     Fields storage for the 3D NSEM solution. Labels polarizations by px and py.
 
diff --git a/SimPEG/NSEM/NSEM.py b/SimPEG/NSEM/NSEM.py
index e9915f5d..d0bd8b77 100644
--- a/SimPEG/NSEM/NSEM.py
+++ b/SimPEG/NSEM/NSEM.py
@@ -1,8 +1,9 @@
+from SimPEG.EM.Utils import omega, mu_0
 from SimPEG import SolverLU as SimpegSolver, PropMaps, Utils, mkvc, sp, np
 from SimPEG.EM.FDEM.FDEM import BaseFDEMProblem
 from SurveyNSEM import Survey, Data
-from FieldsNSEM import BaseNSEMFields
-
+from FieldsNSEM import Fields1D_ePrimSec, Fields3D_ePrimSec
+from SimPEG.NSEM.Utils.MT1Danalytic import getEHfields
 
 class BaseNSEMProblem(BaseFDEMProblem):
     """
@@ -56,9 +57,9 @@ class BaseNSEMProblem(BaseFDEMProblem):
                 # We need fDeriv_m = df/du*du/dm + df/dm
                 # Construct du/dm, it requires a solve
                 # NOTE: need to account for the 2 polarizations in the derivatives.
-                f_src = f[src,:]
+                u_src = f[src,:] # u should be a vector by definition. Need to fix this...
                 # dA_dm and dRHS_dm should be of size nE,2, so that we can multiply by dA_duI. The 2 columns are each of the polarizations.
-                dA_dm = self.getADeriv_m(freq, f_src, v) # Size: nE,2 (u_px,u_py) in the columns.
+                dA_dm = self.getADeriv_m(freq, u_src, v) # Size: nE,2 (u_px,u_py) in the columns.
                 dRHS_dm = self.getRHSDeriv_m(freq, v) # Size: nE,2 (u_px,u_py) in the columns.
                 if dRHS_dm is None:
                     du_dm = dA_duI * ( -dA_dm )
@@ -80,7 +81,7 @@ class BaseNSEMProblem(BaseFDEMProblem):
 
             :param numpy.ndarray m (nC, 1) - conductive model
             :param numpy.ndarray v (nD, 1) - vector
-            :param NSEMfields object u (optional) - NSEM fields object, if not given it is calculated
+            :param NSEMfields object f (optional) - NSEM fields object, if not given it is calculated
             :rtype: NSEMdata object
             :return: Data sensitivities wrt m
         """
@@ -103,7 +104,7 @@ class BaseNSEMProblem(BaseFDEMProblem):
 
             for src in self.survey.getSrcByFreq(freq):
                 ftype = self._fieldType + 'Solution'
-                f_src = f[src, :]
+                f_src = f[src, :] # Need to fix this...
 
                 for rx in src.rxList:
                     # Get the adjoint evalDeriv
@@ -130,3 +131,428 @@ class BaseNSEMProblem(BaseFDEMProblem):
             # Clean the factorization, clear memory.
             ATinv.clean()
         return Jtv
+
+###################################
+## 1D problems
+###################################
+
+class Problem1D_ePrimSec(BaseNSEMProblem):
+    """
+    A NSEM problem soving a e formulation and primary/secondary fields decomposion.
+
+    By eliminating the magnetic flux density using
+
+        .. math ::
+
+            \mathbf{b} = \\frac{1}{i \omega}\\left(-\mathbf{C} \mathbf{e} \\right)
+
+
+    we can write Maxwell's equations as a second order system in \\\(\\\mathbf{e}\\\) only:
+
+    .. math ::
+        \\left(\mathbf{C}^T \mathbf{M^e_{\mu^{-1}}} \mathbf{C} + i \omega \mathbf{M^f_\sigma}] \mathbf{e}_{s} =& i \omega \mathbf{M^f_{\delta \sigma}} \mathbf{e}_{p}
+    which we solve for \\\(\\\mathbf{e_s}\\\). The total field \\\mathbf{e}\\ = \\\mathbf{e_p}\\ + \\\mathbf{e_s}\\.
+
+    The primary field is estimated from a background model (commonly half space ).
+
+
+    """
+
+    # From FDEMproblem: Used to project the fields. Currently not used for NSEMproblem.
+    _solutionType = 'e_1dSolution'
+    _formulation  = 'EF'
+    fieldsPair = Fields1D_ePrimSec
+
+    # Initiate properties
+    _sigmaPrimary = None
+
+
+    def __init__(self, mesh, **kwargs):
+        BaseNSEMProblem.__init__(self, mesh, **kwargs)
+        self.fieldsPair = Fields1D_e
+        # self._sigmaPrimary = sigmaPrimary
+    @property
+    def MeMui(self):
+        """
+            Edge inner product matrix
+        """
+        if getattr(self, '_MeMui', None) is None:
+            self._MeMui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
+        return self._MeMui
+
+    @property
+    def MfSigma(self):
+        """
+            Edge inner product matrix
+        """
+        if getattr(self, '_MfSigma', None) is None:
+            self._MfSigma = self.mesh.getFaceInnerProduct(self.curModel.sigma)
+        return self._MfSigma
+
+    @property
+    def sigmaPrimary(self):
+        """
+        A background model, use for the calculation of the primary fields.
+
+        """
+        return self._sigmaPrimary
+
+    @sigmaPrimary.setter
+    def sigmaPrimary(self, val):
+        # Note: TODO add logic for val, make sure it is the correct size.
+        self._sigmaPrimary = val
+
+    def getA(self, freq):
+        """
+            Function to get the A matrix.
+
+            :param float freq: Frequency
+            :rtype: scipy.sparse.csr_matrix
+            :return: A
+        """
+
+        # Note: need to use the code above since in the 1D problem I want
+        # e to live on Faces(nodes) and h on edges(cells). Might need to rethink this
+        # Possible that _fieldType and _eqLocs can fix this
+        MeMui = self.MeMui
+        MfSigma = self.MfSigma
+        C = self.mesh.nodalGrad
+        # Make A
+        A = C.T*MeMui*C + 1j*omega(freq)*MfSigma
+        # Either return full or only the inner part of A
+        return A
+
+    def getADeriv_m(self, freq, u, v, adjoint=False):
+        """
+        The derivative of A wrt sigma
+        """
+
+        dsig_dm = self.curModel.sigmaDeriv
+        MeMui = self.MeMui
+        #
+        u_src = u['e_1dSolution']
+        dMfSigma_dm = self.mesh.getFaceInnerProductDeriv(self.curModel.sigma)(u_src) * self.curModel.sigmaDeriv
+        if adjoint:
+            return 1j * omega(freq) * (  dMfSigma_dm.T * v )
+        # Note: output has to be nN/nF, not nC/nE.
+        # v should be nC
+        return 1j * omega(freq) * ( dMfSigma_dm * v )
+
+    def getRHS(self, freq):
+        """
+            Function to return the right hand side for the system.
+            :param float freq: Frequency
+            :rtype: numpy.ndarray (nF, 1), numpy.ndarray (nF, 1)
+            :return: RHS for 1 polarizations, primary fields
+        """
+
+        # Get sources for the frequncy(polarizations)
+        Src = self.survey.getSrcByFreq(freq)[0]
+        S_e = Src.S_e(self)
+        return -1j * omega(freq) * S_e
+
+    def getRHSDeriv_m(self, freq, v, adjoint=False):
+        """
+        The derivative of the RHS wrt sigma
+        """
+
+        Src = self.survey.getSrcByFreq(freq)[0]
+        S_eDeriv = Src.S_eDeriv_m(self, v, adjoint)
+        return -1j * omega(freq) * S_eDeriv
+
+    def fields(self, m):
+        '''
+        Function to calculate all the fields for the model m.
+
+        :param np.ndarray (nC,) m: Conductivity model
+        '''
+        # Set the current model
+        self.curModel = m
+        # Make the fields object
+        F = self.fieldsPair(self.mesh, self.survey)
+        # Loop over the frequencies
+        for freq in self.survey.freqs:
+            if self.verbose:
+                startTime = time.time()
+                print 'Starting work for {:.3e}'.format(freq)
+                sys.stdout.flush()
+            A = self.getA(freq)
+            rhs  = self.getRHS(freq)
+            Ainv = self.Solver(A, **self.solverOpts)
+            e_s = Ainv * rhs
+
+            # Store the fields
+            Src = self.survey.getSrcByFreq(freq)[0]
+            # NOTE: only store the e_solution(secondary), all other components calculated in the fields object
+            F[Src, 'e_1dSolution'] = e_s[:,-1] # Only storing the yx polarization as 1d
+
+            # Note curl e = -iwb so b = -curl e /iw
+            # b = -( self.mesh.nodalGrad * e )/( 1j*omega(freq) )
+            # F[Src, 'b_1d'] = b[:,1]
+            if self.verbose:
+                print 'Ran for {:f} seconds'.format(time.time()-startTime)
+                sys.stdout.flush()
+        return F
+
+# Note this is not fully functional.
+# Missing:
+# Fields class corresponding to the fields
+# Update Jvec and Jtvec to include all the derivatives components
+# Other things ...
+class Problem1D_eTotal(BaseNSEMProblem):
+    """
+    A NSEM problem solving a e formulation and a Total bondary domain decompostion.
+
+    Solves the equation:
+
+    Math:
+        Have to do this...
+        Not implement correctly.......
+    """
+
+    # From FDEMproblem: Used to project the fields. Currently not used for NSEMproblem.
+    _solutionType = 'e_1dSolution'
+    _formulation  = 'EF'
+    fieldsPair = Fields1D_eTotal
+
+    def __init__(self, mesh, **kwargs):
+        BaseNSEMProblem.__init__(self, mesh, **kwargs)
+    @property
+    def MeMui(self):
+        """
+            Edge inner product matrix
+        """
+        if getattr(self, '_MeMui', None) is None:
+            self._MeMui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
+        return self._MeMui
+
+    @property
+    def MfSigma(self):
+        """
+            Edge inner product matrix
+        """
+        if getattr(self, '_MfSigma', None) is None:
+            self._MfSigma = self.mesh.getFaceInnerProduct(self.curModel.sigma)
+        return self._MfSigma
+
+    def getA(self, freq, full=False):
+        """
+            Function to get the A matrix.
+
+            :param float freq: Frequency
+            :param logic full: Return full A or the inner part
+            :rtype: scipy.sparse.csr_matrix
+            :return: A
+        """
+
+        MeMui = self.MeMui
+        MfSigma = self.MfSigma
+        # Note: need to use the code above since in the 1D problem I want
+        # e to live on Faces(nodes) and h on edges(cells). Might need to rethink this
+        # Possible that _fieldType and _eqLocs can fix this
+        # MeMui = self.MfMui
+        # MfSigma = self.MfSigma
+        C = self.mesh.nodalGrad
+        # Make A
+        A = C.T*MeMui*C + 1j*omega(freq)*MfSigma
+        # Either return full or only the inner part of A
+        if full:
+            return A
+        else:
+            return A[1:-1,1:-1]
+
+    def getADeriv_m(self, freq, u, v, adjoint=False):
+        raise NotImplementedError('getADeriv is not implemented')
+
+    def getRHS(self, freq):
+        """
+            Function to return the right hand side for the system.
+            :param float freq: Frequency
+            :rtype: numpy.ndarray (nE, 2), numpy.ndarray (nE, 2)
+            :return: RHS for both polarizations, primary fields
+        """
+        # Get sources for the frequency
+        # NOTE: Need to use the source information, doesn't really apply in 1D
+        src = self.survey.getSrcByFreq(freq)
+        # Get the full A
+        A = self.getA(freq,full=True)
+        # Define the outer part of the solution matrix
+        Aio = A[1:-1,[0,-1]]
+        Ed, Eu, Hd, Hu = getEHfields(self.mesh,self.curModel.sigma,freq,self.mesh.vectorNx)
+        Etot = (Ed + Eu)
+        sourceAmp = 1.0
+        Etot = ((Etot/Etot[-1])*sourceAmp) # Scale the fields to be equal to sourceAmp at the top
+        ## Note: The analytic solution is derived with e^iwt
+        eBC = np.r_[Etot[0],Etot[-1]]
+        # The right hand side
+
+        return -Aio*eBC, eBC
+
+    def getRHSderiv_m(self, freq, backSigma, u, v, adjoint=False):
+        raise NotImplementedError('getRHSDeriv not implemented yet')
+        return None
+
+    def fields(self, m):
+        '''
+        Function to calculate all the fields for the model m.
+
+        :param np.ndarray (nC,) m: Conductivity model
+        :param np.ndarray (nC,) m_back: Background conductivity model
+        '''
+        self.curModel = m
+        # RHS, CalcFields = self.getRHS(freq,m_back), self.calcFields
+
+        F = Fields1D_eTotal(self.mesh, self.survey)
+        for freq in self.survey.freqs:
+            if self.verbose:
+                startTime = time.time()
+                print 'Starting work for {:.3e}'.format(freq)
+                sys.stdout.flush()
+            A = self.getA(freq)
+            rhs, e_o = self.getRHS(freq)
+            Ainv = self.Solver(A, **self.solverOpts)
+            e_i = Ainv * rhs
+            e = mkvc(np.r_[e_o[0], e_i, e_o[1]],2)
+            # Store the fields
+            Src = self.survey.getSrcByFreq(freq)
+            # NOTE: only store e fields
+            F[Src, 'e_1dSolution'] = e[:,0]
+            if self.verbose:
+                print 'Ran for {:f} seconds'.format(time.time()-startTime)
+                sys.stdout.flush()
+        return F
+
+
+###################################
+## 3D problems
+###################################
+class Problem3D_ePrimSec(BaseNSEMProblem):
+    """
+    A NSEM problem solving a e formulation and a primary/secondary fields decompostion.
+
+    By eliminating the magnetic flux density using
+
+        .. math ::
+
+            \mathbf{b} = \\frac{1}{i \omega}\\left(-\mathbf{C} \mathbf{e} \\right)
+
+
+    we can write Maxwell's equations as a second order system in \\\(\\\mathbf{e}\\\) only:
+
+    .. math ::
+        \\left(\mathbf{C}^T \mathbf{M^f_{\mu^{-1}}} \mathbf{C} + i \omega \mathbf{M^e_\sigma}] \mathbf{e}_{s} =& i \omega \mathbf{M^e_{\delta \sigma}} \mathbf{e}_{p}
+    which we solve for \\\(\\\mathbf{e_s}\\\). The total field \\\mathbf{e}\\ = \\\mathbf{e_p}\\ + \\\mathbf{e_s}\\.
+
+    The primary field is estimated from a background model (commonly as a 1D model).
+
+    """
+
+    # From FDEMproblem: Used to project the fields. Currently not used for NSEMproblem.
+    _solutionType = [ 'e_pxSolution', 'e_pySolution']  # Forces order on the object
+    _formulation  = 'EB'
+    fieldsPair = Fields3D_ePrimSec
+
+    # Initiate properties
+    _sigmaPrimary = None
+
+    def __init__(self, mesh, **kwargs):
+        BaseNSEMProblem.__init__(self, mesh, **kwargs)
+
+    @property
+    def sigmaPrimary(self):
+        """
+        A background model, use for the calculation of the primary fields.
+
+        """
+        return self._sigmaPrimary
+    @sigmaPrimary.setter
+    def sigmaPrimary(self, val):
+        # Note: TODO add logic for val, make sure it is the correct size.
+        self._sigmaPrimary = val
+
+    def getA(self, freq):
+        """
+            Function to get the A system.
+
+            :param float freq: Frequency
+            :rtype: scipy.sparse.csr_matrix
+            :return: A
+        """
+        Mmui = self.MfMui
+        Msig = self.MeSigma
+        C = self.mesh.edgeCurl
+
+        return C.T*Mmui*C + 1j*omega(freq)*Msig
+
+    def getADeriv_m(self, freq, u, v, adjoint=False):
+        """
+        Calculate the derivative of A wrt m.
+
+        """
+        # Fix u to be a matrix nE,2
+        # This considers both polarizations and returns a nE,2 matrix for each polarization
+        if adjoint:
+            dMe_dsigV = sp.hstack(( self.MeSigmaDeriv( u['e_pxSolution'] ).T, self.MeSigmaDeriv(u['e_pySolution'] ).T ))*v
+        else:
+            # Need a nE,2 matrix to be returned
+            dMe_dsigV = np.hstack(( mkvc(self.MeSigmaDeriv( u['e_pxSolution'] )*v,2), mkvc( self.MeSigmaDeriv(u['e_pySolution'] )*v,2) ))
+        return 1j * omega(freq) * dMe_dsigV
+
+
+    def getRHS(self, freq):
+        """
+        Function to return the right hand side for the system.
+
+            :param float freq: Frequency
+            :rtype: numpy.ndarray (nE, 2), numpy.ndarray (nE, 2)
+            :return: RHS for both polarizations, primary fields
+        """
+
+        # Get sources for the frequncy(polarizations)
+        Src = self.survey.getSrcByFreq(freq)[0]
+        S_e = Src.S_e(self)
+        return -1j * omega(freq) * S_e
+
+    def getRHSDeriv_m(self, freq, v, adjoint=False):
+        """
+        The derivative of the RHS with respect to sigma
+        """
+
+        Src = self.survey.getSrcByFreq(freq)[0]
+        S_eDeriv = Src.S_eDeriv_m(self, v, adjoint)
+        return -1j * omega(freq) * S_eDeriv
+
+    def fields(self, m):
+        '''
+        Function to calculate all the fields for the model m.
+
+        :param np.ndarray (nC,) m: Conductivity model
+        '''
+        # Set the current model
+        self.curModel = m
+
+        F = self.fieldsPair(self.mesh, self.survey)
+        for freq in self.survey.freqs:
+            if self.verbose:
+                startTime = time.time()
+                print 'Starting work for {:.3e}'.format(freq)
+                sys.stdout.flush()
+            A = self.getA(freq)
+            rhs  = self.getRHS(freq)
+            # Solve the system
+            Ainv = self.Solver(A, **self.solverOpts)
+            e_s = Ainv * rhs
+
+            # Store the fields
+            Src = self.survey.getSrcByFreq(freq)[0]
+            # Store the fields
+            # Use self._solutionType
+            F[Src, 'e_pxSolution'] = e_s[:,0]
+            F[Src, 'e_pySolution'] = e_s[:,1]
+            # Note curl e = -iwb so b = -curl/iw
+
+            if self.verbose:
+                print 'Ran for {:f} seconds'.format(time.time()-startTime)
+                sys.stdout.flush()
+            Ainv.clean()
+        return F
diff --git a/SimPEG/NSEM/Problem1D/Probs.py b/SimPEG/NSEM/Problem1D/Probs.py
deleted file mode 100644
index 1ac6771d..00000000
--- a/SimPEG/NSEM/Problem1D/Probs.py
+++ /dev/null
@@ -1,291 +0,0 @@
-from SimPEG.EM.Utils import omega
-from SimPEG import mkvc
-from scipy.constants import mu_0
-from SimPEG.NSEM.NSEM import BaseNSEMProblem
-from SimPEG.NSEM.SurveyNSEM import Survey, Data
-from SimPEG.NSEM.FieldsNSEM import Fields1D_e
-from SimPEG.NSEM.Utils.MT1Danalytic import getEHfields
-import numpy as np
-import multiprocessing, sys, time
-
-
-class eForm_psField(BaseNSEMProblem):
-    """
-    A NSEM problem soving a e formulation and primary/secondary fields decomposion.
-
-    By eliminating the magnetic flux density using
-
-        .. math ::
-
-            \mathbf{b} = \\frac{1}{i \omega}\\left(-\mathbf{C} \mathbf{e} \\right)
-
-
-    we can write Maxwell's equations as a second order system in \\\(\\\mathbf{e}\\\) only:
-
-    .. math ::
-        \\left(\mathbf{C}^T \mathbf{M^e_{\mu^{-1}}} \mathbf{C} + i \omega \mathbf{M^f_\sigma}] \mathbf{e}_{s} =& i \omega \mathbf{M^f_{\delta \sigma}} \mathbf{e}_{p}
-    which we solve for \\\(\\\mathbf{e_s}\\\). The total field \\\mathbf{e}\\ = \\\mathbf{e_p}\\ + \\\mathbf{e_s}\\.
-
-    The primary field is estimated from a background model (commonly half space ).
-
-
-    """
-    # From FDEMproblem: Used to project the fields. Currently not used for NSEMproblem.
-    _fieldType = 'e_1d'
-    _eqLocs    = 'EF'
-    _sigmaPrimary = None
-
-
-    def __init__(self, mesh, **kwargs):
-        BaseNSEMProblem.__init__(self, mesh, **kwargs)
-        self.fieldsPair = Fields1D_e
-        # self._sigmaPrimary = sigmaPrimary
-    @property
-    def MeMui(self):
-        """
-            Edge inner product matrix
-        """
-        if getattr(self, '_MeMui', None) is None:
-            self._MeMui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
-        return self._MeMui
-
-    @property
-    def MfSigma(self):
-        """
-            Edge inner product matrix
-        """
-        if getattr(self, '_MfSigma', None) is None:
-            self._MfSigma = self.mesh.getFaceInnerProduct(self.curModel.sigma)
-        return self._MfSigma
-
-    @property
-    def sigmaPrimary(self):
-        """
-        A background model, use for the calculation of the primary fields.
-
-        """
-        return self._sigmaPrimary
-
-    @sigmaPrimary.setter
-    def sigmaPrimary(self, val):
-        # Note: TODO add logic for val, make sure it is the correct size.
-        self._sigmaPrimary = val
-
-    def getA(self, freq):
-        """
-            Function to get the A matrix.
-
-            :param float freq: Frequency
-            :rtype: scipy.sparse.csr_matrix
-            :return: A
-        """
-
-        # Note: need to use the code above since in the 1D problem I want
-        # e to live on Faces(nodes) and h on edges(cells). Might need to rethink this
-        # Possible that _fieldType and _eqLocs can fix this
-        MeMui = self.MeMui
-        MfSigma = self.MfSigma
-        C = self.mesh.nodalGrad
-        # Make A
-        A = C.T*MeMui*C + 1j*omega(freq)*MfSigma
-        # Either return full or only the inner part of A
-        return A
-
-    def getADeriv_m(self, freq, u, v, adjoint=False):
-        """
-        The derivative of A wrt sigma
-        """
-
-        dsig_dm = self.curModel.sigmaDeriv
-        MeMui = self.MeMui
-        #
-        u_src = u['e_1dSolution']
-        dMfSigma_dm = self.mesh.getFaceInnerProductDeriv(self.curModel.sigma)(u_src) * self.curModel.sigmaDeriv
-        if adjoint:
-            return 1j * omega(freq) * (  dMfSigma_dm.T * v )
-        # Note: output has to be nN/nF, not nC/nE.
-        # v should be nC
-        return 1j * omega(freq) * ( dMfSigma_dm * v )
-
-    def getRHS(self, freq):
-        """
-            Function to return the right hand side for the system.
-            :param float freq: Frequency
-            :rtype: numpy.ndarray (nF, 1), numpy.ndarray (nF, 1)
-            :return: RHS for 1 polarizations, primary fields
-        """
-
-        # Get sources for the frequncy(polarizations)
-        Src = self.survey.getSrcByFreq(freq)[0]
-        S_e = Src.S_e(self)
-        return -1j * omega(freq) * S_e
-
-    def getRHSDeriv_m(self, freq, v, adjoint=False):
-        """
-        The derivative of the RHS wrt sigma
-        """
-
-        Src = self.survey.getSrcByFreq(freq)[0]
-        S_eDeriv = Src.S_eDeriv_m(self, v, adjoint)
-        return -1j * omega(freq) * S_eDeriv
-
-    def fields(self, m):
-        '''
-        Function to calculate all the fields for the model m.
-
-        :param np.ndarray (nC,) m: Conductivity model
-        '''
-        # Set the current model
-        self.curModel = m
-
-        F = Fields1D_e(self.mesh, self.survey)
-        for freq in self.survey.freqs:
-            if self.verbose:
-                startTime = time.time()
-                print 'Starting work for {:.3e}'.format(freq)
-                sys.stdout.flush()
-            A = self.getA(freq)
-            rhs  = self.getRHS(freq)
-            Ainv = self.Solver(A, **self.solverOpts)
-            e_s = Ainv * rhs
-
-            # Store the fields
-            Src = self.survey.getSrcByFreq(freq)[0]
-            # NOTE: only store the e_solution(secondary), all other components calculated in the fields object
-            F[Src, 'e_1dSolution'] = e_s[:,-1] # Only storing the yx polarization as 1d
-
-            # Note curl e = -iwb so b = -curl e /iw
-            # b = -( self.mesh.nodalGrad * e )/( 1j*omega(freq) )
-            # F[Src, 'b_1d'] = b[:,1]
-            if self.verbose:
-                print 'Ran for {:f} seconds'.format(time.time()-startTime)
-                sys.stdout.flush()
-        return F
-
-# Note this is not fully functional.
-# Missing:
-# Fields class corresponding to the fields
-# Update Jvec and Jtvec to include all the derivatives components
-# Other things ...
-class eForm_TotalField(BaseNSEMProblem):
-    """
-    A NSEM problem solving a e formulation and a Total bondary domain decompostion.
-
-    Solves the equation:
-
-    Math:
-
-
-    """
-
-    # From FDEMproblem: Used to project the fields. Currently not used for NSEMproblem.
-    _fieldType = 'e'
-    _eqLocs    = 'EF'
-
-
-    def __init__(self, mesh, **kwargs):
-        BaseNSEMProblem.__init__(self, mesh, **kwargs)
-    @property
-    def MeMui(self):
-        """
-            Edge inner product matrix
-        """
-        if getattr(self, '_MeMui', None) is None:
-            self._MeMui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
-        return self._MeMui
-
-    @property
-    def MfSigma(self):
-        """
-            Edge inner product matrix
-        """
-        if getattr(self, '_MfSigma', None) is None:
-            self._MfSigma = self.mesh.getFaceInnerProduct(self.curModel.sigma)
-        return self._MfSigma
-
-    def getA(self, freq, full=False):
-        """
-            Function to get the A matrix.
-
-            :param float freq: Frequency
-            :param logic full: Return full A or the inner part
-            :rtype: scipy.sparse.csr_matrix
-            :return: A
-        """
-
-        MeMui = self.MeMui
-        MfSigma = self.MfSigma
-        # Note: need to use the code above since in the 1D problem I want
-        # e to live on Faces(nodes) and h on edges(cells). Might need to rethink this
-        # Possible that _fieldType and _eqLocs can fix this
-        # MeMui = self.MfMui
-        # MfSigma = self.MfSigma
-        C = self.mesh.nodalGrad
-        # Make A
-        A = C.T*MeMui*C + 1j*omega(freq)*MfSigma
-        # Either return full or only the inner part of A
-        if full:
-            return A
-        else:
-            return A[1:-1,1:-1]
-
-    def getADeriv_m(self, freq, u, v, adjoint=False):
-        raise NotImplementedError('getADeriv is not implemented')
-
-    def getRHS(self, freq):
-        """
-            Function to return the right hand side for the system.
-            :param float freq: Frequency
-            :rtype: numpy.ndarray (nE, 2), numpy.ndarray (nE, 2)
-            :return: RHS for both polarizations, primary fields
-        """
-        # Get sources for the frequency
-        # NOTE: Need to use the source information, doesn't really apply in 1D
-        src = self.survey.getSrcByFreq(freq)
-        # Get the full A
-        A = self.getA(freq,full=True)
-        # Define the outer part of the solution matrix
-        Aio = A[1:-1,[0,-1]]
-        Ed, Eu, Hd, Hu = getEHfields(self.mesh,self.curModel.sigma,freq,self.mesh.vectorNx)
-        Etot = (Ed + Eu)
-        sourceAmp = 1.0
-        Etot = ((Etot/Etot[-1])*sourceAmp) # Scale the fields to be equal to sourceAmp at the top
-        ## Note: The analytic solution is derived with e^iwt
-        eBC = np.r_[Etot[0],Etot[-1]]
-        # The right hand side
-
-        return -Aio*eBC, eBC
-
-    def getRHSderiv_m(self, freq, backSigma, u, v, adjoint=False):
-        raise NotImplementedError('getRHSDeriv not implemented yet')
-        return None
-
-    def fields(self, m):
-        '''
-        Function to calculate all the fields for the model m.
-
-        :param np.ndarray (nC,) m: Conductivity model
-        :param np.ndarray (nC,) m_back: Background conductivity model
-        '''
-        self.curModel = m
-        # RHS, CalcFields = self.getRHS(freq,m_back), self.calcFields
-
-        F = Fields1D_e(self.mesh, self.survey)
-        for freq in self.survey.freqs:
-            if self.verbose:
-                startTime = time.time()
-                print 'Starting work for {:.3e}'.format(freq)
-                sys.stdout.flush()
-            A = self.getA(freq)
-            rhs, e_o = self.getRHS(freq)
-            Ainv = self.Solver(A, **self.solverOpts)
-            e_i = Ainv * rhs
-            e = mkvc(np.r_[e_o[0], e_i, e_o[1]],2)
-            # Store the fields
-            Src = self.survey.getSrcByFreq(freq)
-            # NOTE: only store e fields
-            F[Src, 'e_1dSolution'] = e[:,0]
-            if self.verbose:
-                print 'Ran for {:f} seconds'.format(time.time()-startTime)
-                sys.stdout.flush()
-        return F
diff --git a/SimPEG/NSEM/Problem1D/__init__.py b/SimPEG/NSEM/Problem1D/__init__.py
deleted file mode 100644
index b0a77250..00000000
--- a/SimPEG/NSEM/Problem1D/__init__.py
+++ /dev/null
@@ -1 +0,0 @@
-from Probs import eForm_TotalField, eForm_psField
\ No newline at end of file
diff --git a/SimPEG/NSEM/Problem2D/Probs.py b/SimPEG/NSEM/Problem2D/Probs.py
deleted file mode 100644
index e69de29b..00000000
diff --git a/SimPEG/NSEM/Problem2D/__init__.py b/SimPEG/NSEM/Problem2D/__init__.py
deleted file mode 100644
index fc80254b..00000000
--- a/SimPEG/NSEM/Problem2D/__init__.py
+++ /dev/null
@@ -1 +0,0 @@
-pass
\ No newline at end of file
diff --git a/SimPEG/NSEM/Problem3D/Probs.py b/SimPEG/NSEM/Problem3D/Probs.py
deleted file mode 100644
index cb217abe..00000000
--- a/SimPEG/NSEM/Problem3D/Probs.py
+++ /dev/null
@@ -1,138 +0,0 @@
-from SimPEG import Survey, Problem, Utils, Models, np, sp, mkvc, SolverLU as SimpegSolver
-from SimPEG.EM.Utils import omega
-from scipy.constants import mu_0
-from SimPEG.NSEM.NSEM import BaseNSEMProblem
-from SimPEG.NSEM.SurveyNSEM import Survey, Data
-from SimPEG.NSEM.FieldsNSEM import Fields3D_e
-import multiprocessing, sys, time
-
-
-
-class eForm_ps(BaseNSEMProblem):
-    """
-    A NSEM problem solving a e formulation and a primary/secondary fields decompostion.
-
-    By eliminating the magnetic flux density using
-
-        .. math ::
-
-            \mathbf{b} = \\frac{1}{i \omega}\\left(-\mathbf{C} \mathbf{e} \\right)
-
-
-    we can write Maxwell's equations as a second order system in \\\(\\\mathbf{e}\\\) only:
-
-    .. math ::
-        \\left(\mathbf{C}^T \mathbf{M^f_{\mu^{-1}}} \mathbf{C} + i \omega \mathbf{M^e_\sigma}] \mathbf{e}_{s} =& i \omega \mathbf{M^e_{\delta \sigma}} \mathbf{e}_{p}
-    which we solve for \\\(\\\mathbf{e_s}\\\). The total field \\\mathbf{e}\\ = \\\mathbf{e_p}\\ + \\\mathbf{e_s}\\.
-
-    The primary field is estimated from a background model (commonly as a 1D model).
-
-    """
-
-    # From FDEMproblem: Used to project the fields. Currently not used for NSEMproblem.
-    _fieldType = 'e'
-    _eqLocs    = 'FE'
-    fieldsPair = Fields3D_e
-    _sigmaPrimary = None
-
-    def __init__(self, mesh, **kwargs):
-        BaseNSEMProblem.__init__(self, mesh, **kwargs)
-
-    @property
-    def sigmaPrimary(self):
-        """
-        A background model, use for the calculation of the primary fields.
-
-        """
-        return self._sigmaPrimary
-    @sigmaPrimary.setter
-    def sigmaPrimary(self, val):
-        # Note: TODO add logic for val, make sure it is the correct size.
-        self._sigmaPrimary = val
-
-    def getA(self, freq):
-        """
-            Function to get the A system.
-
-            :param float freq: Frequency
-            :rtype: scipy.sparse.csr_matrix
-            :return: A
-        """
-        Mmui = self.MfMui
-        Msig = self.MeSigma
-        C = self.mesh.edgeCurl
-
-        return C.T*Mmui*C + 1j*omega(freq)*Msig
-
-    def getADeriv_m(self, freq, u, v, adjoint=False):
-        """
-        Calculate the derivative of A wrt m.
-
-        """
-
-        # This considers both polarizations and returns a nE,2 matrix for each polarization
-        if adjoint:
-            dMe_dsigV = sp.hstack(( self.MeSigmaDeriv( u['e_pxSolution'] ).T, self.MeSigmaDeriv(u['e_pySolution'] ).T ))*v
-        else:
-            # Need a nE,2 matrix to be returned
-            dMe_dsigV = np.hstack(( mkvc(self.MeSigmaDeriv( u['e_pxSolution'] )*v,2), mkvc( self.MeSigmaDeriv(u['e_pySolution'] )*v,2) ))
-        return 1j * omega(freq) * dMe_dsigV
-
-
-    def getRHS(self, freq):
-        """
-        Function to return the right hand side for the system.
-
-            :param float freq: Frequency
-            :rtype: numpy.ndarray (nE, 2), numpy.ndarray (nE, 2)
-            :return: RHS for both polarizations, primary fields
-        """
-
-        # Get sources for the frequncy(polarizations)
-        Src = self.survey.getSrcByFreq(freq)[0]
-        S_e = Src.S_e(self)
-        return -1j * omega(freq) * S_e
-
-    def getRHSDeriv_m(self, freq, v, adjoint=False):
-        """
-        The derivative of the RHS with respect to sigma
-        """
-
-        Src = self.survey.getSrcByFreq(freq)[0]
-        S_eDeriv = Src.S_eDeriv_m(self, v, adjoint)
-        return -1j * omega(freq) * S_eDeriv
-
-    def fields(self, m):
-        '''
-        Function to calculate all the fields for the model m.
-
-        :param np.ndarray (nC,) m: Conductivity model
-        '''
-        # Set the current model
-        self.curModel = m
-
-        F = Fields3D_e(self.mesh, self.survey)
-        for freq in self.survey.freqs:
-            if self.verbose:
-                startTime = time.time()
-                print 'Starting work for {:.3e}'.format(freq)
-                sys.stdout.flush()
-            A = self.getA(freq)
-            rhs  = self.getRHS(freq)
-            # Solve the system
-            Ainv = self.Solver(A, **self.solverOpts)
-            e_s = Ainv * rhs
-
-            # Store the fields
-            Src = self.survey.getSrcByFreq(freq)[0]
-            # Store the fieldss
-            F[Src, 'e_pxSolution'] = e_s[:,0]
-            F[Src, 'e_pySolution'] = e_s[:,1]
-            # Note curl e = -iwb so b = -curl/iw
-
-            if self.verbose:
-                print 'Ran for {:f} seconds'.format(time.time()-startTime)
-                sys.stdout.flush()
-            Ainv.clean()
-        return F
-
diff --git a/SimPEG/NSEM/Problem3D/__init__.py b/SimPEG/NSEM/Problem3D/__init__.py
deleted file mode 100644
index 27e761a9..00000000
--- a/SimPEG/NSEM/Problem3D/__init__.py
+++ /dev/null
@@ -1 +0,0 @@
-from Probs import eForm_ps
\ No newline at end of file

From 70876327016dbf63a8d582c55b433ad4edc0d15d Mon Sep 17 00:00:00 2001
From: GudniRos <gkrosen@gmail.com>
Date: Tue, 31 May 2016 22:06:47 -0700
Subject: [PATCH 06/14] Fixing name space, updating utils

---
 SimPEG/NSEM/{NSEM.py => ProblemNSEM.py} | 10 +++++----
 SimPEG/NSEM/Utils/dataUtils.py          | 29 +++++++++++++------------
 SimPEG/NSEM/__init__.py                 |  4 ++--
 3 files changed, 23 insertions(+), 20 deletions(-)
 rename SimPEG/NSEM/{NSEM.py => ProblemNSEM.py} (99%)

diff --git a/SimPEG/NSEM/NSEM.py b/SimPEG/NSEM/ProblemNSEM.py
similarity index 99%
rename from SimPEG/NSEM/NSEM.py
rename to SimPEG/NSEM/ProblemNSEM.py
index d0bd8b77..18a72981 100644
--- a/SimPEG/NSEM/NSEM.py
+++ b/SimPEG/NSEM/ProblemNSEM.py
@@ -1,9 +1,10 @@
-from SimPEG.EM.Utils import omega, mu_0
+from SimPEG.EM.Utils.EMUtils import omega, mu_0
 from SimPEG import SolverLU as SimpegSolver, PropMaps, Utils, mkvc, sp, np
 from SimPEG.EM.FDEM.FDEM import BaseFDEMProblem
 from SurveyNSEM import Survey, Data
-from FieldsNSEM import Fields1D_ePrimSec, Fields3D_ePrimSec
+from FieldsNSEM import BaseNSEMFields, Fields1D_ePrimSec, Fields3D_ePrimSec
 from SimPEG.NSEM.Utils.MT1Danalytic import getEHfields
+import time, sys
 
 class BaseNSEMProblem(BaseFDEMProblem):
     """
@@ -169,7 +170,6 @@ class Problem1D_ePrimSec(BaseNSEMProblem):
 
     def __init__(self, mesh, **kwargs):
         BaseNSEMProblem.__init__(self, mesh, **kwargs)
-        self.fieldsPair = Fields1D_e
         # self._sigmaPrimary = sigmaPrimary
     @property
     def MeMui(self):
@@ -313,7 +313,7 @@ class Problem1D_eTotal(BaseNSEMProblem):
     # From FDEMproblem: Used to project the fields. Currently not used for NSEMproblem.
     _solutionType = 'e_1dSolution'
     _formulation  = 'EF'
-    fieldsPair = Fields1D_eTotal
+    # fieldsPair = Fields1D_eTotal
 
     def __init__(self, mesh, **kwargs):
         BaseNSEMProblem.__init__(self, mesh, **kwargs)
@@ -399,6 +399,8 @@ class Problem1D_eTotal(BaseNSEMProblem):
         :param np.ndarray (nC,) m: Conductivity model
         :param np.ndarray (nC,) m_back: Background conductivity model
         '''
+
+
         self.curModel = m
         # RHS, CalcFields = self.getRHS(freq,m_back), self.calcFields
 
diff --git a/SimPEG/NSEM/Utils/dataUtils.py b/SimPEG/NSEM/Utils/dataUtils.py
index fa97aa9a..b879c044 100644
--- a/SimPEG/NSEM/Utils/dataUtils.py
+++ b/SimPEG/NSEM/Utils/dataUtils.py
@@ -79,7 +79,7 @@ def plotMT1DModelData(problem,models,symList=None):
     axM = fig.add_axes([0.075,.1,.25,.875])
     axM.set_xlabel('Resistivity [Ohm*m]',fontsize=fontSize)
     axM.set_xlim(1e-1,1e5)
-    axM.set_ylim(-10000,5000)
+    # axM.set_ylim(-10000,5000)
     axM.set_ylabel('Depth [km]',fontsize=fontSize)
     axR = fig.add_axes([0.42,.575,.5,.4])
     axR.set_xscale('log')
@@ -132,24 +132,25 @@ def plotMT1DModelData(problem,models,symList=None):
     freq = simpeg.mkvc(data1D['freq'],2)
     res, phs = appResPhs(freq,allData)
 
-    stdCol = 'gray'
-    axRtw = axR.twinx()
-    axRtw.set_ylabel('Std of log10',color=stdCol)
-    [(t.set_color(stdCol), t.set_rotation(-45)) for t in axRtw.get_yticklabels()]
-    axPtw = axP.twinx()
-    axPtw.set_ylabel('Std ',color=stdCol)
-    [t.set_color(stdCol) for t in axPtw.get_yticklabels()]
-    axRtw.plot(freq, np.std(np.log10(res),1),'--',color=stdCol)
-    axPtw.plot(freq, np.std(phs,1),'--',color=stdCol)
+    if False:
+        stdCol = 'gray'
+        axRtw = axR.twinx()
+        axRtw.set_ylabel('Std of log10',color=stdCol)
+        [(t.set_color(stdCol), t.set_rotation(-45)) for t in axRtw.get_yticklabels()]
+        axPtw = axP.twinx()
+        axPtw.set_ylabel('Std ',color=stdCol)
+        [t.set_color(stdCol) for t in axPtw.get_yticklabels()]
+        axRtw.plot(freq, np.std(np.log10(res),1),'--',color=stdCol)
+        axPtw.plot(freq, np.std(phs,1),'--',color=stdCol)
 
     # Fix labels and ticks
 
-    yMtick = [l/1000 for l in axM.get_yticks().tolist()]
-    axM.set_yticklabels(yMtick)
+    # yMtick = [l/1000 for l in axM.get_yticks().tolist()]
+    # axM.set_yticklabels(yMtick)
     [ l.set_rotation(90) for l in axM.get_yticklabels()]
     [ l.set_rotation(90) for l in axR.get_yticklabels()]
-    [(t.set_color(stdCol), t.set_rotation(-45)) for t in axRtw.get_yticklabels()]
-    [t.set_color(stdCol) for t in axPtw.get_yticklabels()]
+    # [(t.set_color(stdCol), t.set_rotation(-45)) for t in axRtw.get_yticklabels()]
+    # [t.set_color(stdCol) for t in axPtw.get_yticklabels()]
     for ax in [axM,axR,axP]:
         ax.xaxis.set_tick_params(labelsize=fontSize)
         ax.yaxis.set_tick_params(labelsize=fontSize)
diff --git a/SimPEG/NSEM/__init__.py b/SimPEG/NSEM/__init__.py
index 98323827..bf9a8a65 100644
--- a/SimPEG/NSEM/__init__.py
+++ b/SimPEG/NSEM/__init__.py
@@ -1,5 +1,5 @@
 import Utils
 from SurveyNSEM import Rx, Survey, Data
-from FieldsNSEM import Fields1D_e, Fields3D_e
-import Problem1D, Problem2D, Problem3D
+from FieldsNSEM import Fields1D_ePrimSec, Fields3D_ePrimSec
+from ProblemNSEM import Problem1D_ePrimSec, Problem3D_ePrimSec
 import SrcNSEM
\ No newline at end of file

From f34314bba22dddf8afdd9e1d820e43b0605ab58d Mon Sep 17 00:00:00 2001
From: GudniRos <gkrosen@gmail.com>
Date: Tue, 31 May 2016 22:33:57 -0700
Subject: [PATCH 07/14] Fixed tests and Jvec

---
 SimPEG/NSEM/ProblemNSEM.py                          | 2 +-
 tests/mt/test_Problem1D_againstAnalyticHalfspace.py | 4 ++--
 tests/mt/test_Problem3D_againstAnalytic.py          | 6 ++++--
 3 files changed, 7 insertions(+), 5 deletions(-)

diff --git a/SimPEG/NSEM/ProblemNSEM.py b/SimPEG/NSEM/ProblemNSEM.py
index 18a72981..e4f04299 100644
--- a/SimPEG/NSEM/ProblemNSEM.py
+++ b/SimPEG/NSEM/ProblemNSEM.py
@@ -104,7 +104,7 @@ class BaseNSEMProblem(BaseFDEMProblem):
             ATinv = self.Solver(AT, **self.solverOpts)
 
             for src in self.survey.getSrcByFreq(freq):
-                ftype = self._fieldType + 'Solution'
+                ftype = self._solutionType
                 f_src = f[src, :] # Need to fix this...
 
                 for rx in src.rxList:
diff --git a/tests/mt/test_Problem1D_againstAnalyticHalfspace.py b/tests/mt/test_Problem1D_againstAnalyticHalfspace.py
index 1951271b..d848c1ba 100644
--- a/tests/mt/test_Problem1D_againstAnalyticHalfspace.py
+++ b/tests/mt/test_Problem1D_againstAnalyticHalfspace.py
@@ -99,7 +99,7 @@ def appRes_psFieldNorm(sigmaHalf):
 
     # Make the survey
     survey, sigma, mesh = setupSurvey(sigmaHalf,False)
-    problem = NSEM.Problem1D.eForm_psField(mesh, sigmaPrimary = sigma)
+    problem = NSEM.Problem1D_ePrimSec(mesh, sigmaPrimary = sigma)
     problem.pair(survey)
 
     # Get the fields
@@ -117,7 +117,7 @@ def appPhs_psFieldNorm(sigmaHalf):
 
     # Make the survey
     survey, sigma, mesh = setupSurvey(sigmaHalf,False)
-    problem = NSEM.Problem1D.eForm_psField(mesh, sigmaPrimary = sigma)
+    problem = NSEM.Problem1D_ePrimSec(mesh, sigmaPrimary = sigma)
     problem.pair(survey)
 
     # Get the fields
diff --git a/tests/mt/test_Problem3D_againstAnalytic.py b/tests/mt/test_Problem3D_againstAnalytic.py
index 9ad11a33..a267c9d7 100644
--- a/tests/mt/test_Problem3D_againstAnalytic.py
+++ b/tests/mt/test_Problem3D_againstAnalytic.py
@@ -7,6 +7,8 @@ from SimPEG import NSEM
 from SimPEG.Utils import meshTensor
 from scipy.constants import mu_0
 
+np.random.seed(1100)
+
 TOLr = 5e-2
 TOL = 1e-4
 FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order
@@ -127,11 +129,11 @@ def setupSimpegNSEMfwd_eForm_ps(inputSetup,comp='Imp',singleFreq=False,expMap=Tr
     ## Setup the problem object
     sigma1d = M.r(sigBG,'CC','CC','M')[0,0,:]
     if expMap:
-        problem = NSEM.Problem3D.eForm_ps(M,sigmaPrimary= np.log(sigma1d) )
+        problem = NSEM.Problem3D_ePrimSec(M,sigmaPrimary= np.log(sigma1d) )
         problem.mapping = simpeg.Maps.ExpMap(problem.mesh)
         problem.curModel = np.log(sig)
     else:
-        problem = NSEM.Problem3D.eForm_ps(M,sigmaPrimary= sigma1d)
+        problem = NSEM.Problem3D_ePrimSec(M,sigmaPrimary= sigma1d)
         problem.curModel = sig
     problem.pair(survey)
     problem.verbose = False

From fb7f5a53d4d40be2904414285fb9eeed05908bd1 Mon Sep 17 00:00:00 2001
From: GudniRos <gkrosen@gmail.com>
Date: Wed, 1 Jun 2016 01:05:21 -0700
Subject: [PATCH 08/14] Working on examples

---
 SimPEG/Examples/MT_1D_ForwardAndInversion.py | 33 +++++++++++---------
 SimPEG/Examples/MT_3D_Foward.py              | 12 +++----
 2 files changed, 24 insertions(+), 21 deletions(-)

diff --git a/SimPEG/Examples/MT_1D_ForwardAndInversion.py b/SimPEG/Examples/MT_1D_ForwardAndInversion.py
index 69cedd98..59a6c87e 100644
--- a/SimPEG/Examples/MT_1D_ForwardAndInversion.py
+++ b/SimPEG/Examples/MT_1D_ForwardAndInversion.py
@@ -1,9 +1,11 @@
 import SimPEG as simpeg
 import numpy as np
-import SimPEG.MT as MT
+from SimPEG import NSEM
 from scipy.constants import mu_0
 import matplotlib.pyplot as plt
 
+np.random.seed(1983)
+
 def run(plotIt=True):
     """
         MT: 1D: Inversion
@@ -23,7 +25,7 @@ def run(plotIt=True):
     ct = 20
     air = simpeg.Utils.meshTensor([(ct,16,1.4)])
     core = np.concatenate( (  np.kron(simpeg.Utils.meshTensor([(ct,10,-1.3)]),np.ones((5,))) , simpeg.Utils.meshTensor([(ct,5)]) ) )
-    bot = simpeg.Utils.meshTensor([(core[0],10,-1.4)])
+    bot = simpeg.Utils.meshTensor([(core[0],12,-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)
@@ -47,41 +49,43 @@ def run(plotIt=True):
     # Make the background model
     sigma_0 = np.ones(m1d.nCx)*sig_air
     sigma_0[active] = sig_half
+    sigma_0[layer1] = sig_layer1
+    sigma_0[layer2] = .002
     m_0 = np.log(sigma_0[active])
 
     # Set the mapping
-    actMap = simpeg.Maps.ActiveCells(m1d, active, np.log(1e-8), nC=m1d.nCx)
+    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 = []
     for rxType in ['z1dr','z1di']:
-        rxList.append(MT.Rx(simpeg.mkvc(np.array([0.0]),2).T,rxType))
+        rxList.append(NSEM.Rx(simpeg.mkvc(np.array([0.0]),2).T,rxType))
     # Source list
     srcList =[]
     for freq in freqs:
-            srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq))
+            srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,freq))
     # Make the survey
-    survey = MT.Survey(srcList)
+    survey = NSEM.Survey(srcList)
     survey.mtrue = m_true
 
     ## Set the problem
-    problem = MT.Problem1D.eForm_psField(m1d,sigmaPrimary=sigma_0,mapping=mappingExpAct)
+    problem = NSEM.Problem1D_ePrimSec(m1d,sigmaPrimary=sigma_0,mapping=mappingExpAct)
     problem.pair(survey)
 
     ## Forward model data
     # Project the data
     survey.dtrue = survey.dpred(m_true)
-    survey.dobs = survey.dtrue + 0.025*abs(survey.dtrue)*np.random.randn(*survey.dtrue.shape)
+    survey.dobs = survey.dtrue + 0.01*abs(survey.dtrue)*np.random.randn(*survey.dtrue.shape)
 
     if plotIt:
-        fig = MT.Utils.dataUtils.plotMT1DModelData(problem)
+        fig = NSEM.Utils.dataUtils.plotMT1DModelData(problem,[])
         fig.suptitle('Target - smooth true')
 
 
     # Assign uncertainties
-    std = 0.05 # 5% std
+    std = 0.025 # 5% std
     survey.std = np.abs(survey.dobs*std)
     # Assign the data weight
     Wd = 1./survey.std
@@ -92,6 +96,7 @@ def run(plotIt=True):
     # Set the optimization
     opt = simpeg.Optimization.InexactGaussNewton(maxIter = 30)
     opt.counter = C
+    opt.maxStep = m1d.nC * survey.nD
     opt.LSshorten = 0.5
     opt.remember('xc')
     # Data misfit
@@ -99,7 +104,7 @@ def run(plotIt=True):
     dmis.Wd = Wd
     # Regularization - with a regularization mesh
     regMesh = simpeg.Mesh.TensorMesh([m1d.hx[problem.mapping.sigmaMap.maps[-1].indActive]],m1d.x0)
-    reg = simpeg.Regularization.Tikhonov(regMesh)
+    reg = simpeg.Regularization.Tikhonov(m1d,indActive=active)
     reg.mrefInSmooth = True
     reg.alpha_s = 1e-7
     reg.alpha_x = 1.
@@ -109,11 +114,9 @@ def run(plotIt=True):
     # Beta cooling
     beta = simpeg.Directives.BetaSchedule()
     beta.coolingRate = 4
-    betaest = simpeg.Directives.BetaEstimate_ByEig(beta0_ratio=0.75)
+    betaest = simpeg.Directives.BetaEstimate_ByEig(beta0_ratio=1.)
     targmis = simpeg.Directives.TargetMisfit()
     targmis.target = survey.nD
-    saveModel = simpeg.Directives.SaveModelEveryIteration()
-    saveModel.fileName = 'Inversion_TargMisEqnD_smoothTrue'
     # Create an inversion object
     inv = simpeg.Inversion.BaseInversion(invProb, directiveList=[beta,betaest,targmis])
 
@@ -121,7 +124,7 @@ def run(plotIt=True):
     mopt = inv.run(m_0)
 
     if plotIt:
-        fig = MT.Utils.dataUtils.plotMT1DModelData(problem,[mopt])
+        fig = NSEM.Utils.dataUtils.plotMT1DModelData(problem,[mopt])
         fig.suptitle('Target - smooth true')
         plt.show()
 
diff --git a/SimPEG/Examples/MT_3D_Foward.py b/SimPEG/Examples/MT_3D_Foward.py
index da16eeee..01700afe 100644
--- a/SimPEG/Examples/MT_3D_Foward.py
+++ b/SimPEG/Examples/MT_3D_Foward.py
@@ -2,7 +2,7 @@
 
 # Import
 import SimPEG as simpeg
-from SimPEG import MT
+from SimPEG import NSEM
 import numpy as np
 try:
     from pymatsolver import MumpsSolver as Solver
@@ -37,16 +37,16 @@ def run(plotIt=True, nFreq=1):
     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))
+            rxList.append(NSEM.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))
+        srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,freq))
     # Survey MT
-    survey = MT.Survey(srcList)
+    survey = NSEM.Survey(srcList)
 
     ## Setup the problem object
-    problem = MT.Problem3D.eForm_ps(M, sigmaPrimary=sigBG)
+    problem = NSEM.Problem3D_ePrimSec(M, sigmaPrimary=sigBG)
     problem.pair(survey)
     problem.Solver = Solver
 
@@ -55,7 +55,7 @@ def run(plotIt=True, nFreq=1):
     dataVec = survey.eval(fields)
 
     # Make the data
-    mtData = MT.Data(survey,dataVec)
+    mtData = NSEM.Data(survey,dataVec)
     # Add plots
     if plotIt:
         pass

From 8117ee3b062bad3f7eb697bb1b4c64fce373cd03 Mon Sep 17 00:00:00 2001
From: GudniRos <gkrosen@gmail.com>
Date: Wed, 1 Jun 2016 10:22:18 -0700
Subject: [PATCH 09/14] Fixing MT_1D example, runs but inversion results should
 be better.

---
 SimPEG/Examples/MT_1D_ForwardAndInversion.py | 38 +++++++++++---------
 1 file changed, 22 insertions(+), 16 deletions(-)

diff --git a/SimPEG/Examples/MT_1D_ForwardAndInversion.py b/SimPEG/Examples/MT_1D_ForwardAndInversion.py
index 59a6c87e..a1e33d62 100644
--- a/SimPEG/Examples/MT_1D_ForwardAndInversion.py
+++ b/SimPEG/Examples/MT_1D_ForwardAndInversion.py
@@ -19,13 +19,13 @@ def run(plotIt=True):
     ## Setup the forward modeling
     # Setting up 1D mesh and conductivity models to forward model data.
     # Frequency
-    nFreq = 31
-    freqs = np.logspace(3,-3,nFreq)
+    nFreq = 26
+    freqs = np.logspace(2,-3,nFreq)
     # Set mesh parameters
-    ct = 20
-    air = simpeg.Utils.meshTensor([(ct,16,1.4)])
+    ct = 10
+    air = simpeg.Utils.meshTensor([(ct,25,1.4)])
     core = np.concatenate( (  np.kron(simpeg.Utils.meshTensor([(ct,10,-1.3)]),np.ones((5,))) , simpeg.Utils.meshTensor([(ct,5)]) ) )
-    bot = simpeg.Utils.meshTensor([(core[0],12,-1.4)])
+    bot = simpeg.Utils.meshTensor([(core[0],25,-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)
@@ -35,7 +35,7 @@ def run(plotIt=True):
     layer1 = (m1d.vectorCCx<-500.) & (m1d.vectorCCx>=-800.)
     layer2 = (m1d.vectorCCx<-3500.) & (m1d.vectorCCx>=-5000.)
     # Set the conductivity values
-    sig_half = 2e-3
+    sig_half = 1e-2
     sig_air = 1e-8
     sig_layer1 = .2
     sig_layer2 = .2
@@ -49,8 +49,8 @@ def run(plotIt=True):
     # Make the background model
     sigma_0 = np.ones(m1d.nCx)*sig_air
     sigma_0[active] = sig_half
-    sigma_0[layer1] = sig_layer1
-    sigma_0[layer2] = .002
+    # sigma_0[layer1] = sig_layer1
+    # sigma_0[layer2] = .002
     m_0 = np.log(sigma_0[active])
 
     # Set the mapping
@@ -61,7 +61,7 @@ def run(plotIt=True):
     # Receivers
     rxList = []
     for rxType in ['z1dr','z1di']:
-        rxList.append(NSEM.Rx(simpeg.mkvc(np.array([0.0]),2).T,rxType))
+        rxList.append(NSEM.Rx(simpeg.mkvc(np.array([-0.5]),2).T,rxType))
     # Source list
     srcList =[]
     for freq in freqs:
@@ -94,27 +94,32 @@ def run(plotIt=True):
     # Define a counter
     C =  simpeg.Utils.Counter()
     # Set the optimization
-    opt = simpeg.Optimization.InexactGaussNewton(maxIter = 30)
+    opt = simpeg.Optimization.ProjectedGNCG(maxIter = 25)
     opt.counter = C
-    opt.maxStep = m1d.nC * survey.nD
-    opt.LSshorten = 0.5
+    opt.lower = np.log(1e-4)
+    opt.upper = np.log(5)
+    opt.LSshorten = 0.1
     opt.remember('xc')
     # Data misfit
     dmis = simpeg.DataMisfit.l2_DataMisfit(survey)
     dmis.Wd = Wd
     # Regularization - with a regularization mesh
-    regMesh = simpeg.Mesh.TensorMesh([m1d.hx[problem.mapping.sigmaMap.maps[-1].indActive]],m1d.x0)
-    reg = simpeg.Regularization.Tikhonov(m1d,indActive=active)
+    regMesh = simpeg.Mesh.TensorMesh([m1d.hx[active]],m1d.x0)
+    reg = simpeg.Regularization.Tikhonov(regMesh)
+    reg.mref = m_true
     reg.mrefInSmooth = True
-    reg.alpha_s = 1e-7
+    reg.alpha_s = 1e-1
     reg.alpha_x = 1.
+
     # Inversion problem
     invProb = simpeg.InvProblem.BaseInvProblem(dmis, reg, opt)
     invProb.counter = C
     # Beta cooling
     beta = simpeg.Directives.BetaSchedule()
-    beta.coolingRate = 4
+    beta.coolingRate = 4.
+    beta.coolingFactor = 4.
     betaest = simpeg.Directives.BetaEstimate_ByEig(beta0_ratio=1.)
+    # betaest.beta0 = 1.
     targmis = simpeg.Directives.TargetMisfit()
     targmis.target = survey.nD
     # Create an inversion object
@@ -126,6 +131,7 @@ def run(plotIt=True):
     if plotIt:
         fig = NSEM.Utils.dataUtils.plotMT1DModelData(problem,[mopt])
         fig.suptitle('Target - smooth true')
+        fig.axes[0].set_ylim([-10000,500])
         plt.show()
 
 if __name__ == '__main__':

From 8d4581e92f87e74d841987f3e97c915febbe9536 Mon Sep 17 00:00:00 2001
From: GudniRos <gkrosen@gmail.com>
Date: Wed, 1 Jun 2016 11:55:12 -0700
Subject: [PATCH 10/14] Cleaned MT_1D inverison example

---
 SimPEG/Examples/MT_1D_ForwardAndInversion.py | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/SimPEG/Examples/MT_1D_ForwardAndInversion.py b/SimPEG/Examples/MT_1D_ForwardAndInversion.py
index a1e33d62..59e30e01 100644
--- a/SimPEG/Examples/MT_1D_ForwardAndInversion.py
+++ b/SimPEG/Examples/MT_1D_ForwardAndInversion.py
@@ -49,8 +49,6 @@ def run(plotIt=True):
     # Make the background model
     sigma_0 = np.ones(m1d.nCx)*sig_air
     sigma_0[active] = sig_half
-    # sigma_0[layer1] = sig_layer1
-    # sigma_0[layer2] = .002
     m_0 = np.log(sigma_0[active])
 
     # Set the mapping
@@ -119,7 +117,7 @@ def run(plotIt=True):
     beta.coolingRate = 4.
     beta.coolingFactor = 4.
     betaest = simpeg.Directives.BetaEstimate_ByEig(beta0_ratio=1.)
-    # betaest.beta0 = 1.
+    betaest.beta0 = 1.
     targmis = simpeg.Directives.TargetMisfit()
     targmis.target = survey.nD
     # Create an inversion object

From b3029b697b2daa3cb8a42b03139016f69202c0d5 Mon Sep 17 00:00:00 2001
From: GudniRos <gkrosen@gmail.com>
Date: Wed, 1 Jun 2016 12:00:52 -0700
Subject: [PATCH 11/14] Adjust settings for the MT1D inversion example

---
 SimPEG/Examples/MT_1D_ForwardAndInversion.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/SimPEG/Examples/MT_1D_ForwardAndInversion.py b/SimPEG/Examples/MT_1D_ForwardAndInversion.py
index 59e30e01..6af45d0b 100644
--- a/SimPEG/Examples/MT_1D_ForwardAndInversion.py
+++ b/SimPEG/Examples/MT_1D_ForwardAndInversion.py
@@ -104,7 +104,6 @@ def run(plotIt=True):
     # Regularization - with a regularization mesh
     regMesh = simpeg.Mesh.TensorMesh([m1d.hx[active]],m1d.x0)
     reg = simpeg.Regularization.Tikhonov(regMesh)
-    reg.mref = m_true
     reg.mrefInSmooth = True
     reg.alpha_s = 1e-1
     reg.alpha_x = 1.

From 50c4a6caf839f9a508925698d08dd9a0a0dc7e73 Mon Sep 17 00:00:00 2001
From: GudniRos <gkrosen@gmail.com>
Date: Wed, 1 Jun 2016 20:53:12 -0700
Subject: [PATCH 12/14] Refactor and updated tests

---
 SimPEG/Examples/__init__.py                   |  29 +-
 SimPEG/NSEM/Utils/__init__.py                 |   5 +-
 SimPEG/NSEM/Utils/srcUtils.py                 |  46 ---
 SimPEG/NSEM/Utils/testUtils.py                | 198 +++++++++++++
 tests/mt/{ => forward}/__init__.py            |   0
 .../test_AnalyticFunctionVsAppResPhs.py}      |   7 +-
 .../test_Problem1D_AnalyticVsNumeric.py}      |  59 +---
 .../test_Problem1D_VsAnalyticHalfspace.py     | 102 +++++++
 .../test_Problem3D_VsAnalyticSolution.py      |  54 ++++
 tests/mt/inversion/__init__.py                |  12 +
 tests/mt/inversion/test_Problem3D_Adjoint.py  |  58 ++++
 tests/mt/inversion/test_Problem3D_Derivs.py   |  83 ++++++
 ...test_Problem1D_againstAnalyticHalfspace.py | 147 ----------
 tests/mt/test_Problem3D_againstAnalytic.py    | 270 ------------------
 14 files changed, 537 insertions(+), 533 deletions(-)
 delete mode 100644 SimPEG/NSEM/Utils/srcUtils.py
 create mode 100644 SimPEG/NSEM/Utils/testUtils.py
 rename tests/mt/{ => forward}/__init__.py (100%)
 rename tests/mt/{test_ApparentResistivityAnalytic.py => forward/test_AnalyticFunctionVsAppResPhs.py} (85%)
 rename tests/mt/{test_Problem1D_totalDvsPSvsAnalytic.py => forward/test_Problem1D_AnalyticVsNumeric.py} (57%)
 create mode 100644 tests/mt/forward/test_Problem1D_VsAnalyticHalfspace.py
 create mode 100644 tests/mt/forward/test_Problem3D_VsAnalyticSolution.py
 create mode 100644 tests/mt/inversion/__init__.py
 create mode 100644 tests/mt/inversion/test_Problem3D_Adjoint.py
 create mode 100644 tests/mt/inversion/test_Problem3D_Derivs.py
 delete mode 100644 tests/mt/test_Problem1D_againstAnalyticHalfspace.py
 delete mode 100644 tests/mt/test_Problem3D_againstAnalytic.py

diff --git a/SimPEG/Examples/__init__.py b/SimPEG/Examples/__init__.py
index c67652a2..7b2e1f70 100644
--- a/SimPEG/Examples/__init__.py
+++ b/SimPEG/Examples/__init__.py
@@ -1,27 +1,28 @@
 # Run this file to add imports.
 
 ##### AUTOIMPORTS #####
-import DC_Analytic_Dipole
-import DC_Forward_PseudoSection
 import EM_FDEM_1D_Inversion
-import EM_FDEM_Analytic_MagDipoleWholespace
-import EM_Schenkel_Morrison_Casing
+import Mesh_QuadTree_Creation
 import EM_TDEM_1D_Inversion
+import Mesh_QuadTree_FaceDiv
+import Mesh_Tensor_Creation
 import FLOW_Richards_1D_Celia1990
-import Forward_BasicDirectCurrent
+import DC_Forward_PseudoSection
+import Mesh_Operators_CahnHilliard
+import Mesh_Basic_Types
 import Inversion_IRLS
 import Inversion_Linear
-import Mesh_Basic_PlotImage
-import Mesh_Basic_Types
-import Mesh_Operators_CahnHilliard
-import Mesh_QuadTree_Creation
-import Mesh_QuadTree_FaceDiv
-import Mesh_QuadTree_HangingNodes
-import Mesh_Tensor_Creation
-import MT_1D_ForwardAndInversion
+import EM_Schenkel_Morrison_Casing
 import MT_3D_Foward
+import MT_1D_ForwardAndInversion
+import MT_1D_analytic_nlayer_Earth
+import Forward_BasicDirectCurrent
+import EM_FDEM_Analytic_MagDipoleWholespace
+import Mesh_Basic_PlotImage
+import DC_Analytic_Dipole
+import Mesh_QuadTree_HangingNodes
 
-__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_IRLS", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_ForwardAndInversion", "MT_3D_Foward"]
+__examples__ = ["EM_FDEM_1D_Inversion", "Mesh_QuadTree_Creation", "EM_TDEM_1D_Inversion", "Mesh_QuadTree_FaceDiv", "Mesh_Tensor_Creation", "FLOW_Richards_1D_Celia1990", "DC_Forward_PseudoSection", "Mesh_Operators_CahnHilliard", "Mesh_Basic_Types", "Inversion_IRLS", "Inversion_Linear", "EM_Schenkel_Morrison_Casing", "MT_3D_Foward", "MT_1D_ForwardAndInversion", "MT_1D_analytic_nlayer_Earth", "Forward_BasicDirectCurrent", "EM_FDEM_Analytic_MagDipoleWholespace", "Mesh_Basic_PlotImage", "DC_Analytic_Dipole", "Mesh_QuadTree_HangingNodes"]
 
 ##### AUTOIMPORTS #####
 
diff --git a/SimPEG/NSEM/Utils/__init__.py b/SimPEG/NSEM/Utils/__init__.py
index b683f8b4..acfc1a65 100644
--- a/SimPEG/NSEM/Utils/__init__.py
+++ b/SimPEG/NSEM/Utils/__init__.py
@@ -1,4 +1,5 @@
-from MT1Dsolutions import *  # Add the names of the functions
-from MT1Danalytic import *
+from MT1Dsolutions import get1DEfields  # Add the names of the functions
+from MT1Danalytic import getEHfields, getImpedance
 from dataUtils import *
 from ediFilesUtils import *
+from testUtils import *
diff --git a/SimPEG/NSEM/Utils/srcUtils.py b/SimPEG/NSEM/Utils/srcUtils.py
deleted file mode 100644
index 9bf141d8..00000000
--- a/SimPEG/NSEM/Utils/srcUtils.py
+++ /dev/null
@@ -1,46 +0,0 @@
-import SimPEG as simpeg, numpy as np
-
-def homo1DModelSource(mesh,freq,m_back):
-    '''
-        Function that calculates and return background fields for a 3D mesh and model.
-        The calculuations use 1D field solution for a vertical slice throught model (south-western most column),
-        which is assigned at the fields everywhere for the respective polarizations.2
-
-        :param Simpeg mesh object mesh: Holds information on the discretization
-        :param float freq: The frequency to solve at
-        :param np.array m_back: Background model of conductivity to base the calculations on.
-        :rtype: numpy.ndarray (mesh.nE,2)
-        :return: eBG_bp, E fields for the background model at both polarizations.
-
-    '''
-
-    # import
-    from SimPEG.NSEM.Utils import get1DEfields
-    # Get a 1d solution for a halfspace background
-    mesh1d = simpeg.Mesh.TensorMesh([mesh.hz],np.array([mesh.x0[2]]))
-    # Note: Everything is using e^iwt
-    e0_1d = get1DEfields(mesh1d,mesh.r(m_back,'CC','CC','M')[0,0,:],freq)
-    # 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
diff --git a/SimPEG/NSEM/Utils/testUtils.py b/SimPEG/NSEM/Utils/testUtils.py
new file mode 100644
index 00000000..8ac2b940
--- /dev/null
+++ b/SimPEG/NSEM/Utils/testUtils.py
@@ -0,0 +1,198 @@
+import unittest
+import sys
+from scipy.constants import mu_0
+import SimPEG as simpeg
+
+from SimPEG.Utils import meshTensor
+import numpy as np
+
+np.random.seed(1100)
+# Define the tolerances
+TOLr = 5e-2
+TOLp = 5e-2
+
+
+def getAppResPhs(NSEMdata):
+    # Make impedance
+    from SimPEG.NSEM.Utils import appResPhs
+    zList = []
+    for src in NSEMdata.survey.srcList:
+        zc = [src.freq]
+        for rx in src.rxList:
+            if 'i' in rx.rxType:
+                m=1j
+            else:
+                m = 1
+            zc.append(m*NSEMdata[src,rx])
+        zList.append(zc)
+    return [appResPhs(zList[i][0],np.sum(zList[i][1:3])) for i in np.arange(len(zList))]
+
+
+def setup1DSurvey(sigmaHalf,tD=True,structure=False):
+    from SimPEG import NSEM
+    # 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 ['z1dr','z1di']:
+        rxList.append(NSEM.Rx(simpeg.mkvc(np.array([0.0]),2).T,rxType))
+    # Source list
+    srcList =[]
+    if tD:
+        for freq in freqs:
+            srcList.append(NSEM.SrcNSEM.polxy_1DhomotD(rxList,freq))
+    else:
+        for freq in freqs:
+            srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,freq))
+
+    survey = NSEM.Survey(srcList)
+    return survey, sigma, m1d
+
+
+def setupSimpegNSEM_ePrimSec(inputSetup,comp='Imp',singleFreq=False,expMap=True):
+    from SimPEG import NSEM
+
+    M,freqs,sig,sigBG,rx_loc = inputSetup
+    # Make a receiver list
+    rxList = []
+    if comp == 'All':
+        for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi','tzxr','tzxi','tzyr','tzyi']:
+            rxList.append(NSEM.Rx(rx_loc,rxType))
+    elif comp == 'Imp':
+        for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi']:
+            rxList.append(NSEM.Rx(rx_loc,rxType))
+    elif comp == 'Tip':
+        for rxType in ['tzxr','tzxi','tzyr','tzyi']:
+            rxList.append(NSEM.Rx(rx_loc,rxType))
+    else:
+        rxList.append(NSEM.Rx(rx_loc,comp))
+    # Source list
+    srcList =[]
+
+    if singleFreq:
+        srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,singleFreq))
+    else:
+        for freq in freqs:
+            srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,freq))
+    # Survey NSEM
+    survey = NSEM.Survey(srcList)
+
+    ## Setup the problem object
+    sigma1d = M.r(sigBG,'CC','CC','M')[0,0,:]
+    if expMap:
+        problem = NSEM.Problem3D_ePrimSec(M,sigmaPrimary= np.log(sigma1d) )
+        problem.mapping = simpeg.Maps.ExpMap(problem.mesh)
+        problem.curModel = np.log(sig)
+    else:
+        problem = NSEM.Problem3D_ePrimSec(M,sigmaPrimary= sigma1d)
+        problem.curModel = sig
+    problem.pair(survey)
+    problem.verbose = False
+    try:
+        from pymatsolver import MumpsSolver
+        problem.Solver = MumpsSolver
+    except:
+        pass
+
+    return (survey, problem)
+
+def getInputs():
+    """
+    Function that returns Mesh, freqs, rx_loc, elev.
+    """
+    # 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])
+    # M = simpeg.Mesh.TensorMesh([[(1000,6,-1.5),(1000.,6),(1000,6,1.5)],[(1000,6,-1.5),(1000.,2),(1000,6,1.5)],[(1000,6,-1.3),(1000.,6),(1000,6,1.3)]], x0=['C','C','C'])# Setup the model
+    M = simpeg.Mesh.TensorMesh([[(200,6,-1.5),(200.,4),(200,6,1.5)],[(200,6,-1.5),(200.,4),(200,6,1.5)],[(200,8,-1.5),(200.,8),(200,8,1.5)]], x0=['C','C','C'])# Setup the model
+    # Set the frequencies
+    freqs = np.logspace(1,-3,5)
+    elev = 0
+
+    ## Setup the the survey object
+    # Receiver locations
+    rx_x, rx_y = np.meshgrid(np.arange(-350,350,200),np.arange(-350,350,200))
+    rx_loc = np.hstack((simpeg.Utils.mkvc(rx_x,2),simpeg.Utils.mkvc(rx_y,2),elev+np.zeros((np.prod(rx_x.shape),1))))
+
+    return M, freqs, rx_loc, elev
+
+def random(conds):
+    ''' Returns a halfspace model based on the inputs'''
+    M, freqs, rx_loc, elev = getInputs()
+
+    # Backround
+    sigBG = np.ones(M.nC)*conds
+    # Add randomness to the model (10% of the value).
+    sig = np.exp( np.log(sigBG) + np.random.randn(M.nC)*(conds)*1e-1 )
+
+    return (M, freqs, sig, sigBG, rx_loc)
+
+def halfSpace(conds):
+    ''' Returns a halfspace model based on the inputs'''
+    M, freqs, rx_loc, elev = getInputs()
+
+    # Model
+    ccM = M.gridCC
+    # conds = [1e-2]
+    groundInd = ccM[:,2] < elev
+    sig = np.zeros(M.nC) + 1e-8
+    sig[groundInd] = conds
+    # Set the background, not the same as the model
+    sigBG = np.zeros(M.nC) + 1e-8
+    sigBG[groundInd] = conds
+
+    return (M, freqs, sig, sigBG, rx_loc)
+
+def blockInhalfSpace(conds):
+    ''' Returns a halfspace model based on the inputs'''
+    M, freqs, rx_loc, elev = getInputs()
+
+    # Model
+    ccM = M.gridCC
+    # conds = [1e-2]
+    groundInd = ccM[:,2] < elev
+    sig = simpeg.Utils.ModelBuilder.defineBlock(M.gridCC,np.array([-1000,-1000,-1500]),np.array([1000,1000,-1000]),conds)
+    sig[~groundInd] = 1e-8
+    # Set the background, not the same as the model
+    sigBG = np.zeros(M.nC) + 1e-8
+    sigBG[groundInd] = conds[1]
+
+    return (M, freqs, sig, sigBG, rx_loc)
+
+def twoLayer(conds):
+    ''' Returns a 2 layer model based on the conductivity values given'''
+    M, freqs, rx_loc, elev = getInputs()
+
+    # Model
+    ccM = M.gridCC
+    groundInd = ccM[:,2] < elev
+    botInd = ccM[:,2] < -3000
+    sig = np.zeros(M.nC) + 1e-8
+    sig[groundInd] = conds[1]
+    sig[botInd] = conds[0]
+    # Set the background, not the same as the model
+    sigBG = np.zeros(M.nC) + 1e-8
+    sigBG[groundInd] = conds[1]
+
+
+    return (M, freqs, sig, sigBG, rx_loc)
+
diff --git a/tests/mt/__init__.py b/tests/mt/forward/__init__.py
similarity index 100%
rename from tests/mt/__init__.py
rename to tests/mt/forward/__init__.py
diff --git a/tests/mt/test_ApparentResistivityAnalytic.py b/tests/mt/forward/test_AnalyticFunctionVsAppResPhs.py
similarity index 85%
rename from tests/mt/test_ApparentResistivityAnalytic.py
rename to tests/mt/forward/test_AnalyticFunctionVsAppResPhs.py
index 44bc3c23..9a838b4b 100644
--- a/tests/mt/test_ApparentResistivityAnalytic.py
+++ b/tests/mt/forward/test_AnalyticFunctionVsAppResPhs.py
@@ -4,10 +4,7 @@ from SimPEG import NSEM
 
 TOL = 1e-6
 
-def appResPhs(freq,z):
-    app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2
-    app_phs = np.arctan2(-z.imag,z.real)*(180/np.pi)
-    return app_res, app_phs
+
 
 def appResNorm(sigmaHalf):
     nFreq = 26
@@ -25,7 +22,7 @@ def appResNorm(sigmaHalf):
 
     Zarr = np.concatenate(Z)
 
-    app_r, app_p = appResPhs(freqs,Zarr)
+    app_r, app_p = NSEM.Utils.appResPhs(freqs,Zarr)
 
     return np.linalg.norm(np.abs(app_r - np.ones(nFreq)/sigmaHalf)) / np.log10(sigmaHalf)
 
diff --git a/tests/mt/test_Problem1D_totalDvsPSvsAnalytic.py b/tests/mt/forward/test_Problem1D_AnalyticVsNumeric.py
similarity index 57%
rename from tests/mt/test_Problem1D_totalDvsPSvsAnalytic.py
rename to tests/mt/forward/test_Problem1D_AnalyticVsNumeric.py
index 09c68400..90a4b5ca 100644
--- a/tests/mt/test_Problem1D_totalDvsPSvsAnalytic.py
+++ b/tests/mt/forward/test_Problem1D_AnalyticVsNumeric.py
@@ -8,45 +8,6 @@ TOLr = 5e-2
 TOLp = 5e-2
 
 
-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],15,-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
-    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 ['z1dr','z1di']:
-        rxList.append(NSEM.Rx(simpeg.mkvc(np.array([0.0]),2).T,rxType))
-    # Source list
-    srcList =[]
-    if tD:
-        for freq in freqs:
-            srcList.append(NSEM.SrcNSEM.polxy_1DhomotD(rxList,freq))
-    else:
-        for freq in freqs:
-            srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,freq))
-
-    survey = NSEM.Survey(srcList)
-    return survey, sigma, m1d
-
 def getAppResPhs(NSEMdata):
     # Make impedance
     def appResPhs(freq,z):
@@ -70,7 +31,7 @@ def calculateAnalyticSolution(srcList,mesh,model):
     data1D = NSEM.Data(surveyAna)
     for src in surveyAna.srcList:
         elev = src.rxList[0].locs[0]
-        anaEd, anaEu, anaHd, anaHu = NSEM.Utils.NSEM1Danalytic.getEHfields(mesh,model,src.freq,elev)
+        anaEd, anaEu, anaHd, anaHu = NSEM.Utils.MT1Danalytic.getEHfields(mesh,model,src.freq,elev)
         anaE = anaEd+anaEu
         anaH = anaHd+anaHu
         # Scale the solution
@@ -86,8 +47,8 @@ def dataMis_AnalyticTotalDomain(sigmaHalf):
     # Make the survey
 
     # Total domain solution
-    surveyTD, sigma, mesh = setupSurvey(sigmaHalf)
-    problemTD = NSEM.Problem1D.eForm_TotalField(mesh)
+    surveyTD, sigma, mesh = NSEM.Utils.testUtils.setup1DSurvey(sigmaHalf)
+    problemTD = NSEM.Problem1D_eTotal(mesh) # This not fully implemented
     problemTD.pair(surveyTD)
     # Analytic data
     dataAnaObj = calculateAnalyticSolution(surveyTD.srcList,mesh,sigma)
@@ -108,16 +69,16 @@ def dataMis_AnalyticPrimarySecondary(sigmaHalf):
 
     # Make the survey
     # Primary secondary
-    surveyPS, sigmaPS, mesh = setupSurvey(sigmaHalf,tD=False)
-    problemPS = NSEM.Problem1D.eForm_psField(mesh)
-    problemPS.sigmaPrimary = sigmaPS
-    problemPS.pair(surveyPS)
+    survey, sigma, mesh = NSEM.Utils.testUtils.setup1DSurvey(sigmaHalf,False,structure=True)
     # Analytic data
-    dataAnaObj = calculateAnalyticSolution(surveyPS.srcList,mesh,sigmaPS)
+    problem = NSEM.Problem1D_ePrimSec(mesh, sigmaPrimary = sigma)
+    problem.pair(survey)
 
-    dataPS = surveyPS.dpred(sigmaPS)
+    dataAnaObj = calculateAnalyticSolution(survey.srcList,mesh,sigma)
+
+    data = survey.dpred(sigma)
     dataAna = simpeg.mkvc(dataAnaObj)
-    return np.all((dataPS - dataAna)/dataAna < 2.)
+    return np.all((data - dataAna)/dataAna < 2.)
 
 
 
diff --git a/tests/mt/forward/test_Problem1D_VsAnalyticHalfspace.py b/tests/mt/forward/test_Problem1D_VsAnalyticHalfspace.py
new file mode 100644
index 00000000..098ca0ef
--- /dev/null
+++ b/tests/mt/forward/test_Problem1D_VsAnalyticHalfspace.py
@@ -0,0 +1,102 @@
+import unittest
+import SimPEG as simpeg
+from SimPEG import NSEM
+from SimPEG.Utils import meshTensor
+import numpy as np
+# Define the tolerances
+TOLr = 5e-1
+TOLp = 5e-1
+
+
+
+def appRes_TotalFieldNorm(sigmaHalf):
+
+    # Make the survey
+    survey, sigma, mesh = NSEM.Utils.testUtils.setup1DSurvey(sigmaHalf)
+    problem = NSEM.Problem1D_eTotal(mesh)
+    problem.pair(survey)
+
+    # Get the fields
+    fields = problem.fields(sigma)
+
+    # Project the data
+    data = survey.eval(fields)
+
+    # Calculate the app res and phs
+    app_r = np.array(NSEM.Utils.testUtils.getAppResPhs(data))[:,0]
+
+    return np.linalg.norm(np.abs(np.log(app_r) - np.log(np.ones(survey.nFreq)/sigmaHalf))*np.log(sigmaHalf))
+
+def appPhs_TotalFieldNorm(sigmaHalf):
+
+    # Make the survey
+    survey, sigma, mesh = NSEM.Utils.testUtils.setup1DSurvey(sigmaHalf)
+    problem = NSEM.Problem1D_eTotal(mesh)
+    problem.pair(survey)
+
+    # Get the fields
+    fields = problem.fields(sigma)
+
+    # Project the data
+    data = survey.eval(fields)
+
+    # Calculate the app  phs
+    app_p = np.array(NSEM.Utils.testUtils.getAppResPhs(data))[:,1]
+
+    return np.linalg.norm(np.abs(app_p - np.ones(survey.nFreq)*45)/ 45)
+
+def appRes_psFieldNorm(sigmaHalf):
+
+    # Make the survey
+    survey, sigma, mesh = NSEM.Utils.testUtils.setup1DSurvey(sigmaHalf,False)
+    problem = NSEM.Problem1D_ePrimSec(mesh, sigmaPrimary = sigma)
+    problem.pair(survey)
+
+    # Get the fields
+    fields = problem.fields(sigma)
+
+    # Project the data
+    data = survey.eval(fields)
+
+    # Calculate the app res and phs
+    app_r = np.array(NSEM.Utils.testUtils.getAppResPhs(data))[:,0]
+
+    return np.linalg.norm(np.abs(np.log(app_r) - np.log(np.ones(survey.nFreq)/sigmaHalf))*np.log(sigmaHalf))
+
+def appPhs_psFieldNorm(sigmaHalf):
+
+    # Make the survey
+    survey, sigma, mesh = NSEM.Utils.testUtils.setup1DSurvey(sigmaHalf,False)
+    problem = NSEM.Problem1D_ePrimSec(mesh, sigmaPrimary = sigma)
+    problem.pair(survey)
+
+    # Get the fields
+    fields = problem.fields(sigma)
+
+    # Project the data
+    data = survey.eval(fields)
+
+    # Calculate the app  phs
+    app_p = np.array(NSEM.Utils.testUtils.getAppResPhs(data))[:,1]
+
+    return np.linalg.norm(np.abs(app_p - np.ones(survey.nFreq)*45)/ 45)
+
+class TestAnalytics(unittest.TestCase):
+
+    def setUp(self):
+        pass
+    # Total Fields
+    # def test_appRes2en1(self):self.assertLess(appRes_TotalFieldNorm(2e-1), TOLr)
+    # def test_appPhs2en1(self):self.assertLess(appPhs_TotalFieldNorm(2e-1), TOLp)
+
+    # Primary/secondary
+    def test_appRes1en0_ps(self):self.assertLess(appRes_psFieldNorm(1e-0), TOLr)
+    def test_appPhs1en0_ps(self):self.assertLess(appPhs_psFieldNorm(1e-0), TOLp)
+    def test_appRes2en1_ps(self):self.assertLess(appRes_psFieldNorm(2e-1), TOLr)
+    def test_appPhs2en1_ps(self):self.assertLess(appPhs_psFieldNorm(2e-1), TOLp)
+    def test_appRes2en3_ps(self):self.assertLess(appRes_psFieldNorm(2e-3), TOLr)
+    def test_appPhs2en3_ps(self):self.assertLess(appPhs_psFieldNorm(2e-3), TOLp)
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/tests/mt/forward/test_Problem3D_VsAnalyticSolution.py b/tests/mt/forward/test_Problem3D_VsAnalyticSolution.py
new file mode 100644
index 00000000..b736bc17
--- /dev/null
+++ b/tests/mt/forward/test_Problem3D_VsAnalyticSolution.py
@@ -0,0 +1,54 @@
+# Test functions
+from glob import glob
+import numpy as np, sys, os, time, scipy, subprocess
+import SimPEG as simpeg
+import unittest
+from SimPEG import NSEM
+from SimPEG.Utils import meshTensor
+from scipy.constants import mu_0
+
+np.random.seed(1100)
+
+TOLr = 1
+TOLp = 2
+FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order
+CONDUCTIVITY = 1e1
+MU = mu_0
+freq = [1e-1, 2e-1]
+addrandoms = True
+
+def appResPhsHalfspace_eFrom_ps_Norm(sigmaHalf,appR=True,expMap=False):
+    if appR:
+        label = 'resistivity'
+    else:
+        label = 'phase'
+    print 'Apperent {:s} test of eFormulation primary/secondary at {:g}\n\n'.format(label,sigmaHalf)
+
+    # Calculate the app  phs
+    survey, problem = NSEM.Utils.testUtils.setupSimpegNSEM_ePrimSec(NSEM.Utils.testUtils.halfSpace(sigmaHalf),expMap=expMap)
+    data = problem.dataPair(survey,survey.dpred(problem.curModel))
+    recData = data.toRecArray('Complex')
+    app_rpxy, app_rpyx = NSEM.Utils.appResPhs(recData['freq'],recData['zxy'])[0], NSEM.Utils.appResPhs(recData['freq'],recData['zyx'])[0]
+    if appR:
+        return np.linalg.norm( np.abs(np.log10(app_rpxy[0]) - np.log10(1./sigmaHalf)) * np.log10(sigmaHalf ))
+    else:
+        return np.linalg.norm( np.abs(app_rpxy[1] + 135) / 135 )
+
+
+class TestAnalytics(unittest.TestCase):
+
+    def setUp(self):
+        # Make the survey and the problem
+        pass
+
+    # # Test apparent resistivity and phase
+    def test_appRes1en2(self):self.assertLess(appResPhsHalfspace_eFrom_ps_Norm(1e-2),TOLr)
+    def test_appPhs1en2(self):self.assertLess(appResPhsHalfspace_eFrom_ps_Norm(1e-2,False),TOLp)
+
+    def test_appRes1en1(self):self.assertLess(appResPhsHalfspace_eFrom_ps_Norm(1e-1),TOLr)
+    def test_appPhs1en1(self):self.assertLess(appResPhsHalfspace_eFrom_ps_Norm(1e-1,False),TOLp)
+
+
+
+if __name__ == '__main__':
+    unittest.main()
\ No newline at end of file
diff --git a/tests/mt/inversion/__init__.py b/tests/mt/inversion/__init__.py
new file mode 100644
index 00000000..420388ef
--- /dev/null
+++ b/tests/mt/inversion/__init__.py
@@ -0,0 +1,12 @@
+import os
+import glob
+import unittest
+
+if __name__ == '__main__':
+    test_file_strings = glob.glob('test_*.py')
+    module_strings = [str[0:len(str)-3] for str in test_file_strings]
+    suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str
+              in module_strings]
+    testSuite = unittest.TestSuite(suites)
+
+    unittest.TextTestRunner(verbosity=2).run(testSuite)
diff --git a/tests/mt/inversion/test_Problem3D_Adjoint.py b/tests/mt/inversion/test_Problem3D_Adjoint.py
new file mode 100644
index 00000000..2061e1ff
--- /dev/null
+++ b/tests/mt/inversion/test_Problem3D_Adjoint.py
@@ -0,0 +1,58 @@
+# Test functions
+from glob import glob
+import numpy as np, sys, os, time, scipy, subprocess
+import SimPEG as simpeg
+import unittest
+from SimPEG import NSEM
+from SimPEG.Utils import meshTensor
+from scipy.constants import mu_0
+
+
+TOLr = 5e-2
+TOL = 1e-4
+FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order
+CONDUCTIVITY = 1e1
+MU = mu_0
+freq = [1e-1, 2e-1]
+addrandoms = True
+
+
+
+def JvecAdjointTest(inputSetup,comp='All',freq=False):
+    (M, freqs, sig, sigBG, rx_loc) = inputSetup
+    survey, problem = NSEM.Utils.testUtils.setupSimpegNSEM_ePrimSec(inputSetup,comp='All',singleFreq=freq)
+    print 'Adjoint test of eForm primary/secondary for {:s} comp at {:s}\n'.format(comp,str(survey.freqs))
+
+    m  = sig
+    u = problem.fields(m)
+
+    v = np.random.rand(survey.nD,)
+    # print problem.PropMap.PropModel.nP
+    w = np.random.rand(problem.mesh.nC,)
+
+    vJw = v.ravel().dot(problem.Jvec(m, w, u))
+    wJtv = w.ravel().dot(problem.Jtvec(m, v, u))
+    tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR])
+    print ' vJw   wJtv  vJw - wJtv     tol    abs(vJw - wJtv) < tol'
+    print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol
+    return np.abs(vJw - wJtv) < tol
+
+
+class NSEM_AdjointTests(unittest.TestCase):
+
+    def setUp(self):
+        pass
+
+    # Test the adjoint of Jvec and Jtvec
+    # def test_JvecAdjoint_zxxr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxxr',.1))
+    # def test_JvecAdjoint_zxxi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxxi',.1))
+    # def test_JvecAdjoint_zxyr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxyr',.1))
+    # def test_JvecAdjoint_zxyi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxyi',.1))
+    # def test_JvecAdjoint_zyxr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyxr',.1))
+    # def test_JvecAdjoint_zyxi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyxi',.1))
+    # def test_JvecAdjoint_zyyr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyyr',.1))
+    # def test_JvecAdjoint_zyyi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyyi',.1))
+    def test_JvecAdjoint_All(self):self.assertTrue(JvecAdjointTest(NSEM.Utils.testUtils.random(1e-2),'All',.1))
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/tests/mt/inversion/test_Problem3D_Derivs.py b/tests/mt/inversion/test_Problem3D_Derivs.py
new file mode 100644
index 00000000..9e365ec9
--- /dev/null
+++ b/tests/mt/inversion/test_Problem3D_Derivs.py
@@ -0,0 +1,83 @@
+# Test functions
+from glob import glob
+import numpy as np, sys, os, time, scipy, subprocess
+import SimPEG as simpeg
+import unittest
+from SimPEG import NSEM
+from SimPEG.Utils import meshTensor
+from scipy.constants import mu_0
+
+np.random.seed(1100)
+
+TOLr = 5e-2
+TOL = 1e-4
+FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order
+CONDUCTIVITY = 1e1
+MU = mu_0
+freq = [1e-1, 2e-1]
+addrandoms = True
+
+
+# Test the Jvec derivative
+def DerivJvecTest(inputSetup,comp='All',freq=False,expMap=True):
+    (M, freqs, sig, sigBG, rx_loc) = inputSetup
+    survey, problem = NSEM.Utils.testUtils.setupSimpegNSEM_ePrimSec(inputSetup,comp=comp,singleFreq=freq,expMap=expMap)
+    print 'Derivative test of Jvec for eForm primary/secondary for {:s} comp at {:s}\n'.format(comp,survey.freqs)
+    # problem.mapping = simpeg.Maps.ExpMap(problem.mesh)
+    # problem.sigmaPrimary = np.log(sigBG)
+    x0 = np.log(sigBG)
+    # cond = sig[0]
+    # x0 = np.log(np.ones(problem.mesh.nC)*cond)
+    # problem.sigmaPrimary = x0
+    # if True:
+    #     x0  = x0 + np.random.randn(problem.mesh.nC)*cond*1e-1
+    survey = problem.survey
+    def fun(x):
+        return survey.dpred(x), lambda x: problem.Jvec(x0, x)
+    return simpeg.Tests.checkDerivative(fun, x0, num=3, plotIt=False, eps=FLR)
+
+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\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.eval(src,survey.mesh,f), lambda t: rx.evalDeriv(src,survey.mesh,f0,simpeg.mkvc(t,2))
+
+    return simpeg.Tests.checkDerivative(fun, u0, num=3, plotIt=False, eps=FLR)
+
+
+
+class NSEM_DerivTests(unittest.TestCase):
+
+    def setUp(self):
+        pass
+
+    # Do a derivative test of Jvec
+    # def test_derivJvec_zxxr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxxr',.1))
+    # def test_derivJvec_zxxi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxxi',.1))
+    # def test_derivJvec_zxyr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxyr',.1))
+    # def test_derivJvec_zxyi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxyi',.1))
+    # def test_derivJvec_zyxr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyxr',.1))
+    # def test_derivJvec_zyxi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyxi',.1))
+    # def test_derivJvec_zyyr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyyr',.1))
+    # def test_derivJvec_zyyi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyyi',.1))
+    def test_derivJvec_All(self):self.assertTrue(DerivJvecTest(NSEM.Utils.testUtils.random(1e-2),'All',.1))
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/tests/mt/test_Problem1D_againstAnalyticHalfspace.py b/tests/mt/test_Problem1D_againstAnalyticHalfspace.py
deleted file mode 100644
index d848c1ba..00000000
--- a/tests/mt/test_Problem1D_againstAnalyticHalfspace.py
+++ /dev/null
@@ -1,147 +0,0 @@
-import unittest
-import SimPEG as simpeg
-from SimPEG import NSEM
-from SimPEG.Utils import meshTensor
-import numpy as np
-# Define the tolerances
-TOLr = 5e-2
-TOLp = 5e-2
-
-
-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(NSEM.Rx(simpeg.mkvc(np.array([0.0]),2).T,rxType))
-    # Source list
-    srcList =[]
-    if tD:
-        for freq in freqs:
-            srcList.append(NSEM.SrcNSEM.polxy_1DhomotD(rxList,freq))
-    else:
-        for freq in freqs:
-            srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,freq))
-
-    survey = NSEM.Survey(srcList)
-    return survey, sigma, m1d
-
-def getAppResPhs(NSEMdata):
-    # Make impedance
-    def appResPhs(freq,z):
-        app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2
-        app_phs = np.arctan2(z.imag,z.real)*(180/np.pi)
-        return app_res, app_phs
-    zList = []
-    for src in NSEMdata.survey.srcList:
-        zc = [src.freq]
-        for rx in src.rxList:
-            if 'i' in rx.rxType:
-                m=1j
-            else:
-                m = 1
-            zc.append(m*NSEMdata[src,rx])
-        zList.append(zc)
-    return [appResPhs(zList[i][0],np.sum(zList[i][1:3])) for i in np.arange(len(zList))]
-
-def appRes_TotalFieldNorm(sigmaHalf):
-
-    # Make the survey
-    survey, sigma, mesh = setupSurvey(sigmaHalf)
-    problem = NSEM.Problem1D.eForm_TotalField(mesh)
-    problem.pair(survey)
-
-    # Get the fields
-    fields = problem.fields(sigma)
-
-    # Project the data
-    data = survey.eval(fields)
-
-    # Calculate the app res and phs
-    app_r = np.array(getAppResPhs(data))[:,0]
-
-    return np.linalg.norm(np.abs(app_r - np.ones(survey.nFreq)/sigmaHalf)*sigmaHalf)
-
-def appPhs_TotalFieldNorm(sigmaHalf):
-
-    # Make the survey
-    survey, sigma, mesh = setupSurvey(sigmaHalf)
-    problem = NSEM.Problem1D.eForm_TotalField(mesh)
-    problem.pair(survey)
-
-    # Get the fields
-    fields = problem.fields(sigma)
-
-    # Project the data
-    data = survey.eval(fields)
-
-    # Calculate the app  phs
-    app_p = np.array(getAppResPhs(data))[:,1]
-
-    return np.linalg.norm(np.abs(app_p - np.ones(survey.nFreq)*45)/ 45)
-
-def appRes_psFieldNorm(sigmaHalf):
-
-    # Make the survey
-    survey, sigma, mesh = setupSurvey(sigmaHalf,False)
-    problem = NSEM.Problem1D_ePrimSec(mesh, sigmaPrimary = sigma)
-    problem.pair(survey)
-
-    # Get the fields
-    fields = problem.fields(sigma)
-
-    # Project the data
-    data = survey.eval(fields)
-
-    # Calculate the app res and phs
-    app_r = np.array(getAppResPhs(data))[:,0]
-
-    return np.linalg.norm(np.abs(app_r - np.ones(survey.nFreq)/sigmaHalf)*sigmaHalf)
-
-def appPhs_psFieldNorm(sigmaHalf):
-
-    # Make the survey
-    survey, sigma, mesh = setupSurvey(sigmaHalf,False)
-    problem = NSEM.Problem1D_ePrimSec(mesh, sigmaPrimary = sigma)
-    problem.pair(survey)
-
-    # Get the fields
-    fields = problem.fields(sigma)
-
-    # Project the data
-    data = survey.eval(fields)
-
-    # Calculate the app  phs
-    app_p = np.array(getAppResPhs(data))[:,1]
-
-    return np.linalg.norm(np.abs(app_p - np.ones(survey.nFreq)*45)/ 45)
-
-class TestAnalytics(unittest.TestCase):
-
-    def setUp(self):
-        pass
-    # Total Fields
-    # def test_appRes2en1(self):self.assertLess(appRes_TotalFieldNorm(2e-1), TOLr)
-    # def test_appPhs2en1(self):self.assertLess(appPhs_TotalFieldNorm(2e-1), TOLp)
-
-    # Primary/secondary
-    def test_appRes2en2_ps(self):self.assertLess(appRes_psFieldNorm(2e-2), TOLr)
-    def test_appPhs2en2_ps(self):self.assertLess(appPhs_psFieldNorm(2e-2), TOLp)
-
-if __name__ == '__main__':
-    unittest.main()
diff --git a/tests/mt/test_Problem3D_againstAnalytic.py b/tests/mt/test_Problem3D_againstAnalytic.py
deleted file mode 100644
index a267c9d7..00000000
--- a/tests/mt/test_Problem3D_againstAnalytic.py
+++ /dev/null
@@ -1,270 +0,0 @@
-# Test functions
-from glob import glob
-import numpy as np, sys, os, time, scipy, subprocess
-import SimPEG as simpeg
-import unittest
-from SimPEG import NSEM
-from SimPEG.Utils import meshTensor
-from scipy.constants import mu_0
-
-np.random.seed(1100)
-
-TOLr = 5e-2
-TOL = 1e-4
-FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order
-CONDUCTIVITY = 1e1
-MU = mu_0
-freq = [1e-1, 2e-1]
-addrandoms = True
-
-
-def getInputs():
-    """
-    Function that returns Mesh, freqs, rx_loc, elev.
-    """
-    # 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])
-    # M = simpeg.Mesh.TensorMesh([[(1000,6,-1.5),(1000.,6),(1000,6,1.5)],[(1000,6,-1.5),(1000.,2),(1000,6,1.5)],[(1000,6,-1.3),(1000.,6),(1000,6,1.3)]], x0=['C','C','C'])# Setup the model
-    M = simpeg.Mesh.TensorMesh([[(1000,6,-1.5),(1000.,4),(1000,6,1.5)],[(1000,6,-1.5),(1000.,4),(1000,6,1.5)],[(500,8,-1.3),(500.,8),(500,8,1.3)]], x0=['C','C','C'])# Setup the model
-    # Set the frequencies
-    freqs = np.logspace(1,-3,5)
-    elev = 0
-
-    ## Setup the the survey object
-    # Receiver locations
-    rx_x, rx_y = np.meshgrid(np.arange(-1000,1001,500),np.arange(-1000,1001,500))
-    rx_loc = np.hstack((simpeg.Utils.mkvc(rx_x,2),simpeg.Utils.mkvc(rx_y,2),elev+np.zeros((np.prod(rx_x.shape),1))))
-
-    return M, freqs, rx_loc, elev
-
-def random(conds):
-    ''' Returns a halfspace model based on the inputs'''
-    M, freqs, rx_loc, elev = getInputs()
-
-    # Backround
-    sigBG = np.ones(M.nC)*conds
-    # Add randomness to the model (10% of the value).
-    sig = np.exp( np.log(sigBG) + np.random.randn(M.nC)*(conds)*1e-1 )
-
-    return (M, freqs, sig, sigBG, rx_loc)
-
-def halfSpace(conds):
-    ''' Returns a halfspace model based on the inputs'''
-    M, freqs, rx_loc, elev = getInputs()
-
-    # Model
-    ccM = M.gridCC
-    # conds = [1e-2]
-    groundInd = ccM[:,2] < elev
-    sig = np.zeros(M.nC) + 1e-8
-    sig[groundInd] = conds
-    # Set the background, not the same as the model
-    sigBG = np.zeros(M.nC) + 1e-8
-    sigBG[groundInd] = conds
-
-    return (M, freqs, sig, sigBG, rx_loc)
-
-def blockInhalfSpace(conds):
-    ''' Returns a halfspace model based on the inputs'''
-    M, freqs, rx_loc, elev = getInputs()
-
-    # Model
-    ccM = M.gridCC
-    # conds = [1e-2]
-    groundInd = ccM[:,2] < elev
-    sig = simpeg.Utils.ModelBuilder.defineBlock(M.gridCC,np.array([-1000,-1000,-1500]),np.array([1000,1000,-1000]),conds)
-    sig[~groundInd] = 1e-8
-    # Set the background, not the same as the model
-    sigBG = np.zeros(M.nC) + 1e-8
-    sigBG[groundInd] = conds[1]
-
-    return (M, freqs, sig, sigBG, rx_loc)
-
-def twoLayer(conds):
-    ''' Returns a 2 layer model based on the conductivity values given'''
-    M, freqs, rx_loc, elev = getInputs()
-
-    # Model
-    ccM = M.gridCC
-    groundInd = ccM[:,2] < elev
-    botInd = ccM[:,2] < -3000
-    sig = np.zeros(M.nC) + 1e-8
-    sig[groundInd] = conds[1]
-    sig[botInd] = conds[0]
-    # Set the background, not the same as the model
-    sigBG = np.zeros(M.nC) + 1e-8
-    sigBG[groundInd] = conds[1]
-
-
-    return (M, freqs, sig, sigBG, rx_loc)
-
-
-
-def setupSimpegNSEMfwd_eForm_ps(inputSetup,comp='Imp',singleFreq=False,expMap=True):
-    M,freqs,sig,sigBG,rx_loc = inputSetup
-    # Make a receiver list
-    rxList = []
-    if comp == 'All':
-        for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi','tzxr','tzxi','tzyr','tzyi']:
-            rxList.append(NSEM.Rx(rx_loc,rxType))
-    elif comp == 'Imp':
-        for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi']:
-            rxList.append(NSEM.Rx(rx_loc,rxType))
-    elif comp == 'Tip':
-        for rxType in ['tzxr','tzxi','tzyr','tzyi']:
-            rxList.append(NSEM.Rx(rx_loc,rxType))
-    else:
-        rxList.append(NSEM.Rx(rx_loc,comp))
-    # Source list
-    srcList =[]
-
-    if singleFreq:
-        srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,singleFreq))
-    else:
-        for freq in freqs:
-            srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,freq))
-    # Survey NSEM
-    survey = NSEM.Survey(srcList)
-
-    ## Setup the problem object
-    sigma1d = M.r(sigBG,'CC','CC','M')[0,0,:]
-    if expMap:
-        problem = NSEM.Problem3D_ePrimSec(M,sigmaPrimary= np.log(sigma1d) )
-        problem.mapping = simpeg.Maps.ExpMap(problem.mesh)
-        problem.curModel = np.log(sig)
-    else:
-        problem = NSEM.Problem3D_ePrimSec(M,sigmaPrimary= sigma1d)
-        problem.curModel = sig
-    problem.pair(survey)
-    problem.verbose = False
-    try:
-        from pymatsolver import MumpsSolver
-        problem.Solver = MumpsSolver
-    except:
-        pass
-
-    return (survey, problem)
-
-def getAppResPhs(NSEMdata):
-    # Make impedance
-    def appResPhs(freq,z):
-        app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2
-        app_phs = np.arctan2(z.imag,z.real)*(180/np.pi)
-        return app_res, app_phs
-    recData = NSEMdata.toRecArray('Complex')
-    return appResPhs(recData['freq'],recData['zxy']), appResPhs(recData['freq'],recData['zyx'])
-
-def JvecAdjointTest(inputSetup,comp='All',freq=False):
-    (M, freqs, sig, sigBG, rx_loc) = inputSetup
-    survey, problem = setupSimpegNSEMfwd_eForm_ps(inputSetup,comp='All',singleFreq=freq)
-    print 'Adjoint test of eForm primary/secondary for {:s} comp at {:s}\n'.format(comp,str(survey.freqs))
-
-    m  = sig
-    u = problem.fields(m)
-
-    v = np.random.rand(survey.nD,)
-    # print problem.PropMap.PropModel.nP
-    w = np.random.rand(problem.mesh.nC,)
-
-    vJw = v.ravel().dot(problem.Jvec(m, w, u))
-    wJtv = w.ravel().dot(problem.Jtvec(m, v, u))
-    tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR])
-    print ' vJw   wJtv  vJw - wJtv     tol    abs(vJw - wJtv) < tol'
-    print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol
-    return np.abs(vJw - wJtv) < tol
-
-# Test the Jvec derivative
-def DerivJvecTest(inputSetup,comp='All',freq=False,expMap=True):
-    (M, freqs, sig, sigBG, rx_loc) = inputSetup
-    survey, problem = setupSimpegNSEMfwd_eForm_ps(inputSetup,comp=comp,singleFreq=freq,expMap=expMap)
-    print 'Derivative test of Jvec for eForm primary/secondary for {:s} comp at {:s}\n'.format(comp,survey.freqs)
-    # problem.mapping = simpeg.Maps.ExpMap(problem.mesh)
-    # problem.sigmaPrimary = np.log(sigBG)
-    x0 = np.log(sigBG)
-    # cond = sig[0]
-    # x0 = np.log(np.ones(problem.mesh.nC)*cond)
-    # problem.sigmaPrimary = x0
-    # if True:
-    #     x0  = x0 + np.random.randn(problem.mesh.nC)*cond*1e-1
-    survey = problem.survey
-    def fun(x):
-        return survey.dpred(x), lambda x: problem.Jvec(x0, x)
-    return simpeg.Tests.checkDerivative(fun, x0, num=3, plotIt=False, eps=FLR)
-
-def DerivProjfieldsTest(inputSetup,comp='All',freq=False):
-
-    survey, problem = setupSimpegNSEMfwd_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.eval(src,survey.mesh,f), lambda t: rx.evalDeriv(src,survey.mesh,f0,simpeg.mkvc(t,2))
-
-    return simpeg.Tests.checkDerivative(fun, u0, num=3, plotIt=False, eps=FLR)
-
-def appResPhsHalfspace_eFrom_ps_Norm(sigmaHalf,appR=True,expMap=False):
-    if appR:
-        label = 'resistivity'
-    else:
-        label = 'phase'
-    # Make the survey and the problem
-    survey, problem = setupSimpegNSEMfwd_eForm_ps(halfSpace(sigmaHalf),expMap=expMap)
-    print 'Apperent {:s} test of eFormulation primary/secondary at {:g}\n\n'.format(label,sigmaHalf)
-
-    data = problem.dataPair(survey,survey.dpred(problem.curModel))
-    # Calculate the app  phs
-    app_rpxy, app_rpyx = np.array(getAppResPhs(data))
-    if appR:
-        return np.all(np.abs(app_rpxy[0,:] - 1./sigmaHalf) * sigmaHalf < .4)
-    else:
-        return np.all(np.abs(app_rpxy[1,:] + 135) / 135 < .4)
-
-class TestAnalytics(unittest.TestCase):
-
-    def setUp(self):
-        pass
-    # # Test apparent resistivity and phase
-    def test_appRes1en2(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-2))
-    def test_appPhs1en2(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-2,False))
-
-    def test_appRes1en3(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-3))
-    def test_appPhs1en3(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-3,False))
-
-    # Do a derivative test of Jvec
-    # def test_derivJvec_zxxr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxxr',.1))
-    # def test_derivJvec_zxxi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxxi',.1))
-    # def test_derivJvec_zxyr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxyr',.1))
-    # def test_derivJvec_zxyi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxyi',.1))
-    # def test_derivJvec_zyxr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyxr',.1))
-    # def test_derivJvec_zyxi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyxi',.1))
-    # def test_derivJvec_zyyr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyyr',.1))
-    # def test_derivJvec_zyyi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyyi',.1))
-    def test_derivJvec_All(self):self.assertTrue(DerivJvecTest(random(1e-2),'All',.1))
-
-    # Test the adjoint of Jvec and Jtvec
-    # def test_JvecAdjoint_zxxr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxxr',.1))
-    # def test_JvecAdjoint_zxxi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxxi',.1))
-    # def test_JvecAdjoint_zxyr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxyr',.1))
-    # def test_JvecAdjoint_zxyi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxyi',.1))
-    # def test_JvecAdjoint_zyxr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyxr',.1))
-    # def test_JvecAdjoint_zyxi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyxi',.1))
-    # def test_JvecAdjoint_zyyr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyyr',.1))
-    # def test_JvecAdjoint_zyyi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyyi',.1))
-    def test_JvecAdjoint_All(self):self.assertTrue(JvecAdjointTest(random(1e-2),'All',.1))
-
-if __name__ == '__main__':
-    unittest.main()

From d1f7d0c37a18f2c5e44fa8bbdc2126f94228b172 Mon Sep 17 00:00:00 2001
From: GudniRos <gkrosen@gmail.com>
Date: Wed, 1 Jun 2016 21:54:19 -0700
Subject: [PATCH 13/14] Adding ipywidgets to the conda install list for travis

---
 .travis.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.travis.yml b/.travis.yml
index a4614a6b..e91ae24c 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -35,7 +35,7 @@ before_install:
 
 # Install packages
 install:
-  - conda install --yes pip python=$TRAVIS_PYTHON_VERSION numpy scipy matplotlib cython ipython nose vtk
+  - conda install --yes pip python=$TRAVIS_PYTHON_VERSION numpy scipy matplotlib cython ipython ipywidgets nose vtk
   - pip install nose-cov python-coveralls
 
   - git clone https://github.com/rowanc1/pymatsolver.git

From cf5181016fe261c1c2d3abdfedfe70073e76b5b6 Mon Sep 17 00:00:00 2001
From: GudniRos <gkrosen@gmail.com>
Date: Tue, 7 Jun 2016 09:45:02 -0700
Subject: [PATCH 14/14] Made MT_1d_analytic_nLayer_Earth.run have default n
 layers of 3. Trying to fix example testing.

---
 SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py b/SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py
index 62ab2a68..b8edd9bb 100644
--- a/SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py
+++ b/SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py
@@ -399,7 +399,7 @@ def PlotAppRes3LayersInteract(h1,h2,sigl1,sigl2,sigl3,mul1,mul2,mul3,epsl1,epsl2
     PlotAppRes(frangn,thick3,sig3,chg3_0,taux3,c3,mu3,eps3,3,F_Envelope,PlotEnvelope)
 
 
-def run(n,plotIt=True):
+def run(n=3,plotIt=True):
     # something to make a plot
 
     F = frange(-5.,5.,20)
@@ -420,7 +420,7 @@ def run(n,plotIt=True):
     return Res, Phase
 
 if __name__ == '__main__':
-    run(3)
+    run()