diff --git a/example/.gmtdefaults4 b/example/.gmtdefaults4 index 3d0aae9..2ecdffe 100644 --- a/example/.gmtdefaults4 +++ b/example/.gmtdefaults4 @@ -1,5 +1,5 @@ # -# GMT-SYSTEM 4.5.15 [64-bit] Defaults file +# GMT-SYSTEM 4.5.16 [64-bit] Defaults file # #-------- Plot Media Parameters ------------- PAGE_COLOR = white diff --git a/example/DSurfTomo.in b/example/DSurfTomo.in index 7cafb5a..a7b6c89 100644 --- a/example/DSurfTomo.in +++ b/example/DSurfTomo.in @@ -1,16 +1,16 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c INPUT PARAMETERS cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc -surfdataTB.dat c: data file +surfdata.dat c: data file 18 18 9 c: nx ny nz (grid number in lat lon and depth direction) 25.2 121.35 c: goxd gozd (upper left point,[lat,lon]) 0.015 0.017 c: dvxd dvzd (grid interval in lat and lon direction) -25 c: nsrc*maxf +20 c: nsrc*maxf 4.0 0.0 c: weight damp 3 c: nsublayer (numbers of sublayers for each grid interval:grid --> layer) 0.5 2.8 c: minimum velocity, maximum velocity 10 c: maxiter (iteration number) -0.2 c: sparsity fraction +0.1 c: sparsity fraction 26 c: kmaxRc (followed by periods) 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0 0 c: kmaxRg @@ -18,4 +18,4 @@ surfdataTB.dat c: data file 0 c: kmaxLg 0 c: synthetic flag(0:real data,1:synthetic) 0.02 c: noiselevel -2.5 c: threshold +0.5 c: threshold diff --git a/scripts/GenerateIniMOD.py b/scripts/GenerateIniMOD.py index af79196..2fc776a 100644 --- a/scripts/GenerateIniMOD.py +++ b/scripts/GenerateIniMOD.py @@ -6,20 +6,19 @@ import numpy as np #parameters need to be changed #start -nx=75 -ny=96 -nz=17 +nx=18 +ny=18 +nz=8 minvel=0.8 velgrad=0.5 -dep1=1.5+np.array([-1.5, 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0,10.0,11.0,13.0,16.0,20.0,30.0]) -vel1=np.loadtxt('mod.1d') +dep1=np.array([0,0.2,0.4,0.6,0.8,1.1,1.4,1.8,2.5]) #end vs1=np.zeros(nz) mod=np.zeros((nz*ny,nx)) for k in range(nz): for j in range(ny): for i in range(nx): - mod[k*ny+j,i]= vel1[k]/1.75#minvel+dep1[k]*velgrad + mod[k*ny+j,i]= minvel+dep1[k]*velgrad with open('MOD','w') as fp: for i in range(nz): fp.write('%5.1f' % dep1[i]) diff --git a/scripts/GenerateTrueMOD.py b/scripts/GenerateTrueMOD.py index a0c0ea1..0ccb769 100644 --- a/scripts/GenerateTrueMOD.py +++ b/scripts/GenerateTrueMOD.py @@ -12,7 +12,7 @@ nz=9 minvel=0.9 velgrad=0.6 dep1=[0,0.2,0.4,0.6,0.8,1.1,1.4,1.8,2.5] -anosize=1.0 +anosize=0.5 amplitude=0.4 #end x=range(1,nx+1) @@ -35,7 +35,7 @@ for k in range(nz): for i in range(nx): mod[(k)*ny+j,i]=bg[k,i,j]+pxy[k,i,j]*amplitude k=5 -plt.imshow(mod[k*ny:(k+1)*ny,:],cmap='jet_r',interpolation='bicubic') +plt.imshow(mod[k*ny:(k+1)*ny,:],cmap='jet_r') np.savetxt('MOD.true',mod,fmt='%4.4f') plt.colorbar() plt.show() diff --git a/src/DSurfTomo b/src/DSurfTomo index 00a895b..c22dd26 100755 Binary files a/src/DSurfTomo and b/src/DSurfTomo differ diff --git a/src/main.f90 b/src/main.f90 index 7c7e848..f0ac900 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -382,12 +382,36 @@ program SurfTomo cbst(i) = obst(i) - dsyn(i) enddo - threshold = threshold0+(maxiter/2-iter)/3*0.5 + ! write out rw, iw and cbst for testing + !open(44,file='iw.dat') + !write(44,*) nar + !do i = 2,nar+1 + !write(44,*) iw(i) + !enddo + !do i = 1,nar + !write(44,*) col(i) + !enddo + !close(44) + + !open(44,file='rw.dat') + !do i = 1,nar + !write(44,*) rw(i) + !enddo + !close(44) + + !open(44,file='residal.dat') + !do i = 1,dall + !write(44,*) cbst(i) + !enddo + !close(44) + + !threshold = threshold0+(maxiter/2-iter)/3*0.5 do i = 1,dall - datweight(i) = 1.0 - if(abs(cbst(i)) > threshold) then - datweight(i) = exp(-(abs(cbst(i))-threshold)) - endif +! datweight(i) = 1.0 +! if(abs(cbst(i)) > threshold) then +! datweight(i) = exp(-(abs(cbst(i))-threshold)) +! endif + datweight(i) = 1.0/(1+0.05*exp(cbst(i)**2*threshold0)) cbst(i) = cbst(i)*datweight(i) enddo