diff --git a/src/VoroTomo_2D_MPI.py b/src/VoroTomo_2D_MPI.py
index c5ef4f9..fc9b4e9 100644
--- a/src/VoroTomo_2D_MPI.py
+++ b/src/VoroTomo_2D_MPI.py
@@ -39,14 +39,14 @@
 ##########################################
 
 # define the absolute path of data and project
-rootpath = '/Users/chengxinjiang/Documents/Github/VoroTomo_SW/example/real_data'
-data = glob.glob(os.path.join(rootpath,'dispersion/LVC_snr8_2wl_8s.dat'))
+rootpath = '/Users/chengxin/Documents/ANU/St_Helens/Voro_Tomo/synthetic/dispersion/Love'
+data = glob.glob(os.path.join(rootpath,'syn_selection_*.dat'))
 
 # useful parameters for location and inversion
-ncell = 300                                                                 # number of Voronoi cells for the target region
-nsets = 300                                                                 # number of realizations
-latmin,dlat,nlat = 35.5,0.05,110                                            # latitude range of the target region
-lonmin,dlon,nlon = -122.5,0.05,100                                          # longitude range of the target region
+ncell = 350                                                                 # number of Voronoi cells for the target region
+nsets = 50                                                                 # number of realizations
+latmin,dlat,nlat = 44.5,0.05,76                                             # latitude range of the target region
+lonmin,dlon,nlon = -124.8,0.05,124                                          # longitude range of the target region
 velall           = np.zeros(((nlon+1)*(nlat+1),nsets))                      # vector of the final model averaged over (nsets) times of realizations
 geogrid = {'latmin':latmin,'lonmin':lonmin,                                 # assemble geographic info into a dict
             'dlat':dlat,'dlon':dlon,
@@ -60,8 +60,6 @@
 
 if rank==0:
     # make sub-directory to save intermediate files
-    if not os.path.isdir(os.path.join(rootpath,'figures')):
-        os.mkdir(os.path.join(rootpath,'figures'))
     if not os.path.isdir(os.path.join(rootpath,'iterations')):
         os.mkdir(os.path.join(rootpath,'iterations'))
 
@@ -86,10 +84,11 @@
 
         # define the size of subsampling
         obs_all = len(dfile)
-        obs_sub = int(np.floor(obs_all*0.7))
+        obs_sub = int(np.floor(obs_all*0.8))
 
         # read traveltime info
-        ave     = float("{:.3f}".format(np.mean(dfile['vel'])))                     # averaged velocity from the observation
+        #ave     = float("{:.3f}".format(np.mean(dfile['vel'])))                     # averaged velocity from the observation
+        ave     = 3
         sphgrid = subf.construct_grid(geogrid,ave)                                  # spherical grid data and initial velocity model
 
         # generate parameters
@@ -167,7 +166,7 @@
         np.savez_compressed(tname,vel=vel)
 
         # plot each inversion results
-        subf.plot_tomo(geogrid,vel,rootpath,iset)
+        #subf.plot_tomo(geogrid,vel,rootpath,iset)
 
 ############################################################
 ######### RUN construct_final_model_residuals NEXT #########
@@ -176,4 +175,4 @@
 # syn time to exit
 comm.barrier()
 if rank == 0:
-    sys.exit()
\ No newline at end of file
+    sys.exit()
diff --git a/src/__pycache__/subfunction.cpython-37.pyc b/src/__pycache__/subfunction.cpython-37.pyc
index f546d41..a8a5ff0 100644
Binary files a/src/__pycache__/subfunction.cpython-37.pyc and b/src/__pycache__/subfunction.cpython-37.pyc differ
diff --git a/src/construct_final_model_residula.py b/src/construct_final_model_residula.py
index a6151d1..4be4299 100644
--- a/src/construct_final_model_residula.py
+++ b/src/construct_final_model_residula.py
@@ -1,5 +1,7 @@
 import os 
 import glob 
+import sys
+import scipy
 import numpy as np 
 import pandas as pd 
 import subfunction as subf 
@@ -19,17 +21,20 @@
 ###################################
 
 # define absolute path
-rootpath = '/Users/chengxinjiang/Documents/ANU/Voro_Tomo/synthetic'
+#rootpath = '/Users/chengxin/Documents/ANU/St_Helens/Voro_Tomo/synthetic/Oct4_WithNoise_0.5s_spacing2/Love'
+rootpath = "/Users/chengxin/Documents/ANU/St_Helens/Voro_Tomo/real_data_inversion_Cascadia/Love_phase/roundII"
 
 # control parameter
-cal_resid  = True                                                           # calculate the residual of initial/final model
-plot_final = False                                                          # plot the average model
-plot_each  = False                                                          # plot tomo at each iteration 
-real_data  = False                                                          # whether this is for real data or not
+cal_raypath = True                                                           # calculate the ray path of the averaged model
+cal_resid  = False                                                           # calculate the residual of initial/final model
+plot_final = False                                                           # plot the average model
+plot_each  = False                                                           # plot tomo at each iteration 
+write_each = False                                                           # write each model into txt file
+real_data  = False                                                           # whether this is for real data or not
 
 # several useful parameters
-pper    = [2,3,4,5,6,8,10,12,14,16,18,20,24,28]                                   # all period range
-nset    = 300                                                               # number of inversions
+pper    = [3,4,5,6,8,10,12,14,16,18,20,24,28,32,36,40]                                    # all period range
+nset    = 70                                                                 # number of inversions
 
 # check whether dir exists or not
 if not os.path.isdir(os.path.join(rootpath,'models')):
@@ -40,12 +45,19 @@
     os.mkdir(os.path.join(rootpath,'figures'))
 
 # define geographic information
-latmin,dlat,nlat = 35.5,0.02,100                                            # latitude range of the target region
-lonmin,dlon,nlon = -120.0,0.02,100                                          # longitude range of the target region
+latmin,dlat,nlat = 24.5,0.2,123                                            # latitude range of the target region
+lonmin,dlon,nlon = -124.8,0.2,289                                          # longitude range of the target region
+
+latmin,dlat,nlat = 44.5,0.02,190                                            # latitude range of the target region
+lonmin,dlon,nlon = -124.8,0.02,310                                          # longitude range of the target region
 geogrid = {'latmin':latmin,'lonmin':lonmin,                                 # assemble geographic info into a dict
             'dlat':dlat,'dlon':dlon,
             'nlat':nlat,'nlon':nlon}
 
+# grid for output files
+latgrid = np.linspace(latmin+nlat*dlat,latmin,nlat+1) 
+longrid = np.linspace(lonmin,lonmin+nlon*dlon,nlon+1)
+lon,lat = np.meshgrid(longrid,latgrid)
 
 #######################################
 ######## LOOPS THROUGH PERIODS ########
@@ -53,11 +65,18 @@
 
 # loop through all periods
 for iper in range(len(pper)):
+    tdlat,tnlat = 0.1,38
+    tdlon,tnlon = 0.1,62
+    tlatgrid = np.linspace(latmin+tnlat*tdlat,latmin,tnlat+1)
+    tlongrid = np.linspace(lonmin,lonmin+tnlon*tdlon,tnlon+1)
+    tlon,tlat = np.meshgrid(tlongrid,tlatgrid)
+    tcount = np.zeros(tlon.size,dtype=np.int16)
     per = pper[iper]
     velall = np.zeros(((nlon+1)*(nlat+1),nset))                             # vector of the final model averaged over (nsets) times of realizations
 
     # find the observational data
-    dfile = os.path.join(rootpath,'dispersion/data_snr8_'+str(per)+'s.dat')
+    #dfile = os.path.join(rootpath,'dispersion/Love_phase/new_L{0:02d}_USANT15.txt'.format(per))
+    dfile = os.path.join(rootpath,'../dispersion/selection_'+str(per)+'s.dat')
     if not os.path.isfile(dfile):
         raise ValueError('double check! cannot find %s'%dfile)
 
@@ -69,65 +88,112 @@
         ave = 3.0
 
     # load the velocity model
-    allfiles = glob.glob(os.path.join(rootpath,'iterations/'+str(per)+'s_*.npz'))
+    allfiles = sorted(glob.glob(os.path.join(rootpath,'iterations/{0:d}s_*.npz'.format(per))))
     if not len(allfiles):
         raise ValueError('no model generated! double check')
 
     # load all models into a big matrix
     ii = 0
-    for ifile in allfiles:
+    for ifile in allfiles[:nset]:
         try:
             tvel = np.load(ifile)['vel']
             if plot_each:
                 subf.plot_tomo(geogrid,-tvel*ave*ave+ave,rootpath,ii)
+            if write_each:
+                tvel1 = -tvel*ave*ave+ave
+                # format the data column
+                fout0 = open(rootpath+'/models/new_tomo_'+str(per)+'s_'+str(ii)+'m.txt','w')
+                for jj in range(lon.size):
+                    fout0.write('%6d %8.3f %8.3f %6.3f %6.3f 0.0\n'%(jj,lon.flatten()[jj],lat.flatten()[jj],tvel[jj],tvel1[jj]))
+                fout0.close()
             velall[:,ii]  = tvel
             ii+=1
         except Exception as err:
             print('cannot load %s because %s'%(ifile,err))
 
+    print('total # of models %d'%ii)
     # remove empty traces
     velall = velall[:,:ii]
 
     # averaged model and get stad
-    tvel = np.mean(velall,axis=1)
-    vel_std = np.std(velall,axis=1)
-
-    # calculate the absolute vel and vel pertb
-    vel_abs = -tvel*ave*ave+ave
+    tvel = -velall*ave*ave+ave
+    vel_abs = np.mean(tvel,axis=1)
+    vel_std = np.std(tvel,axis=1)
     vel_per = (vel_abs-ave)/ave*100
-    #print('old and new averages are %6.3f %6.3f'%(ave,nave))
 
-    # output into txt file
-    latgrid = np.linspace(latmin+nlat*dlat,latmin,nlat+1) 
-    longrid = np.linspace(lonmin,lonmin+nlon*dlon,nlon+1)
-    lon,lat = np.meshgrid(longrid,latgrid)
+    #######################################
+    ######## CALCULATE RAYPATPH ###########
+    #######################################
 
+    if cal_raypath:
+        sphgrid = subf.construct_grid(geogrid,ave)                                  # spherical grid data and initial velocity model
+        xygrid = np.zeros(shape=(tlon.size,2),dtype=np.float32)
+        ttheta,tphi = subf.geo2sph(tlat.flatten(),tlon.flatten())
+        xygrid[:,0],xygrid[:,1] = subf.sph2xyz(ttheta,tphi)
+        sphgrid['vel1'] = np.reshape(vel_abs,(1,sphgrid['ntheta'],sphgrid['nphi']))
+
+        #######################################
+        ######## CALCULATE RESIDUALS ##########
+        #######################################
+
+        if cal_resid:
+            fname = os.path.join(rootpath,'residuals/tomo_'+str(per)+'s.txt')                     
+            res_initial,res_final = subf.calculate_residuals(data,geogrid,sphgrid,fname)
+            
+            # plot the histogram
+            nbin = 30
+            plt.close('all')
+            plt.hist(res_initial,nbin,density=True,facecolor='g',alpha=0.6,label='initial')
+            plt.hist(res_final,nbin,density=True,facecolor='b',alpha=0.6,label='final')
+            plt.legend(loc='upper right')
+            plt.savefig(rootpath+'/residuals/hist'+str(per)+'s.pdf',format='pdf',mpi=300)
+        
+        srcs    = data.evt_id.unique()
+        # loop through each virtual source
+        for isc in srcs:
+            print(len(srcs),isc)
+
+            # source location to initialize the traveltime field
+            arr4rc = data[data.evt_id.values==isc]
+            lat_s,lon_s   = arr4rc['evt_lat'].iloc[0],arr4rc['evt_lon'].iloc[0]
+
+            # generate parameters
+            solver = subf.solve_source_region(geogrid,sphgrid,lat_s,lon_s,1)
+
+            # loop through each receiver 
+            for irc in range(len(arr4rc)):
+                try:
+
+                    # receiver location
+                    lat_r,lon_r = arr4rc['sta_lat'].iloc[irc],arr4rc['sta_lon'].iloc[irc]
+                    theta_r,phi_r = subf.geo2sph(lat_r,lon_r)
+                    src_pos = (RR,theta_r,phi_r)
+
+                    # do ray tracing here
+                    ray = solver.trace_ray(np.array(src_pos))
+                    ray = pd.DataFrame(ray,columns=('rho','theta','phi'))
+                    ray['x'] = ray['rho']*np.sin(ray['theta'])*np.cos(ray['phi'])
+                    ray['y'] = ray['rho']*np.sin(ray['theta'])*np.sin(ray['phi'])
+                    dist = scipy.spatial.distance.cdist(ray[['x','y']].values,xygrid)
+                    argmin = np.argmin(dist,axis=1)
+                    uargmin = np.unique(argmin)
+                    tcount[uargmin] += 1
+
+                except Exception as err:
+                    print('ERROR for',irc,isc,err)               
+
+    fout1 = open('{0:s}/ray_{1:02d}s.txt'.format(rootpath,per),'w')
+    ttlon,ttlat = tlon.flatten(),tlat.flatten()
+    for ii in range(tlon.size):
+       fout1.write('%8.3f %8.3f %d\n'%(ttlon[ii],ttlat[ii],tcount[ii]))
+    fout1.close()
     # format the data column
-    fout = open(rootpath+'/models/tomo_'+str(per)+'s.txt','w')
+    fout = open('{0:s}/models/L_tomo_{1:02d}s.txt'.format(rootpath,per),'w')
     for ii in range(lon.size):
         fout.write('%6d %8.3f %8.3f %6.3f %6.3f %8.4f\n'%(ii,lon.flatten()[ii],lat.flatten()[ii],vel_per[ii],\
             vel_abs[ii],vel_std[ii]))
     fout.close()
 
-
-    #######################################
-    ######## CALCULATE RESIDUALS ##########
-    #######################################
-
-    if cal_resid:
-        sphgrid = subf.construct_grid(geogrid,ave)  
-        sphgrid['vel1'] = np.reshape(vel_abs,(1,sphgrid['ntheta'],sphgrid['nphi']))
-        fname = os.path.join(rootpath,'residuals/tomo_'+str(per)+'s.txt')                     
-        res_initial,res_final = subf.calculate_residuals(data,geogrid,sphgrid,fname)
-        
-        # plot the histogram
-        nbin = 30
-        plt.close('all')
-        plt.hist(res_initial,nbin,density=True,facecolor='g',alpha=0.6,label='initial')
-        plt.hist(res_final,nbin,density=True,facecolor='b',alpha=0.6,label='final')
-        plt.legend(loc='upper right')
-        plt.savefig(rootpath+'/residuals/hist'+str(per)+'s.pdf',format='pdf',mpi=300)
-
     # plot the model
     if plot_final:
-        subf.plot_tomo(geogrid,vel_abs,rootpath,-1)
\ No newline at end of file
+        subf.plot_tomo(geogrid,vel_abs,rootpath,-1)
diff --git a/src/make_syn_2D.py b/src/make_syn_2D.py
index 00614e2..813a7b7 100644
--- a/src/make_syn_2D.py
+++ b/src/make_syn_2D.py
@@ -18,12 +18,12 @@
 syn_option = 'checkerboard'
 
 # absoluate path for traveltime data and for output
-rootpath = '/Users/chengxinjiang/Documents/Github/VoroTomo_SW/example/synthetic'
-dfiles   = glob.glob(os.path.join(rootpath,'dispersion/LVC_snr8_2wl_8s.dat'))
+rootpath = '/Users/chengxin/Documents/ANU/St_Helens/Voro_Tomo/synthetic/Oct4_WithNoise_0.5s_spacing2/Love'
+dfiles   = glob.glob(os.path.join(rootpath,'selection_*s.dat'))
 
 # full 2D research area
-latmin,dlat,nlat = 35.5,0.05,110                                            # latitude range of the target region
-lonmin,dlon,nlon = -122.5,0.05,100                                          # longitude range of the target region
+latmin,dlat,nlat = 44.5,0.02,190                                            # latitude range of the target region
+lonmin,dlon,nlon = -124.8,0.02,310                                          # longitude range of the target region
 geogrid = {'latmin':latmin,'lonmin':lonmin,                                 # assemble geographic info into a dict 
             'dlat':dlat,'dlon':dlon,
             'nlat':nlat,'nlon':nlon}
@@ -31,8 +31,9 @@
 if syn_option=='checkerboard':
     # target anomaly information
     anomaly = {'avs':3,'amp':0.3,                                           # amplitude in percentage
-                'size':15,'spacing': True,                                  # size of anomaly and whether of spacing
-                'noise':0.3}
+                'size':15,'spacing': 2,                                      # size of anomaly and whether of spacing
+                'noise_ave':0,
+                'noise_std':0.5}
 
     # make sphgrid for traveltime prediction
     sphgrid = subf.construct_grid(geogrid,anomaly['avs'])                       # spherical grid data and initial velocity model
@@ -46,7 +47,7 @@
     anomaly = {'avs':3,'amp':-0.3,                                               # amplitude in percentage
                 'latmin':46.1,'latmax':46.4,                                     # size of anomaly and whether of spacing
                 'lonmin':-122.4,'lonmax':-122.1,
-                'noise':0.3}
+                'noise':0}
 
     # make sphgrid for traveltime prediction
     sphgrid = subf.construct_grid(geogrid,anomaly['avs'])                       # spherical grid data and initial velocity model
@@ -97,8 +98,9 @@
                 theta_r,phi_r = subf.geo2sph(lat_r,lon_r)
                 src_pos  = (RR,theta_r,phi_r)
                 syn_data = solver.traveltime.value(np.array(src_pos))
-                if anomaly['noise']:
-                    syn_data += np.random.normal(loc=anomaly['noise'],scale=0.05)
+                if anomaly['noise_std']:
+                    syn_data += np.random.normal(loc=anomaly['noise_ave'],scale=anomaly['noise_std'])
+
                 fout.write('%6d,%6d,%8.3f,%8.3f,%6d,%8.3f,%8.3f,%8.3f,%8.3f\n'%(arr4rc['index'].iloc[irc],\
                     isc,lon_s,lat_s,arr4rc['sta_id'].iloc[irc],lon_r,lat_r,arr4rc['dist'].iloc[irc]/syn_data,\
                     arr4rc['dist'].iloc[irc]))
@@ -107,4 +109,4 @@
                 print('ERROR for',irc,isc,err)
 
     fout.close()
-    print('synthetic done for %s'%dfiles[iper])
\ No newline at end of file
+    print('synthetic done for %s'%dfiles[iper])
diff --git a/src/remove_data_residuals.py b/src/remove_data_residuals.py
index ed318ae..c8c7cac 100644
--- a/src/remove_data_residuals.py
+++ b/src/remove_data_residuals.py
@@ -16,8 +16,8 @@
 '''
 
 # define rootpath and find all files
-rootpath = '/Users/chengxinjiang/Documents/ANU/Voro_Tomo/real_data'
-pper = [3,6,8,12,14,18,20,24,28]
+rootpath = '/Users/chengxin/Documents/ANU/St_Helens/Voro_Tomo/real_data_inversion_Cascadia/Love_group'
+pper = [2,3,4,5,6,8,10,12,14,16,18,20,24,28,32,36,40]
 
 # loop through each period
 for ii in range(len(pper)):
@@ -34,8 +34,8 @@
     dist = np.array(data['dist'])
     vel  = np.array(data['vel'])
     tt   = np.array(dist/vel)
-    res1 = np.array(data['res_initial'])
-    res2 = np.array(data['res_final'])
+    res1 = np.array(data['res_initial']).astype(np.float32)
+    res2 = np.array(data['res_final']).astype(np.float32)
 
     # get statistics
     ave1 = np.mean(res1)
@@ -45,7 +45,7 @@
 
     # make selections
     indx = np.where((res2>=ave2-2*std2)&(res2<=ave2+2*std2))[0]
-    fout = open(rootpath+'/dispersion/selection2_'+str(per)+'s.dat','w')
+    fout = open(rootpath+'/dispersion/selection_'+str(per)+'s.dat','w')
     fout.write('index,evt_id,evt_lon,evt_lat,sta_id,sta_lon,sta_lat,vel,dist\n')
     for jj in range(len(indx)):
         tindx = indx[jj]
@@ -81,6 +81,6 @@
     plt.text(15+np.min(dist),0.85*np.max(tt),text)
     plt.title(str(per)+'s final model')
     plt.colorbar()
-    outname = rootpath+'/dispersion/obs2_residual_'+str(per)+'s.pdf'
+    outname = rootpath+'/dispersion/obs_residual_'+str(per)+'s.pdf'
     plt.savefig(outname,format='pdf',dpi=300)
     #plt.show()
diff --git a/src/select_iterations.py b/src/select_iterations.py
new file mode 100644
index 0000000..f0400ba
--- /dev/null
+++ b/src/select_iterations.py
@@ -0,0 +1,124 @@
+import os 
+import glob 
+import numpy as np 
+import pandas as pd 
+import subfunction as subf
+import matplotlib.pyplot as plt 
+
+'''
+approximate the evolution of MSE (mean squared error) with the number 
+of iterations for selection of optimal # of iterations
+
+by Chengxin @ANU (Mar2020)
+'''
+
+###################################
+######## PARAMETER SECTION ########
+###################################
+
+# define absolute path
+rootpath = '/Users/chengxin/Documents/ANU/St_Helens/Voro_Tomo/real_data/phase/iter2'
+
+# control parameter
+plot_each  = False                                                          # plot tomo at each iteration 
+real_data  = True                                                           # whether this is for real data or not
+write_modle= True                                                           # write the final model into file
+
+# several parameters
+nset  = 300                                                                 # number of inversions
+per   = 4                                                                   # target period range
+nnmod = [2,4,6,8,10,14,18,25,40,60,80,120,160,200,250,300]                      # number of models for sub-sampling
+niter = len(nnmod)                                                          # number of final models to output
+
+# check whether dir exists or not
+if not os.path.isdir(os.path.join(rootpath,'select_iterations')):
+    os.mkdir(os.path.join(rootpath,'select_iterations'))
+
+# define geographic information
+latmin,dlat,nlat = 45.0,0.02,150                                            # latitude range of the target region
+lonmin,dlon,nlon = -124.3,0.02,250                                          # longitude range of the target region
+geogrid = {'latmin':latmin,'lonmin':lonmin,                                 # assemble geographic info into a dict
+            'dlat':dlat,'dlon':dlon,
+            'nlat':nlat,'nlon':nlon}
+
+# find the observational data
+dfile = os.path.join(rootpath,'dispersion_selected/Rayleigh/selection2_'+str(per)+'s.dat')
+if not os.path.isfile(dfile):
+    raise ValueError('double check! cannot find %s'%dfile)
+
+# get averaged velocity
+data  = pd.read_csv(dfile)
+if real_data:
+    ave = np.mean(data['vel'])
+else:
+    ave = 4.0
+
+###################################
+######## DO SUB-SAMPLING ##########
+###################################
+
+# load the velocity model
+allfiles = glob.glob(os.path.join(rootpath,'Rayleigh/iterations/'+str(per)+'s_*.npz'))
+if not len(allfiles):
+    raise ValueError('no model generated! double check')
+
+'''
+# loop through each subsampling
+for iii in range(niter):
+    nmod   = nnmod[iii]
+    velall = np.zeros(((nlon+1)*(nlat+1),nmod))                       # vector of the final model averaged over (nsets) times of realizations
+    indx = np.random.randint(0,nset-1,size=nmod)
+
+    # load all models into a big matrix
+    ii = 0
+    for tindx in indx:
+        ifile = allfiles[tindx]
+        print(ifile,tindx)
+        try:
+            velall[:,ii] = np.load(ifile)['vel']
+            ii+=1
+        except Exception as err:
+            print('cannot load %s because %s'%(ifile,err))
+
+    # remove empty traces
+    velall = velall[:,:ii]
+
+    # averaged model and get stad
+    tvel = -velall*ave*ave+ave
+    vel_abs = np.mean(tvel,axis=1)
+    vel_std = np.std(tvel,axis=1)
+    vel_per = (vel_abs-ave)/ave*100
+
+    # output into txt file
+    latgrid = np.linspace(latmin+nlat*dlat,latmin,nlat+1) 
+    longrid = np.linspace(lonmin,lonmin+nlon*dlon,nlon+1)
+    lon,lat = np.meshgrid(longrid,latgrid)
+
+    # calculate residuals 
+    sphgrid = subf.construct_grid(geogrid,ave)  
+    sphgrid['vel1'] = np.reshape(vel_abs,(1,sphgrid['ntheta'],sphgrid['nphi']))
+    fname = os.path.join(rootpath,'select_iterations/residual_'+str(per)+'s_'+str(nmod)+'iter.dat')                     
+    res_initial,res_final = subf.calculate_residuals(data,geogrid,sphgrid,fname)
+
+    # format the data column
+    if write_modle:
+        fout = open(rootpath+'/select_iterations/tomo_'+str(per)+'s'+str(nmod)+'iter.dat','w')
+        for ii in range(lon.size):
+            fout.write('%6d %8.3f %8.3f %6.3f %6.3f %8.4f\n'%(ii,lon.flatten()[ii],lat.flatten()[ii],vel_per[ii],\
+                vel_abs[ii],vel_std[ii]))
+        fout.close()
+'''
+
+# summarize the MSE and plot the results
+mse = np.zeros(len(nnmod),dtype=np.float32)
+for ii,nmod in enumerate(nnmod):
+    tfile = os.path.join(rootpath,'select_iterations/residual_'+str(per)+'s_'+str(nmod)+'iter.dat')
+    tdata = pd.read_csv(tfile)
+    res   = np.array(tdata['res_final'])
+    mse[ii] = np.sum(res*res)/len(res)
+
+# plot the results
+plt.figure(figsize=(6,4))
+plt.plot(nnmod,mse,'o-')
+plt.savefig(rootpath+'/select_iterations/evolution.pdf', format='pdf', dpi=400)
+plt.close()
\ No newline at end of file
diff --git a/src/subfunction.py b/src/subfunction.py
index 3d2faeb..114ed8c 100644
--- a/src/subfunction.py
+++ b/src/subfunction.py
@@ -8,6 +8,7 @@
 from scipy.sparse import coo_matrix
 from collections import defaultdict
 from matplotlib import pyplot as plt
+from scipy.ndimage import gaussian_filter
 
 '''
 subfunctions used in VoroTomo_2D
@@ -288,10 +289,8 @@ def make_checkerboard(geogrid,anomaly_info):
     '''
 
     # check the basic info of the targeted checkerboard
-    checker_vs   = anomaly_info['avs']
     checker_size = anomaly_info['size']
     checker_spa  = anomaly_info['spacing']
-    random_noise = anomaly_info['noise']
 
     # load parameters from the geogrid dict
     latmin,dlat,nlat = geogrid['latmin'],geogrid['dlat'],geogrid['nlat']+1
@@ -300,10 +299,10 @@ def make_checkerboard(geogrid,anomaly_info):
     # construct one basic anomaly
     if checker_spa:
         # spike input
-        anomaly1 = np.ones(4*checker_size,dtype=np.float32)
+        anomaly1 = np.ones(2*(checker_spa+1)*checker_size,dtype=np.float32)
         anomaly1[:checker_size] = -1
-        anomaly1[checker_size:2*checker_size] = 0
-        anomaly1[3*checker_size:] = 0
+        anomaly1[checker_size:(checker_spa+1)*checker_size] = 0
+        anomaly1[(checker_spa+2)*checker_size:] = 0
         anomaly2 = anomaly1*-1
         #anomaly3 = np.concatenate((anomaly1[checker_size:],anomaly1[:checker_size]))
     else:
@@ -325,6 +324,9 @@ def make_checkerboard(geogrid,anomaly_info):
         else:
             amp.append(np.zeros(nlon))
 
+    # smooth using a gaussian filter
+    amp= gaussian_filter(amp, sigma=1.5)
+
     return amp
 
 def make_syn(geogrid,anomaly):
@@ -451,4 +453,4 @@ def plot_tomo(geogrid,vel,rootpath,iset=-1):
     else:
         outfname = rootpath+'/figures/tomo_'+str(iset)+'.pdf'
     plt.savefig(outfname, format='pdf', dpi=400)
-    plt.close()
\ No newline at end of file
+    plt.close()
diff --git a/src/test/convert_to_fmst.py b/src/test/convert_to_fmst.py
new file mode 100644
index 0000000..4950b3c
--- /dev/null
+++ b/src/test/convert_to_fmst.py
@@ -0,0 +1,45 @@
+import os 
+import glob 
+import numpy as np 
+import pandas as pd 
+
+'''
+this script reads the traveltime for VoroTomo_2D and convert it
+for fmst code
+'''
+
+# absolute path for the data
+rootpath = '/Users/chengxin/Documents/ANU/St_Helens/Voro_Tomo/synthetic/fmst_tomo/ttime/Rayleigh'
+dfiles   = glob.glob(os.path.join(rootpath,'selection*.dat'))
+
+# total number of stations
+nsta = 190
+std  = 0.2
+
+# loop through each file
+for ii in range(len(dfiles)):
+    dfile = dfiles[ii]
+    per   = dfile.split('/')[-1].split('.')[0].split('_')[-1]
+
+    # load data
+    data  = pd.read_csv(dfile)
+    indxs = np.array(data['evt_id'])
+    indxr = np.array(data['sta_id'])
+    vel   = np.array(data['vel'])
+    dist  = np.array(data['dist'])
+
+    # ready for conversion
+    fout = open(rootpath+'/otimes_'+str(per)+'.dat','w')
+    for iss in range(nsta):
+        tindx1 = np.where(indxs==iss)[0]
+        for irr in range(nsta):
+            tindx2 = np.where(indxr==irr)[0]
+            tindx  = np.intersect1d(tindx1,tindx2)
+            print('found %d for %d and %d pairs'%(len(tindx),iss,irr))
+
+            if len(tindx)==1:
+                fout.write('1 %10.6f %3.1f\n'%(dist[tindx]/vel[tindx],std))
+            else:
+                fout.write('0 0 10\n')
+    
+    fout.close()
\ No newline at end of file
diff --git a/src/test/make_inputs.py b/src/test/make_inputs.py
index d969681..738a61d 100644
--- a/src/test/make_inputs.py
+++ b/src/test/make_inputs.py
@@ -10,12 +10,12 @@
 '''
 
 # absolute path for the data
-rootpath = '/Users/chengxinjiang/Documents/Github/VoroTomo_SW/example/real_data'
-dfiles   = glob.glob(os.path.join(rootpath,'dispersion/LVC_snr8_*s.dat'))
+rootpath = '/Users/chengxin/Documents/ANU/St_Helens/CCFs/merge/Love_phase'
+dfiles   = glob.glob(os.path.join(rootpath,'data_snr8_*s.dat'))
 
 # loop through the file
 for ifile in dfiles:
-    data = pd.read_csv(ifile)
+    data = pd.read_csv(ifile,sep='\s+')
 
     # find the unique stations
     slon = np.array(data['lons'])
@@ -45,7 +45,7 @@
             sta_list.append([rlon[ii],rlat[ii]])
 
     # output to a file
-    tfile = rootpath+'/dispersion/new_'+ifile.split('/')[-1]
+    tfile = rootpath+'/new_'+ifile.split('/')[-1]
     fout = open(tfile,'w')
     fout.write('index,evt_id,evt_lon,evt_lat,sta_id,sta_lon,sta_lat,vel,dist\n')