#! /bin/sh
# PSDemo --- example script for sufdmod2
# finite-difference modeling: display model and movie
# Author: John Stockwell
# Modified from various XDemos by Chris Liner 2/12/02
# Big changes...
#  1. deleted (n1,n2) params to unif2 (unused)
#     note on use of (dx,dz) in unif2 hardwired model
#     added 2x2 smoothing to velocity model
#  2. uses ps2pdf to generate pdf output
#  3. cleaned up unit labeling inconsistency
#  4. moved shot to x=50 m to show more far offset action
#     moved shot to z=15 m depth (more realistic)
#     moved horiz recs to 5 m depth (more realistic)
#     moved vsp to x=150 m (somewhat realistic)
#  5. changed top boundary to free surface (ghost,multiples)
#  6. output of sufdmod2 saved to file movie.su
#  7. added display of snapshot 56
#  8. added ps displays of hseis and vseis
#     including a double flip of vseis for standard vsp plot
#  9. added 4-panel display of frame 56, vels, hseis, vseis
# 10. added pscube display of finite difference data
# 11. for clarity, su format files now have .su extension
# 12. added lots of comments

# parameters used in programs below
n1=100 d1=5 f1=0.0 label1="Depth (m)"
n2=100 d2=5 f2=0.0 label2="Distance (m)"
xs=50 zs=15 hsz=5 vsx=150 verbose=2
vsfile="vseis.su" ssfile="sseis.su" hsfile="hseis.su"
tmax=.4 mt=5 fpeak=35 fmax=40

# make 5-interface demonstration velocity file
# 1. the model.txt file is hardwired to
#    (nx,dx,nz,dz)=(100,10,100,10)
# 2. model.txt is a text file that can be edited
#    to change the interface geometry
unif2 tfile=model.out

# generate the velocity model v(x,z)
# Note: with the default model.txt, dx and dz are
#       in effect dummy parameters that determine
#       how much of the complete model to use.
#       Setting (dx,dz)=(10,10) uses the entire model,
#       while smaller (dx,dz) gives some portion of
#       the upper left corner.  Using dx not equal
#       dz changes the aspect ratio of the model
#       (e.g., try dx=10 dy=10 to see whole model)
unif2 < model.out nx=$n2 nz=$n1 \
      dx=7 dz=7 \
      method=spline \
  | smooth2 n1=$n1 n2=$n2 r1=2 r2=2 \
  > vel.out

# make image of the velocity model
psimage < vel.out \
  title="Wavespeed Profile" \
  n1=$n1 n2=$n2 d1=$d1 d2=$d2 \
  legend=1 lnice=1 units="m/s" \
  label1="$label1" label2="$label2" \
  > vel.ps
ps2pdf vel.ps

# 1. run finite difference code
# 2. set dummy dt header word for supsimage plotting
# 3. set tracl header word from 1->ntr for suwind trace selection
sufdmod2 <vel.out \
    nz=$n1 dz=$d1 nx=$n2 dx=$d2 verbose=1 \
    fpeak=$fpeak fmax=$fmax \
    xs=$xs zs=$zs hsz=$hsz vsx=$vsx hsfile=$hsfile \
    vsfile=$vsfile ssfile=$ssfile verbose=$verbose \
    tmax=$tmax abs=0,1,1,1 mt=$mt \
  | sushw key=dt,tracl a=4000,1 b=0,1 \
  > movie.su

# Make image of snapshot 56.
# Note: Snapshot 56 corresponds to time step 280
#       because mt=5, five compute steps for every
#       output time sample.  The actual time value
#       associated with time step 280 can be found
#       from the verbose output
suwind < movie.su key=tracl min=5600 max=5700 \
  | supsimage perc=98 \
      title="Wavefield (t=0.28 s; source at x=50 m, z=5 m)" \
      n1=$n1 n2=$n2 d1=$d1 d2=$d2 f2=0 \
      label1="$label1" label2="$label2" \
  > snap.ps
ps2pdf snap.ps

# make image of horizontal data
supsimage < hseis.su perc=98 \
  title="Surface data (dx=5 m, z=5 m)" \
  d2=$d2 f2=0 \
  label1="Time (s)" label2="$label2" \
> hseis.ps
ps2pdf hseis.ps

# make image of vertical (vsp) data
# 1. flip the data to have depth vertical and time horizontal
suflip < vseis.su flip=-1 \
  | suflip flip=3 \
  > tmp.su
# make the image
supsimage < tmp.su perc=98 \
  title="VSP data (x=150 m, dz=5 m)" \
  d1=$d1 f1=0 f2=0 d2=.001 \
  label1="Depth (m)" label2="Time (s)" \
> vseis.ps
ps2pdf vseis.ps

# make composite image of vel model, a snapshot,
# streamer seismic, and vsp seismic
merge4 snap.ps vel.ps hseis.ps vseis.ps \
  > all.ps
ps2pdf all.ps

# make ps movie of all the snap shots
# Note: fixed clip for absolute scaling between frames
supsmovie < movie.su clip=1 \
  title="Acoustic Finite-Differencing" \
  title2="Frame" \
  label1="$label1" label2="$label2" \
  n1=$n1 d1=$d1 f1=$f1 n2=$n2 d2=$d2 f2=$f2 \
  > movie.ps
ps2pdf movie.ps

# cube view of data subset
# 1. the cube is (nx,ny,nt)=(n1,n2,n3)
# 2. skip first 44 snapshots of 100 traces each
# 3. chop off top half of the cube (pscube thinks dir1 is time)
suwind < movie.su \
  key=tracl min=4400 \
  itmin=50 \
  > cube.su
# make the image
# 1. n1 = 100 - 50 = 50 ; f1 = 50*5 = 250 m
# 2. n3 = 81 - 44 = 37 ; f3 = 44 * .004 = 0.176 s
supscube < cube.su \
  n1=50 n2=100 \
  f1=250 d1=5 f2=0 d2=5 f3=0.176 d3=.004 \
  f1num=250 f2num=0 f3num=0.2 \
  label1="$label1" label2="$label2" \
  label3="Time (0.176 - 0.2 s)" \
  title="Cube display of finite difference data" \
  > cube.ps
ps2pdf cube.ps

supscubecontour < cube.su \
  nc=20 \
  n1=50 n2=100 \
  f1=250 d1=5 f2=0 d2=5 f3=0.176 d3=.004 \
  f1num=250 f2num=0 f3num=0.2 \
  label1="$label1" label2="$label2" \
  label3="Time (0.176 - 0.2 s)" \
  title="Cube contour display of finite difference data" \
  > contourcube.ps
ps2pdf contourcube.ps


exit 0


