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