include ${YANGLIB}/makefileY.top

#Test and compare with Yaxun's example


#Marmousi Model
#Do fwd and ajoint operation to Test the WET operators.

#Background velocity always uniform 1000, reflector for born data always at reflDp.H
#Four vel model, constant up 1100, constant down 900
#High Vel anomaly 1100, low vel anomaly 900
R = ../fig
B = ../bin
P = ../par
Y = ~yang/from-other/yaxun_DSO
strip = ~yang/Research/codelib/bin/StripOutDummyAxis.x
tool = ~yang/Research/codelib/bin

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

marm = ../tomo/marm



##Problematic,asymmetric with big jump at t=0


#lowcut
dat.%.lc.H: dat.%.H 
	<$< Bandpass flo=2 fhi=60 > $@
dat.%.bp.H: dat.%.H 
	<$< Bandpass flo=1 fhi=60 > $@

wavY = ../NR/wavY.HH

gid = 1
setGID := gpuID=${gid}
pad := Pad beg1=50 end1=50 
shotsg = T=3.0 nshots=1 shot1_ix=250 shot_inc_x=4 dt_snap=0.05

wemBIN = ~yang/from-other/yaxun_DSO/Bin

MARM_MOD_PAR = ogx=6000 dgx=20 ngx=301 osx=6000 dsx=120 nsx=51
MigParMarm = pdip=45 offset=n ntpx=25 nws=64 oper="foufd" nxpad_beg=50 nxpad_end=50 #image_xmin=-2000 image_xmax=2000
screenratio=0.3 screenwd=12.0 screenht=3.6

BSC_PARA = typbsc=2 nopbsc=3
OnewayParMarm = bUseOneway=1 pdip=45 nws=64 csou=marm.csou.H ${BSC_PARA}

##Oneway data
#First use a Hamming Window shape the src signature, so it is not exactly the same as in time domain
FREQ = f1=5 f2=10 f3=35 f4=40  


#circular dependence, only one target can be enabled at the same time.
../NR/%.HH: %.H
	Cp $< $@ out=stdout

v.T.marm.H: ../NR/v.T.marm.HH
	Cp $< $@
v.0.marm.H: ../NR/v.0.marm.HH
	Cp $< $@
v.z.marm.H: ../NR/v.z.marm.HH
	Cp $< $@
v.sm.marm.H: ../NR/v.sm.marm.HH
	Cp $< $@
ref.marm.H: ../NR/ref.marm.HH
	Cp $< $@
marm.csou.H: ../NR/marm.csou.HH
	Cp $< $@



npadend = 20
#v.T.marm.H: ${marm}/bwi.marm.vmod.H
#	<$< StripOutDummyAxis.x output=a.H
#	<a.H Transp plane=12 | Pad beg1=${npad} end1=${npadend} extend=1 >$@
#v.0.marm.H: ${marm}/bwi.marm.bvel.H
#	<$< StripOutDummyAxis.x output=a.H
#	<a.H Transp plane=12 | Pad beg1=${npad} end1=${npadend} extend=1 >$@
#
#v.z.marm.H: v.T.marm.H
#	<$< StackSpread.x >$@
#
#v.sm.marm.H: v.T.marm.H
#	<$< Smooth rect1=2 rect2=4 repeat=1 >$@
#


#Graduate increase from 1600 to 3200
v.zg.marm.H: v.T.marm.H
	Vel n1=101 n2=301 n3=1 o3=0 d3=1 o1=600 o2=0 vc=1600 d1=20 d2=20 \
			z1=600 vr1=1600 smooth1=3200 >t.H
	< t.H Pad beg1=50 beg2=100 end2=100 extend=1 >$@
	echo o2=4000 >> $@

v.c.marm.H: v.T.marm.H
	Vel n1=101 n2=301 n3=1 o3=0 d3=1 o1=600 o2=0 vc=3000 d1=20 d2=20 > t.H
	< t.H Pad beg1=50 beg2=100 end2=100 extend=1 >$@
	echo o2=4000 >> $@

v.zg_${obj}.marm.H: 
	< Smarm1w.zg.${obj}.sgE0.semb1.wd.recip.H Pad beg1=50 end1=20 beg2=100 end2=100 extend=1 >$@

#The mask is from 7000 -> 11000, cannot set o2=7000 in Vel. Bug in Vel
marm.gradmask.H: 
	Vel n1=101 n2=301 n3=1 o3=0 d3=1 o1=600 o2=0 vc=1 d1=20 d2=20 >t.H
	<t.H Pad beg1=50 beg2=100 end2=100 extend=0 >$@
	echo o2=4000 >>$@

dS.%.marm.H:v.%.marm.H
	Math file1=$< file2=v.sm.marm.H exp='1/file2-1/file1' >$@

dv.%.marm.H:v.%.marm.H
	Math file1=$< file2=v.sm.marm.H exp='file2-file1' >$@

dvi.%.${obj}.marm.H: Smarm1w.%.${obj}.sgE0.semb${semb}.recip.H
	Math file1=$< file2=v.sm.marm.H exp='file2-file1' >$@

%.sm.H:%.H
	<$< Smooth rect2=10 repeat=2 >$@

#Window slowness result
%.wd.H:%.H
	< $< Window f1=50 n1=81 f2=100 n2=301 >$@


#from (z,x) to (hx,hy,x,y,z) and truncate top

v1w.%.marm.H: v.%.marm.H
	<$< Window f1=${npad} | Transp reshape=1,3 plane=12 >$@


#the survey size 16km, 45deg is the limit for oneway, thus Cable length 4km is enough
dat1w.marm.H: v1w.T.marm.H  marm.csou.H ref.marm.H
	${wemBIN}/wem3d.x csou=marm.csou.H vmod=$< imag=ref.marm.H crec=$@ \
        adj=n ${MigParMarm} ${MARM_MOD_PAR}


Test_img1w.%.marm.H: dat1w.marm.H v1w.%.marm.H ${wavY}
	${wemBIN}/wem3d.x csou=marm.csou.H vmod=v1w.$*.marm.H crec=$< imag=$@ \
        adj=y ${MigParMarm}

dat21wsg.%.${ref}.H: velmarm1w.%.H  csou.H
	${wemBIN}/wem3d.x csou=csou.H vmod=$< imag=${ref}.1w.H crec=$@ \
        adj=n ngx=201 ogx=0 dgx=20 nsx=1 offset=y dsx=80 osx=0 ${MigPar2}

img21wsg.%.${ref}.H: dat21wsg.%.${ref}.H velmarm1w.%.H csou.H
	${wemBIN}/wem3d.x csou=csou.H vmod=velmarm1w.$*.H crec=$< imag=$@ \
        adj=y ${MigPar2} offset=y



######Oneway migration
img1w.%.marm.H: dat1w.marm.H v.%.marm.H
	$B/CudaRTM2d.gpux par=$P/marm.P vel=v.$*.marm.H \
	wavelet=${wavY} seis_record=$< ${OnewayParMarm} ${dump} \
	do_offset=1 final_off_image=$@ final_image=imgz0.$*.H 

####One way Tomadjoint
dI1w.%.marm.fwd.H: dS.%.marm.H
	$B/CudaTomFwd.gpux par=$P/marm.P vel=v.$*.marm.H \
	 wavelet=${wavY} data_obs=dat1w.marm.H ${OnewayParMarm} ${dump} \
	do_offset=1 >$@ deltaI=$@ deltaS=$<

dS1w.%.marm.adj.H: dI1w.%.marm.fwd.H
	$B/CudaTomAdj.gpux par=$P/marm.P vel=v.$*.marm.H \
	 wavelet=${wavY} data_obs=dat1w.marm.H ${dump} ${OnewayParMarm} \
	do_offset=1 >$@ deltaIh=$< deltaSs=t.H

bUseSingleMoveoutPar = 0
obj		= zref1
objFunc 	= 1
sgEvent 	= 0
bUWLS 	= 0
thresMax = 0.005

mEventPar_marm = beam_z_interv=100 beam_z_window=400 beam_z_o=790 beam_z_n=16
ifeq (0,${sgEvent})
	mEventPar = ${mEventPar_marm}
endif


### Do inversion using oneway
ifneq (,$(findstring bf2, ${obj})) #Use second param
	bUseSingleMoveoutPar = 0
	objFunc = 5
endif

ifeq (${obj},bf1) #Use second param
	bUseSingleMoveoutPar = 1
	objFunc = 5
endif

bPick = 0
#if use bf method
ifneq (,$(findstring bf, ${obj}))
ifneq (,$(findstring g, ${obj}))
	bUseGauss = 1
endif
ifneq (,$(findstring p, ${obj}))
	bPick = 1
endif
endif

ifneq (,$(findstring bf1m, ${obj}))
	objFunc = 55
	bUseSingleMoveoutPar = 1
	bUWLS = 1
endif

ifneq (,$(findstring bf1mn,${obj}))
	objFunc = 55
	bUseSingleMoveoutPar = 1
	bUWLS = 0
	bNoLSfit = 1
endif

ifeq (${obj},1) #Maximize stack power
	objFunc = 1
endif
ifeq (${obj},dso)
	objFunc = 6
endif

#If use zref method
ifneq (,$(findstring zref, ${obj}))
	objFunc = 4
	ifeq (zref0,${obj})
		corrUsage = 0
	endif
	ifeq (zref1,${obj})
		corrUsage = 1
	endif
	ifeq (zref2,${obj})
		corrUsage = 2
	endif
endif

semb = 1
#sembPar = bUseSemblance=${semb} ax_semb_d=5.0 ax_semb_n=201 fSembSmoothWindow=400 bUseGaussianDerivForDJDU=${bUseGauss}
#For constant velocity

sembPar = bUseSemblance=${semb} ax_semb_d=5.0 ax_semb_n=201 fSembSmoothWindow=500 bUseGaussianDerivForDJDU=${bUseGauss} bNormSemb=${bNormSemb}
ifeq (${semb},1)
ifneq (,$(findstring nm, ${obj})) #Use second param
  bNormSemb = 1
else
  bNormSemb = 0
endif
endif


aa = 1

ifeq (${aa} , 1) #Notice the space before the $ is crucial, if there is one space, then ' ${aa}' != 1, !!!! Makefile Headache Grammar.
 #$(error specify aa = 1)
else
 #$(error specify aa = ${aa} sj)
endif

test.marm:
	echo ${aa}
	echo obj=${objFunc}

aa = 2

otherpar = corrUsage=${corrUsage} bPickMoveoutPar=${bPick}  bNoLSfit=${bNoLSfit} bDivideUseMask=1 \
					 bUseWeightedLSfit=${bUWLS} gradMask=marm.gradmask.H bMaskGradient=1 \
					 ax_ang_n=61  



 

sens_dImarm1w.%.${obj}.H:
	$B/DimgWEMVA.x <img1w.$*.marm.H output=$@ par=$P/marm.P objFunc=${objFunc}  \
		dJdu=sens.dJdu.$*.${obj}.sgE${sgEvent}.semb${semb}.H Ecube=Ecubemarm1w.$*.${obj}.sgE${sgEvent}.semb${semb}.H Fcube=Fcubemarm1w.$*.${obj}.sgE${sgEvent}.semb${semb}.H \
		Emask=Emaskmarm1w.$*.${obj}.sgE${sgEvent}.semb${semb}.H finalWt=fWtsensmarm1w.$*.${obj}.H thresMax=${thresMax} bUseSingleMoveoutParam=${bUseSingleMoveoutPar} bSensKernel=1 bSingleEvent=0 \
		${mEventPar} ${sembPar} ${otherpar} \
		${dump} 


sens_dSmarm1w.%.${obj}.H: sens_dImarm1w.%.${obj}.H
	$B/CudaTomAdj.gpux par=$P/marm.P vel=v.$*.marm.H \
	 wavelet=${wavY} data_obs=dat1w.marm.H ${OnewayParMarm} \
	do_offset=1 >$@ deltaIh=$< 
	${setGID} ${dump}

#img1w.%.marm.H
dImarm1w.%.${obj}.sgE${sgEvent}.semb${semb}.H: 
	$B/DimgWEMVA.x <img1w.$*.marm.H output=$@ par=$P/marm.P objFunc=${objFunc}  \
		angCube=angCubemarm1w.$*.H angCubeBp=angCubeBpmarm1w.$*.${obj}.sgE${sgEvent}.semb${semb}.H Acube=Acubemarm1w.$*.${obj}.sgE${sgEvent}.semb${semb}.H Bcube=Bcubemarm1w.$*.${obj}.sgE${sgEvent}.semb${semb}.H ABcube=ABcubemarm1w.$*.${obj}.sgE${sgEvent}.semb${semb}.H \
		dJdu=dImarm1w.dJdu.$*.${obj}.sgE${sgEvent}.semb${semb}.H Ecube=Ecubemarm1w.$*.${obj}.sgE${sgEvent}.semb${semb}.H Fcube=Fcubemarm1w.$*.${obj}.sgE${sgEvent}.semb${semb}.H \
		Emask=Emaskmarm1w.$*.${obj}.sgE${sgEvent}.semb${semb}.H finalWt=fWtmarm1w.$*.${obj}.sgE${sgEvent}.semb${semb}.H thresMax=${thresMax} bUseSingleMoveoutParam=${bUseSingleMoveoutPar} bSensKernel=0 bSingleEvent=${sgEvent} \
		${mEventPar} ${sembPar} ${otherpar} \
		sembS=semb.$*.${obj}.sgE${sgEvent}.H corr=corr.$*.${obj}.sgE${sgEvent}.H wtAll=wtmarm1w.$*.${obj}.semb${semb}.H \
		${dump} 


dSmarm1w.%.${obj}.sgE${sgEvent}.semb${semb}.H: dImarm1w.%.${obj}.sgE${sgEvent}.semb${semb}.H 
	$B/CudaTomAdj.gpux par=$P/marm.P vel=v.$*.marm.H \
	 wavelet=${wavY} data_obs=dat1w.marm.H ${OnewayParMarm} \
	do_offset=1 >$@ deltaIh=$< 
	${setGID} ${dump}


#The target region is 81 X 301 pts, the vel model is 151X501 pts, 20m spacing inital smooth scale is 1.5wl * 10wl = 300m X 1000m = 10pts X 7pts
solvPar =	maxNonLinIter=20 maxPerturbScale=0.1 bDoSmooth=1 n1spl=500 n2spl=15 \
					bVaryStepLength=1 	bVarySmoothing=1 fmaxVel=4800 fminVel=1500 \
				  maxPerturbBase=0.3 maxPerturbDecay=0.9

#We need to do inversion now

#Test inversion using single shots
sgOffpar = img_dhx=1 img_nhx=11 img_ohx=-5
S1wsg.%.${obj}.sgE${sgEvent}.semb${semb}.H: dat1wsg.%.H
	$B/CudaBeamWETrefl.gpux >$@ par=$P/marm.P objFunc=${objFunc} data_obs=dat1wsg.$*.H vel=vel2.c0.H \
		wavelet=${wavY} grads=gradssg.$*.${obj}.sgE${sgEvent}.semb${semb}.H mods=modssg.$*.${obj}.sgE${sgEvent}.semb${semb}.H dats=datssg.$*.${obj}.sgE${sgEvent}.semb${semb}.H eerror=errsg.$*.${obj}.sgE${sgEvent}.semb${semb}.H fImage=fImgsg.$*.${obj}.sgE${sgEvent}.semb${semb}.H ${sembPar} ${OnewayParMarm} ${solvPar} \
		thresMax=${thresMax} bUseSingleMoveoutParam=${bUseSingleMoveoutPar} bSensKernel=0 bSingleEvent=${sgEvent} ${mEventPar} ${otherpar} bUseWeightedLSfit=${bUWLS} \
		angCube=SangCube1wsg.$*.${obj}.H angCubeBp=SangCubeBp1wsg.$*.${obj}.H ABcube=SABcube1wsg.$*.${obj}.H dJdu=Smarm1w.dJdusg.$*.${obj}.H Ecube=SEcube1wsg.$*.${obj}.H Fcube=SFcube1wsg.$*.${obj}.H \
		Emask=SEmask1wsg.$*.${obj}.H corr=Scorr.$*.${obj}.sgE${sgEvent}.H finalWt=SfWt1wsg.$*.${obj}.H ${sgOffpar} ${dump} 



Smarm1w.%.${obj}.sgE${sgEvent}.semb${semb}.H: v.%.marm.H dat1w.marm.H ${wayY} marm.gradmask.H
	$B/beamWETreflOneway.x >$@ par=$P/marm.P objFunc=${objFunc} data_obs=dat1w.marm.H vel=v.$*.marm.H \
		wavelet=${wavY} grads=grads.$*.${obj}.sgE${sgEvent}.semb${semb}.H mods=mods.$*.${obj}.sgE${sgEvent}.semb${semb}.H dats=dats.$*.${obj}.sgE${sgEvent}.semb${semb}.H error=err.$*.${obj}.sgE${sgEvent}.semb${semb}.H fImage=fImg.$*.${obj}.sgE${sgEvent}.semb${semb}.H ${sembPar} ${OnewayParMarm} ${solvPar} \
		thresMax=${thresMax} bUseSingleMoveoutParam=${bUseSingleMoveoutPar} bSensKernel=0 bSingleEvent=${sgEvent} ${mEventPar} ${otherpar} bUseWeightedLSfit=${bUWLS} \
		angCube=SangCubemarm1w.$*.${obj}.H angCubeBp=SangCubeBpmarm1w.$*.${obj}.H ABcube=SABcubemarm1w.$*.${obj}.H dJdu=Smarm1w.dJdu.$*.${obj}.H Ecube=SEcubemarm1w.$*.${obj}.H Fcube=SFcubemarm1w.$*.${obj}.H Acube=SAcubemarm1w.$*.${obj}.H Bcube=SBcubemarm1w.$*.${obj}.H \
		Emask=SEmaskmarm1w.$*.${obj}.H finalWt=SfWtmarm1w.$*.${obj}.H ${dump} 

##To compare the final image of marm result
%.slopem.H: %.H
	<$< Transp plane=23 | Off2slope.x amax=50 n=61 ${dump} f_AA=1.0 b_AA=1| Transp plane=23 >$@

#Convert to angle domain
fImgwd.%.dso.H:
	make fImg.$*.dso.sgE0.semb1.slopem.H
	make fImg.$*.dso.sgE0.semb1.slopem.wd.H 
	Cp fImg.$*.dso.sgE0.semb1.slopem.wd.H $@
fImgwd.%.bf2g.H:
	< fImg.$*.bf2g.sgE0.semb1.H Transp plane=23 > t.H
	make t.wd.H
	Cp t.wd.H $@
#Migrate using the inverted slowness
fImgOff.zg.bf2g.H:
	Cp  Smarm1w.zg.bf2g.sgE0.semb1.recip.H v.zgbf2g.marm.H
	make img1w.zgbf2g.marm.H
	make img1w.zgbf2g.marm.slopem.H
	cp img1w.zgbf2g.marm.H $@
	
imgwd.T.H:
	make img1w.T.marm.slopem.wd.H
	Cp 	img1w.T.marm.slopem.wd.H $@
imgwd.zg.H:
	make img1w.zg.marm.slopem.wd.H
	Cp 	img1w.zg.marm.slopem.wd.H $@
viewImg:
	${Icube} imgwd.zg.H imgwd.T.H fImgwd.zg.dso.H fImgwd.zg.bf2g.H nviews=4

#The image stack for four images, true img, initial img, dso img, and RMO image
img.stks.zg.H:
	< img1w.T.marm.CT3d.wd.H  Window f1=5 > t1.H
	< img1w.zg.marm.CT3d.wd.H  Window f1=5 > t2.H
	< fImg.zg.dso.sgE0.semb1.CT3d.wd.H  Window f1=5 >t3.H
	< fImgOff.zg.bf2g.CT3d.wd.H  Window f1=5 > t4.H
	Cat t1.H t2.H t3.H t4.H >$@

# 0,1,2,3
$R/img.stk.zg.%.v: img.stks.zg.H
	<$< Window f3=$* n3=1 | ${gra} clip=0.65 labelfat=3 titlefat=3 label1='z(m)' label2='x(m)' \
		screenratio=0.4 screenwd=14 screenht=5.6 title=' ' out=$@ ${nul} ${bar0}

graph_label = labelfat=3 titlefat=3 label1='x(m)' plotfat=3
grey_label =  labelfat=3 titlefat=3 label1='z(m)' label2='x(m)' screenratio=0.4 screenwd=12.0 screenht=4.8
blabel  = barlabel="velocity(m/s)"
velclip = newclip=1 bclip=1600 eclip=3000 


# For SEP REPORT
$R/v.%.marm.v:
	< v.$*.marm.wd.H ${gcn} out=$@ ${nul} ${grey_label} ${blabel} ${velclip}
$R/v.%.marm.ic.v:
	< v.$*.marm.wd.H ${gcn} out=$@ ${nul} ${grey_label} ${blabel} ${velclip}
$R/img.%.marm.0off.v:
	< img1w.$*.marm.CT3d.wd.H ${gre}  out=$@ ${nul} ${grey_label} >$@ ${bar0}

#Plot the inverted velocity model
$R/vinv.%.dso.marm.v:
	< Smarm1w.$*.dso.sgE0.semb1.wd.recip.H ${gcn}  out=$@ ${nul} ${grey_label} ${blabel} ${velclip} >$@ 
$R/vinv.%.bf2g.marm.v:
	< Smarm1w.$*.bf2g.sgE0.semb1.wd.recip.H ${gcn}  out=$@ ${nul} ${grey_label} ${blabel} ${velclip} >$@ 

$R/vinv.%.dso.marm.ic.v:
	< Smarm1w.$*.dso.sgE0.semb1.wd.recip.H  ${gcn}  out=$@ ${nul} ${grey_label} ${blabel} ${velclip} >$@ 
$R/vinv.%.bf2g.marm.ic.v:
	< Smarm1w.$*.bf2g.sgE0.semb1.wd.recip.H ${gcn}  out=$@ ${nul} ${grey_label} ${blabel} ${velclip} >$@ 



$R/grad.0.marm.bf2.v:
	< dSmarm1w.0.bf2.sgE0.semb0.wd.H ${gce} out=$@ ${nul} ${grey_label} title='slow grad dJ/du'
$R/grad.0.marm.dso.v:
	< dSmarm1w.0.dso.sgE0.wd.H ${gce} out=$@ ${nul} ${grey_label} title='slow grad dso'
$R/grad.0.marm.zref0.v:
	< dSmarm1w.0.zref0.sgE0.semb0.wd.H ${gce} out=$@ ${nul} ${grey_label} title='slow grad zref0'


### An example that correlation method is error prone when multiple events runs into together.
corr_error:
	< corr.zg.zref0.sgE0.H Window n4=1 min4=10620 n3=1 min3=1190 | ${gce} | Tube& #The corr panel
	< dImarm1w.dJdu.zg.zref0.sgE0.semb1.H Window n3=1 min3=10620 | Transp plane=12 | ${gce} | Tube& #The errrouneous picked max lag
	< angCubemarm1w.zg.H Window n3=1 min3=10620 > ang.H
	< ang.H ${gce} | Tube& #The angle gather
	< ang.H StackSpread.x > angf.H
	< ang.H Window min1=1000 max1=1380 	| ${gce} title='ang'  | Tube &
	< angf.H Window min1=1000 max1=1380 | ${gce} title='angf' | Tube &
corr2:
	#Corr.x f1=ang.H f2=ang.H > cc.H startLag=1 endLag=10
	< ang.H StackSpread.x > angf.H
###Plot something
###Semblance panel


#Taper the data on the edge to control edge effects
dat.%.tp.H: dat.%.H
	<$< TaperEdgeData.x taperLength=200 >$@



