# . . Test mig for EAGE Salt BODY
#
RCUT = n1=250 f1=1 
nt=4096
dt=0.004
ot=0.
rwepar2  = nref=6 verbose=1 forward=1 kinematic=1. min_region_pct=.2 del_dist=.1 verb=0 \
	niter_lloyd=20 ntap=10 nloop=3
Crwepar2 = nref=3 verbose=1 forward=1 kinematic=1. min_region_pct=.2 del_dist=.1 verb=0 \
	niter_lloyd=20 ntap=10 nloop=3
nx3=1400
ox3=2000
dx3=10

ny3=1
oy3=0
dy3=50

nz3=550
dz3=10
oz3=0
CCout3 = azn=$(nz3) azo=$(oz3) azd=$(dz3)\
	 axn=$(nx3) axo=$(ox3) axd=$(dx3)\
	 ayn=$(ny3) ayo=$(oy3) ayd=$(dy3)\
	 nsx=2 nsz=2 nsy=0 adj=1 fold=1
###############
#
# . . Cart ray params

CPnz=550
CPoz=0
CPdz=10
CPnx=512
CPox=3000
CPdx=20
HWTgraph = Graph crowd=${crowd} wantaxis=n title=" " label1=" " label2=" "

CPcartpar= a1n=$(CPnz) a1d=$(CPdz) a1o=$(CPoz) \
	   a2n=$(CPnx) a2d=$(CPdx) a2o=$(CPox) 
GREYCUT4   = min2=3500 max2=13000 min1=0 max1=5000
GREYCUT   = min2=3000 max2=13000 min1=0 max1=5000
GRAPHCUT  = min1=2500 max1=13500 min2=0 max2=5500
GRAPHCUT2 = min1=2500 max1=13500 min2=0 max2=5500
GREYCUT2  = min2=4000 max2=12500 min1=0 max1=5500
GREYCUT3  = min2=4500 max2=13250 min1=0 max1=5500
c4=0.85
XB = ./Bin/
##########################################
#
# . . Velocity and Rays
#

%.x:
	@-cd ./Src; make default; cd ..
PSeikonal.H:
	< PS.vel1.H Interp d1out=10 d2out=10 > t1.H
	< t1.H FMeikonal order=3 zshot=0 yshot=9000 xshot=0 |\
		Smooth rect1=3 rect2=3 repeat=5 > $@

PS.rays.HH PS.rays2.HH:
	matlab < ./matlab/Cnt2.m 

# . . The velocity model slice we're interested in
PS.vel1.H: 
	Math file1="Eslow.HH" exp="1./file1" > t$@  
	< t$@ Pad beg2=200 end2=200 extend=1 end1=175 > $@
	echo "o3=0." >> $@
	\rm -f t$@ 

# . . Rays for Riemannian
PS.rays.H: 
	< Imesh.20.HH Real | Interp d1out=4.95 d2out=9.95 n2out=512 n1out=512 type=1 > xx.H
	< Imesh.20.HH Imag | Interp d1out=4.95 d2out=9.95 n2out=512 n1out=512 type=1 > zz.H
	Cmplx xx.H zz.H > $@
# . . velocity on Riemannian Rays
PS.vel.H: PS.vel1.H $(XBIN)/CC2RC.x            PS.rays.H
	< PS.vel1.H Smooth rect2=2 rect1=2 repeat=1 | $(XBIN)/CC2RC.x adj=0 rays=PS.rays.H >$@ 

# . . Rays for Cartesian
CPS.ray.H: $(XBIN)/Lin2.x
	$B/Lin2.x axis=2 $(CPcartpar) >xx.H
	$B/Lin2.x axis=1 $(CPcartpar) >zz.H
	Cmplx xx.H zz.H > CPS.ray2.H

# . . Velocity on Cart Rays
CPS.vel.H: PS.vel1.H $(XBIN)/CC2RC.x            CPS.ray2.H
	<  PS.vel1.H $(XBIN)/CC2RC.x adj=0 rays=CPS.ray2.H > t1.H
	< t1.H	Smooth rect1=2 rect2=2 repeat=2 tridiag=1 >$@

###################################################
#
# . . Data
#
# . . Riemannian Migration
PS.data.H Data.H:               PS.vel.H
	Get n2 n3 o2 o3 d2 d3 < PS.vel.H >par.P
	Spike n1=$(nt) o1=$(ot) d1=$(dt) par=par.P \
		nsp=7 k1=50,125,200,275,350,425,500 \
		k2=0,0,0,0,0,0,0 k3=0,0,0,0,0,0,0 \
		mag=.001,.001,.001,.001,.001,.001,.001 > Data.H
	<  Data.H Transf f_min=1 f_min1=2 f_max1=25 f_max=30 wei=y verb=y >t$@
	<  t$@ Transp plane=35 > PS.data.H
	\rm t$@ 
# . . Cartesian Migration
CPS.data.H CData.H:             CPS.vel.H
	Get n2 n3 o2 o3 d2 d3 < CPS.vel.H >par.P
	Spike n1=$(nt) o1=$(ot) d1=$(dt) par=par.P \
		nsp=7 k1=63,138,213,288,363,438,513\
		 k2=301,301,301,301,301,301,301 k3=0,0,0,0,0,0,0\
		 mag=.001,.001,.001,.001,.001,.001,.001> CData.H
	<  CData.H Transf f_min=1 f_min1=2 f_max1=25 f_max=30 wei=y verb=y >t$@
	<  t$@ Transp plane=35 > CPS.data.H
	\rm t$@ 

#####################################################
#
# . . Migration
#
rweparx2 = nref=10 verbose=1 forward=1 kinematic=1. min_region_pct=.1\
		 del_dist=.05 verb=0 niter_lloyd=20
PS.mig.H PS.CC.H: PS.data.H PS.vel.H PS.rays.H $(B)/RWE3D_source.x              
	< PS.rays.H Real > xx.H
	< PS.rays.H Imag > zz.H
	Math file1=xx.H exp='0.*file1' > yy.H
	Cat xx.H yy.H zz.H axis=4 > rayray.H
	<   PS.data.H $(XBIN)/RWE3D_source.x rays=rayray.H \
		vel=PS.vel.H swf=PS.data.H $(rweparx2)> PS.mig.H
	< PS.mig.H  $(XBIN)/C2Rsinc3D.x $(CCout3) rays=rayray.H fold=0> PS.CC.H

CPS.mig.H CPS.CC.H: CPS.data.H  CPS.vel.H CPS.ray2.H $(XBIN)/RWE3D_source.x                   
	< CPS.ray2.H Real > xx.H
	< CPS.ray2.H Imag > zz.H
	Math file1=xx.H exp='0.*file1' > yy.H
	Cat xx.H yy.H zz.H axis=4 > rayray.H
	< CPS.data.H $(XBIN)/RWE3D_source.x rays=rayray.H \
		vel=CPS.vel.H swf=CPS.data.H $(Crwepar2)> CPS.mig.H
	< CPS.mig.H  $(XBIN)/C2Rsinc3D.x $(CCout3) rays=rayray.H > CPS.CC.H 


# . . Finite-difference comparison
sou.H:
	Wavelet n1=500 d1=0.001 \
		wavelet=ricker0         \
		domain=time tdelay=0.05  \
		fund=9.0 fhigh=15.    | \
		Scale rscale=1000 >$@
	echo "label1=t_"  >>$@


TW = $(XB)/aWFvNB.x
ACMOD = $(TW) equord=2 verb=y nt=3200 dt=0.000666 jsnap=150 \
		nsz=1 dsz=1. osz=0 \
		nsx=1 dsx=1. osx=9000
ro.H:     CPS.vel.H
	< CPS.vel.H Scale rscale=0.0 | Add add=1000.0 >$@

Twoway.H:      sou.H $(XB)/aWFvNB.x  ro.H    CPS.vel.H 
	time < sou.H $(ACMOD) ro=ro.H vp=CPS.vel.H wfld=$@ > shot.H

Good.H:   Twoway.H
	< Twoway.H Window3d n3=1 f3=2| Scale rscale=1.0 > t1.H
	< Twoway.H Window3d n3=1 f3=5| Scale rscale=1.25> t2.H
	< Twoway.H Window3d n3=1 f3=8| Scale rscale=1.5 > t3.H
	< Twoway.H Window3d n3=1 f3=11| Scale rscale=1.75> t4.H
	< Twoway.H Window3d n3=1 f3=14| Scale rscale=2.0 > t5.H
	< Twoway.H Window3d n3=1 f3=17| Scale rscale=2.25> t6.H
	< Twoway.H Window3d n3=1 f3=20| Scale rscale=2.5 > t7.H
	Add t1.H t2.H t3.H t4.H t5.H t6.H t7.H > $@

########################################################
#
# . . Output Figures
#
SPexamp.v: PS.vel1.H PS.rays.H CPS.rays.H PS.mig.H PS.CC.H
	< PS.vel1.H Scale rscale=0.001 > t1.H
	< t1.H Grey newclip=1 bpclip=0 epclip=100 out=t1.v $(dn) color=j\
		label1=" " label2=" " $(GREYCUT) \
		title=" " min1=0 wantaxis=n crowd=$(crowd) crowd=$(c4)
	< PS.rays.H Window3d n4=2 f4=0 j4=2 > tt.H
	< tt.H Real | Window3d n3=1 f3=0 > xx.H 
	< tt.H Imag | Window3d n3=1 f3=1 > zz.H
	Cmplx xx.H zz.H > r.H
	<  r.H Window3d j1=1 j2=15 > kk1.H
	< kk1.H ${HWTgraph} $(dn) out=jr.v plotcol=2 wantaxis=n label1=" " \
		label2=" " $(GRAPHCUT) crowd=$(c4)
	<  r.H Window3d j1=8 j2=1 | Transp > kk2.H
	< kk2.H ${HWTgraph} label2=" " wantaxis=n $(GRAPHCUT) \
		label1=" " labelrot=n $(dn) out=jw.v plotcol=1 crowd=$(c4)
	vp_Overlay jr.v jw.v > t2.v
	vp_Overlay t1.v t2.v > tl.v
	< PS.vel1.H Scale rscale=0.001 > t1.H
	echo "d1=1 d1=1 o1=1 o2=1" >> t1.H
	< t1.H Grey newclip=1 bpclip=0 epclip=100 out=t1.v $(dn) color=j\
		label1=" " label2=" "  \
		title=" " min1=0 wantaxis=n crowd=$(crowd) crowd=$(c4)
	< CPS.rays.H Window3d n4=2 f4=0 j4=2 > tt.H
	< tt.H Window3d n3=1 f3=0 > xx.H
	< tt.H Window3d n3=1 f3=1 > zz.H
	Cmplx xx.H zz.H > r.H
	<  r.H Window3d j1=1 j2=15 > kk1.H
	< kk1.H ${HWTgraph} $(dn) out=jr.v plotcol=2 wantaxis=n label1=" " \
		label2=" " $(GRAPHCUT2) crowd=$(c4)
	<  r.H Window3d j1=8 j2=1 | Transp > kk2.H
	< kk2.H ${HWTgraph} label2=" " wantaxis=n $(GRAPHCUT2) \
		label1=" " labelrot=n $(dn) out=jw.v plotcol=1 crowd=$(c4)
	vp_Overlay jr.v jw.v > t2.v
	vp_Overlay t1.v t2.v > tr.v
	< PS.mig.H Grey out=br.v wantaxis=n pclip=99.5 $(dn) title=" " crowd=$(c4)
	<  PS.CC.H Grey out=bl.v wantaxis=n pclip=99.5 $(dn) $(GREYCUT) title=" " crowd=$(c4)
	vp_SideBySideAniso tl.v tr.v > top.v
	vp_SideBySideAniso bl.v br.v > bottom.v
	vp_OverUnderAniso top.v bottom.v > $@
	\rm -f t1.H t1.v tt.H xx.H zz.H r.H kk1.H kk2.H jw.v jr.vt2.v tl.v t2.v \
		tr.v bl.v br.v top.v bottom.v


SPcompare.v: PS.CC.H PS.vel1.H CPS.CC.H Good.H CPS.vel.H
# . . Riemannian
	< PS.CC.H   Window3d min2=3000 max2=13000 max1=5000 > t1.H
	< PS.vel1.H Window3d min2=3000 max2=13000 max1=5000 > tt.H
	< tt.H Interp d1out=10 d2out=10 n1out=500 n2out=1000> t2.H
	Math file1="t1.H" file2="t2.H" exp="80000*file1+file2" > tt.H
	< tt.H Grey newclip=1 bpclip=0.025 epclip=99 out=a.v $(dn) \
		label1=" " label2=" " $(GREYCUT) title=" " polarity=-1 \
		wantaxis=n crowd=0.85
# . . Cartesian
	< CPS.CC.H   Window3d min2=3000 max2=13000 max1=5000 > t1.H
	<  PS.vel1.H Window3d min2=3000 max2=13000 max1=5000 > tt.H
	< tt.H Interp d1out=10 d2out=10 n1out=500 n2out=1000> t2.H
	Math file1="t1.H" file2="t2.H" exp="1000000*file1+file2" > tt.H
	< tt.H Grey newclip=1 bpclip=0.025 epclip=99 out=b.v $(dn) \
		label1=" " label2=" " $(GREYCUT) title=" " polarity=-1 \
		wantaxis=n crowd=0.85
# . . Finite-difference
	< Good.H    Window3d min2=3500 max2=13000 max1=5000 > t1.H
	< CPS.vel.H Window3d min2=3500 max2=13000 max1=5000 > t2.H
	Math file1="t1.H" file2="t2.H" exp="500000*file1+file2" > tt.H
	< tt.H Grey newclip=1 bpclip=0 epclip=99.95 out=c.v $(dn) \
		label1="Depth (m) " label2="Distance (m)" $(GREYCUT4)\
		 title="Finite-Difference Generated Wavefield" polarity=-1 \
		wantaxis=y crowd=0.8
	vp_OverUnderAniso c.v b.v b.v> $@	
	\rm -f t1.H t2.H tt.H a.v b.v c.v