mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-30 20:36:55 +08:00
Move sub functions to PF.Magnetics
Add Mag Integral forward operator Begin Inversion script for magnetic integral
This commit is contained in:
D Fournier
+54
-1
@@ -187,4 +187,57 @@ def readUBCmagObs(obs_file):
|
||||
dobs[ii,:] = np.array(line.split(),dtype=float)
|
||||
line = fid.readline()
|
||||
|
||||
return B, M, dobs
|
||||
return B, M, dobs
|
||||
|
||||
def read_MAGfwr_inp(input_file):
|
||||
"""Read input files for forward modeling MAG data with integral form
|
||||
INPUT:
|
||||
input_file: File name containing the forward parameter
|
||||
|
||||
OUTPUT:
|
||||
mshfile
|
||||
obsfile
|
||||
modfile
|
||||
magfile
|
||||
topofile
|
||||
# All files should be in the working directory, otherwise the path must
|
||||
# be specified.
|
||||
|
||||
Created on Jul 17, 2013
|
||||
|
||||
@author: dominiquef
|
||||
"""
|
||||
|
||||
|
||||
fid = open(input_file,'r')
|
||||
|
||||
line = fid.readline()
|
||||
l_input = line.split('!')
|
||||
mshfile = l_input[0].rstrip()
|
||||
|
||||
line = fid.readline()
|
||||
l_input = line.split('!')
|
||||
obsfile = l_input[0].rstrip()
|
||||
|
||||
line = fid.readline()
|
||||
l_input = line.split('!')
|
||||
modfile = l_input[0].rstrip()
|
||||
|
||||
line = fid.readline()
|
||||
l_input = line.split('!')
|
||||
if l_input=='null':
|
||||
magfile = []
|
||||
|
||||
else:
|
||||
magfile = l_input[0].rstrip()
|
||||
|
||||
|
||||
line = fid.readline()
|
||||
l_input = line.split('!')
|
||||
if l_input=='null':
|
||||
topofile = []
|
||||
|
||||
else:
|
||||
topofile = l_input[0].rstrip()
|
||||
|
||||
return mshfile, obsfile, modfile, magfile, topofile
|
||||
@@ -57,7 +57,7 @@ def fwr_MAG_data(mesh,B,M,rxLoc,model,flag):
|
||||
d = np.zeros(ndata)
|
||||
|
||||
elif flag=='xyz':
|
||||
d = np.zeros((3,ndata))
|
||||
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..."
|
||||
@@ -70,7 +70,7 @@ def fwr_MAG_data(mesh,B,M,rxLoc,model,flag):
|
||||
Gxyz = np.vstack((tx,ty,tz))*Mxyz
|
||||
|
||||
if flag=='xyz':
|
||||
d[:,ii] = Gxyz.dot(model)
|
||||
d[ii:ndata:] = mkvc(Gxyz.dot(model))
|
||||
|
||||
elif flag=='tmi':
|
||||
d[ii] = Ptmi.dot(Gxyz.dot(model))
|
||||
@@ -79,13 +79,13 @@ def fwr_MAG_F(mesh,B,M,rxLoc,flag):
|
||||
Mxyz = sp.spdiags(M,0,mesh.nC * 3,mesh.nC * 3)
|
||||
|
||||
if flag == 'tmi':
|
||||
F = np.zeros(ndata, mesh.nC)
|
||||
F = np.zeros((ndata, mesh.nC))
|
||||
|
||||
elif flag == 'xyz':
|
||||
F = np.zeros(int(3*ndata), mesh.nC)
|
||||
F = np.zeros((int(3*ndata), mesh.nC))
|
||||
|
||||
elif flag == 'full':
|
||||
F = np.zeros(int(3*ndata), int(3*mesh.nC))
|
||||
F = np.zeros((int(3*ndata), int(3*mesh.nC)))
|
||||
|
||||
else:
|
||||
print """Flag must be either 'tmi' | 'xyz' | 'full', please revised"""
|
||||
@@ -8,15 +8,15 @@ os.chdir(home_dir)
|
||||
|
||||
#%%
|
||||
from SimPEG import np, Utils
|
||||
from simpegPF import BaseMag
|
||||
import simpegPF as PF
|
||||
|
||||
## New scripts to be added to basecode
|
||||
from fwr_MAG_data import fwr_MAG_data
|
||||
from read_MAGfwr_inp import read_MAGfwr_inp
|
||||
#from fwr_MAG_data import fwr_MAG_data
|
||||
#from read_MAGfwr_inp import read_MAGfwr_inp
|
||||
|
||||
#%%
|
||||
# Read input file
|
||||
[mshfile, obsfile, modfile, magfile, topofile] = read_MAGfwr_inp(inpfile)
|
||||
[mshfile, obsfile, modfile, magfile, topofile] = PF.BaseMag.read_MAGfwr_inp(inpfile)
|
||||
|
||||
# Load mesh file
|
||||
mesh = Utils.meshutils.readUBCTensorMesh(mshfile)
|
||||
@@ -25,13 +25,13 @@ mesh = Utils.meshutils.readUBCTensorMesh(mshfile)
|
||||
model = Utils.meshutils.readUBCTensorModel(modfile,mesh)
|
||||
|
||||
# Load in observation file
|
||||
[B,M,dobs] = BaseMag.readUBCmagObs(obsfile)
|
||||
[B,M,dobs] = PF.BaseMag.readUBCmagObs(obsfile)
|
||||
|
||||
rxLoc = dobs[:,0:3]
|
||||
ndata = rxLoc.shape[0]
|
||||
|
||||
# Compute forward model using integral equation
|
||||
d = fwr_MAG_data(mesh,B,M,rxLoc,model,'tmi')
|
||||
d = PF.Magnetics.Intgrl_Fwr_Data(mesh,B,M,rxLoc,model,'tmi')
|
||||
|
||||
# Form data object with coordinates and write to file
|
||||
data = np.c_[rxLoc , d , np.zeros((ndata,1))]
|
||||
@@ -43,6 +43,7 @@ with file('FWR_data.dat','w') as fid:
|
||||
fid.write('%i\n' %(ndata) )
|
||||
np.savetxt(fid, data, fmt='%e',delimiter=' ',newline='\n')
|
||||
|
||||
|
||||
print "Observation file saved to " + home_dir + '\FWR_data.dat'
|
||||
|
||||
|
||||
@@ -8,10 +8,8 @@ os.chdir(home_dir)
|
||||
from SimPEG import *
|
||||
import matplotlib.pyplot as plt
|
||||
import simpegPF as PF
|
||||
from simpegPF import BaseMag
|
||||
import matplotlib
|
||||
|
||||
from fwr_MAG_obs import fwr_MAG_obs
|
||||
#from fwr_MAG_data import fwr_MAG_data
|
||||
|
||||
plt.close('all')
|
||||
|
||||
@@ -39,7 +37,7 @@ 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)]
|
||||
@@ -54,6 +52,8 @@ for ii in range(d_iter):
|
||||
|
||||
mcell = mesh.nC
|
||||
|
||||
print 'Mesh size: ' + str(mcell)
|
||||
|
||||
sph_ind = PF.MagAnalytics.spheremodel(mesh, 0, 0, 0, R)
|
||||
|
||||
chibkg = 0.
|
||||
@@ -62,10 +62,10 @@ for ii in range(d_iter):
|
||||
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,:])
|
||||
d = PF.Magnetics.Intgrl_Fwr_Data(mesh,B,M,rxLoc,model,'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')
|
||||
@@ -90,30 +90,32 @@ for ii in range(d_iter):
|
||||
|
||||
|
||||
#%% Plot results
|
||||
print 'Residual between analytical sphere and integral forward'
|
||||
for ii in range(d_iter):
|
||||
nc = 3**(ii+1)
|
||||
print "Residual= " + str(lrl[ii]) + "\t dx= " + str(1./nc)
|
||||
|
||||
print "||r||= " + str(lrl[ii]) + "\t dx= " + str(1./nc)
|
||||
|
||||
#%% Plot fields
|
||||
|
||||
plt.figure(1)
|
||||
ax = plt.subplot(131)
|
||||
ax = plt.subplot(221)
|
||||
plt.imshow(np.reshape(bxa,X.shape).T, interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()], origin = 'lower')
|
||||
plt.colorbar()
|
||||
plt.colorbar(fraction=0.04)
|
||||
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)
|
||||
ax = plt.subplot(222)
|
||||
plt.imshow(np.reshape(bya,X.shape).T, interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()], origin = 'lower')
|
||||
plt.colorbar()
|
||||
plt.colorbar(fraction=0.04)
|
||||
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)
|
||||
ax = plt.subplot(212)
|
||||
plt.imshow(np.reshape(bza,X.shape).T, interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()], origin = 'lower')
|
||||
plt.colorbar()
|
||||
plt.colorbar(fraction=0.04)
|
||||
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')
|
||||
@@ -121,23 +123,23 @@ ax.set_title('Sphere Ana Bz')
|
||||
#%% Plot the forward solution from integral
|
||||
|
||||
plt.figure(2)
|
||||
ax = plt.subplot(131)
|
||||
ax = plt.subplot(221)
|
||||
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.colorbar(fraction=0.04)
|
||||
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')
|
||||
|
||||
ax = plt.subplot(132)
|
||||
ax = plt.subplot(222)
|
||||
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.colorbar(fraction=0.04)
|
||||
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')
|
||||
|
||||
ax = plt.subplot(133)
|
||||
ax = plt.subplot(212)
|
||||
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.colorbar(fraction=0.04)
|
||||
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')
|
||||
@@ -145,23 +147,23 @@ ax.set_title('Sphere Ana Bz')
|
||||
|
||||
#%% Plot foward data
|
||||
plt.figure(3)
|
||||
ax = plt.subplot(131)
|
||||
ax = plt.subplot(221)
|
||||
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.colorbar(fraction=0.04)
|
||||
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)
|
||||
ax = plt.subplot(222)
|
||||
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.colorbar(fraction=0.04)
|
||||
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)
|
||||
ax = plt.subplot(212)
|
||||
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.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')
|
||||
@@ -0,0 +1,37 @@
|
||||
import os
|
||||
|
||||
home_dir = 'C:\Users\dominiquef.MIRAGEOSCIENCE\Documents\GIT\SimPEG\simpegpf\simpegPF\Dev'
|
||||
|
||||
inpfile = 'MAG3Cinv.inp'
|
||||
|
||||
os.chdir(home_dir)
|
||||
|
||||
#%%
|
||||
from SimPEG import np, Utils
|
||||
import simpegPF as PF
|
||||
|
||||
## New scripts to be added to basecode
|
||||
#from fwr_MAG_data import fwr_MAG_data
|
||||
#from read_MAGfwr_inp import read_MAGfwr_inp
|
||||
|
||||
#%%
|
||||
# Read input file
|
||||
[mshfile, obsfile, modfile, magfile, topofile] = PF.BaseMag.read_MAGfwr_inp(inpfile)
|
||||
|
||||
# Load mesh file
|
||||
mesh = Utils.meshutils.readUBCTensorMesh(mshfile)
|
||||
|
||||
# Load model file
|
||||
model = Utils.meshutils.readUBCTensorModel(modfile,mesh)
|
||||
|
||||
# Load in observation file
|
||||
[B,M,dobs] = PF.BaseMag.readUBCmagObs(obsfile)
|
||||
|
||||
rxLoc = dobs[:,0:3]
|
||||
ndata = rxLoc.shape[0]
|
||||
|
||||
# Create forward operator
|
||||
F = PF.Magnetics.Intrgl_Fwr_Op(mesh,B,M,rxLoc,'tmi')
|
||||
|
||||
|
||||
|
||||
+341
-4
@@ -409,10 +409,347 @@ if __name__ == '__main__':
|
||||
# plt.show()
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
def Intgrl_Fwr_Data(mesh,B,M,rxLoc,model,flag):
|
||||
"""
|
||||
Forward model magnetic data using integral equation
|
||||
|
||||
INPUT:
|
||||
xn, yn, zn = Mesh nodes location
|
||||
B = Inducing field parameter [Binc, Bdecl, B0]
|
||||
M = Magnetization matrix [Minc, Mdecl]
|
||||
rxLox = Observation location informat [obsx, obsy, obsz]
|
||||
model = Model associated with mesh
|
||||
flag = Data type "tmi" | "xyz"
|
||||
|
||||
OUTPUT:
|
||||
dobs =Observation array in format [obsx, obsy, obsz, data]
|
||||
|
||||
Created on Oct 7, 2015
|
||||
|
||||
@author: dominiquef
|
||||
"""
|
||||
|
||||
xn = mesh.vectorNx;
|
||||
yn = mesh.vectorNy;
|
||||
zn = mesh.vectorNz;
|
||||
|
||||
mcell = (len(xn)-1) * (len(yn)-1) * (len(zn)-1)
|
||||
|
||||
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
|
||||
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':
|
||||
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,:])
|
||||
Gxyz = np.vstack((tx,ty,tz))*Mxyz
|
||||
|
||||
if flag=='xyz':
|
||||
d[ii::ndata] = mkvc(Gxyz.dot(model))
|
||||
|
||||
elif flag=='tmi':
|
||||
d[ii] = Ptmi.dot(Gxyz.dot(model))
|
||||
|
||||
# Display progress
|
||||
count = progress(ii,count,ndata)
|
||||
|
||||
|
||||
print "Done 100% ...forward modeling completed!!\n"
|
||||
|
||||
return d
|
||||
|
||||
def Intrgl_Fwr_Op(mesh,B,M,rxLoc,flag):
|
||||
"""
|
||||
Magnetic forward operator in integral form
|
||||
|
||||
INPUT:
|
||||
mesh = Mesh in SimPEG format
|
||||
B = Inducing field parameter [Binc, Bdecl, B0]
|
||||
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
|
||||
"""
|
||||
|
||||
xn = mesh.vectorNx;
|
||||
yn = mesh.vectorNy;
|
||||
zn = mesh.vectorNz;
|
||||
|
||||
mcell = (len(xn)-1) * (len(yn)-1) * (len(zn)-1)
|
||||
|
||||
ndata = rxLoc.shape[0]
|
||||
|
||||
|
||||
# Convert Bdecination from north to cartesian
|
||||
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
|
||||
if len(M) == 3:
|
||||
|
||||
# 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])
|
||||
|
||||
Mxyz = sp.vstack((Mx,My,Mz));
|
||||
|
||||
# Otherwise if given a vector 3*ncells
|
||||
elif len(M) == mesh.nC * 3:
|
||||
|
||||
Mxyz = sp.spdiags(M,0,mesh.nC * 3,mesh.nC * 3)
|
||||
|
||||
if flag == 'tmi':
|
||||
F = np.zeros((ndata, mesh.nC))
|
||||
|
||||
# Projection matrix
|
||||
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':
|
||||
F = np.zeros((int(3*ndata), mesh.nC))
|
||||
|
||||
elif flag == 'full':
|
||||
F = np.zeros((int(3*ndata), int(3*mesh.nC)))
|
||||
|
||||
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):
|
||||
|
||||
tx, ty, tz = get_T_mat(xn,yn,zn,rxLoc[ii,:])
|
||||
|
||||
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':
|
||||
F[ii,:] = tx
|
||||
F[ii+ndata,:] = ty
|
||||
F[ii+2*ndata,:] = tz
|
||||
|
||||
|
||||
# Display progress
|
||||
count = progress(ii,count,ndata)
|
||||
|
||||
|
||||
print "Done 100% ...forward modeling completed!!\n"
|
||||
|
||||
return F
|
||||
|
||||
def get_T_mat(xn,yn,zn,rxLoc):
|
||||
"""
|
||||
Load in the nodes of a tensor mesh and computes the magnetic tensor
|
||||
for a given observation location rxLoc[obsx, obsy, obsz]
|
||||
OUTPUT:
|
||||
Tx = [Txx Txy Txz]
|
||||
Ty = [Tyx Tyy Tyz]
|
||||
Tz = [Tzx Tzy Tzz]
|
||||
|
||||
where each elements have dimension 1-by-mcell.
|
||||
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.
|
||||
"""
|
||||
|
||||
ncx = len(xn)-1
|
||||
ncy = len(yn)-1
|
||||
ncz = len(zn)-1
|
||||
|
||||
mcell = ncx*ncy*ncz
|
||||
|
||||
# Pre-allocate space for 1D array
|
||||
Tx = np.zeros((1,3*mcell))
|
||||
Ty = np.zeros((1,3*mcell))
|
||||
Tz = np.zeros((1,3*mcell))
|
||||
|
||||
yn2,xn2,zn2 = np.meshgrid(yn[1:], xn[1:], zn[1:])
|
||||
yn1,xn1,zn1 = np.meshgrid(yn[0:ncy], xn[0:ncx], zn[0:ncz])
|
||||
|
||||
yn2 = mkvc(yn2)
|
||||
yn1 = mkvc(yn1)
|
||||
|
||||
zn2 = mkvc(zn2)
|
||||
zn1 = mkvc(zn1)
|
||||
|
||||
xn2 = mkvc(xn2)
|
||||
xn1 = mkvc(xn1)
|
||||
#%%
|
||||
#==============================================================================
|
||||
|
||||
|
||||
dz2 = rxLoc[2] - zn1;
|
||||
dz1 = rxLoc[2] - zn2;
|
||||
|
||||
|
||||
dy2 = yn2 - rxLoc[1];
|
||||
dy1 = yn1 - rxLoc[1];
|
||||
|
||||
|
||||
dx2 = xn2 - rxLoc[0];
|
||||
dx1 = xn1 - rxLoc[0];
|
||||
|
||||
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 );
|
||||
arg4 = np.sqrt( dz1**2 + R2 );
|
||||
arg5 = np.sqrt( dz2**2 + R3 );
|
||||
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 ) ) +\
|
||||
- np.arctan2( dy1 * dz1 , ( dx2 * arg8 ) ) +\
|
||||
np.arctan2( dy2 * dz2 , ( dx1 * arg1 ) ) +\
|
||||
- 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 ) ) +\
|
||||
- np.arctan2( dx1 * dz1 , ( dy2 * arg4 ) ) +\
|
||||
np.arctan2( dx2 * dz2 , ( dy1 * arg5 ) ) +\
|
||||
- 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);
|
||||
|
||||
|
||||
return Tx,Ty,Tz
|
||||
|
||||
def progress(iter,prog,final):
|
||||
|
||||
arg = np.floor(float(iter)/float(final)*10.);
|
||||
|
||||
if arg > prog:
|
||||
|
||||
strg = "Done " + str(arg*10) + " %"
|
||||
print strg
|
||||
prog = arg;
|
||||
|
||||
return prog
|
||||
|
||||
|
||||
|
||||
|
||||
Reference in New Issue
Block a user