#! /bin/sh

## Andreas Rueger March 02, 1997
## Compute seimic signatures of reflected 
## wavefields in anisotropic media
## Type the name of the executable to get
## on-line documentation

### output files
### Change these paths or define your own scratch 
### directory


#####  Path to executables
#B=${HOME}/Release/Rays/bin

file=file

## Plotting parameters
hbox=2.2 wbox=5
titlesize=18

## Model dimension
xmax=5.0 zmax=1.2

## Number of phase angles, anglerange 
nangle=201 fangle=0.0 langle=70.0 nplotr=$nangle

## Source coordinates
xs=0.0 zs=0.0

## number of points per plot, need to be increased
## for more complex raypaths
nxz=20

## Wavefront time. If you want to plot more than 
## one wavefront, you need to use the nt-parameter
tw=1.0

## wave mode (0=P 1=SV 2=SH)
mode=0


## elastic parameters
vp1=2.0 vs1=0.5 e1=.1 d1=0.0 g1=0.0 a1=0.0
vp2=4.5 vs2=1.0 e2=0  a2=0.0 dens=2.7

## title=""

## Built the model. This model builder is very similar
## to the one described in CWP-U26
##
elamodel xmin=0 zmin=0 xmax=$xmax zmax=$zmax maxangle=5 \
1	xedge=0,5 \
	zedge=0,0 \
2	xedge=0,5 \
	zedge=1,1 \
3	xedge=0,5 \
	zedge=$zmax,$zmax \
	fill=0.5,0.1,$vp1,$vs1,$e1,$d1,$g1,$a1,$dens \
	fill=0.5,1.1,$vp2,$vs2,$dens \
	kedge=1,2,3 \
	>$file.model

# generate PostScript output
elaps <$file.model  gtri=-1 \
	hbox=$hbox wbox=$wbox gmin=0.1 gmax=0.7 \
	title=$title \
	>$file.model.eps

####   shoot rays. Rays can only be reflected or
####   converted at internal boundaries 
elaray <$file.model >$file.rayend \
        rayfile=$file.ray  wavefile=$file.wave \
        nangle=$nangle mode=$mode \
        nxz=$nxz xs=$xs zs=$zs infofile=info1\
        fangle=$fangle langle=$langle \
	tw=$tw refseq=2,1,1 refseq=1,-1

subset <$file.ray >$file.temp n1=2 n2=$nxz n3=$nangle id3s=5

####   plot the rays with CWP graphics
psgraph <$file.temp \
        n=$nxz nplot=41  \
	hbox=$hbox wbox=$wbox linegray=0.3 \
	title=$title \
	x1beg=0 x1end=$zmax x2beg=0 x2end=$xmax \
        style=seismic linegray=1 >$file.ray.eps


transp <$file.wave >$file.twave n1=1 n2=$nangle nbpe=8

####    plot the wavefronts
psgraph <$file.twave >$file.wave.eps  \
	linewidth=1 mark=8 marksize=2 \
        nplot=1 n=$nangle hbox=$hbox wbox=$wbox \
	title=$title \
        x1beg=0 x1end=$zmax x2beg=0 x2end=$xmax \
        style=seismic

### Here you merge the model/rays/wavefronts 
psmerge in=$file.model.eps in=$file.ray.eps in=$file.wave.eps > $file.all.eps

### You may want to use <gs> instead
echo "You may view the file $file.all.eps with your PostScript Previewer"

## extract traveltime information
raydata<$file.rayend ascci=0 t=1 kend=1

## This tells you how many rays provide traveltimes
n=`cat outpar`
echo 'number of points '$n


psgraph <x_t.data n=$n x1beg=0.9 x1end=2.5 style=seismic\
	  f1num=1 d1num=1 grid1=dot \
          hbox=$hbox wbox=$wbox \
	  label1="time (s)" label2="offset (km)" \
	  titlesize=$titlesize \
          title="P vp=$vp1,vs=$vs1,e=$e1,d=$d1,g=$g1,a=$a1"\
         >$file.t.eps

echo "You may view the file $file.t.eps with your PostScript Previewer"

# Plot seismogram
compon=0  # both components
ng=101    # number of receivers
nt=1001   # number of trace samples
ft=0      # first time
dt=0.004  # sample interval

krecord=1   # recording interface
nameref=2   # reflecting interface
lscale=$ng  # max extrapolation length for ray information
inter=1     # parabolic interpolation
      

elasyn xg=0,5 zg=0,0  ng=$ng <$file.rayend \
    xfile=x_file.bin  zfile=z_file.bin compon=$compon \
    ft=$ft  nt=$nt  dt=$dt \
    krecord=$krecord lscale=$lscale nameref=$nameref \
    inter=$inter

xwigb < x_file.bin n1=$nt  d2=.05 \
      label1="Time (s)" label2="x-position (km)" \
      title="x-component" \
      perc=99.9 &

xwigb <z_file.bin n1=$nt d2=.05 xbox=600 \
      label1="Time (s)" label2="x-position (km)" \
      title="z-component" \
      perc=99.9 &

exit 0
