def trav(thk,v,dist,ns): # n=0 means surface n is the bottom of the layer
v = 1/(v)**2
lmax = ns
jo = len(thk)
t0,p0 = find2(dist,lmax,v,ns,thk)
for lmax in range(ns+1,jo):
t,p = find2(dist,lmax,v,ns,thk)
if (t<t0):
t0 = t
p0 = p
return t0.real,p0.real
def find2(dist,lmax,v,ns,thk):
p1 = 1e-20j
pm = 9999
for i in range(lmax):
pm = np.min((pm,v[lmax-1]))
p2 = np.sqrt(pm)+p1
p0 = 0.5*(p1+p2)
while((p2-p1)>1e-6):
dtdp0 = dtdp(dist,p0,ns,lmax,v,thk)
if np.abs(dtdp0)<1e-10:
pm =np.min((pm,v[lmax]))
if dtdp0>0:
p1 = p0
else:
p2 = p0
p0 = 0.5*(p1+p2)
pm = np.sqrt(pm)
if (lmax<len(v)-1):
if (lmax>ns) and pm<np.real(p0):
p0 = pm
t0 = taup(p0,dist,ns,lmax,v,thk)
return t0,p0
def taup(p,dist,ns,lmax,v,thk):
tau = p*dist
pp = p*p
for i in range(ns):
tau = tau+np.sqrt(v[i]-pp)*thk[i]
for i in range(ns,lmax):
tau = tau+np.sqrt(v[i]-pp)*thk[i]*2
return tau
def dtdp(dist,p,ns,lmax,v,thk):
pp = p*p
dtdp=0
for i in range(ns):
dtdp = dtdp-thk[i]/np.sqrt(v[i]-pp)
for i in range(ns,lmax):
dtdp = dtdp-2*thk[i]/np.sqrt(v[i]-pp)
dtdp = dist+p*dtdp
return dtdp
# print('fine')
model = np.array(([5.43, 0.00],
[5.80, 5.00],
[5.90, 8.00],
[5.95,11.00],
[6.00,14.00],
[6.09,20.00],
[6.27,30.00],
[6.50,45.00],
[8.00,55.00]))
thick = np.hstack((0,np.diff(model[:,1])))
vp = model[:,0]
vs = vp/1.73
tp = np.array([trav(thick,vp,dist,2)[0] for j,dist in enumerate(dists)])
ts = np.array([trav(thick,vs,dist,2)[0] for j,dist in enumerate(dists)])
trav(thick,vp,dist,2)中的2代表第三个depth
def trav_dep(model,mode,dist,dep):
if dep in model[:,1]:
print('yes')
thick = np.hstack((0,np.diff(model[:,1])))
vp = model[:,0]
vs = vp/1.73
ind_0 = np.where(model[:,1]==dep)[0][0]
ind_1 = ind_0+1
t, p = trav(thick, vp, dist, ind_1)
else:
ind_0 = np.where(model[:,1]>dep)[0][0]
ind_1 = ind_0+1
model_ = np.zeros((len(model)+1,2))
model_[:ind_0] = model[:ind_0]
model_[ind_0] = np.array([(model[ind_0-1][0]+model[ind_1-1][0])/2, dep])
model_[ind_0] = np.array([(model[ind_0-1][0]+model[ind_0-1][0])/2, dep])
model_[ind_0+1:] = model[ind_0:]
# model = np.concatenate((model[ind_0],[(model[ind_0][0]+model[ind_1][0])/2, dep],model[ind_1:]))
thick = np.hstack((0,np.diff(model_[:,1])))
vp = model_[:,0]
vs = vp/1.73
if mode=='p':
t, p = trav(thick, vp, dist, ind_1)
elif mode=='s':
t, p = trav(thick, vs, dist, ind_1)
return t, p
Updated
