include ${YANGLIB}/makefileY.top

#Test and compare with Yaxun's example



#Four vel model, constant up 1100, constant down 900
#High Vel anomaly 1100, low vel anomaly 900

B = ../bin
P = ../par
Y = ~yang/from-other/yaxun_DSO
strip = ~yang/Research/codelib/bin/StripOutDummyAxis.x
tool = ~yang/Research/codelib/bin

toolTags =   tool1=tool1.H tool2=tool2.H 
dump = dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
		dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H \
				dump9=dump9.H dump10=dump10.H dump11=dump11.H dump12=dump12.H 

modExtend = Pad beg1=50 end1=50 
npad = 50

########################## Try beam focusing ####################, in this case, it is
## is transmission case, make bigger difference, 25 percent difference
##Create Data pad=50
n1 := 200
n2 := 501

vel.t0.H: 
	Vel n1=${n1} n2=${n2} n3=1 o3=0 d3=1 o1=-500 o2=-2500 vc=2000 d1=10 d2=10 > $@
vel.tc1.H: 
	Vel n1=${n1} n2=${n2} n3=1 o3=0 d3=1 o1=-500 o2=-2500 vc=2500 d1=10 d2=10 > $@
vel.tc2.H: 
	Vel n1=${n1} n2=${n2} n3=1 o3=0 d3=1 o1=-500 o2=-2500 vc=1600 d1=10 d2=10 > $@
#small perturbation %1
vel.tc3.H: 
	Vel n1=${n1} n2=${n2} n3=1 o3=0 d3=1 o1=-500 o2=-2500 vc=1980 d1=10 d2=10 > $@
vel.tc4.H: 
	Vel n1=${n1} n2=${n2} n3=1 o3=0 d3=1 o1=-500 o2=-2500 vc=2020 d1=10 d2=10 > $@

#middle perturbation %10
vel.tc5.H: 
	Vel n1=${n1} n2=${n2} n3=1 o3=0 d3=1 o1=-500 o2=-2500 vc=2200 d1=10 d2=10 > $@
vel.tc6.H: 
	Vel n1=${n1} n2=${n2} n3=1 o3=0 d3=1 o1=-500 o2=-2500 vc=1800 d1=10 d2=10 > $@

#anomaly up/down
vel.ta1.H: $B/lens_twin.x 
	$< nx=401 nz=200 ox=-2000 dx=10 oz=-500 dz=10 x11=0 x13=495 x21=4000 x23=445 \
		width1=300 height1=150 \
		alpha1=-0.25 alpha2=0 zref=800 v0=2000	\
		refl=tref0.H bvel=t.H | Pad beg2=50 end2=50 extend=1 > $@
#small error anomaly
vel.ta5.H: $B/lens_twin.x 
	$< nx=401 nz=200 ox=-2000 dx=10 oz=-500 dz=10 x11=0 x13=495 x21=4000 x23=445 \
		width1=400 height1=200 \
		alpha1=-0.05 alpha2=0 zref=800 v0=2000	\
		refl=tref0.H bvel=t.H | Pad beg2=50 end2=50 extend=1 > $@

#one anomaly, at asymmetric location
vel.ta6.H: $B/lens_twin.x 
	$< nx=401 nz=200 ox=-2000 dx=10 oz=-500 dz=10 x11=-1000 x13=250 x21=4000 x23=445 \
		width1=500 height1=300 \
		alpha1=-0.15 alpha2=0 zref=800 v0=2000	\
		refl=tref0.H bvel=t.H | Pad beg2=50 end2=50 extend=1 > $@
		
vel.ta2.H: $B/lens_twin.x
	$< nx=401 nz=200 ox=-2000 dx=10 oz=-500 dz=10 x11=0 x13=495 x21=4000 x23=845 \
		width1=300 height1=150 \
		alpha1=+0.20 alpha2=0 zref=800 v0=2000	\
		refl=tref0.H bvel=t.H | Pad beg2=50 end2=50 extend=1 > $@
#with anomaly and wrong background const velocity
vel.ta4.H: $B/lens_twin.x
	$< nx=401 nz=200 ox=-2000 dx=10 oz=-500 dz=10 x11=200 x13=300 x21=4000 x23=845 \
		width1=500 height1=250 \
		alpha1=+0.20 alpha2=0 zref=800 v0=2100	\
		refl=tref0.H bvel=t.H | Pad beg2=50 end2=50 extend=1 > $@

#twin anomaly model
vel.ta3.H: $B/lens_twin.x
	$< nx=401 nz=200 ox=-2000 dx=10 oz=-500 dz=10 x11=-1200 x13=300 x21=1000 x23=500 \
		width1=300 height1=150 width2=1000 height=300 \
		alpha1=+0.20 alpha2=-0.15 zref=800 v0=2000	\
		refl=tref0.H bvel=t.H | Pad beg2=50 end2=50 extend=1 > $@

#Bigger shape
vel.ta7.H: $B/lens_twin.x
	$< nx=401 nz=200 ox=-2000 dx=10 oz=-500 dz=10 x11=0 x13=400 x21=1000 x23=500 \
		width1=1000 height1=300 width2=1000 height=300 \
		alpha1=-0.50 alpha2=0 zref=800 v0=2000	\
		refl=tref0.H bvel=t.H | Pad beg2=50 end2=50 extend=1 > $@
vel.ta8.H: $B/lens_twin.x
	$< nx=401 nz=200 ox=-2000 dx=10 oz=-500 dz=10 x11=0 x13=400 x21=1000 x23=500 \
		width1=1000 height1=300 width2=1000 height=300 \
		alpha1=+0.33 alpha2=0 zref=800 v0=2000	\
		refl=tref0.H bvel=t.H | Pad beg2=50 end2=50 extend=1 > $@

dv.t%.H: vel.t%.H
	Math file1=$< file2=vel.t0.H exp='file1-file2' >$@
wavY = wavY.tr.H
pad := Pad beg1=50 end1=50 

wavY.tr.H:
	Wavelet n1=250 o1=0 d1=0.004 domain=time fund=20 wavelet=ricker2 tdelay=0.12 | Add scale=-1 > $@
	echo "peakFreq=20" >> $@

#data seperately
#for debug
dummy:
	echo ${setGID}
datnl.t%.sg.H: vel.t%.H
	$B/CudaModeling2d.gpux par=$P/tr.P vel=$< \
	wavelet=${wavY} refl_model=ref.H T=3.0 nshots=1 shot1_ix=250 \
	do_offset=0 seis_record=$@ bLinearModeling=0 npad=50 \
	take2ndTimeDerivRTM=1 ${dump}  j_snap=1 \
	src_snap=c1.sgmod.src.H src_snap2=c1.sgmod.src2.H

## %=c1,c2,a1,a2
datnl.t%.H: vel.t%.H
	$B/CudaModeling2d.gpux par=$P/tr.P vel=$< \
	wavelet=${wavY} \
	do_offset=0 seis_record=$@ \
	take2ndTimeDerivRTM=0 bLinearModeling=0 ${dump} \
	src_snap=src.snap.t$*.mod.H \
	src_snap2=src.snap2.t$*.mod.H ${setGID}
#Taper the data on the edge to control edge effects
dat.t%.tp.H: dat.t%.H
	<$< TaperEdgeData.x taperLength=200 >$@


#slowness perturbation plot
$R/vels.v:
	<dS.c1.H ${gcn} out=1.v title="const neg"
	<dS.c2.H ${gcn} out=2.v title="const pos"
	<dS.a1.H ${gcn} out=3.v title="gauss neg"
	<dS.a2.H ${gcn} out=4.v title="gauss pos"
	vppen 1.v 2.v 3.v 4.v gridnum=2.2 >$@
R = ../fig
$R/grad.51s.v:
	<grad.tc1.51.H ${gcspe} out=1.v title="const neg"
	<grad.tc2.51.H ${gcspe} out=2.v title="const pos"
	<grad.ta1.51.H ${gcspe} out=3.v title="gauss neg"
	<grad.ta2.51.H ${gcspe} out=4.v title="gauss pos"
	vppen 1.v 2.v 3.v 4.v gridnum=2.2 >$@

obj = 1

#Test
grad.t%.sg.obj${obj}.H: datnl.t%.sg.H datnl.t0.sg.H
	$B/CudaBeamWET_grad.gpux par=$P/tr.P nshots=1 shot1_ix=250 vel=vel.t0.H data_obs=$< dcal=datnl.t0.sg.H \
	wavelet=${wavY} src_snap=snap1A.H recv_snap=snap1B.H ${toolTags} \
	grad=$@ grads=grads.t$*.sg.H  corr=corr.t$*.sg.H \
	backProjData=bdb.t$*.sg.obj${obj}.H ${dump} aCube=a.$*.sg.H bCube=b.$*.sg.H aCube0=a0.H bCube0=b0.H \
	v0=2000 vTrue=2500 \
	DJDU=DJDU.t$*.sg.H  H=H.t$*.sg.H E=E.t$*.sg.H objFunc=${obj} \
	trcA=trcA.t$*.sg.H trcA1=trcA1.t$*.sg.H trcB1=trcB1.t$*.sg.H \
	A2=A2.t$*.sg.H B2=B2.t$*.sg.H H12=H12.t$*.sg.H \
	Emat=Emat.t$*.sg.H Fmat=Fmat.t$*.sg.H


#RTM using the backprojected data
gradsg.bpj.%.H: %.H
	$B/CudaRTM2d.gpux par=$P/tr.P nshots=1 shot1_ix=250 vel=vel.t0.H	\
	wavelet=${wavY} seis_record=$<  \
	do_offset=0 final_image=$@ ${dump}

%.trc.H: %.H
	make $*.CT2d.H
	cp $*.CT2d.H $@
#Plot sensitivity kernel Take a trace and then Pad to the proper size 
sens.bpj.%.H: %.trc.H
	<$< Pad extend=0 beg2=200 end2=200 >tdat.H; sleep 0.1s
	$B/CudaRTM2d.gpux par=$P/tr.P nshots=1 shot1_ix=250 vel=vel.t0.H	\
	wavelet=${wavY} seis_record=tdat.H  \
	do_offset=0 final_image=$@ ${dump}
	
mdir = /homes/sep/yang/Research/matlab

sens.t%.sg.obj${obj}.H: ${mdir}/bdb.t%.sg.obj${obj}.bp.H
	echo \n peakFreq=20 >> ${mdir}/wavY.bp.H
	echo '\n' d2=10 o2=0 d3=40 o3=0 n3=1 >> $<
	<$< Pad extend=0 beg2=200 end2=200 >tdat.H; sleep 0.1s
	$B/CudaRTM2d.gpux par=$P/tr.P nshots=1 shot1_ix=250 vel=vel.t0.H	\
	wavelet=${mdir}/wavY.bp.H seis_record=tdat.H  \
	do_offset=0 final_image=$@ ${dump}

gradbp.t%.sg.obj${obj}.H: ${mdir}/bdb.t$*.sg.obj${obj}.bp.H 
	$B/CudaRTM2d.gpux par=$P/tr.P nshots=1 shot1_ix=250 vel=vel.t0.H	\
	wavelet=${mdir}/wavY.bp.H seis_record=$<  \
	do_offset=0 final_image=$@ ${dump}

#Do multiple shots
grad.t%.2.H: datnl.t%.H datnl.t0.H
	$B/CudaBeamWET_grad.gpux par=$P/tr.P vel=vel.t0.H data_obs=$< dcal=datnl.t0.H \
	wavelet=${wavY} src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ grads=grads.t$*.2.H objFunc=2 corr=corr.t$*.H \
	backProjData=bdb.t$*.2.H ${dump} ${setGID}
#gradient, Yi Luo
grad.t%.7.H: datnl.t%.H datnl.t0.H
	$B/CudaBeamWET_grad.gpux par=$P/tr.P vel=vel.t0.H data_obs=$< dcal=datnl.t0.H \
	wavelet=${wavY} src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ grads=grads.t$*.7.H objFunc=7 corr=corr.t$*.H \
	backProjData=bdb.t$*.7.H dCdx2=dCdx2.t$*.H ${dump} ${setGID}

#Beam focus result
grad.t%.3.H: datnl.t%.H datnl.t0.H
	$B/CudaBeamWET_grad.gpux par=$P/tr.P vel=vel.t0.H data_obs=$< dcal=datnl.t0.H \
	wavelet=${wavY} src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ grads=grads.t$*.3.H objFunc=3 corr=corr.t$*.H \
	backProjData=bdb.t$*.3.H ${dump} ${setGID} 

#Use Yang, mimimize the slope
grad.t%.5.H: datnl.t%.H datnl.t0.H
	$B/CudaBeamWET_grad.gpux par=$P/tr.P vel=vel.t0.H data_obs=$< dcal=datnl.t0.H \
	wavelet=${wavY} src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ grads=grads.t$*.5.H objFunc=5 corr=corr.t$*.H \
	backProjData=bdb.t$*.5.H ${dump} ${setGID} ${toolTags} \
	bCube=b.t$*.5.H bCube0=b0.t$*.5.H 
#Use Yang's global shift obj function
grad.t%.51.H: datnl.t%.H datnl.t0.H
	$B/CudaBeamWET_grad.gpux par=$P/tr.P vel=vel.t0.H data_obs=$< dcal=datnl.t0.H \
	wavelet=${wavY} src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ grads=grads.t$*.51.H objFunc=51 corr=corr.t$*.H \
	backProjData=bdb.t$*.51.H ${dump} ${setGID} ${toolTags} \
	bCube=b.t$*.51.H bCube0=b0.t$*.51.H 

#Use two param fitting: slope and shift
grad.t%.41.H: datnl.t%.H datnl.t0.H
	$B/CudaBeamWET_grad.gpux par=$P/tr.P vel=vel.t0.H data_obs=$< dcal=datnl.t0.H \
	wavelet=${wavY} src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ grads=grads.t$*.41.H objFunc=41 corr=corr.t$*.H \
	backProjData=bdb.t$*.41.H ${dump} ${setGID} ${toolTags} \
	trcA=trcA.t$*.41.H trcA1=trcA1.t$*.41.H trcB1=trcB1.t$*.41.H \
	A2=A2.t$*.41.H B2=B2.t$*.41.H \
	Emat=Emat.t$*.41.H Fmat=Fmat.t$*.41.H

#Use global obj func, but only use shift moveout params 
grad.t%.32.H: datnl.t%.H datnl.t0.H
	$B/CudaBeamWET_grad.gpux par=$P/tr.P vel=vel.t0.H data_obs=$< dcal=datnl.t0.H \
	wavelet=${wavY} src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ grads=grads.t$*.32.H objFunc=32 corr=corr.t$*.H \
	backProjData=bdb.t$*.32.H ${dump} ${setGID} ${toolTags} \
	trcA=trcA.t$*.32.H trcA1=trcA1.t$*.32.H trcB1=trcB1.t$*.32.H \
	A2=A2.t$*.32.H B2=B2.t$*.32.H \
	Emat=Emat.t$*.32.H Fmat=Fmat.t$*.32.H

grad.t%.35.H: datnl.t%.H datnl.t0.H
	$B/CudaBeamWET_grad.gpux par=$P/tr.P vel=vel.t0.H data_obs=$< dcal=datnl.t0.H \
	wavelet=${wavY} \
	grad=$@ grads=grads.t$*.35.H objFunc=35 corr=corr.t$*.H \
	backProjData=bdb.t$*.35.H ${dump} ${setGID} ${toolTags} \
	trcA=trcA.t$*.35.H trcA1=trcA1.t$*.35.H A2=A2.t$*.35.H \
	Emat=Emat.t$*.35.H Fmat=Fmat.t$*.35.H

#Smaller lag for a1,a2 case to avoid artifacts.
grad.ta%.7.sl.H: datnl.ta%.H datnl.t0.H
	$B/CudaBeamWET_grad.gpux par=$P/tr.P vel=vel.t0.H data_obs=$< dcal=datnl.t0.H \
	start_lag=-25 end_lag=25 \
	wavelet=${wavY} src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ grads=grads.ta$*.7.sl.H objFunc=7 corr=corr.ta$*.sl.H \
	backProjData=bdb.ta$*.7.sl.H dCdx2=dCdx2.ta$*.sl.H ${dump} ${setGID}

#gradient, Yang's guassian Deriv Obj func
grad.t%.9.H: datnl.t%.H datnl.t0.H
	$B/CudaBeamWET_grad.gpux par=$P/tr.P vel=vel.t0.H data_obs=$< dcal=datnl.t0.H \
	wavelet=${wavY} \
	grad=$@ grads=grads.t$*.9.H objFunc=9 corr=corr.t$*.H \
	backProjData=bdb.t$*.9.H widthGaussDeriv=0.8 wt=wt.t$*.9.H \
	start_lag=-500 end_lag=500 \
	${dump} ${setGID}

#Show the velocity of tc1,tc2,ta7,ta8
velshow.9:
	${gce} ${bar0} <vel.tc1.H bias=2000 | Tube &
	${gce} ${bar0} <vel.tc2.H bias=2000 | Tube &
	${gce} ${bar0} <vel.ta7.H bias=2000 | Tube &
	${gce} ${bar0} <vel.ta8.H bias=2000 | Tube &
gradshow.9:
#	${gce} ${bar0} <grad.tc1.9.H  | Tube &
#	${gce} ${bar0} <grad.tc2.9.H  | Tube &
#	${gce} ${bar0} <grad.ta7.9.H  | Tube &
	${gce} ${bar0} <grad.ta8.9.H  | Tube &

bpjdatshow.9:
	${gce} ${bar0} <bdb.tc1.1.CT3d.H  | Tube &
	${gce} ${bar0} <bdb.tc1.9.CT3d.H  | Tube &

#Do inversion
#Test single shot
S.t%.sg.H: datnl.t%.sg.H 
	$B/CudaBeamWET.gpux par=$P/tr.P vel=vel.t0.H data_obs=$< \
	wavelet=${wavY}  >$@ T=3.0 nshots=1 shot1_ix=250 do_offset=0 \
	grads=grads.t$*.sg.2.H models=ms.t$*.sg.2.H fData=fd.t$*.sg.2.H \
	objFunc=2 error=err.t$*.sg.H ${dump} ${setGID} \
	maxNonLinIter=5 maxPerturbScale=0.05 bDoSrcInt=1 bDoRecvInt=0 \
	bDoSmooth=1 n1spl=5 n2spl=20 bVaryStepLength=1
	

S.t%.2.H: datnl.t%.H 
	$B/CudaBeamWET.gpux par=$P/tr.P vel=vel.t0.H data_obs=$< \
	wavelet=${wavY}  >$@ do_offset=0 \
	grads=grads.t$*.2.H models=ms.t$*.2.H fData=fd.t$*.2.H \
	objFunc=2 error=err.t$*.sg.H ${dump} ${setGID} \
	maxNonLinIter=3 maxPerturbScale=0.05 bDoSrcInt=0 bDoRecvInt=0 \
	bDoSmooth=1 n1spl=5 n2spl=20 bVaryStepLength=1
#gradient, Yi Luo



#Test correlation
datnl.t.cr.H: datnl.t0.c.H
	Corr.x startLag=-200 endLag=200 f1=$< f2=datnl.tc1.c.H > c01.H 
	Corr.x startLag=-200 endLag=200 f1=$< f2=datnl.tc2.c.H > c02.H 
	Corr.x startLag=0 endLag=500 f1=c01.H f2=datnl.tc1.c.H > t01.H 
	Corr.x startLag=0 endLag=500 f1=c02.H f2=datnl.tc2.c.H > t02.H 

#Test Bspline smoothing
bTest: grad.ta1.sg.H
	<$< Window n1=1 min1=10 >tmp.H
	<tmp.H TestOpBspline.x n1=50 > a.H    dump1=dump1.H dump2=dump2.H
	<$< TestOpBspline.x n1=16 n2=50 > b.H dump1=dump1.H dump2=dump2.H


