#! /bin/sh
set -x
# Shell to generate Model 2 in MULLER's paper. Journal of Geophysics 1985.

# output files
wfp=./pressure1
wfr=./radial1
wfz=./vertical1
outf=./info1

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

#source parameters
stype=1			# source type.=1 for explosive, 2 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.75,1.75,1.75,1.75,2.575,2.6,2.8,3.0,3.1,3.35,3.8,4.35,3.7,2.75,4.65,4.2,3.8,4.65,2.5,4.65,3.65,4.7,4.15,4.55,3.2,4.55,3.2,4.35,3.65,4.0,3.35,3.7,3.15,4.4,2.5,4.85,3.65,4.8,3.3,4.5,4.2,2.8,4.5,3.15,4.4,3.9,3.0,4.7,2.85,4.45,3.15,4.5

ct=1.01,1.01,1.01,1.01,1.5,1.55,1.65,1.7,1.73,1.85,2.2,2.65,2.2,1.6,2.7,2.3,1.85,2.75,1.4,2.75,2.1,2.75,2.2,2.7,1.7,2.75,1.7,2.6,2.0,2.25,1.85,2.1,1.75,2.5,1.4,2.75,2.1,2.73,1.8,2.55,2.4,1.65,2.6,1.8,2.02,2.25,1.75,2.65,1.65,2.02,1.8,2.55

ql=100,100,100,100,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000

qt=50,50,50,50,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500

rho=2.05,2.05,2.05,2.05,2.22,2.24,2.28,2.29,2.30,2.35,2.47,2.63,2.26,2.26,2.66,2.5,2.35,2.68,2.19,2.68,2.43,2.07,2.47,2.64,2.29,2.64,2.29,2.61,2.4,2.49,2.35,2.43,2.31,2.57,2.19,2.66,2.43,2.65,2.33,2.59,2.54,2.28,2.61,2.33,2.41,2.49,2.31,2.63,2.28,2.41,2.33,2.59

t=0.001,0.001,0.028,0.010,0.06,0.1,0.1,0.1,0.02,0.025,0.025,0.14,0.01,0.02,0.035,0.07,0.01,0.02,0.01,0.015,0.01,0.05,0.005,0.005,0.005,0.035,0.005,0.055,0.005,0.005,0.005,0.005,0.005,0.01,0.005,0.04,0.005,0.005,0.005,0.03,0.015,0.005,0.01,0.005,0.015,0.01,0.005,0.02,0.005,0.01,0.005,1.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=0		# reducing velocity, equal to maxcimum cl if set to 0
fpeak=35		# peak frequency for ricker or akb wavelet
nf=512			# number of frequencies in the output data
nt=512			# number of time samples in output traces
tlag=0.0		# time lag in seconds
nx=20			# number of output traces


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

# plotting parameters:
agc=0			# apply agc to the plot 
wagc=0.5		# agc window length (sec) 
pbal=0			# trace balancing by rms value
qbal=0			# trace balancing by clip value
qclip=0.9		# clip by quantile on absolute value of trace
perc=98			# perc to apply to plot
xcur=0.6
va=0
wt=1
fill=0


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

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 \
	wt=$wt va=$va xcur=$xcur fill=$fill label1="time (s)" \
	label2="Trace Number" &

suxwigb < $wfr title="PSV radial field perc=$perc f=$fpeak" perc=$perc \
	wt=$wt va=$va xcur=$xcur fill=$fill label1="time (s)" \
	label2="Trace Number" xbox=300 &

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

exit
