change weight scheme

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

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

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])

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()

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