wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-06-28 13:20:13 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'Dom' of https://github.com/simpeg/simpegpf into Dom

Browse Source 
Conflicts:
	simpegPF/Magnetics.py
This commit is contained in:
D Fournier
2016-01-15 00:03:37 -08:00
parent c0445e5db5 dbcd57bc4e
commit e4e19ceac1
2 changed files with 239 additions and 192 deletions
Download Patch File Download Diff File Expand all files Collapse all files
simpegPF/Dev/Intgrl_MAG_Fwr_SphereTest.py
+33 -31
View File
@@ -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')
ax.set_title('Sphere Ana Bz')
plt.show()
simpegPF/Magnetics.py
+206 -161
View File
@@ -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