rpsv4

su 的源代码库

#! /bin/sh

set -x
# Computes a PSV seismogram for Mallick-Frazier model 2. Geophysics, 1987

# output files
wfp=./pressure4
wfr=./radial4
wfz=./vertical4
outf=./info4

# input parameters
tsec=14.336		# length of computed trace (seconds)
fs=0.07			# sampling paramter 0.07<fs<0.12
decay=50.0		# decay factor to avoid time series wraparound
flt=0			# earth flattening flag
lsource=3		# layer on top of which the source is located
nlayers=19		# number of layers to process
nw=256			# number of frequencies to process for ref. modeling
nor=1			# number of receivers
dt=0.028		# time sampling interval (secs)
p2w=1.0			# maximum ray parameter (s/km)
lobs=2			# layers on top of which the receivers are located
bx=1.0			# beginning range (in km)
fx=14			# final range (offset) (in km)
dx=1.0			# range increment (in km)
pw1=0.0			# begin of low end Hanning window for ray parameter s/km
pw2=0.0			# end of low end Hanning window for ray parameter s/km
pw3=0.0			# begin of high end Hanning window for ray parameter 
pw4=0.0			# end of high end Hanning window for ray parameter s/km

#source parameters
stype=2			# source type.=2 for explosive, 1 for fault plane
h1=0.0			# vertical linear part of the source
h2=0.0			# horizontal linear part of the source
m0=1.0			# seismic moment
m1=1.0			# [1][1] component of the moment tensor (stpe=1)
m2=0.0			# [1][2]=[2][1] component of the moment tensor (stype=1)
m3=1.0			# [2][2] component of the moment tensor (stype=1)
delta=0.0		# dip in degrees (stype=2)
lambda=0.0		# rake in degrees (stype=2)
phis=0.0		# fault plane azimuth in degrees (stype=2)
phi=0.0			# receiver location azimuth in degrees (stype=2)

#parameters needed only in layer interpolation is required
nlint=0			# number of times layer interpolation is required
nintlayers=		# array[nlint] of number of layers to interpolate
intlayers=		# array[nlint] of layers on top of which to interpolate
intlayth=		# array[nlint] of layer thicknesses to interpolate

# parameters required only if random velocity layers are requested
rand=0			# =1 to request random velocity layers
nrand_layers=0		# maximum number of random layers allowed
layer=0			# layer on top of which the random velocity layers
			# are inserted
zlayer=0		# thickness of random layers
sdcl=0			# standard deviation for compressional velocities
sdct=0			# standard deviation for shear velocities

# parameters required only if qoption is requested
qopt=0			# =1 to request the q-option
layern=0		# layer on top of which the q-option is invoked	
wrefp=1.0		# reference frequency for compressional velocities
wrefs=1.0		# reference frequency for shear velocities
epsp=0.001		# reference amplitude for compressional velocities
epss=0.001		# reference amplitude for shear velocities
sigp=0.1		# xxxx for comporessional velocities
sigs=0.1		# xxxx for shear velocities

# key input parameters, (can be input from files, clfile, ctfile, etc) 
cl=1.5,2.2,3.0,3.8,4.2,4.4,4.6,4.85,5.05,5.2,5.5,5.75,6.0,6.25,6.5,6.75,6.9,6.5,8.2
ct=0.35,0.35,0.35,2.2,2.4,2.6,2.75,2.85,3.0,3.1,3.2,3.3,3.4,3.5,3.6,3.7,2.75,3.4,4.7
ql=250,250,250,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,600
qt=700,700,700,700,700,700,700,750,750,750,750,750,750,750,750,750,750,750,850
rho=1.4,1.55,1.7,2.0,2.1,2.2,2.25,2.3,2.35,2.4,2.5,2.6,2.7,2.7,2.8,2.9,2.9,2.6,3.2
t=0.1,0.1,0.1,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.2,0.2,0.2,0.2,0.2,0.25,0.25,0.1,2.0

# output parameters
wtype=1			# =1 for PSV, =2 for SH
wfield=2		# =1 for displacement, 2 for velocity, 3 for accelera
vsp=0			# =1 for vsp, 0 for surface seismogram
win=0			# =1 for windowed frequencies
wavelet_type=2		# =1 for spike, =2 for ricker1, =3 for ricker2,
			# =4 for an akb wavelet
verbose=3		# =0 for no processing information, =1 for processing
			# information in the screen, =2 for processing 
			# information to a file (outf) =3 for both
red_vel=8.2		# reducing velocity, equal to maxcimum cl if set to 0
fpeak=25		# peak frequency for ricker or akb wavelet
nf=256			# number of frequencies in the output data
nt=512			# number of time samples in output traces
tlag=3.0		# time lag in seconds
nx=40			# number of output traces


# output filter parameters (can be input from files, filtypefile,dbpofile, etc)
nfilters=		# number of filters to apply to output traces
filters_type=		# =1 for high cut filter, =2 for low-cut filter
			# =3 for notch filter
filters_phase=		# =0 for zero phase filter, =1 for minimum phase filter
dbpo=			# DB/octave for filter slopes
f1=			# frequency to start filter action
f2=			# high end frequency filter

# plotting parameters:
agc=0			# apply agc to the plot 
wagc=0.5		# agc window length (sec) 
pbal=1			# trace balancing by rms value
qbal=0			# trace balancing by clip value
qclip=1.0		# clip by quantile on absolute value of trace
perc=99			# perc to apply to plot


##############################################################################

sureflpsvsh fs=$fs decay=$decay flt=$flt lsource=$lsource nw=$nw nor=$nor \
	tsec=$tsec dt=$dt lobs=$lobs bx=$bx fx=$fx dx=$dx nx=$nx \
	pw1=$pw1 pw2=$pw2 pw3=$pw3 pw4=$pw4 p2w=$p2w \
	stype=$stype h1=$h1 h2=$h2 m0=$m0 m1=$m1 m2=$m2 m3=$m3 delta=$delta \
	lambda=$lambda phis=$phis phi=$phi \
	rand=$rand nrand_layers=$nrand_layers layer=$layer zlayer=$zlayer \
	sdcl=$sdcl sdct=$sdct \
	qopt=$qopt layern=$layern wrefp=$wrefp wrefs=$wrefs epsp=$epsp \
	epss=$epss sigp=$sigp sigs=$sigs \
	cl=$cl ct=$ct ql=$ql qt=$qt rho=$rho t=$t nlayers=$nlayers \
	wtype=$wtype wfield=$wfield red_vel=$red_vel nf=$nf win=$win vsp=$vsp \
	wavelet_type=$wavelet_type verbose=$verbose nfilters=$nfilters \
	filters_type=$filters_type dbpo=$dbpo f1=$f1 f2=$f2 nt=$nt tlag=$tlag \
	fpeak=$fpeak filters_phase=$filters_phase \
	wfp=$wfp wfr=$wfr wfz=$wfz outf=$outf

# plot the results
suxwigb < $wfp title="####PSV pressure field perc=$perc f=$fpeak" perc=$perc \
	label1="time (s)" label2="Trace Number" &

suxwigb < $wfr title="####PSV radial field perc=$perc f=$fpeak" perc=$perc \
	label1="time (s)" label2="Trace Number" xbox=600 &

suxwigb < $wfz title="####PSV vertical field perc=$perc f=$fpeak" perc=$perc \
	label1="time (s)" label2="Trace Number" xbox=300 ybox=300 &

exit 0