diff --git a/simpegPF/BaseMag.py b/simpegPF/BaseMag.py index bd78009d..df75013b 100644 --- a/simpegPF/BaseMag.py +++ b/simpegPF/BaseMag.py @@ -187,4 +187,57 @@ def readUBCmagObs(obs_file): dobs[ii,:] = np.array(line.split(),dtype=float) line = fid.readline() - return B, M, dobs \ No newline at end of file + 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 \ No newline at end of file diff --git a/simpegPF/Dev/functions/fwr_MAG_data.py b/simpegPF/Dev/Archives/fwr_MAG_data.py similarity index 96% rename from simpegPF/Dev/functions/fwr_MAG_data.py rename to simpegPF/Dev/Archives/fwr_MAG_data.py index 88ba50c2..6bd3dbe1 100644 --- a/simpegPF/Dev/functions/fwr_MAG_data.py +++ b/simpegPF/Dev/Archives/fwr_MAG_data.py @@ -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)) diff --git a/simpegPF/Dev/functions/get_MAG_F.py b/simpegPF/Dev/Archives/get_MAG_F.py similarity index 96% rename from simpegPF/Dev/functions/get_MAG_F.py rename to simpegPF/Dev/Archives/get_MAG_F.py index d9431b16..54172d71 100644 --- a/simpegPF/Dev/functions/get_MAG_F.py +++ b/simpegPF/Dev/Archives/get_MAG_F.py @@ -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""" diff --git a/simpegPF/Dev/functions/get_T_mat.py b/simpegPF/Dev/Archives/get_T_mat.py similarity index 100% rename from simpegPF/Dev/functions/get_T_mat.py rename to simpegPF/Dev/Archives/get_T_mat.py diff --git a/simpegPF/Dev/functions/get_UBC_mesh.py b/simpegPF/Dev/Archives/get_UBC_mesh.py similarity index 100% rename from simpegPF/Dev/functions/get_UBC_mesh.py rename to simpegPF/Dev/Archives/get_UBC_mesh.py diff --git a/simpegPF/Dev/functions/read_MAG_obs.py b/simpegPF/Dev/Archives/read_MAG_obs.py similarity index 100% rename from simpegPF/Dev/functions/read_MAG_obs.py rename to simpegPF/Dev/Archives/read_MAG_obs.py diff --git a/simpegPF/Dev/functions/read_MAGfwr_inp.py b/simpegPF/Dev/Archives/read_MAGfwr_inp.py similarity index 100% rename from simpegPF/Dev/functions/read_MAGfwr_inp.py rename to simpegPF/Dev/Archives/read_MAGfwr_inp.py diff --git a/simpegPF/Dev/MAG3D_fwr_Driver.py b/simpegPF/Dev/Intgrl_MAG_Fwr_Driver.py similarity index 76% rename from simpegPF/Dev/MAG3D_fwr_Driver.py rename to simpegPF/Dev/Intgrl_MAG_Fwr_Driver.py index 4870159e..448e6b83 100644 --- a/simpegPF/Dev/MAG3D_fwr_Driver.py +++ b/simpegPF/Dev/Intgrl_MAG_Fwr_Driver.py @@ -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' diff --git a/simpegPF/Dev/MAG3D_fwr_Test.py b/simpegPF/Dev/Intgrl_MAG_Fwr_SphereTest.py similarity index 85% rename from simpegPF/Dev/MAG3D_fwr_Test.py rename to simpegPF/Dev/Intgrl_MAG_Fwr_SphereTest.py index 1d30dbb6..ddea3760 100644 --- a/simpegPF/Dev/MAG3D_fwr_Test.py +++ b/simpegPF/Dev/Intgrl_MAG_Fwr_SphereTest.py @@ -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') \ No newline at end of file diff --git a/simpegPF/Dev/Intgrl_MAG_Inv_Driver.py b/simpegPF/Dev/Intgrl_MAG_Inv_Driver.py new file mode 100644 index 00000000..296d0a1b --- /dev/null +++ b/simpegPF/Dev/Intgrl_MAG_Inv_Driver.py @@ -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') + + + diff --git a/simpegPF/Magnetics.py b/simpegPF/Magnetics.py index dc96baa8..8b8e88cd 100644 --- a/simpegPF/Magnetics.py +++ b/simpegPF/Magnetics.py @@ -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