diff --git a/simpegPF/Dev/Intgrl_MAG_Fwr_SphereTest.py b/simpegPF/Dev/Intgrl_MAG_Fwr_SphereTest.py index fbf48ae7..4aed811b 100644 --- a/simpegPF/Dev/Intgrl_MAG_Fwr_SphereTest.py +++ b/simpegPF/Dev/Intgrl_MAG_Fwr_SphereTest.py @@ -1,8 +1,8 @@ import os -home_dir = 'C:\Users\dominiquef.MIRAGEOSCIENCE\Documents\GIT\SimPEG\simpegpf\simpegPF\Dev' +# home_dir = 'C:\Users\dominiquef.MIRAGEOSCIENCE\Documents\GIT\SimPEG\simpegpf\simpegPF\Dev' -os.chdir(home_dir) +# os.chdir(home_dir) #%% from SimPEG import * @@ -16,12 +16,12 @@ plt.close('all') #%% Create survey B = np.array(([-45.,315.,50000.])) - + M = np.array(([-45.,315.])) - + # Sphere radius R = 0.25 - + # # Or create juste a plane grid xr = np.linspace(-2., 2., 5) yr = np.linspace(-2., 2., 5) @@ -37,67 +37,67 @@ lrl = np.zeros(d_iter) # Create mesh using simpeg and write out in GIF format 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 - + print 'Mesh size: ' + str(mcell) - + sph_ind = PF.MagAnalytics.spheremodel(mesh, 0, 0, 0, R) - + chibkg = 0. chiblk = 0.01 model = np.ones(mcell)*chibkg model[sph_ind] = chiblk - + actv = np.ones(mcell) - + #%% Forward mode ldata d = PF.Magnetics.Intgrl_Fwr_Data(mesh,B,M,rxLoc,model,actv,'xyz') fwr_x = d[0:ndata] fwr_y = d[ndata:2*ndata] fwr_z = d[2*ndata:] - + #%% 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] - + 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 - - + lrl[ii] = sum( r_Bx**2 + r_By**2 + r_Bz**2 ) **0.5 + + + #%% Plot results -print 'Residual between analytical sphere and integral forward' +print 'Residual between analytical sphere and integral forward' for ii in range(d_iter): nc = 3**(ii+1) print "||r||= " + str(lrl[ii]) + "\t dx= " + str(1./nc) - + #%% Plot fields plt.figure(1) @@ -168,4 +168,6 @@ plt.imshow(np.reshape(r_Bz,X.shape).T, interpolation="bicubic", extent=[xr.min() plt.colorbar(fraction=0.04) 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 +ax.set_title('Sphere Ana Bz') + +plt.show() diff --git a/simpegPF/Magnetics.py b/simpegPF/Magnetics.py index 0a167674..92996c3b 100644 --- a/simpegPF/Magnetics.py +++ b/simpegPF/Magnetics.py @@ -59,7 +59,7 @@ class MagneticsDiffSecondary(Problem.BaseProblem): .. math :: \mathbf{rhs} = \Div(\MfMui)^{-1}\mathbf{M}^f_{\mu_0^{-1}}\mathbf{B}_0 - \Div\mathbf{B}_0+\diag(v)\mathbf{D} \mathbf{P}_{out}^T \mathbf{B}_{sBC} - + """ B0 = self.getB0() Dface = self.mesh.faceDiv @@ -68,13 +68,13 @@ class MagneticsDiffSecondary(Problem.BaseProblem): mu = self.mapping*m chi = mu/mu_0-1 - + #temporary fix Bbc, Bbc_const = CongruousMagBC(self.mesh, self.survey.B0, chi) self.Bbc = Bbc self.Bbc_const = Bbc_const # return self._Div*self.MfMuI*self.MfMu0*B0 - self._Div*B0 + Mc*Dface*self._Pout.T*Bbc - return self._Div*self.MfMuI*self.MfMu0*B0 - self._Div*B0 + return self._Div*self.MfMuI*self.MfMu0*B0 - self._Div*B0 def getA(self, m): """ @@ -410,13 +410,13 @@ if __name__ == '__main__': def Intgrl_Fwr_Data(mesh,B,M,rxLoc,model,actv,flag): - """ + """ Forward model magnetic data using integral equation - + INPUT: - xn, yn, zn = Mesh nodes location + mesh = SimPEG.TensorMesh B = Inducing field parameter [Binc, Bdecl, B0] - M = Magnetization matrix [Minc, Mdecl] + M = Magnetization matrix [Minc, Mdecl] -90:90, 0:360 rxLox = Observation location informat [obsx, obsy, obsz] model = Model associated with mesh actv = Active cells from topo (from getActiveTopo) @@ -428,111 +428,118 @@ def Intgrl_Fwr_Data(mesh,B,M,rxLoc,model,actv,flag): Created on Oct 7, 2015 @author: dominiquef - """ - + """ + inds = np.nonzero(actv)[0] P = sp.csr_matrix((np.ones(inds.size),(inds, range(inds.size))), shape=(mesh.nC, len(inds))) xn = mesh.vectorNx; yn = mesh.vectorNy; zn = mesh.vectorNz; - + yn2,xn2,zn2 = np.meshgrid(yn[1:], xn[1:], zn[1:]) yn1,xn1,zn1 = np.meshgrid(yn[0:-1], xn[0:-1], zn[0:-1]) - + Yn = P.T*np.c_[mkvc(yn1), mkvc(yn2)] Xn = P.T*np.c_[mkvc(xn1), mkvc(xn2)] Zn = P.T*np.c_[mkvc(zn1), mkvc(zn2)] - + mcell = len(inds) - - ndata = rxLoc.shape[0] - + + ndata = rxLoc.shape[0] + # Convert declination from north to cartesian Md = (450.-float(M[1]))%360. - + # Create magnetization matrix mx = np.cos(np.deg2rad(M[0])) * np.cos(np.deg2rad(Md)) my = np.cos(np.deg2rad(M[0])) * np.sin(np.deg2rad(Md)) mz = np.sin(np.deg2rad(M[0])) - + Mx = Utils.sdiag(np.ones([mcell])*mx*B[2]) My = Utils.sdiag(np.ones([mcell])*my*B[2]) Mz = Utils.sdiag(np.ones([mcell])*mz*B[2]) - + #matplotlib.pyplot.spy(scipy.sparse.csr_matrix(Mx)) #plt.show() Mxyz = sp.vstack((Mx,My,Mz)); - + #%% Create TMI projector - - # Convert Bdecination from north to cartesian + + # Convert Bdecination from north to cartesian D = (450.-float(B[1]))%360. - - 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; - + + if flag=='tmi': + 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) - + elif flag=='xyz': d = np.zeros(int(3*ndata)) - + # Loop through all observations and create forward operator (ndata-by-mcell) print "Begin forward modeling " +str(int(ndata)) + " data points..." - + # Add counter to dsiplay progress. count = -1 for ii in range(ndata): - - tx, ty, tz = get_T_mat(Xn,Yn,Zn,rxLoc[ii,:]) + + tx, ty, tz = get_T_mat(Xn,Yn,Zn,rxLoc[ii,:]) Gxyz = np.vstack((tx,ty,tz))*Mxyz - + + # Remove non-active cells if flag=='xyz': d[ii::ndata] = mkvc(Gxyz.dot(P.T*model)) - + elif flag=='tmi': d[ii] = Ptmi.dot(Gxyz.dot(P.T*model)) - + # Display progress count = progress(ii,count,ndata) - - + + print "Done 100% ...forward modeling completed!!\n" - + return d +<<<<<<< HEAD def Intrgl_Fwr_Op(mesh,B,M,rxLoc,actv,flag): """ +======= +def Intrgl_Fwr_Op(mesh,B,M,rxLoc,flag): + """ +>>>>>>> dbcd57bc4e5d0d6258a10b760966b58419e0a9d2 Magnetic forward operator in integral form - + INPUT: mesh = Mesh in SimPEG format B = Inducing field parameter [Binc, Bdecl, B0] - M = Magnetization information + M = Magnetization information [OPTIONS] 1- [Minc, Mdecl] : Assumes uniform magnetization orientation 2- [mx1,mx2,..., my1,...,mz1] : cell-based defined magnetization direction 3- diag(M): Block diagonal matrix with [Mx, My, Mz] along the diagonal - + rxLox = Observation location informat [obsx, obsy, obsz] - + flag = 'tmi' | 'xyz' | 'full' [OPTIONS] 1- tmi : Magnetization direction used and data are projected onto the inducing field direction F.shape([ndata, nc]) - + 2- xyz : Magnetization direction used and data are given in 3-components F.shape([3*ndata, nc]) - + 3- full: Full tensor matrix stored with shape([3*ndata, 3*nc]) - + OUTPUT: F = Linear forward modeling operation Created on Dec, 20th 2015 @author: dominiquef +<<<<<<< HEAD """ # Find non-zero cells inds = np.nonzero(actv)[0] @@ -560,9 +567,21 @@ def Intrgl_Fwr_Op(mesh,B,M,rxLoc,actv,flag): # Convert Bdecination from north to cartesian +======= + """ + + xn = mesh.vectorNx; + yn = mesh.vectorNy; + zn = mesh.vectorNz; + + ndata = rxLoc.shape[0] + + + # Convert Bdecination from north to cartesian +>>>>>>> dbcd57bc4e5d0d6258a10b760966b58419e0a9d2 D = (450.-float(B[1]))%360. - - + + # Pre-allocate space and create magnetization matrix if required if (flag=='tmi') | (flag == 'xyz'): # If assumes uniform magnetization direction @@ -570,16 +589,16 @@ def Intrgl_Fwr_Op(mesh,B,M,rxLoc,actv,flag): print 'Magnetization vector must be Nc x 3' return - - + + Mx = Utils.sdiag(M[:,0]*B[2]) My = Utils.sdiag(M[:,1]*B[2]) Mz = Utils.sdiag(M[:,2]*B[2]) - - Mxyz = sp.vstack((Mx,My,Mz)) - - + Mxyz = sp.vstack((Mx,My,Mz)) + + + if flag == 'tmi': F = np.zeros((ndata, mcell)) @@ -587,56 +606,77 @@ def Intrgl_Fwr_Op(mesh,B,M,rxLoc,actv,flag): 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; - + elif flag == 'xyz': +<<<<<<< HEAD F = np.zeros((int(3*ndata), mcell)) elif flag == 'full': F = np.zeros((int(3*ndata), int(3*mcell))) +======= + F = np.zeros((int(3*ndata), mesh.nC)) + + elif flag == 'full': + F = np.zeros((int(3*ndata), int(3*mesh.nC))) + +>>>>>>> dbcd57bc4e5d0d6258a10b760966b58419e0a9d2 else: print """Flag must be either 'tmi' | 'xyz' | 'full', please revised""" return - - + + # Loop through all observations and create forward operator (ndata-by-mcell) print "Begin calculation of forward operator: " + flag - + # Add counter to dsiplay progress. Good for large problems count = -1; for ii in range(ndata): +<<<<<<< HEAD tx, ty, tz = get_T_mat(Xn,Yn,Zn,rxLoc[ii,:]) +======= + + tx, ty, tz = get_T_mat(xn,yn,zn,rxLoc[ii,:]) + +>>>>>>> dbcd57bc4e5d0d6258a10b760966b58419e0a9d2 if flag=='tmi': F[ii,:] = Ptmi.dot(np.vstack((tx,ty,tz)))*Mxyz - + elif flag == 'xyz': F[ii,:] = tx*Mxyz F[ii+ndata,:] = ty*Mxyz F[ii+2*ndata,:] = tz*Mxyz - - elif flag == 'full': + + elif flag == 'full': F[ii,:] = tx F[ii+ndata,:] = ty F[ii+2*ndata,:] = tz - - + + # Display progress count = progress(ii,count,ndata) +<<<<<<< HEAD print "Done 100% ...forward operator completed!!\n" +======= + + + print "Done 100% ...forward modeling completed!!\n" + +>>>>>>> dbcd57bc4e5d0d6258a10b760966b58419e0a9d2 return F def get_T_mat(Xn,Yn,Zn,rxLoc): - """ - Load in the active nodes of a tensor mesh and computes the magnetic tensor + """ + Load in the active nodes of a tensor mesh and computes the magnetic tensor for a given observation location rxLoc[obsx, obsy, obsz] - + INPUT: - Xn, Yn, Zn: Node location matrix for the lower and upper most corners of + Xn, Yn, Zn: Node location matrix for the lower and upper most corners of all cells in the mesh shape[nC,2] M OUTPUT: @@ -648,20 +688,27 @@ def get_T_mat(Xn,Yn,Zn,rxLoc): Only the upper half 5 elements have to be computed since symetric. Currently done as for-loops but will eventually be changed to vector indexing, once the topography has been figured out. - + Created on Oct, 20th 2015 @author: dominiquef +<<<<<<< HEAD """ eps = 1e-10 # add a small value to the locations to avoid /0 +======= + """ + + +>>>>>>> dbcd57bc4e5d0d6258a10b760966b58419e0a9d2 mcell = Xn.shape[0] - + # Pre-allocate space for 1D array Tx = np.zeros((1,3*mcell)) Ty = np.zeros((1,3*mcell)) Tz = np.zeros((1,3*mcell)) +<<<<<<< HEAD dz2 = rxLoc[2] - Zn[:,0] + eps dz1 = rxLoc[2] - Zn[:,1] + eps @@ -672,12 +719,24 @@ def get_T_mat(Xn,Yn,Zn,rxLoc): dx2 = Xn[:,1] - rxLoc[0] + eps dx1 = Xn[:,0] - rxLoc[0] + eps +======= + + dz2 = rxLoc[2] - Zn[:,0] + dz1 = rxLoc[2] - Zn[:,1] + + dy2 = Yn[:,1] - rxLoc[1] + dy1 = Yn[:,0] - rxLoc[1] + + dx2 = Xn[:,1] - rxLoc[0] + dx1 = Xn[:,0] - rxLoc[0] + +>>>>>>> dbcd57bc4e5d0d6258a10b760966b58419e0a9d2 R1 = ( dy2**2 + dx2**2 ) R2 = ( dy2**2 + dx1**2 ) R3 = ( dy1**2 + dx2**2 ) R4 = ( dy1**2 + dx1**2 ) - - + + arg1 = np.sqrt( dz2**2 + R2 ) arg2 = np.sqrt( dz2**2 + R1 ) arg3 = np.sqrt( dz1**2 + R1 ) @@ -686,9 +745,9 @@ def get_T_mat(Xn,Yn,Zn,rxLoc): arg6 = np.sqrt( dz2**2 + R4 ) arg7 = np.sqrt( dz1**2 + R4 ) arg8 = np.sqrt( dz1**2 + R3 ) - - - + + + Tx[0,0:mcell] = np.arctan2( dy1 * dz2 , ( dx2 * arg5 ) ) +\ - np.arctan2( dy2 * dz2 , ( dx2 * arg2 ) ) +\ np.arctan2( dy2 * dz1 , ( dx2 * arg3 ) ) +\ @@ -697,13 +756,13 @@ def get_T_mat(Xn,Yn,Zn,rxLoc): - np.arctan2( dy1 * dz2 , ( dx1 * arg6 ) ) +\ np.arctan2( dy1 * dz1 , ( dx1 * arg7 ) ) +\ - np.arctan2( dy2 * dz1 , ( dx1 * arg4 ) ); - - + + Ty[0,0:mcell] = np.log( ( dz2 + arg2 ) / (dz1 + arg3 ) ) +\ -np.log( ( dz2 + arg1 ) / (dz1 + arg4 ) ) +\ np.log( ( dz2 + arg6 ) / (dz1 + arg7 ) ) +\ -np.log( ( dz2 + arg5 ) / (dz1 + arg8 ) ); - + Ty[0,mcell:2*mcell] = np.arctan2( dx1 * dz2 , ( dy2 * arg1 ) ) +\ - np.arctan2( dx2 * dz2 , ( dy2 * arg2 ) ) +\ np.arctan2( dx2 * dz1 , ( dy2 * arg3 ) ) +\ @@ -712,57 +771,57 @@ def get_T_mat(Xn,Yn,Zn,rxLoc): - np.arctan2( dx1 * dz2 , ( dy1 * arg6 ) ) +\ np.arctan2( dx1 * dz1 , ( dy1 * arg7 ) ) +\ - np.arctan2( dx2 * dz1 , ( dy1 * arg8 ) ); - + R1 = (dy2**2 + dz1**2); R2 = (dy2**2 + dz2**2); R3 = (dy1**2 + dz1**2); R4 = (dy1**2 + dz2**2); - + Ty[0,2*mcell:] = np.log( ( dx1 + np.sqrt( dx1**2 + R1 ) ) / (dx2 + np.sqrt( dx2**2 + R1 ) ) ) +\ -np.log( ( dx1 + np.sqrt( dx1**2 + R2 ) ) / (dx2 + np.sqrt( dx2**2 + R2 ) ) ) +\ np.log( ( dx1 + np.sqrt( dx1**2 + R4 ) ) / (dx2 + np.sqrt( dx2**2 + R4 ) ) ) +\ -np.log( ( dx1 + np.sqrt( dx1**2 + R3 ) ) / (dx2 + np.sqrt( dx2**2 + R3 ) ) ); - + R1 = (dx2**2 + dz1**2); R2 = (dx2**2 + dz2**2); R3 = (dx1**2 + dz1**2); R4 = (dx1**2 + dz2**2); - + Tx[0,2*mcell:] = np.log( ( dy1 + np.sqrt( dy1**2 + R1 ) ) / (dy2 + np.sqrt( dy2**2 + R1 ) ) ) +\ -np.log( ( dy1 + np.sqrt( dy1**2 + R2 ) ) / (dy2 + np.sqrt( dy2**2 + R2 ) ) ) +\ np.log( ( dy1 + np.sqrt( dy1**2 + R4 ) ) / (dy2 + np.sqrt( dy2**2 + R4 ) ) ) +\ -np.log( ( dy1 + np.sqrt( dy1**2 + R3 ) ) / (dy2 + np.sqrt( dy2**2 + R3 ) ) ); - + Tz[0,2*mcell:] = -( Ty[0,mcell:2*mcell] + Tx[0,0:mcell] ); Tz[0,mcell:2*mcell] = Ty[0,2*mcell:]; Tx[0,mcell:2*mcell] = Ty[0,0:mcell]; Tz[0,0:mcell] = Tx[0,2*mcell:]; - - - + + + Tx = Tx/(4*np.pi); Ty = Ty/(4*np.pi); - Tz = Tz/(4*np.pi); - - + Tz = Tz/(4*np.pi); + + return Tx,Ty,Tz def progress(iter,prog,final): - """ + """ progress(iter,prog,final) - + Function measuring the progress of a process and print to screen the %. Useful to estimate the remaining runtime of a large problem. - + Created on Dec, 20th 2015 @author: dominiquef """ arg = np.floor(float(iter)/float(final)*10.); - + if arg > prog: - - strg = "Done " + str(arg*10) + " %" + + strg = "Done " + str(arg*10) + " %" print strg prog = arg; @@ -771,73 +830,73 @@ def progress(iter,prog,final): def dipazm_2_xyz(dip,azm_N): """ dipazm_2_xyz(dip,azm_N) - + Function converting degree angles for dip and azimuth from north to a 3-components in cartesian coordinates. - + INPUT dip : Value or vector of dip from horizontal in DEGREE azm_N : Value or vector of azimuth from north in DEGREE - + OUTPUT M : [n-by-3] Array of xyz components of a unit vector in cartesian - + Created on Dec, 20th 2015 @author: dominiquef """ mcell = len(azm_N) - + M = np.zeros((mcell,3)) - + # Modify azimuth from North to Cartesian-X azm_X = (450.- azm_N) % 360. - + D = np.deg2rad(dip) I = np.deg2rad(azm_X) - + M[:,0] = np.cos(D) * np.cos(I) ; M[:,1] = np.cos(D) * np.sin(I) ; M[:,2] = np.sin(D) ; - + return M def get_dist_wgt(mesh,rxLoc,R,R0): """ get_dist_wgt(xn,yn,zn,rxLoc,R,R0) - + Function creating a distance weighting function required for the magnetic inverse problem. - + INPUT xn, yn, zn : Node location rxLoc : Observation locations [obsx, obsy, obsz] R : Decay factor (mag=3, grav =2) R0 : Small factor added (default=dx/4) - + OUTPUT wr : [mcell] Vector of distance weighting - + Created on Dec, 20th 2015 @author: dominiquef """ - + p = 1/np.sqrt(3); - + # Create cell center location Ym,Xm,Zm = np.meshgrid(mesh.vectorCCy, mesh.vectorCCx, mesh.vectorCCz) - hY,hX,hZ = np.meshgrid(mesh.hy, mesh.hx, mesh.hz) - + hY,hX,hZ = np.meshgrid(mesh.hy, mesh.hx, mesh.hz) + V = np.reshape(mesh.vol,hY.shape) wr = np.zeros(hY.shape) - + ndata = rxLoc.shape[0] count = -1; print "Begin calculation of distance weighting for R= " + str(R) - + for dd in range(ndata): - + nx1 = (Xm - hX * p - rxLoc[dd,0])**2 nx2 = (Xm + hX * p - rxLoc[dd,0])**2 @@ -855,26 +914,26 @@ def get_dist_wgt(mesh,rxLoc,R,R0): R6 = np.sqrt(nx1 + ny2 + nz2) R7 = np.sqrt(nx2 + ny2 + nz1) R8 = np.sqrt(nx2 + ny2 + nz2) - + temp = (R1 + R0)**-R + (R2 + R0)**-R + (R3 + R0)**-R + (R4 + R0)**-R + (R5 + R0)**-R + (R6 + R0)**-R + (R7 + R0)**-R + (R8 + R0)**-R wr = wr + (V*temp/8.)**2 - + count = progress(dd,count,ndata) - - + + wr = np.sqrt(wr)/V wr = mkvc(wr) wr = np.sqrt(wr/(np.max(wr))) - + return wr - + def writeUBCobs(filename,B,M,rxLoc,d,wd): """ writeUBCobs(filename,B,M,rxLoc,d,wd) - + Function writing an observation file in UBC-MAG3D format. - + INPUT filename : Name of out file including directory B : Inducing field parameters [Inc, Decl, Intensity] @@ -882,43 +941,43 @@ def writeUBCobs(filename,B,M,rxLoc,d,wd): rxLoc : Observation locations [obsx, obsy, obsz] d : Data vector wd : Uncertainty vector - + OUTPUT - Obsfile - + Obsfile + Created on Dec, 27th 2015 @author: dominiquef """ - + data = np.c_[rxLoc , d , wd] with file(filename,'w') as fid: fid.write('%6.2f %6.2f %6.2f\n' %(B[0], B[1], B[2]) ) - fid.write('%6.2f %6.2f %6.2f\n' %(M[0], M[1], 1) ) - fid.write('%i\n' %len(d) ) + fid.write('%6.2f %6.2f %6.2f\n' %(M[0], M[1], 1) ) + fid.write('%i\n' %len(d) ) np.savetxt(fid, data, fmt='%e',delimiter=' ',newline='\n') - - + + print "Observation file saved to: " + filename - + def getActiveTopo(mesh,topo,flag): """ getActiveTopo(mesh,topo) - + Function creates an active cell model from topography - + INPUT mesh : Mesh in SimPEG format topo : Scatter points defining topography [x,y,z] - + OUTPUT - actv : Active cell model - + actv : Active cell model + Created on Dec, 27th 2015 @author: dominiquef - """ + """ import scipy.interpolate as interpolation if (flag=='N'): @@ -926,42 +985,28 @@ def getActiveTopo(mesh,topo,flag): # wght = np.zeros((mesh.nNx,mesh.nNy)) cx = mesh.vectorNx cy = mesh.vectorNy - + F = interpolation.NearestNDInterpolator(topo[:,0:2],topo[:,2]) [Y,X] = np.meshgrid(cy,cx) - + Zn = F(X,Y) -#============================================================================== -# for ii in range(topo.shape[0]): -# dx = topo[ii,0] - cx + 1e-8 -# dy = topo[ii,1] - cy + 1e-8 -# -# [Wx,Wy] = np.meshgrid(dx**2.,dy**2.) -# -# W = np.sqrt(Wx + Wy)**-1. -# -# wght = wght + W -# Zn = Zn + topo[ii,2]*W -# -# Zn = Zn/wght -#============================================================================== - + actv = np.zeros((mesh.nCx, mesh.nCy, mesh.nCz)) - + if (flag=='N'): Nz = mesh.vectorNz[1:] - - - for jj in range(mesh.nCy): - + + + for jj in range(mesh.nCy): + for ii in range(mesh.nCx): - + temp = [kk for kk in range(len(Nz)) if np.all(Zn[ii:(ii+2),jj:(jj+2)] > Nz[kk]) ] actv[ii,jj,temp] = 1 - - + + actv = mkvc(actv) - + return actv - - \ No newline at end of file + +