#! /bin/sh

# Author:	Werner M. Heigl
# Usage:	./WedgeModel

# We would like to create a reflectivity model of a wedge, which is
# often used to model subcropping layers against an unconformity.
# The purpose of this is to study how the wavelets interfere when
# the thickness of the subcropping layer decreases. This gives an
# idea about the frequency/bandwidth that is needed if one wants
# to resolve or detect the subcropping layer for a given thickness
# and velocities.

# It is less clear what resolution and detection means for seismic
# images. It is clear though what the two concepts depend on, namely
# anything that affects the interference, i.e. functional form of
# the wavelet and its parameters, the medium geometry and properties

# The strategy now is to break the problem into simple parts, solve
# each individually, and put the pieces together at the end.
# A wedge model has two reflectors; they differ only in their
# reflectivity and inclination. So, we create a flat reflector first,
# then an inclined second, and add the two panels of traces together.

#
#                     velocity v1
# 0                                                xmax
# *------------------------------------------------+
#     *                                            + 
#		   *                velocity v2            +          
#			    *                                  +
#				     *                             +
#						  *                        h
#							   *                   +
#									*              +
#		velocity v3  					 *         +
#											  *    +
#												   *


# set parameters
xmax=20000			# model dimension in feet
h=800				# thickness of subcropping layer at xmax
ntr=200
dx=$(echo "$xmax/$ntr" | bc)
dh=$(echo "scale=3; ($h/$xmax)*$dx" | bc)

v1=6000				# velocities in ft/s
v2=8000				
v3=6000
# vlmo reflects h=200 at 5000 feet offset and v2=8000
vlmo=$(echo "scale=3; $xmax/($h/$v2)" | bc)
echo $vlmo

r1=$(echo "scale=2; ($v2-$v1)/($v2+$v1)" | bc)
r2=$(echo "scale=2; ($v3-$v2)/($v3+$v2)" | bc)
echo $r1
echo $r2

nt=500
dt=0.001
tmax=$(echo "scale=3; $dt*$nt" | bc)
t0=0.2
waveform=ricker1
fpeak=50	# anti-alias limit is 0.2/dt=200

verbose=0


# 1. Make a wavelet
suwaveform type=$waveform dt=$dt fpeak=$fpeak verbose=$verbose |
sunormalize norm=max verbose=$verbose >wavelet.su

# 2. Make a flat reflector
sunull nt=$nt dt=$dt ntr=$ntr verbose=$verbose |
sushw key=offset a=0 b=$dx verbose=$verbose |
suaddevent vel=$v1 t0=$t0 type=flat amp=$r1 verbose=$verbose |
suconv sufile=wavelet.su verbose=$verbose |
suwind tmin=0 tmax=$tmax verbose=$verbose >flat.su
sunormalize <flat.su norm=max | suxwigb clip=1 key=offset & 

# 3. Make an inclined reflector
sunull nt=$nt dt=$dt ntr=$ntr verbose=$verbose |
sushw key=offset a=0 b=$dx verbose=$verbose |
suaddevent vel=$vlmo t0=$t0 type=lmo amp=$r2 verbose=$verbose |
suconv sufile=wavelet.su verbose=$verbose |
suwind tmin=0 tmax=$tmax verbose=$verbose >inclined.su
sunormalize <inclined.su norm=max | suxwigb clip=1 key=offset & 

# 4. Add the two panels
suop2 flat.su inclined.su op=sum verbose=$verbose >sum.su
#sunormalize <sum.su norm=max verbose=$verbose |
suximage <sum.su perc=100 f2=0 d2=$dh \
         label2="Thickness (ft)" label1="Time (s)" \
		 title="$fpeak Hz $waveform Wavelet" \
		 verbose=$verbose & 


exit 0
