#!/bin/sh

# This script demonstrates prestack PSPI migration of shot records from
# FD modeling. In the shot records the direct arrivals are muted and the 
# seismograms are filtered to get rid of low frequency content due to
# the tail of the 2D Green's function.
# Finally the shot records are prestack PSPI migrated where the cross-
# correlation is normalized by the incident wave field and then stacked.

WBOX=945
HBOX=360

unif2 < unif2.ascii par=parfile1.par > vfile       # gridded velocities

NX=`grep nx parfile1.par | sed 's/nx=//'`
NZ=`grep nz parfile1.par | sed 's/nz=//'`
DX=`grep dx parfile1.par | sed 's/dx=//'`
DZ=`grep dz parfile1.par | sed 's/dz=//'`

suaddhead < vfile ns=$NZ | \
suximage title="velocity model for FD modeling for Prestack Migration" \
         label1="Depth [m]" label2="Distance [m]" \
         d1=$DZ d2=$DX wbox=$WBOX hbox=$HBOX legend=1 &

NSHOTS=20                  # number of shots
FSHOT=600                  # first horizontal shot position
DSHOT=$((2000 / $NSHOTS))  # distance between shots

VEL=2200     # velocity used for muting direct arrival
T2=0.12      # mute time at shot position

XS=$FSHOT
ISHOT=1

rm -f sections.su
touch sections.su

echo
echo "  generating $NSHOTS shot records, this may take a few minutes ..."
while [ $ISHOT -le $NSHOTS ]; do

  T1=`echo "10k $XS $VEL / $T2 + p" | dc`          # mute time (left)
  T3=`echo "10k 3245 $XS - $VEL / $T2 + p" | dc`   # mute time (right)

  sufdmod2 < vfile par=parfile2.par fmax=50. hsz=10. hsfile=tmp tmax=1. xs=$XS zs=10. > /dev/null

# mute and filter the shot records
  cat tmp | \
  sushw key=fldr,scalco a=$ISHOT,1  | \
  sumute key=gx xmute=0,$XS,3145 tmute=$T1,$T2,$T3 ntaper=100 | \
  sutaper ntr=$NX tr1=20 | \
  subfilt hicut=0 fpasslo=10 fstoplo=5  >> sections.su

  ISHOT=`expr $ISHOT + 1`
  XS=`expr $XS + $DSHOT`

done
echo
echo "  shotrecords done."

suwind < sections.su tmin=0.03 | sushw key=tracl,tracr b=1,1 d=1,1 > data.su
suximage < data.su title="$NSHOTS muted shot records from FD modeling" \
           label1="time [s]" label2="tracl" legend=1 &

# remove upper 2 velocity grid lines corresponding to datum at 10 m depth
NZ1=`expr $NZ - 2`
suaddhead < vfile ns=$NZ | sushw key=dt a=1 | suwind itmin=2 | sustrip > tmp

# transpose the velocity grid
transp < tmp n1=$NZ1 > vel_tr
