#############################################
#
# . . Post-stack example
#

nz=801
oz=0
dz=37.5
nx=6400
ox=-30000
dx=37.5
Cartpar=a1n=$(nz) a1d=$(dz) a1o=$(oz) a2n=$(nx) a2o=$(ox) a2d=$(dx) 
weipar3= max5=25 squeeze=n 

################
# 
# Receiver Wavefield
#	On the left side
#
Vel.HH: # given in model

rwf2.H: # zerodd.mute.HH included in directory
	< zerodd.mute.HH Window3d max2=5 j2=2 > t0.H
	< t0.H Pad beg2=50 end2=50 end1=717> t1.H
	< t1.H Interp n2out=1024 o2out=-20.1016 d2out=0.0257 maxsize=400 > tt1.H
	< tt1.H Transf wei=y verb=y > t2.H
	< t2.H Window3d $(weipar3) | Window3d> t3.H
	< t3.H Window3d j2=1 min2=0.2 max2=24.5 > t4.H
	< t4.H Pad beg1=3072 | Window3d n1=4046 > t5.H
	< t5.H Pad end1=50 > $@
	rm -f t0.H t1.H tt1.H t2.H t3.H t4.H t5.H

PM.H rays.H PM.RWE.H AAA.H:  rwf2.H Vel.HH $(XBIN)/RWE2D_ZERO.x
	time $(XBIN)/RWE2D_ZERO.x rwf="rwf2.H" vel="Vel.HH"\
		image="$@" Rimage="Rimage.H" verbose=1 \
		verb=0 ntap=50 norm=1 nsx=2 nsz=2 minang=.5 maxang=179.5 \
		ot=0. dt=0.00065 nt=1201 
	Cp Rimage.H PM.RWE.H

CartparA = a1n=21 a1d=60 a1o=0 a2n=26 a2d=1 a2o=-20 
RWE.VEL.v: $(XBIN)/Lin2.x $(XBIN)/CC2RC.x Vel.HH
	$(XBIN)/Lin2.x axis=2 $(CartparA) > xx.H
	$(XBIN)/Lin2.x axis=1 $(CartparA) > zz.H
	Cmplx xx.H zz.H > Rays.H
	< Rays.H Graph    plotcol=1 out=t1.v $(dn) wantaxis=n title=" " \
		min1=-20 max1=5 min2=0 max2=1200
	< Rays.H Transp | Graph plotcol=1 out=t2.v $(dn) wantaxis=n title=" " \
		min1=-20 max1=5 min2=0 max2=1200
	vp_Overlay t1.v t2.v > test.v
	< Vel.HH $(XBIN)/CC2RC.x rays=Rays.HH | Window3d j2=4 j1=2 min2=1.5 > t1.H
	echo "d1=2 o1=0 o2=-20 d2=0.047" >> t1.H
	< AAA.H Window3d n2=536 n1=301 f1=400 | Interp n1out=601 d1out=0.00125 > t2.H
	Math file1="t1.H" file2="t2.H" exp="file2/file1" >> t3.H
	< t3.H Grey newclip=1 bpclip=0 epclip=100 out=t$@ $(dn)  title="b)" polarity=-1 \
	 	label1="Extrapolation Step" \
		label2="Surface take-off location (km)" min2=-20 max2=5 min1=0 max1=1200 \
		labelrot=n
	vp_Overlay t$@ test.v > $@

Single.H Single.RWE.H: rwf2.H $(XBIN)/RWE2D_ZERO.x Vel.HH 
	< rwf2.H Window3d n2=1 min2=2.5 > R.H
	< Vel.HH Pad extend=1 beg2=500 end1=500 > vvv.H
	time $(XBIN)/RWE2D_ZERO.x rwf="R.H" vel="vvv.H"\
		image="$@" Rimage="Rimage.H" verbose=1 \
		verb=0 ntap=75 norm=1 nsx=2 nsz=2 \
		minang=90 maxang=179 \
		ot=0. dt=0.000425 nt=1801  # was 0.00065
	Cp Rimage.H Single.RWE.H

Cartpar = a1n=100 a1d=.025 a1o=0 a2n=100 a2d=0.025 a2o=.8 
Figure4.v: Single.H PM.H Vel.HH PM.RWE.H rays.H Single.RWE.H $(XBIN)/Lin2.x
	< rays.H Window3d j1=40 |Transp | Graph min1=5 max1=21 max2=20 min2=0 out=r1.v\
		 $(dn) plotcol=1 yreverse=y min2=0 wantaxis=n title=" "
	< rays.H Window3d j2=40 f1=23   | Graph min1=5 max1=21 max2=20 min2=0 out=r2.v\
		 $(dn) plotcol=2 yreverse=y min2=0 wantaxis=n title=" "
	vp_Overlay r2.v r1.v > ray.v
	< Vel.HH   Window3d min2=5 max2=21 max1=20 j1=2 j2=2 > t0.H
	< Single.H Window3d min2=5 max2=21 max1=20 j2=2 j1=2 > t1.H
	Math file1=t1.H file2=t0.H exp=file1*300000+file2 |\
		 Grey out=n1.v title=" " label1="Depth (km)"\
		label2="Distance (km)" $(dn)
	vp_Overlay n1.v ray.v > UL.v
	$(XBIN)/Lin2.x axis=2 $(Cartpar) > xx.H
	$(XBIN)/Lin2.x axis=1 $(Cartpar) > zz.H
	Cmplx xx.H zz.H > Rays.H
	< Rays.H Window3d j1=5 |Transp | Graph min1=2.3 max1=3.1 max2=0.78 min2=0 out=r1.v\
		 $(dn) plotcol=1 yreverse=y min2=0 wantaxis=n title=" "
	< Rays.H Window3d j2=5         | Graph min1=2.3 max1=3.1 max2=0.78 min2=0 out=r2.v\
		 $(dn) plotcol=2 yreverse=y min2=0 wantaxis=n title=" "
	vp_Overlay r2.v r1.v > Ray.v
	< Single.RWE.H Window3d  min2=2.3 max2=3.1 max1=0.78 min1=0 j2=2 j1=2 |\
		Grey pclip=99.2 title=" "  label1="Extrapolation Direction" label2="Angle (rad)"\
		out=x1.v $(dn)
	vp_Overlay x1.v Ray.v > UR.v
	< PM.H Window3d min2=5 max2=21 max1=20 j1=2 j2=2 > t1.H
	Math file1=t1.H file2=t0.H exp=file1*5000+file2 |\
		Grey out=BL.v title=" " label1="Depth (km)"\
		label2="Distance (km)" $(dn) pclip=99.95
	< PM.RWE.H Window3d  min2=2.3 max2=3.1 max1=0.78 min1=0  j2=2 j1=2|\
		Grey pclip=99.9 title=" "  label1="Extrapolation Direction" label2="Angle (rad)"\
		out=BR.v $(dn)
	vp_SideBySideAniso UL.v UR.v > top.v
	vp_SideBySideAniso BL.v BR.v > bot.v
	vp_OverUnderAniso top.v bot.v > $@
#	rm -f top.v bot.v UL.v UR.v BL.v BR.v Rays.H xx.H zz.H ray.v

Figure2.v Figure3.v: $(XBIN)/Lin2.x $(XBIN)/CC2RC.x Vel.HH AAA.H
	$(XBIN)/Lin2.x axis=2 $(CartparA) > xx.H
	$(XBIN)/Lin2.x axis=1 $(CartparA) > zz.H
	Cmplx xx.H zz.H > Rays.H
	< Rays.H Graph plotcol=1 out=t1.v $(dn) wantaxis=n title=" " \
		min1=-20 max1=5 min2=0 max2=1200
	< Rays.H Transp | Graph plotcol=1 out=t2.v $(dn) wantaxis=n title=" " \
		min1=-20 max1=5 min2=0 max2=1200
	vp_Overlay t1.v t2.v > test.v
	< Vel.HH $(XBIN)/CC2RC.x rays=Rays.HH > zzz.H
	< zzz.H  Window j2=4 j1=2 > t1.H
	< AAA.H Window j2=4 j1=2 > t2.H
	echo "d1=2 o1=0 " >> t1.H
	echo "d1=2 o1=0 " >> t2.H
	Math file1="t1.H" file2="t2.H" exp="file2*file1" >> t3.H
	< t3.H Grey newclip=1 bpclip=0 epclip=100 out=tt.v $(dn)  title="b)" polarity=-1 \
	 	label1="Extrapolation Step" label2="Surface take-off location (rad)"\
		 min2=1.5707963 max2=3.1415926 min1=0 max1=1200 labelrot=n
	vp_Overlay tt.v test.v > RWE.VEL.vv
	< zerodd.HH Window j2=6 j1=3 min2=-19.5 max2=21.5 > zz.H
	< zz.H Grey pclip=99.75 out=Figure2.v $(dn) labelrot=n\
		label1="Time (s)" label2="Distance (km)" title=" " 
	< Vel.HH Window3d         j2=4 j1=4 min2=-19.5 max2=21.5 max1=35   |\
		Grey newclip=1 bpclip=0 epclip=100 \
		out=R.v $(dn) labelrot=n \
		label1="Depth (km)" label2="Distance (km)" title="a)"
	< Rays.HH Window j2=64 j1=64 |\
		Graph yreverse=y plotcol=1 min1=-19.5 max1=21.5 max2=35 min2=0 out=r1.v \
		$(dn label2=" " label2=" " wantaxis=n title=" "
	< Rays.HH Window j2=64 j1=64 | Transp | \
		Graph yreverse=y plotcol=1 min1=-19.5 max1=21.5 max2=35 min2=0 out=r2.v \
		$(dn)  label2=" " label2=" " wantaxis=n title=" "
	vp_Overlay r1.v r2.v > r.v
	vp_Overlay R.v r.v .>  RR.v
	vp_SideBySideAniso RR.v RWE.VEL.vv > Figure3.v
	rm -f L.v R.v r1.v r2.l RR.v

Xmin=10
Xmax=40
Zmin=0
Zmax=12

loc1=0
loc2=5


Figure5.ps:
	matlab < RC2.m

#Fig/RC2.v:
#	Window3d n3=1 f3=0 j1=10 j2=10 < rays.H > r.H
#	Window3d n3=1 f3=1 j1=10 j2=10 < rays.H > i.H
#	Cmplx r.H i.H > BPrays.H
#	< rays.H Window3d n3=1 f3=0 > x.H
#	< rays.H Window3d n3=1 f3=1 > z.H
#	Cmplx x.H z.H > rayray.H	
#	< BPAIT-vel.HH Window3d min j1=4 j2=2 min2=0 max2=50000 > CV.H
#	echo "d1=0.025 d2=0.025 o1=0 o2=0" >> CV.H
#	< CV.H Grey newclip=1 bpclip=5 epclip=95 out=tt.v $(dn) max1=$(Zmax)\
#		 min2=$(Xmin) max2=$(Xmax) color=j \
#		 labelrot=n label1="Depth (km)" label2="Distance (km)" title=" " \
#		polarity=1 title="a)" labelsz=14 titlesz=18 geophysics=y
#	Math file1="BPrays.H" exp_real="0.001*file1.re+${loc1}" exp_imag="0.001*file1.im"  > t1.H
#	< t1.H  Window3d j2=10 j1=1 f2=3  > ray11.H
#	< ray11.H Graph yreverse=y plotcol=0 min1=$(Xmin) max1=$(Xmax) max2=$(Zmax) min2=$(Zmin) \
#		out=r1.v $(dn)\
#		 title=" " label2=" " label2=" " wantaxis=n plotfat=3
#	Math file1="BPrays.H" exp_real="0.001*file1.re+${loc1}" exp_imag="0.001*file1.im" |\
#			Window3d j1=5 f2=0 | Transp > ray12.H
#	< ray12.H  Graph yreverse=y plotcol=0 min1=$(Xmin) max1=$(Xmax) max2=$(Zmax) min2=$(Zmin)\
#		 out=r2.v $(dn) \
#		title=" " label2=" " label2=" " wantaxis=n plotfat=3
#	vp_Overlay r1.v r2.v > rr.v
#	vp_Overlay tt.v rr.v > a.v
#	< CV.H Grey newclip=1 bpclip=5 epclip=95 out=tt.v $(dn) max1=$(Zmax) min2=$(Xmin) max2=$(Xmax) color=j\
#		 labelrot=n label1="Depth (km)" label2="Distance (km)" title=" " polarity=1 title="c)"\
#		 labelsz=14 titlesz=18 geophysics=y
#	Math file1="BPrays.H" exp_real="0.001*file1.re+${loc2}" exp_imag="0.001*file1.im"  > t1.H
#	< t1.H  Window3d j2=10 j1=1 f2=3  > ray21.H
#	< ray21.H Graph yreverse=y plotcol=0 min1=$(Xmin) max1=$(Xmax) max2=$(Zmax) min2=$(Zmin) out=r1.v $(dn)\
#		 title=" " label2=" " label2=" " wantaxis=n plotfat=3
#	Math file1="BPrays.H" exp_real="0.001*file1.re+${loc2}" exp_imag="0.001*file1.im" |\
#		Window3d j1=5 f2=0 | Transp > ray22.H
#	< ray22.H Graph yreverse=y plotcol=0 min1=$(Xmin) \
#		max1=$(Xmax) max2=$(Zmax) min2=$(Zmin) out=r2.v $(dn) \
#		title=" " label2=" " label2=" " wantaxis=n plotfat=3
#	vp_Overlay r1.v r2.v > rr.v
#	vp_Overlay tt.v rr.v >  c.v
#	$(XBIN)/Lin2.x axis=2 $(CartparA) > xx.H
#	$(XBIN)/Lin2.x axis=1 $(CartparA) > zz.H
#	Cmplx xx.H zz.H > Rays.H
#	< Rays.H Graph    plotcol=0 out=t1.v $(dn) wantaxis=n title=" " \
#		min1=-20 max1=5 min2=0 max2=1200 plotfat=3
#	< Rays.H Transp | Graph plotcol=0 out=t2.v $(dn) wantaxis=n title=" " \
#		min1=-20 max1=5 min2=0 max2=1200 plotfat=3
#	vp_Overlay t1.v t2.v > test.v
#	Math file1="rayray.H" exp_real="0.001*file1.re+${loc1}" exp_imag="0.001*file1.im" >rayshift1.H
#	< CV.H $(XBIN)/CC2RC.x rays=rayshift1.H   > TT1.H
#	echo "d2=1 o1=20 d1=0.004545" >> t1.H
#	< TT1.H Transp plane=12 | Grey newclip=1 bpclip=0 epclip=100 \
#		out=tx.v $(dn) title="b)" color=j polarity=1\
#	 	label1="Extrapolation Step" label2="Surface take-off location (km)"\
#		 labelrot=n labelsz=14 titlesz=18 geophysics=y
#	 vp_Overlay tx.v test.v > b.v
#	Math file1="rayray.H" exp_real="0.001*file1.re+${loc2}" \
#		exp_imag="0.001*file1.im" >rayshift2.H
#	< CV.H $(XBIN)/CC2RC.x rays=rayshift2.H   > TT2.H
#	echo "d2=1 o1=25 d1=0.004545" >> t1.H
#	< TT2.H Transp plane=12 | Grey newclip=1 bpclip=0 \
#		epclip=100 out=tx.v $(dn) title="d)" color=j polarity=1\
#	 	label1="Extrapolation Step" label2="Surface take-off location (km)"\
#		 labelrot=n labelsz=14 titlesz=18 geophysics=y 
#	vp_Overlay tx.v test.v > d.v
#	vp_SideBySideAniso a.v b.v > a1.v
#	vp_SideBySideAniso c.v d.v > a2.v
#	vp_OverUnderAniso a1.v a2.v > RC2.v
#	pstexpen RC2.v RC2.ps invras=y color=y
#	epstopdf RC2.ps
#	mv RC2.pdf ~/Thesis/radcig/Fig/RC2.pdf


clean: jclean