diff --git a/seispy/ccp3d.py b/seispy/ccp3d.py index d005908a..3b9ad196 100644 --- a/seispy/ccp3d.py +++ b/seispy/ccp3d.py @@ -1,7 +1,7 @@ import numpy as np from seispy.geo import km2deg, extrema, skm2srad, rad2deg from seispy import distaz -from seispy.rfcorrect import DepModel +from seispy.core.depmodel import DepModel from seispy.setuplog import setuplog from scikits.bootstrap import ci from seispy.ccppara import ccppara, CCPPara diff --git a/seispy/ccpprofile.py b/seispy/ccpprofile.py index 91f3895c..3be8b2fa 100644 --- a/seispy/ccpprofile.py +++ b/seispy/ccpprofile.py @@ -4,7 +4,7 @@ geo2sph, sph2geo from seispy.setuplog import setuplog from seispy.distaz import distaz -from seispy.rfcorrect import DepModel +from seispy.core.depmodel import DepModel from seispy.rf2depth_makedata import Station from seispy.ccppara import ccppara, CCPPara from scikits.bootstrap import ci diff --git a/seispy/core/depmodel.py b/seispy/core/depmodel.py index bc9f5c37..1d3d8f6a 100644 --- a/seispy/core/depmodel.py +++ b/seispy/core/depmodel.py @@ -40,26 +40,26 @@ def _search_vel_file(mode_name): try: raw_model = np.loadtxt(filename) except: - raise IOError("failed while reading velocity model {}".format(filename)) + raise IOError # check cols of file - if raw_model.shape[1]<3 or raw_model.shape[1]>4: - raise ValueError("failed while reading velocity model {}, plz check your format".format(filename)) + if raw_model.shape[1] < 3 or raw_model.shape[1] > 4: + raise ValueError # cal rho if rho is not given in vel files. - if raw_model.shape[1]==3: - model = np.zeros((raw_model.shape[0],4)) - model[:3,:] = raw_model[:,:] - _p, rho = vs2vprho(raw_model[:,2]) - model[3,:] = rho + if raw_model.shape[1] == 3: + model = np.zeros((raw_model.shape[0], 4)) + model[:3, :] = raw_model[:, :] + _p, rho = vs2vprho(raw_model[:, 2]) + model[3, :] = rho return model else: return raw_model + def _layer2grid(dep_range, model): """ - trans model from layer_model 2 layer_grid - grid in between discontinuities are set to 0 + trans model from layer_model to layer_grid leave stuffing and interp 2 discretize @@ -72,94 +72,26 @@ def _layer2grid(dep_range, model): """ - neo_model = np.zeros((len(dep_range),4)).astype(float) - # vp_dep = np.zeros_like(dep_range).astype(float) - # vs_dep = np.zeros_like(dep_range).astype(float) - # rho_dep = np.zeros_like(dep_range).astype(float) + neo_model = np.zeros((len(dep_range), 4)).astype(float) + # vp_dep = np.zeros_like(dep_range).astype(float) + # vs_dep = np.zeros_like(dep_range).astype(float) + # rho_dep = np.zeros_like(dep_range).astype(float) - picks = np.searchsorted(model[:,0],dep_range, side="left") - for _i,_j in enumerate(picks): - neo_model[_i,:]=model[_j,:] + picks = np.searchsorted(model[:, 0], dep_range, side="left") + for _i, _j in enumerate(picks): + neo_model[_i, :] = model[_j, :] # return vp_dep, vs_dep, rho_dep - print(neo_model[:,1]) - return neo_model[:,1], neo_model[:,2], neo_model[:,3] - - -def _descretize(dep_range, model): - """ - - - Parameters - ---------- - dep_range - model - - Returns - ------- - - """ - - - -def discretize(self,raw_depth, raw_vp, raw_vs, raw_rho): - """ - full init Depth model by - interp 1d grid value - - - """ - if self.elevation == 0: - self.depths_elev = self.depths_layer - self.depths_extend = self.depths_layer - else: - dep_append = np.arange(self.depths_layer[-1] + self.dep_val, - self.depths_layer[-1] + self.dep_val + np.floor(self.elevation / self.dep_val + 1), self.dep_val) - self.depths_extend = np.append(self.depths_layer, dep_append) - self.depths_elev = np.append(self.depths_layer, dep_append) - self.elevation - - self.dz = np.append(0, np.diff(self.depths_extend)) - self.thickness = np.append(np.diff(self.depths_extend), 0.) - - self.vp = interp1d(raw_depth, raw_vp, bounds_error=False, - fill_value=raw_vp[0])(self.depths_elev) - self.vs = interp1d(raw_depth, raw_vs, bounds_error=False, - fill_value=raw_vs[0])(self.depths_elev) - self.rho = interp1d(raw_depth, raw_vs, bounds_error=False, - fill_value=raw_rho[0])(self.depths_elev) - - self.R = 6371.0 - self.depths_elev - - -def _diff_elev(elevatioon=0, dep_range): - """ - elevatioon - dep_range - ------------------ - gens: - depths_elev - depths_extend - dz - thickness - vp, vs, rho,R - - - """ - if elevatioon == 0: - depths_elev = dep_range - depths_extend = dep_range - else: - dep_append = np.arange(self.depths[-1] + self.dep_val, self.elevation + self.dep_val, self.dep_val) - self.depths_elev = np.concatenate((self.depths - self.elevation, dep_append)) - - self.depths_extend = np.concatenate((self.depths, dep_append)) - self.dz = np.diff(self.depths_extend, prepend=0) - self.thickness = np.diff(self.depths_extend, append=0) + return neo_model[:, 1], neo_model[:, 2], neo_model[:, 3] +def _intep_mod(model, depths_elev): + vp = interp1d(model[:,0], model[:,1], bounds_error=False, + fill_value=model[0,1])(depths_elev) + vs = interp1d(model[:,0], model[:,2], bounds_error=False, + fill_value=model[0,2])(depths_elev) + rho = interp1d(model[:,0], model[:,3], bounds_error=False, + fill_value=model[0,3])(depths_elev) - depths_extend = np.concatenate((depths, dep_append)) - dz = np.diff(depths_extend, prepend=0) - thickness = np.diff(depths_extend, append=0) - + return vp, vs, rho class DepModel(object): @@ -167,7 +99,16 @@ class DepModel(object): radiu_s is used to call piercing point tps,tpsps or so are used to cal displacement + examples: + >>> model = DepModel(np.array([0, 20.1, 35.1, 100])) + >>> print(model.dz) + [ 0. 20.1 15. 64.9] + >>> print(model.vp) + [5.8 6.5 8.04001059 8.04764706] + >>> print(model.vs) + [3.36 3.75 4.47003177 4.49294118] """ + def __init__(self, dep_range, velmod='iasp91', elevation=0., layer_mod=False): """Class for computing back projection of Ps Ray paths. @@ -185,89 +126,133 @@ def __init__(self, dep_range, velmod='iasp91', elevation=0., layer_mod=False): self.isrho = False self.elevation = elevation self.layer_mod = layer_mod - self.depths_layer = dep_range.astype(float) - self.dep_val = np.average(np.diff(self.depths_layer)) + # dep layer for CCP or other purpose, relative to sea level, dont contain any depth infomation + self.depths = dep_range.astype(float) + self.dep_val = np.average(np.diff(self.depths)) try: self.model_array = _search_vel_file(velmod) - except (IOError,ValueError): - raise IOError + except (IOError, ValueError): + raise IOError(" failed while loading vel model {}".format(velmod)) else: + self._elevation() if layer_mod: self.vp, self.vs, self.rho = \ - _layer2grid(self.depths_layer, self.model_array) + _layer2grid(self.depths, self.model_array) else: self.vp, self.vs, self.rho = \ - _discretize(self.depths_layer, self.model_array) - + _intep_mod(self.model_array, self.depths_elev) @classmethod def read_layer_model(cls, dep_range, h, vp, vs, rho=None, elevation=0): mod = cls(dep_range, velmod=None, layer_mod=True, elevation=elevation) return mod + def _elevation(self): + """ + set all depth related values: + 1. depths_elev: depth array contains layer above sea level(represent by value lower than 0 + 2. depths_extend: depth array contains layer above sea level( evaluate from 0 to original depth + 3. dz: 0, thick_1, thick_2.... + 4. thickness: thick_1, thick_2, ... , 0 + requires elevation, depths_range or other component + >>> model = DepModel(np.array([0, 20.1, 35.1, 100])) + >>> model.depths_elev + array([ 0. , 20.1, 35.1, 100. ]) + >>> model.depths_extend + array([ 0. , 20.1, 35.1, 100. ]) + >>> model.dz + array([ 0. , 20.1, 15. , 64.9]) + >>> model.thickness + array([20.1, 15. , 64.9, 0. ]) + + >>> model = DepModel(np.array([0, 20.1, 35.1, 100]),elevation=10.) + >>> print(model.depths_elev) + [-10. 10.1 25.1 90. 123.33333333] + >>> print(model.depths_extend) + [ 0. 20.1 35.1 100. 133.33333333] + >>> print(model.dz) + [ 0. 20.1 15. 64.9 33.33333333] + >>> print(model.thickness) + [20.1 15. 64.9 33.33333333 0. ] + """ + if self.elevation == 0: + self.depths_elev = self.depths + self.depths_extend = self.depths + else: + + + depths_append = np.arange(self.depths[-1]+self.dep_val, + self.depths[-1]+self.dep_val+np.floor(self.elevation/self.dep_val+1), self.dep_val) + + self.depths_extend = np.append(self.depths, depths_append) + self.depths_elev = np.append(self.depths, depths_append) - self.elevation + + + self.dz = np.append(0, np.diff(self.depths_extend)) + self.thickness = np.append(np.diff(self.depths_extend), 0.) + + self.R = 6371.0 - self.depths_elev + + + + def plot_model(self, show=True): plt.style.use('bmh') if self.isrho: - self.model_fig = plt.figure(figsize=(6,6)) + self.model_fig = plt.figure(figsize=(6, 6)) fignum = 2 else: - self.model_fig = plt.figure(figsize=(4,6)) + self.model_fig = plt.figure(figsize=(4, 6)) fignum = 1 - self.model_ax = self.model_fig.add_subplot(1,fignum,1) - self.model_ax.step(self.vp, self.depths_layer, where='pre', label='Vp') - self.model_ax.step(self.vs, self.depths_layer, where='pre', label='Vs') + self.model_ax = self.model_fig.add_subplot(1, fignum, 1) + self.model_ax.step(self.vp, self.depths, where='pre', label='Vp') + self.model_ax.step(self.vs, self.depths, where='pre', label='Vs') self.model_ax.legend() self.model_ax.set_xlabel('Velocity (km/s)') self.model_ax.set_ylabel('Depth (km)') - self.model_ax.set_ylim([self.depths_layer[0], self.depths_layer[-1]]) + self.model_ax.set_ylim([self.depths[0], self.depths[-1]]) self.model_ax.invert_yaxis() if self.isrho: - self.rho_ax = self.model_fig.add_subplot(1,fignum,2) - self.rho_ax.step(self.rho, self.depths_layer, where='pre', color='C2', label='Density') + self.rho_ax = self.model_fig.add_subplot(1, fignum, 2) + self.rho_ax.step(self.rho, self.depths, where='pre', color='C2', label='Density') self.rho_ax.legend() self.rho_ax.set_xlabel('Density (km/s)') - self.rho_ax.set_ylim([self.depths_layer[0], self.depths_layer[-1]]) + self.rho_ax.set_ylim([self.depths[0], self.depths[-1]]) self.rho_ax.invert_yaxis() if show: plt.show() - - - def tpds(self, rayps, raypp, sphere=True): if sphere: radius = self.R else: radius = 6371. tps = np.cumsum((np.sqrt((radius / self.vs) ** 2 - rayps ** 2) - - np.sqrt((radius / self.vp) ** 2 - raypp ** 2)) * + np.sqrt((radius / self.vp) ** 2 - raypp ** 2)) * (self.dz / radius)) return tps - def tpppds(self, rayps, raypp, sphere=True): if sphere: radius = self.R else: radius = 6371. tps = np.cumsum((np.sqrt((radius / self.vs) ** 2 - rayps ** 2) + - np.sqrt((radius / self.vp) ** 2 - raypp ** 2)) * + np.sqrt((radius / self.vp) ** 2 - raypp ** 2)) * (self.dz / radius)) return tps - def tpspds(self, rayps, sphere=True): if sphere: radius = self.R else: radius = 6371. - tps = np.cumsum(2*np.sqrt((radius / self.vs) ** 2 - rayps ** 2)* + tps = np.cumsum(2 * np.sqrt((radius / self.vs) ** 2 - rayps ** 2) * (self.dz / radius)) return tps - def radius_s(self, rayp, phase='P', sphere=True): """ calculate piercing point, P for Sp and S for Ps @@ -279,6 +264,14 @@ def radius_s(self, rayp, phase='P', sphere=True): Returns ------- + >>> model = DepModel(np.array([0, 20.1, 35.1, 100])) + >>> model.dz + array([ 0. , 20.1, 15. , 64.9]) + >>> model.R + array([6371. , 6350.9, 6335.9, 6271. ]) + >>> model.radius_s(1.2,phase="S", sphere=False)*111.2 + array([0. , 0.0002478 , 0.00046823, 0.00142685]) + """ if phase == 'P': @@ -290,9 +283,9 @@ def radius_s(self, rayp, phase='P', sphere=True): else: radius = 6371. hor_dis = np.cumsum((self.dz / radius) / np.sqrt((1. / (rayp ** 2. * (radius / vel) ** -2)) - 1)) + #hor_dis = np.sqrt((1. / (rayp ** 2. * (radius / vel) ** -2)) - 1) return hor_dis - def raylength(self, rayp, phase='P', sphere=True): if phase == 'P': vel = self.vp @@ -305,8 +298,7 @@ def raylength(self, rayp, phase='P', sphere=True): raylen = (self.dz * radius) / (np.sqrt(((radius / self.vs) ** 2) - (rayp ** 2)) * vel) return raylen + if __name__ == "__main__": - depth = np.array([0, 20.1, 35.1, 100]) - model = DepModel(depth) - rayp = np.arange(0.04, 0.09, 0.01) - print(model.vp) + import doctest + doctest.testmod() diff --git a/seispy/modcreator.py b/seispy/modcreator.py index 27b605f1..dcab38d8 100644 --- a/seispy/modcreator.py +++ b/seispy/modcreator.py @@ -1,6 +1,6 @@ import numpy as np from scipy.interpolate import griddata, interp1d -from seispy.rfcorrect import DepModel +from seispy.core.depmodel import DepModel from seispy.setuplog import setuplog import argparse diff --git a/seispy/rfcorrect.py b/seispy/rfcorrect.py index 78ef09f9..24a4473d 100644 --- a/seispy/rfcorrect.py +++ b/seispy/rfcorrect.py @@ -3,6 +3,8 @@ from scipy.interpolate import interp1d, interpn from scipy.signal import resample from os.path import dirname, join, exists, isfile, abspath + +import seispy.core.depmodel from seispy.geo import skm2srad, sdeg2skm, rad2deg, latlon_from, \ asind, tand, srad2skm, km2deg from seispy.psrayp import get_psrayp @@ -11,7 +13,8 @@ from seispy.harmonics import Harmonics from seispy.plotRT import plotrt as _plotrt from seispy.plotR import plotr as _plotr -from seispy.utils import DepModel, Mod3DPerturbation, scalar_instance +from seispy.utils import Mod3DPerturbation, scalar_instance +from seispy.core.depmodel import DepModel import warnings import glob @@ -557,7 +560,9 @@ def psrf2depth(stadatar, YAxisRange, velmod='iasp91', srayp=None, normalize='sin return ps_rfdepth, endindex, x_s, x_p -def xps_tps_map(dep_mod, srayp, prayp, is_raylen=False, sphere=True, phase=1): +def xps_tps_map(dep_mod: seispy.core.depmodel.DepModel + , srayp, prayp, + is_raylen=False, sphere=True, phase=1): """Calculate horizontal distance and time difference at depths :param dep_mod: 1D velocity model class @@ -611,7 +616,7 @@ def xps_tps_map(dep_mod, srayp, prayp, is_raylen=False, sphere=True, phase=1): tps = interp1d(dep_mod.depths_elev, tps, bounds_error=False, fill_value=(np.nan, tps[-1]))(dep_mod.depths) if is_raylen: raylength_s = interp1d(dep_mod.depths_elev, raylength_s, bounds_error=False, fill_value=(np.nan, raylength_s[-1]))(dep_mod.depths) - raylength_p = interp1d(dep_mod.depths_elev, raylength_p, bounds_error=False, fill_value=(np.nan, raylength_p[-1]))(dep_mod.depths) + raylength_p = interp1d(dep_mod.depths_elev, raylength_p, bounds_error=False, fill_value=(np.nan, raylength_p[-1]))(dep_mod.depths) if is_raylen: return tps, x_s, x_p, raylength_s, raylength_p else: diff --git a/seispy/utils.py b/seispy/utils.py index d5a82e19..38f312ea 100644 --- a/seispy/utils.py +++ b/seispy/utils.py @@ -1,11 +1,13 @@ from os.path import join, dirname, exists +import array + from scipy.io import loadmat from matplotlib.colors import ListedColormap -import matplotlib.pyplot as plt import numpy as np from scipy.interpolate import interp1d, interpn import pandas as pd -import array + +#from seispy.core.depmodel import DepModel def vs2vprho(vs): @@ -44,201 +46,6 @@ def read_rfdep(path): raise FileNotFoundError('Cannot open file of {}'.format(path)) -def _from_layer_model(dep_range, h, vp, vs, rho=None): - dep = 0 - vp_dep = np.zeros_like(dep_range).astype(float) - vs_dep = np.zeros_like(dep_range).astype(float) - rho_dep = np.zeros_like(dep_range).astype(float) - for i, layer in enumerate(h): - if (layer == 0 and i == len(h)-1) or \ - (dep+layer < dep_range[-1] and i == len(h)-1): - idx = np.where(dep_range>=dep)[0] - else: - idx = np.where((dep_range >= dep) & (dep_range < dep+layer))[0] - if idx.size == 0: - raise ValueError('The thickness of layer {} less than the depth interval'.format(i+1)) - vp_dep[idx] = vp[i] - vs_dep[idx] = vs[i] - if rho is not None: - rho_dep[idx] = rho[i] - dep += layer - if dep > dep_range[-1]: - break - return vp_dep, vs_dep, rho_dep - - -class DepModel(object): - def __init__(self, dep_range, velmod='iasp91', elevation=0., layer_mod=False): - """Class for computing back projection of Ps Ray paths. - - Parameters - ---------- - dep_range : numpy.ndarray - Depth range for conversion - velmod : str, optional - Text file of 1D velocity model with first 3 columns of depth/thickness, Vp and Vs, - by default 'iasp91' - elevation : float, optional - Elevation in km, by default 0.0 - """ - self.isrho = False - self.elevation = elevation - self.layer_mod = layer_mod - self.depths = dep_range.astype(float) - try: - self.model_array = np.loadtxt(self.from_file(velmod)) - except (ValueError, TypeError): - return - else: - self.read_model_file() - self.discretize() - - def read_model_file(self): - self.dep_val = np.average(np.diff(self.depths)) - if self.layer_mod: - self.depthsraw = self.depths - if self.model_array.shape[1] == 4: - self.isrho = True - rho_array = self.model_array[:, 3] - else: - rho_array = None - self.vpraw, self.vsraw, self.rhoraw = _from_layer_model(self.depths, - self.model_array[:, 0], - self.model_array[:, 1], - self.model_array[:, 2], - rho=rho_array - ) - else: - self.depthsraw = self.model_array[:, 0] - self.vpraw = self.model_array[:, 1] - self.vsraw = self.model_array[:, 2] - if self.model_array.shape[1] == 4: - self.isrho = True - self.rhoraw = self.model_array[:, 3] - - @classmethod - def read_layer_model(cls, dep_range, h, vp, vs, rho=None, elevation=0): - mod = cls(dep_range, velmod=None, layer_mod=True, elevation=elevation) - if rho is not None: - mod.isrho = True - mod.depthsraw = mod.depths - mod.vpraw, mod.vsraw, mod.rhoraw = _from_layer_model(mod.depths, h, vp, vs, rho=rho) - mod.discretize() - return mod - - def plot_model(self, show=True): - plt.style.use('bmh') - if self.isrho: - self.model_fig = plt.figure(figsize=(6,6)) - fignum = 2 - else: - self.model_fig = plt.figure(figsize=(4,6)) - fignum = 1 - self.model_ax = self.model_fig.add_subplot(1,fignum,1) - self.model_ax.step(self.vp, self.depths, where='pre', label='Vp') - self.model_ax.step(self.vs, self.depths, where='pre', label='Vs') - self.model_ax.legend() - self.model_ax.set_xlabel('Velocity (km/s)') - self.model_ax.set_ylabel('Depth (km)') - self.model_ax.set_ylim([self.depths[0], self.depths[-1]]) - self.model_ax.invert_yaxis() - if self.isrho: - self.rho_ax = self.model_fig.add_subplot(1,fignum,2) - self.rho_ax.step(self.rho, self.depths, where='pre', color='C2', label='Density') - self.rho_ax.legend() - self.rho_ax.set_xlabel('Density (km/s)') - self.rho_ax.set_ylim([self.depths[0], self.depths[-1]]) - self.rho_ax.invert_yaxis() - if show: - plt.show() - - def discretize(self): - if self.elevation == 0: - self.depths_elev = self.depths - self.depths_extend = self.depths - else: - dep_append = np.arange(self.depths[-1]+self.dep_val, - self.depths[-1]+self.dep_val+np.floor(self.elevation/self.dep_val+1), self.dep_val) - self.depths_extend = np.append(self.depths, dep_append) - self.depths_elev = np.append(self.depths, dep_append) - self.elevation - self.dz = np.append(0, np.diff(self.depths_extend)) - self.thickness = np.append(np.diff(self.depths_extend), 0.) - self.vp = interp1d(self.depthsraw, self.vpraw, bounds_error=False, - fill_value=self.vpraw[0])(self.depths_elev) - self.vs = interp1d(self.depthsraw, self.vsraw, bounds_error=False, - fill_value=self.vsraw[0])(self.depths_elev) - if self.isrho: - self.rho = interp1d(self.depthsraw, self.rhoraw, bounds_error=False, - fill_value=self.rhoraw[0])(self.depths_elev) - else: - _, self.rho = vs2vprho(self.vs) - self.R = 6371.0 - self.depths_elev - - def from_file(self, mode_name): - if not isinstance(mode_name, str): - raise TypeError('velmod should be in str type') - if exists(mode_name): - filename = mode_name - elif exists(join(dirname(__file__), 'data', mode_name.lower()+'.vel')): - filename = join(dirname(__file__), 'data', mode_name.lower()+'.vel') - else: - raise ValueError('No such file of velocity model') - return filename - - def tpds(self, rayps, raypp, sphere=True): - if sphere: - radius = self.R - else: - radius = 6371. - tps = np.cumsum((np.sqrt((radius / self.vs) ** 2 - rayps ** 2) - - np.sqrt((radius / self.vp) ** 2 - raypp ** 2)) * - (self.dz / radius)) - return tps - - def tpppds(self, rayps, raypp, sphere=True): - if sphere: - radius = self.R - else: - radius = 6371. - tps = np.cumsum((np.sqrt((radius / self.vs) ** 2 - rayps ** 2) + - np.sqrt((radius / self.vp) ** 2 - raypp ** 2)) * - (self.dz / radius)) - return tps - - def tpspds(self, rayps, sphere=True): - if sphere: - radius = self.R - else: - radius = 6371. - tps = np.cumsum(2*np.sqrt((radius / self.vs) ** 2 - rayps ** 2)* - (self.dz / radius)) - return tps - - def radius_s(self, rayp, phase='P', sphere=True): - if phase == 'P': - vel = self.vp - else: - vel = self.vs - if sphere: - radius = self.R - else: - radius = 6371. - hor_dis = np.cumsum((self.dz / radius) / np.sqrt((1. / (rayp ** 2. * (radius / vel) ** -2)) - 1)) - return hor_dis - - def raylength(self, rayp, phase='P', sphere=True): - if phase == 'P': - vel = self.vp - else: - vel = self.vs - if sphere: - radius = self.R - else: - radius = 6371. - raylen = (self.dz * radius) / (np.sqrt(((radius / self.vs) ** 2) - (rayp ** 2)) * vel) - return raylen - - class Mod3DPerturbation: def __init__(self, modpath, YAxisRange, velmod='iasp91'): dep_mod = DepModel(YAxisRange, velmod=velmod) diff --git a/test/test_case05.py b/test/test_case05.py index 7c050d62..f4f3744c 100644 --- a/test/test_case05.py +++ b/test/test_case05.py @@ -1,4 +1,4 @@ -from seispy.utils import DepModel +from seispy.core.depmodel import DepModel from seispy.seisfwd import SynSeis from seispy.rfcorrect import RFStation import pytest @@ -8,6 +8,8 @@ def test_sub01(): depth = np.array([0, 20.1, 35.1, 100]) model = DepModel(depth) + print(model.vp) + print(model.dz) rayp = np.arange(0.04, 0.09, 0.01) ss = SynSeis(model, rayp, 0.1, 2400) ss.run_fwd()