correcting the readme

.gitignore

@@ -4,4 +4,4 @@ Applications
 example_*
 .gitignore
 temp
-
+DSurfTomo

README.md

@@ -1,11 +1,8 @@
-DSurfTomo is a surface wave inversion program which can irectly invert surface wave dispersion data directly to 3D shear wavespeed models, without the intermediate step of constructing the phase or group velocity maps.
+DSurfTomo is a surface wave tomography program which inverts surface wave dispersion data directly to 3D shear wavespeed models without the intermediate step of constructing the phase or group velocity maps.
 The fast marching method (FMM) (Rawlinson et al., 2004) is used to compute, at each period, surface wave travel times and ray paths between sources and receivers. This avoids the assumption of great-circle propagation that is used in most surface wave tomographic studies, but which is not appropriate in complex media. 
-Please refer to Fang et al. (2015 , GJI) for the detail description of the method.
+To show its usage, an example of Taipei Basin tomography is provided, including scripts that are used for data reformating and plotting results. Interested users are recommended to refer to Fang et al. (2015 , GJI) for the detail description of the method.
 
 Fang, H., Yao, H., Zhang, H., Huang, Y. C., & van der Hilst, R. D. (2015). Direct inversion of surface wave dispersion for three-dimensional shallow crustal structure based on ray tracing: methodology and application. Geophysical Journal International, 201(3), 1251-1263. 
 
 Rawlinson, N. & Sambridge, M., 2004. Wave front evolution in strongly heterogeneous layered media using the fast marching method, Geophys. J. Int., 156(3), 631–647 
 
-For V2.0, I incorporated the random projections based inversion using Poisson Voronoi cells. However, it has not been fully tested, so any feedbacks will be helpful. For more information about the inversion method, please refer to:
-
-Fang, H., van der Hilst, R. D., de Hoop, M. V., Kothari, K., Gupta, S., & Dokmanić, I. (2020). Parsimonious Seismic Tomography with Poisson Voronoi Projections: Methodology and Validation. Seismological Research Letters, 91(1), 343-355.

scripts/GenerateDSurfTomoInputFile.py

@@ -0,0 +1,146 @@
+#/* -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.
+#
+# File Name : GenerateDSurfTomoInputFile.py
+#
+# Purpose : Generate the Input File for DSurfTomo using the data themselves
+#
+# Creation Date : 05-07-2023
+#
+# Last Modified : Tue 18 Jul 2023 08:10:22 PM CST
+#
+# Created By : Hongjian Fang: fanghj1990@gmail.com
+#
+#_._._._._._._._._._._._._._._._._._._._._.*/
+
+import os
+import numpy as np
+import glob
+import pandas as pd
+
+#These are for generating the right format
+DataPath = 'TaipeiRawData'
+DataNames = 'CDisp*.dat'
+wavetype = 2
+phasetype = 0
+
+# These following can be changed in DSurfTomo.in
+nz = 9
+sublayers = 3
+minvel = 0.5
+maxvel = 2.8
+noiselevel = 0.02
+datathreshold = 1.5
+DataFiles = glob.glob(DataPath+'/'+DataNames)
+
+sta1latall = []
+sta1lonall = []
+sta2latall = []
+sta2lonall = []
+periodsall = []
+dispersionall = []
+stationpairID = []
+
+for ifile in DataFiles:
+    dispersion = np.loadtxt(ifile)
+    pairID = ifile.split('/')[-1]
+    sta1lon = dispersion[0,0]
+    sta1lat = dispersion[0,1]
+    sta2lon = dispersion[1,0]
+    sta2lat = dispersion[1,1]
+    nperiods,_ = dispersion.shape
+    nperiods = nperiods - 2
+    periods = np.zeros(nperiods,)
+    disper = np.zeros(nperiods,)
+    for ii in range(nperiods):
+        periods[ii] = dispersion[ii+2,0]
+        disper[ii] = dispersion[ii+2,1]
+   
+    dispersionall = np.hstack([dispersionall,disper])
+    periodsall = np.hstack([periodsall,periods])
+    sta1latall = np.hstack([sta1latall,np.ones(nperiods,)*sta1lat])
+    sta1lonall = np.hstack([sta1lonall,np.ones(nperiods,)*sta1lon])
+    sta2latall = np.hstack([sta2latall,np.ones(nperiods,)*sta2lat])
+    sta2lonall = np.hstack([sta2lonall,np.ones(nperiods,)*sta2lon])
+    stationpairID = stationpairID + [pairID] * nperiods
+
+dataall = pd.DataFrame({'sta1lat':sta1latall, 'sta1lon':sta1lonall, \
+                        'sta2lat':sta2latall, 'sta2lon':sta2lonall, \
+                        'periods':periodsall, 'dispersion': dispersionall, \
+                        'pairid':stationpairID})
+
+dataall = dataall.sort_values(by = ['periods', 'pairid'])
+dataall = dataall.set_index('periods')
+with open('surfdata.dat','w') as fout:
+    UniqPeriods = dataall.index.unique()
+    for iperiod,period in enumerate(UniqPeriods):
+        datasubperiod = dataall.loc[period]
+        if isinstance(datasubperiod,pd.Series):
+            continue
+        datasubperiod = datasubperiod.reset_index().set_index('sta1lat')
+        sta1lat = datasubperiod.index.unique()
+        for ista1 in sta1lat:
+            datasubstation = datasubperiod.loc[ista1]
+            if isinstance(datasubstation,pd.DataFrame):
+                fout.write(f'# {datasubstation.index[0]} {datasubstation["sta1lon"].iloc[0]} {iperiod+1} {wavetype} {phasetype}\n')
+                for ista2 in range(len(datasubstation)):
+                    fout.write(f'{datasubstation["sta2lat"].iloc[ista2]} {datasubstation["sta2lon"].iloc[ista2]} {datasubstation["dispersion"].iloc[ista2]}\n')
+            else:
+                fout.write(f'# {datasubstation.name} {datasubstation["sta1lon"]} {iperiod+1} {wavetype} {phasetype}\n')
+                fout.write(f'{datasubstation["sta2lat"]} {datasubstation["sta2lon"]} {datasubstation["dispersion"]}\n')
+
+
+print('Finish reformatting disperison data')
+
+# Determin the grid interval, 1/3 of the smallest station distance
+
+# Weight and Damp, you will have to run the code a couple of time to decide
+# The default values is 5.0 and 1.0, respectively
+weight = 1.0
+damp = 1.0
+
+# Sparsity fraction, this parameter represent the sparsity of the sensitivity matrix and is set to be 0.2 for safe keeping, but could be as low as 0.02.it is not important as long as your code does not throw you memory erros
+sparsityfraction = 0.1
+
+# Iteration number, the default is 10 time, most often you get stable results after 4-5 iterations.
+maxiteration = 10
+
+
+distall = np.sqrt((dataall['sta2lat']-dataall['sta1lat'])**2+(dataall['sta2lon']-dataall['sta1lon'])**2)
+mindist = distall.min()
+gridintval = int(np.ceil(mindist/3*1000))/1000
+# Determine the origin, 1 grids
+originLat = np.max([np.max(dataall.sta1lat.max()),np.max(dataall.sta2lat.max())]) + 1*gridintval
+originLon = np.min([np.min(dataall.sta1lon.min()),np.min(dataall.sta2lon.min())]) - 1*gridintval
+largestLat = np.max([dataall['sta1lat'].max()-dataall['sta1lat'].min(),dataall['sta2lat'].max()-dataall['sta2lat'].min()])
+largestLon = np.max([dataall['sta1lon'].max()-dataall['sta1lon'].min(),dataall['sta2lon'].max()-dataall['sta2lon'].min()])
+inputfile = open('DSurfTomo.in','w')
+dummy = 'cccc'
+inputfile.write(f'{dummy*10}\n')
+inputfile.write(f'{dummy*10}\n')
+inputfile.write(f'{dummy*10}\n')
+inputfile.write('surfdata.dat\n')
+nx = int(np.ceil(largestLat/gridintval))+5
+ny = int(np.ceil(largestLon/gridintval))+5
+nreceivers = len(dataall.sta1lat.unique())+1
+inputfile.write(f'{nx} {ny} {nz}\n')
+inputfile.write(f'{originLat:<7.3f} {originLon:7.3f}\n')
+inputfile.write(f'{gridintval:<7.3f} {gridintval:7.3f}\n')
+inputfile.write(f'{nreceivers}\n')
+inputfile.write(f'{weight} {damp}\n')
+inputfile.write(f'{sublayers}\n')
+inputfile.write(f'{minvel} {maxvel}\n')
+inputfile.write(f'{maxiteration}\n')
+inputfile.write(f'{sparsityfraction}\n')
+nperiods = len(UniqPeriods)
+inputfile.write(f'{nperiods}\n')
+#inputfile.write(f'{*UniqPeriods}\n')
+inputfile.write(' '.join([str(iperiod) for iperiod in UniqPeriods])+'\n')
+inputfile.write(f'0\n')
+inputfile.write(f'0\n')
+inputfile.write(f'0\n')
+inputfile.write(f'0\n')
+inputfile.write(f'{noiselevel}\n')
+inputfile.write(f'{datathreshold}\n')
+inputfile.close()
+
+print('Finishing generating input file\n')

src/CalSurfG.f90

@@ -1378,13 +1378,11 @@ subroutine CalSurfG(nx,ny,nz,nparpi,vels,iw,rw,col,dsurf, &
         !  if required.
         !
         if (igrt(srcnum,knumi) == 0 .or. (ig == 2 .and. igrt(srcnum,knumi) == 1)) then
-        ! a little stupid, remember to change latter
-        !if (igrt(srcnum,knumi) == 1) then
-        !call gridder(velf0)
-        !endif
         count11=count11+1
         CALL rpaths(x,z,fdm,rcxf(istep,srcnum,knumi),rczf(istep,srcnum,knumi))
         row(1:nparpi)=0.0
+        ! change the Brocher relationship if depth is larger than 35 km         
+        if (depz(nz-1)<35.0) then
         do jj=1,nvz
           do kk=1,nvx
             if(abs(fdm(jj,kk)).ge.ftol) then
@@ -1402,6 +1400,28 @@ subroutine CalSurfG(nx,ny,nz,nparpi,vels,iw,rw,col,dsurf, &
             endif
           enddo
         enddo
+      ! use a different formulation between Vp, Vs, and Rho, according to Luosong's
+      else
+        do jj=1,nvz
+          do kk=1,nvx
+            if(abs(fdm(jj,kk)).ge.ftol) then
+              coe_a=(2.2110-0.8984*2*vels(kk+1,jj+1,1:nz-1)+&
+              0.2786*3*vels(kk+1,jj+1,1:nz-1)**2-0.02412*4*vels(kk+1,jj+1,1:nz-1)**3)
+              vpft=0.9098 + 2.2110*vels(kk+1,jj+1,1:nz-1) - 0.8984*vels(kk+1,jj+1,1:nz-1)**2+ &
+              0.2786*vels(kk+1,jj+1,1:nz-1)**3 - 0.02412*vels(kk+1,jj+1,1:nz-1)**4
+              coe_rho=coe_a*(1.6612-0.4721*2*vpft+&
+              0.0671*3*vpft**2-0.0043*4*vpft**3+&
+              0.000106*5*vpft**4)
+
+              row((jj-1)*nvx+kk:(nz-2)*nvz*nvx+(jj-1)*nvx+kk:nvx*nvz)=&
+                (sen_vp(jj*(nvx+2)+kk+1,knumi,1:nz-1)*coe_a+&
+                sen_rho(jj*(nvx+2)+kk+1,knumi,1:nz-1)*coe_rho+&
+                sen_vs(jj*(nvx+2)+kk+1,knumi,1:nz-1))*fdm(jj,kk)
+            endif
+          enddo
+        enddo
+      endif
+
         do nn=1,nparpi
           if(abs(row(nn)).gt.ftol) then
             nar=nar+1