#############################################
#
# . . BP Velocity Benchmark
#
nz=801
oz=0
dz=37.5
nx=8000
ox=-80000
dx=37.5
Cartpar=a1n=$(nz) a1d=$(dz) a1o=$(oz) a2n=$(nx) a2o=$(ox) a2d=$(dx) 
VPAR = n1=1201 d1=10 o1=0 n2=2401 d2=10 o2=-12000 vc=2000 n3=1 d3=20 o3=0
bppar = flo=1.0 fhi=32.5
BPweipar= min5=5 n5=283 j5=1 squeeze=n
sidecut=50
nfreq=250
XBIN = ./Bin/
dn = >/dev/null
RESDIR=./Fig/


###########
#
# . . Source Codes
haha: 
	@-cd Src/; ${MAKE} default ; cd ..

%.x: haha

######################
#
# . . Ray Coords
MARMwindow.H:
	< /data/marmousi/marmvel.H Window3d j1=2 j2=2 > $@

MARMsmooth.H: MARMwindow.H
	< $< Smooth rect1=10 rect2=10 repeat=5 > $@

MARMrays.H: MARMsmooth.H
	Hwt2d  < $< xsou=4000 zsou=0 oT=0. nT=700 dT=0.0023 oG=-90. nG=360 dG=0.5  >$@ 

MARMrays.v: MARMsmooth.H MARMrays.H
	< MARMsmooth.H Grey color=j newclip=1 bpclip=0 epclip=100 min1=0 max1=3000 \
		min2=0 max2=8000 out=t1.v $(dn) title=" " label1=" " label2=" " wantaxis=n
	< MARMrays.H Window3d j1=1 j2=1 | Graph min2=0 max2=3000 min1=0 max1=8000 plotcol=1 out=t2.v \
		label1=" " $(dn) wantaxis=n yreverse=y title=" " 
	< MARMrays.H Transp | Window3d j2=24 j1=1 | Graph min2=0 max2=3000 min1=0 max1=8000 plotcol=2 \
		out=t3.v $(dn) wantaxis=n yreverse=y title=" "
	vp_Overlay t1.v t2.v t3.v > $@
	rm -f t1.v t2.v t3.v

#################################################
#
# . . BP velocity model
BPAIT = /data/bpait/
MVEL.HH: 
	< $(BPAIT)/vel.H Window3d j1=2> $@

# . . Smooth velocity model for ray tracing
SmoothMVel.H: MVEL.HH
	< $< Smooth rect1=25 rect2=25 repeat=25 tridiag=y maxmem=100000000000> $@

# . . Ray tracing
Mrays.H tMrays.H:  SmoothMVel.H
	< SmoothMVel.H Hwt2d xsou=$(xshot) zsou=0 nT=900 dT=0.005 oT=0.01 \
		oG=-85 nG=2048 dG=0.08304836345872007816 > t$@
	< t$@ Real > rr.H
	< t$@ Imag > ii.H
	Cat axis=3 space=n rr.H ii.H > $@
	echo "o2=-85 d2=0.08304836345872007816" >> $@

RVel.HH: MVEL.HH tMrays.H
	< MVEL.HH $(XBIN)/CC2RC.x rays=tMrays.H > $@

rwf.H:
	Wavelet n1=2048 n2=1 n3=1 wavelet=ricker2 fund=25 dt=0.004 ot=-0.12 tdelay=1.12 > t1a.H
	Wavelet n1=2048 n2=1 n3=1 wavelet=ricker2 fund=25 dt=0.004 ot=-0.12 tdelay=2.12 > t1b.H
	Wavelet n1=2048 n2=1 n3=1 wavelet=ricker2 fund=25 dt=0.004 ot=-0.12 tdelay=3.12 > t1c.H
	Wavelet n1=2048 n2=1 n3=1 wavelet=ricker2 fund=25 dt=0.004 ot=-0.12 tdelay=4.12 > t1d.H
	Add t1a.H t1b.H t1c.H t1d.H > t1.H
	< t1.H Transf f_min=2 f_min1=5. f_max1=45 f_max=50 > t2.H
	< t2.H Pad beg1=1023 end1=1024 | Window3d n2=1 > $@
	echo "o1=-85 d1=0.0830078125 o3=xshot" >> $@
	rm -f t1a.H t1b.H t1c.H t1d.H t1.H t2.H 

vcut.H:   MVEL.HH
	< MVEL.HH Window3d min2=$(xmin) max2=$(xmax) max3=$(zmax)> $@

rho.H: MVEL.HH  vcut.H
	Math file1="vcut.H" exp="0.*file1+1." > $@

TwoWay.H: vcut.H rho.H
	Wavelet n1=10000 n2=1 n3=1 wavelet=ricker2 fund=10 dt=0.004 ot=-0.05 tdelay=0.05 > t1.H
	< t1.H $(XBIN)/aWFv.x vp="vcut.H" verb=1 jsnap=4000 nt=180000 dt=0.000025 \
		nsx=1 dsx=1. osx=$(xshot) nsz=1 dsz=1. osz=50. ro="rho.H" wfld="wfld.H" bndz=10 bndx=10> $@

wfldsnap%.v: TwoWay.H
	< wfld.H Window3d n3=1 f3=$* | Grey pclip=100 out=$@ $(dn)

wfldsnap.H:
	< wfld.H Window3d n3=1 f3=10 > t1.H
	< wfld.H Window3d n3=1 f3=20 > t2.H
	< wfld.H Window3d n3=1 f3=30 > t3.H
	< wfld.H Window3d n3=1 f3=40 > t4.H
	Add t1.H t2.H t3.H t4.H > $@
	rm -f t1.H t2.H t3.H t4.H

TwoWay2.H: vcut.H rho.H
	Fdmod < vcut.H nxs=1 oxs=$(xshot) dxs=1. nzs=1. ozs=0. dzs=1. tmax=4 mint=0 \
		dfile=rho.H  nx=1 xs=$(xshot) nz=1 zs=0. > $@ 

xshot=27500	
xmin=19000
xmax=36000
zmax=12000

EllipticImage.H: rwf.H $(XBIN)/RWE2DMODEL.x MVEL.HH
	Cp rwf.H rwf2.H
	echo "o1=-5110 d1=5 o3=$(xshot)" >> rwf2.H
	$(XBIN)/RWE2D_MODEL.x rwf="rwf2.H" vel="MVEL.HH" image="$@" \
		verbose=1 verb=1 eps=0.0001 nsx=3 nsz=2 nh=1 method=2 \
		oe3=0. de3=0.002 ne3=1000 minang=0.5 maxang=179.5 nxtap=50
	< $@ Grey min2=$(xmin) max2=$(xmax) max1=$(zmax) min1=0 pclip=99.75 out=EllipticImage.v $(dn)

CartesianImage.H:	rwf.H $(XBIN)/RWE2D_MODEL.x MVEL.HH
	Cp rwf.H rwf2.H
	echo "o1=-10220 d1=10 o3=$(xshot)" >> rwf2.H
	$(XBIN)/RWE2D_MODEL.x rwf="rwf2.H" vel="MVEL.HH" image="$@" \
		verbose=1 verb=1 eps=0.0001 nsx=3 nsz=2 nh=1 method=1 \
		oe3=0. de3=20. ne3=600 minang=0.5 maxang=179.5 nxtap=50
	< $@ Grey min2=$(xmin) max2=$(xmax) max1=$(zmax) min1=0  pclip=99.75 out=CartesianImage.v $(dn)


# . . RWE image
CartImage.H: rwf.H RVel.HH Mrays.H $(XBIN)/RWE2D_source.x 
	Cp rwf.H rwf2.H
	echo "o3=$(xshot)" >> rwf2.H
	$(XBIN)/RWE2D_source.x swf="rwf2.H" rwf="rwf2.H" vel="RVel.HH" \
		rays="Mrays.H" Rimage="RWEimage.H"  \
		min_region_pct=1. del_dist=1. verb=0 niter_lloyd=25\
		ntap=10 nloop=4 nref=256 forward=1 norm=1 \
		xmin=0 xmax=67425 zmin=0 zmax=11937.5 dxx=12.5 dzz=12.5 nsx=3 nsz=3 image="$@"
	< $@ Grey min2=$(xmin) max2=$(xmax) max1=$(zmax) min1=0  pclip=99.75 out=CartImage.v $(dn)

# . . Four panel comparison
Fig/MCOMP.v: wfldsnap.H vcut.H CartesianImage.H CartImage.H EllipticImage.H MVEL.HH
# . . Two way
	Math file1="wfldsnap.H" file2="vcut.H" exp="file1/0.26+0.1*(file2-3000)/1500" > vv1.H
	< vv1.H  Grey min2=$(xmin) max2=$(xmax) max1=$(zmax) min1=0 newclip=1 bpclip=0 \
		epclip=99.75 out=aa.v $(dn) title="a) Two-way Modeling" label1="Depth (m)" label2="Distance (m)"\
		labelsz=13  titlesz=16
# . . Cartesian
	Math file1="CartesianImage.H" file2="MVEL.HH" exp="file1/206+0.1*(file2-3000)/1500" > vv2.H
	< vv2.H Grey min2=$(xmin) max2=$(xmax) max1=$(zmax) min1=0 newclip=1 bpclip=0 \
		epclip=99.75 out=bb.v $(dn) title="b) Cartesian" label1="Depth (m)" label2="Distance (m)"\
		labelsz=13  titlesz=16
# . . Ray-based image
	Math file1="CartImage.H" file2="MVEL.HH" exp="file1/1300000+0.1*(file2-3000)/1500" > vv3.H
	< vv3.H      Grey min2=$(xmin) max2=$(xmax) max1=$(zmax) min1=0 newclip=1 bpclip=0 \
		epclip=99.75 out=cc.v $(dn) title="c) Ray-based RWE" label1="Depth (m)" label2="Distance (m)"\
		labelsz=13  titlesz=16
# . . Elliptic Image
	Math file1="EllipticImage.H" file2="MVEL.HH" exp="file1/1110+0.1*(file2-3000)/1500" > vv4.H
	< vv4.H  Grey min2=$(xmin) max2=$(xmax) max1=$(zmax) min1=0 newclip=1 bpclip=0 \
		epclip=99.75 out=dd.v $(dn) title="d) Analytic RWE" label1="Depth (m)" label2="Distance (m)"\
		labelsz=13  titlesz=16
	vp_SideBySideAniso aa.v bb.v > top.v
	vp_SideBySideAniso cc.v dd.v > bot.v
	vp_OverUnderAniso top.v bot.v > $@
	\rm -f aa.v bb.v cc.v dd.v top.v bot.v
	pstexpen ./Fig/MCOMP.v ./Fig/MCOMP.ps fat=1 fatmult=1.5
	epstopdf ./Fig/MCOMP.ps

# . . Images
UL.v: MVEL.HH tMrays.H
	< tMrays.H          Window3d j2=20 | Graph plotcol=1 min1=$(xmin) max1=$(xmax) max2=$(zmax) min2=0 \
		wantaxis=n out=t1.v $(dn) yreverse=y title=" " wantaxis=n min2=0
	< tMrays.H Transp | Window3d j2=30 | Graph plotcol=2 min1=$(xmin) max1=$(xmax) max2=$(zmax) min2=0  \
		wantaxis=n out=t2.v $(dn) yreverse=y title=" " wantaxis=n min2=0
	< MVEL.HH Grey min2=$(xmin) max2=$(xmax) max1=$(zmax) min1=0 \
			out=t3.v $(dn) title="a)" label1="Depth (m)" label2="Distance (m)"\
		labelsz=13  titlesz=16
	vp_Overlay t3.v t1.v t2.v > $@
	rm -f t3.v t1.v t2.v

ULC.v: MVEL.HH tMrays.H
	< tMrays.H          Window3d j2=20 | Graph plotcol=1 min1=$(xmin) max1=$(xmax) max2=$(zmax) min2=0 \
		wantaxis=n out=t1.v $(dn) yreverse=y title=" " wantaxis=n min2=0
	< tMrays.H Transp | Window3d j2=30 | Graph plotcol=2 min1=$(xmin) max1=$(xmax) max2=$(zmax) min2=0  \
		wantaxis=n out=t2.v $(dn) yreverse=y title=" " wantaxis=n min2=0
	< MVEL.HH Grey color=j min2=$(xmin) max2=$(xmax) max1=$(zmax) min1=0 \
			out=t3.v $(dn) title="a)" label1="Depth (m)" label2="Distance (m)"\
		labelsz=13  titlesz=16 newclip=1 bpclip=0 epclip=100
	vp_Overlay t3.v t1.v t2.v > $@
	rm -f t3.v t1.v t2.v

LL.v: CartImage.H  tMrays.H
	< CartImage.H Grey pclip=99 min1=0  min2=$(xmin) max2=$(xmax) max1=$(zmax) out=t1.v \
		label1="Depth (m)" label2="Distance (m)" title="d)" \
		labelsz=13  titlesz=16 $(dn)
	< tMrays.H Window3d f2=6 j2=100 | Graph plotcol=1 out=t2.v $(dn) \
		wantaxis=n title=" " yreverse=y min1=$(xmin) max1=$(xmax) max2=$(zmax) min2=0 
	< tMrays.H Transp | Window3d j2=200 | Graph plotcol=2 out=t3.v $(dn) \
		wantaxis=n title=" " yreverse=y min2=0. min1=$(xmin) max1=$(xmax) max2=$(zmax) 
	vp_Overlay t1.v t2.v t3.v > $@
	rm -f t1.v t2.v t3.v

CPnxx=51
CPoxx=-100
CPdxx=4
CPnzx=26
CPozx=0
CPdzx=0.2 

CPcartparx= a1n=$(CPnzx) a1d=$(CPdzx) a1o=$(CPozx) \
            a2n=$(CPnxx) a2d=$(CPdxx) a2o=$(CPoxx) 
AMAX=85
UR.v: RVel.HH $(XBIN)/Lin2.x 
	$(XBIN)/Lin2.x axis=2 $(CPcartparx) > xx.H
	$(XBIN)/Lin2.x axis=1 $(CPcartparx) > zz.H
	Cmplx xx.H zz.H > r.H
	< RVel.HH Grey min2=-$(AMAX) max2=$(AMAX) min1=0 max1=4.5 out=t1.v $(dn) \
		label1="Travel Time (s)" \
		label2="Shooting Angle (deg)" title="b)"\
		labelsz=13 titlesz=16
	< r.H Graph min1=-$(AMAX) max1=$(AMAX) min2=0 max2=4.5 out=t2.v $(dn) \
		wantaxis=n title=" " plotcol=1 yreverse=n
	< r.H Transp | Graph plotcol=2 min1=-$(AMAX) max1=$(AMAX) min2=0 max2=4.5 out=t3.v $(dn) \
		wantaxis=n title=" " yreverse=n
	vp_Overlay t1.v t2.v t3.v > $@
	rm -f t1.v t2.v t3.v xx.H zz.H r.H 

URC.v: RVel.HH $(XBIN)/Lin2.x 
	$(XBIN)/Lin2.x axis=2 $(CPcartparx) > xx.H
	$(XBIN)/Lin2.x axis=1 $(CPcartparx) > zz.H
	Cmplx xx.H zz.H > r.H
	< RVel.HH Grey min2=-$(AMAX) max2=$(AMAX) min1=0 max1=4.5 out=t1.v $(dn) \
		label1="Travel Time (s)" \
		label2="Shooting Angle (deg)" title="b)"\
		labelsz=13 titlesz=16 color=j newclip=1 bpclip=0 epclip=100
	< r.H Graph min1=-$(AMAX) max1=$(AMAX) min2=0 max2=4.5 out=t2.v $(dn) \
		wantaxis=n title=" " plotcol=1 yreverse=n
	< r.H Transp | Graph plotcol=2 min1=-$(AMAX) max1=$(AMAX) min2=0 max2=4.5 out=t3.v $(dn) \
		wantaxis=n title=" " yreverse=n
	vp_Overlay t1.v t2.v t3.v > $@
	rm -f xx.H zz.H r.H t1.v t2.v t3.v

CPnx=11
CPox=-100
CPdx=20
CPnz=11
CPoz=-0.5
CPdz=0.5

CPcartpar= 	a1n=$(CPnz) a1d=$(CPdz) a1o=$(CPoz) \
		a2n=$(CPnx) a2d=$(CPdx) a2o=$(CPox) 

LR.v: RVel.HH  $(XBIN)/Lin2.x RWEimage.H
	$(XBIN)/Lin2.x axis=2 $(CPcartpar) > xx.H
	$(XBIN)/Lin2.x axis=1 $(CPcartpar) > zz.H
	Cmplx xx.H zz.H > r.H
	< RWEimage.H Grey min2=-$(AMAX) max2=$(AMAX) min1=0. max1=4.5 out=t1.v $(dn)\
		label1="Travel Time (s)" label2="Shooting Angle (deg)" title="c)" \
		labelsz=13  titlesz=16
	< r.H Graph min1=-$(AMAX) max1=$(AMAX) min2=0.0 max2=4.5 out=t2.v $(dn) \
		wantaxis=n plotcol=1 yreverse=n title=" "
	< r.H Transp | Graph min1=-$(AMAX) max1=$(AMAX) min2=0.0 max2=4.5 out=t3.v $(dn)\
		wantaxis=n plotcol=2 yreverse=n title=" "
	vp_Overlay t1.v t2.v t3.v > $@
	rm -f xx.H zz.H t1.v t2.v t3.v r.H

Fig/RWEexample.v: UL.v UR.v LL.v LR.v
	vp_SideBySideAniso UL.v UR.v > top.v
	vp_SideBySideAniso LL.v LR.v > bot.v
	vp_OverUnderAniso top.v bot.v > $@
	pstexpen $@ ./Fig/RWEexample.ps fat=1 fatmult=1.5
	epstopdf ./Fig/RWEexample.ps
	rm -f top.v bot.v UL.v UR.v LL.v LR.v

RWEexamplecolor.v: ULC.v URC.v LL.v LR.v
	vp_SideBySideAniso ULC.v URC.v > top.v
	vp_SideBySideAniso LL.v LR.v > bot.v
	vp_OverUnderAniso top.v bot.v > $@