diff --git a/simpegPF/Dev/MAG3D_fwr_Test.py b/simpegPF/Dev/MAG3D_fwr_Test.py
index aa97e332..1d30dbb6 100644
--- a/simpegPF/Dev/MAG3D_fwr_Test.py
+++ b/simpegPF/Dev/MAG3D_fwr_Test.py
@@ -13,86 +13,155 @@ import matplotlib
 
 from fwr_MAG_obs import fwr_MAG_obs
 
+plt.close('all')
+
 #%% Create survey
 
-B = np.array(([90.,0.,50000.]))
+B = np.array(([-45.,315.,50000.]))
   
-M = np.array(([90.,0.]))
+M = np.array(([-45.,315.]))
  
+# Sphere radius
+R = 0.25
        
 # # Or create juste a plane grid
-xr = np.linspace(-1./2., 1./2., 10)
-yr = np.linspace(-1./2., 1./2., 10)
+xr = np.linspace(-2., 2., 5)
+yr = np.linspace(-2., 2., 5)
 X, Y = np.meshgrid(xr, yr)
-Z = np.ones((xr.size, yr.size))*.75
+Z = np.ones((xr.size, yr.size)) * 2.5
 rxLoc = np.c_[Utils.mkvc(X), Utils.mkvc(Y), Utils.mkvc(Z)]
 
 ndata = rxLoc.shape[0]
 
-#%%
+d_iter = 4
+lrl = np.zeros(d_iter)
+#%% Loop through decreasing meshes and measure the residual
 # Create mesh using simpeg and write out in GIF format
-nc = 30.
 
-hxind = [(1./nc, nc)]
-hyind = [(1./nc, nc)]
-hzind = [(1./nc, nc)]
+for ii in range(d_iter):
+    
+    nc = 3**(ii+1)
+    
+    hxind = [(1./nc, nc)]
+    hyind = [(1./nc, nc)]
+    hzind = [(1./nc, nc)]
+    
+    mesh = Mesh.TensorMesh([hxind, hyind, hzind], 'CCC')
+    
+    xn = mesh.vectorNx
+    yn = mesh.vectorNy
+    zn = mesh.vectorNz
+    
+    mcell = mesh.nC
+    
+    sph_ind = PF.MagAnalytics.spheremodel(mesh, 0, 0, 0, R)
+    
+    chibkg = 0.
+    chiblk = 0.01
+    model = np.ones(mcell)*chibkg
+    model[sph_ind] = chiblk
+    
+    #%% Forward mode ldata
+    d = fwr_MAG_obs(mesh,B,M,rxLoc,model)
+    fwr_x = mkvc(d[0,:])
+    fwr_y = mkvc(d[1,:])
+    fwr_z = mkvc(d[2,:])
+    
+    #%% Get the analystical answer and compute the residual
+    bxa,bya,bza = PF.MagAnalytics.MagSphereAnaFunA(rxLoc[:,0],rxLoc[:,1],rxLoc[:,2],R,0.,0.,0.,chiblk, np.array(([0.,0.,B[2]])),'secondary')
+    Bd = (450.-float(B[1]))%360.
+    Bi = B[0]; # Convert dip to horizontal to cartesian 
+    
+    Bx = np.cos(np.deg2rad(Bi)) * np.cos(np.deg2rad(Bd)) * B[2] 
+    By = np.cos(np.deg2rad(Bi)) * np.sin(np.deg2rad(Bd)) * B[2] 
+    Bz = np.sin(np.deg2rad(Bi)) * B[2] 
+    
+    Bo = np.c_[Bx, By, Bz]
+        
+    bxa,bya,bza = PF.MagAnalytics.MagSphereFreeSpace(rxLoc[:,0],rxLoc[:,1],rxLoc[:,2],R,0.,0.,0.,chiblk, Bo)
+    #bxa,bya,bza = PF.MagAnalytics.MagSphereAnaFunA(rxLoc[:,0],rxLoc[:,1],rxLoc[:,2],R,0.,0.,0.,chiblk, np.array(([0.,0.,B[2]])),'secondary')
+    
+    r_Bx = fwr_x - bxa
+    r_By = fwr_y - bya
+    r_Bz = fwr_z - bza
+    
+    lrl[ii] = sum( r_Bx**2 + r_By**2 + r_Bz**2 ) **0.5
+    
 
-mesh = Mesh.TensorMesh([hxind, hyind, hzind], 'CCC')
-
-xn = mesh.vectorNx
-yn = mesh.vectorNy
-zn = mesh.vectorNz
-
-mcell = mesh.nC
-
-sph_ind = PF.MagAnalytics.spheremodel(mesh, 0, 0, 0, 0.25)
-
-chibkg = 0.
-chiblk = 0.01
-model = np.ones(mcell)*chibkg
-model[sph_ind] = chiblk
-
-#%% Forward mode ldata
-d = fwr_MAG_obs(mesh,B,M,rxLoc,model)
-
-#%% Get the analystical answer and compute the residual
-bxa,bya,bza = PF.MagAnalytics.MagSphereAnaFunA(rxLoc[:,0],rxLoc[:,1],rxLoc[:,2],.25,0.,0.,0.,chiblk, np.array(([0.,0.,B[2]])),'secondary')
-
-r_Bz = mkvc(d) - bza
-lrl = sum( r_Bz**2 ) **0.5
-
-print "Residual between analytical and integral= " + str(lrl)
+    
+#%% Plot results
+for ii in range(d_iter):
+    nc = 3**(ii+1)
+    print "Residual= " + str(lrl[ii]) + "\t dx= " + str(1./nc)
+    
 #%% Plot fields
 
 plt.figure(1)
 ax = plt.subplot(131)
-plt.imshow(np.reshape(bxa,X.shape), interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()])
-plt.contour(X,Y, np.reshape(bxa,X.shape),10)
-plt.scatter(X,Y, c='k', s=5)
+plt.imshow(np.reshape(bxa,X.shape).T, interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()], origin = 'lower')
+plt.colorbar()
+plt.contour(X,Y, np.reshape(bxa,X.shape).T,10)
+plt.scatter(X,Y, c=np.reshape(bxa,X.shape).T, s=20)
 ax.set_title('Sphere Ana Bx')
 
 ax = plt.subplot(132)
-plt.imshow(np.reshape(bya,X.shape), interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()])
-plt.contour(X,Y, np.reshape(bya,X.shape),10)
-plt.scatter(X,Y, c='k', s=5)
+plt.imshow(np.reshape(bya,X.shape).T, interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()], origin = 'lower')
+plt.colorbar()
+plt.contour(X,Y, np.reshape(bya,X.shape).T,10)
+plt.scatter(X,Y, c=np.reshape(bya,X.shape).T, s=20)
 ax.set_title('Sphere Ana By')
 
 ax = plt.subplot(133)
-plt.imshow(np.reshape(bza,X.shape), interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()])
-plt.contour(X,Y, np.reshape(bza,X.shape),10)
-plt.scatter(X,Y, c='k', s=5)
+plt.imshow(np.reshape(bza,X.shape).T, interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()], origin = 'lower')
+plt.colorbar()
+plt.contour(X,Y, np.reshape(bza,X.shape).T,10)
+plt.scatter(X,Y, c=np.reshape(bza,X.shape).T, s=20)
 ax.set_title('Sphere Ana Bz')
-#%% Plot foward data
+
+#%% Plot the forward solution from integral
+
 plt.figure(2)
-plt.subplot(121)
-d2D = np.reshape(d,X.shape)
-plt.imshow(d2D, interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()])
-plt.contour(X,Y, d2D,10)
-plt.scatter(X,Y, c='k', s=5)
+ax = plt.subplot(131)
+plt.imshow(np.reshape(fwr_x,X.shape).T, interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max() ], origin = 'lower')
+plt.colorbar()
+plt.contour(X,Y, np.reshape(fwr_x,X.shape).T,10)
+plt.scatter(X,Y, c=np.reshape(fwr_x,X.shape).T, s=20)
+ax.set_title('Sphere Ana Bx')
 
-#%% Compare fields
+ax = plt.subplot(132)
+plt.imshow(np.reshape(fwr_y,X.shape).T, interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()], origin = 'lower')
+plt.colorbar()
+plt.contour(X,Y, np.reshape(fwr_y,X.shape).T,10)
+plt.scatter(X,Y, c=np.reshape(fwr_y,X.shape).T, s=20)
+ax.set_title('Sphere Ana By')
 
-plt.subplot(122)
-plt.imshow(np.reshape(r_Bz,X.shape), interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()])
-plt.contour(X,Y, np.reshape(dBz,X.shape),10)
-plt.scatter(X,Y, c='k', s=5)
\ No newline at end of file
+ax = plt.subplot(133)
+plt.imshow(np.reshape(fwr_z,X.shape).T, interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()], origin = 'lower')
+plt.colorbar()
+plt.contour(X,Y, np.reshape(fwr_z,X.shape).T,10)
+plt.scatter(X,Y, c=np.reshape(fwr_z,X.shape).T, s=20)
+ax.set_title('Sphere Ana Bz')
+
+
+#%% Plot foward data
+plt.figure(3)
+ax = plt.subplot(131)
+plt.imshow(np.reshape(r_Bx,X.shape).T, interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()], origin = 'lower')
+plt.colorbar()
+plt.contour(X,Y, np.reshape(r_Bx,X.shape).T,10)
+plt.scatter(X,Y, c=np.reshape(r_Bx,X.shape).T, s=20)
+ax.set_title('Sphere Ana Bx')
+
+ax = plt.subplot(132)
+plt.imshow(np.reshape(r_By,X.shape).T, interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()], origin = 'lower')
+plt.colorbar()
+plt.contour(X,Y, np.reshape(r_By,X.shape).T,10)
+plt.scatter(X,Y, c=np.reshape(r_By,X.shape).T, s=20)
+ax.set_title('Sphere Ana By')
+
+ax = plt.subplot(133)
+plt.imshow(np.reshape(r_Bz,X.shape).T, interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()], origin = 'lower')
+plt.colorbar()
+plt.contour(X,Y, np.reshape(r_Bz,X.shape).T,10)
+plt.scatter(X,Y, c=np.reshape(r_Bz,X.shape).T, s=20)
+ax.set_title('Sphere Ana Bz')
\ No newline at end of file
diff --git a/simpegPF/Dev/fwr_MAG_obs.py b/simpegPF/Dev/fwr_MAG_obs.py
index 623ea89a..72fb8eff 100644
--- a/simpegPF/Dev/fwr_MAG_obs.py
+++ b/simpegPF/Dev/fwr_MAG_obs.py
@@ -53,7 +53,7 @@ def fwr_MAG_obs(mesh,B,M,rxLoc,model):
     
     Ptmi = mkvc(np.r_[np.cos(np.deg2rad(B[0]))*np.cos(np.deg2rad(D)),np.cos(np.deg2rad(B[0]))*np.sin(np.deg2rad(D)),np.sin(np.deg2rad(B[0]))],2).T;
     
-    d = np.zeros((ndata,1))
+    d = np.zeros((3,ndata))
     
     # Loop through all observations and create forward operator (ndata-by-mcell)
     print "Begin forward modeling " +str(int(ndata)) + " data points..."
@@ -64,10 +64,11 @@ def fwr_MAG_obs(mesh,B,M,rxLoc,model):
     
         tx, ty, tz = get_T_mat(xn,yn,zn,rxLoc[ii,:])  
         
-        G = Ptmi.dot(np.vstack((tx,ty,tz)))*Mxyz
+        # G = Ptmi.dot(np.vstack((tx,ty,tz)))*Mxyz
+        Gxyz = np.vstack((tx,ty,tz))*Mxyz
         #%%
         # Forward operator
-        d[ii,0] = G.dot(model)
+        d[:,ii] = Gxyz.dot(model)
     
         d_iter = np.floor(float(ii)/float(ndata)*10.);
         
diff --git a/simpegPF/MagAnalytics.py b/simpegPF/MagAnalytics.py
index 34034bae..bea7d508 100644
--- a/simpegPF/MagAnalytics.py
+++ b/simpegPF/MagAnalytics.py
@@ -192,6 +192,53 @@ def IDTtoxyz(Inc, Dec, Btot):
 
     return np.r_[Bx, By, Bz]
 
+def MagSphereFreeSpace(x, y, z, R, xc, yc, zc, chi, Bo):
+    """
+        Computing boundary condition using Congrous sphere method.
+        This is designed for secondary field formulation.
+        >> Input
+        mesh:   Mesh class
+        Bo:     np.array([Box, Boy, Boz]): Primary magnetic flux
+        Chi:    susceptibility at cell volume
+
+        .. math::
+
+            \\vec{B}(r) = \\frac{\mu_0}{4\pi}\\frac{m}{\| \\vec{r}-\\vec{r}_0\|^3}[3\hat{m}\cdot\hat{r}-\hat{m}]
+
+    """
+    if (~np.size(x)==np.size(y)==np.size(z)):
+        print "Specify same size of x, y, z"
+        return
+
+    x = Utils.mkvc(x)
+    y = Utils.mkvc(y)
+    z = Utils.mkvc(z)
+
+    nobs = len(x)
+    
+    Bot = np.sqrt(sum(Bo**2))
+    
+    mx = np.ones([nobs]) * Bo[0,0] *  R**3 / 3. * chi
+    my = np.ones([nobs]) * Bo[0,1] *  R**3 / 3. * chi
+    mz = np.ones([nobs]) * Bo[0,2] *  R**3 / 3. * chi
+
+    M = np.c_[mx, my, mz]
+    
+    rx = (x - xc)
+    ry = (y - yc)
+    rz = (zc - z)
+    
+    rvec = np.c_[rx, ry, rz]
+    r = np.sqrt((rx)**2+(ry)**2+(rz)**2 )
+
+    B = -Utils.sdiag(1./r**3)*M + Utils.sdiag((3 * np.sum(M*rvec,axis=1))/r**5)*rvec
+
+    Bx = B[:,0]
+    By = B[:,1]
+    Bz = B[:,2]
+
+    return Bx, By, Bz
+    
 if __name__ == '__main__':
 
     hxind = [(0,25,1.3),(21, 12.5),(0,25,1.3)]