include $(SEPINC)/SEP.top


F90LIBS = -lsepimagef90  -lsupersetf90 -lsepgeef90 -lsuperset -lsepauxf90 -lsepmathf90 -lsep2df90 -lsep3df90 -lsep3d -lsepf90 -lsep -L/opt/FFTW/lib -lfftw3f 


RESULTSER =

RESULTCR = 
RESULTSNR = 

#UF90FLAGS = -O0 -g -traceback -fpe:0 -check all

BINDIR = ./bin
SRCDIR = ./src
RESDIR = ./fig
OBJDIR = ./obj

B = ${BINDIR}
S = ${SRCDIR}
R = ${RESDIR}
O = ${OBJDIR}

HTMLDIR = .

default1:

default:
	make -k -i default1; make -k -i default1

dn = >/dev/null
grey = Grey $(dn)

#clip = pclip=99 wantscalebar=1
clip = newclip=1 crowd=0.80 # wantscalebar=1
FAT = 3

$B/Simpvelmodel.x: ${OBJDIR}/simpvelmod.o

#################### Finite Difference Modeling 4th order in time and 10th order in space
seis.H: $B/acoustic2d_4t8x.x
	$B/acoustic2d_4t8x.x par=par.P vel=vel-homo.H snap=snap.H

#### Window and subsample the marmousi model to dx=10m dz=4m ######
vel/vel-smarm.H: ~xukai/huo/fortran/Vmodel/marmousi/marmvel.H
	< $< Window3d min2=3000 max2=7000 > vel\vel1.H
	< vel\vel1.H Interp d1out=8 d2out=10 >$@

#### Generate the simple two layer velocity model with a dip velocity boundary using c
n1 = 200
n2 = 250
d1 = 10
d2 = 10

vel/simpmod-flat.H: $B/SimpVelModel2d.x
	$< d1=${d1} d2=${d2} n1=${n1} n2=${n2}  > $@

vel/simpmod-dipping.H: $B/SimpVelModel2d.x
	$< d1=${d1} d2=${d2} n1=${n1} n2=${n2} slope=0.3  > $@

vel/simpmod-steep.H: $B/SimpVelModel2d.x
	$< d1=${d1} d2=${d2} n1=${n1} n2=${n2} slope=-2.0  > $@

#### Create the data file from modeling program

### Use model Marmousi
#seis-smarm.H: vel/vel-smarm.H  $B/acoustic2d_4t10x.x
#	$B/acoustic2d_4t10x.x par=par/marm1.P < $< Snap=snap-smarm.H >$@


### Use model simple
#seis-simp-flat.H: vel/simpmod-flat.H  $B/acoustic2d_4t10x.x
#	$B/acoustic2d_4t10x.x par=par/simp.P < $< Snap=snap-simp-flat.H srcwav=wavlet.H >$@

centerFreq = 15
maxFreq = 25
#generate the wavelet
#since peak freq=30Hz. the length should be at least 1/30*2/dt = 16pts

#ricker.H:
#	Wavelet n1=256 d1=0.001 domain=time wavelet=ricker2 fund=${centerFreq}\
#		 tdelay=0.128 phase=none \
#	  >tmp.H
#		Math file1=tmp.H exp='file1+0.0' >$@
#	< tmp.H Add scale=-1 | Bandpass flo=2 nplo=5 nphi=5 fhi=${maxFreq} > $@
#	< tmp2.H Bandpass flo=2 fhi={maxFreq} > $@  || this will make the zero-freq worse
#	echo peakFreq=${maxFreq} >> $@

rk.H:
	Wavelet n1=1024 d1=0.001 domain=time wavelet=ricker2 fund=30 tdelay=0.512 >$@

hf = 3
rk.%.${hf}.H: rk.H
	<rk.H Bandpass flo=$* fhi=${hf} > $@

#	< $< Bandpass flo=10 fhi=60 phase=0 >$@
delta.H: 
	Spike nsp=1 n1=64 k1=1 mag=1.0 d1=0.004 >$@
	echo freq=40 >> $@

#Impulse response of wave modeling
seis-simp-flat.H: vel/simpmod-flat.H ricker.H $B/Modeling2d.x 
	$B/Modeling2d.x par=par/simp.flat.P vel=$< \
	src_snap=snap-simp-flat.H wavelet=delta.H \
	wavelet_interped=wavInterped.H recv_snap=dummy.H \
	seis_record=$@ 

#actually 20Hz is the central frequency.
ricker.H:
	cp ~/Research/matlab/ricker.20.H $@
	echo '' >> $@
	echo  'peakFreq=30' >> $@

#bench mark Fdmod by SEPlib
seis-std.H: vel/simpmod-flat.H
	Fdmod intag=$< outtag=snap-std.H hsfile=$@ mt=28 no_stdout=0 tmax=2.0 fpeak=30 xs=625 zs=500

#Do the modeling using new RTM code
seis.simp.flat.H: vel/simpmod-flat.H ricker.H $B/Modeling2d.x 
	$B/Modeling2d.x par=par/simp.flat.P vel=$< \
	src_snap=snap.s.f.H wavelet=ricker.H \
	recv_snap=dummy.H \
	seis_record=$@ dump=dump.H dump2=dump2.H

#Do the RTM
vel/simpmod-flat-smooth.H: vel/simpmod-flat.H
	< $< Smooth rect1=20 rect2=10 > $@

#Adj RTM
refl.s.f.H: vel/simpmod-flat.H vel/simpmod-flat-smooth.H
	Math file1=vel/simpmod-flat.H file2=vel/simpmod-flat-smooth.H exp=file1-file2 >$@ 


# RTM
img.s.H imgs.s.H: vel/simpmod-flat-smooth.H $B/RTM2d.x par/RTMSimpVel.P
	$B/RTM2d.x par=par/RTMSimpVel.P vel=$< \
	src_snap=snap.s.H wavelet=ricker.H \
	recv_snap=recv.snap.s.H \
	seis_record=seis.lin.s.H RTM_image=imgs.s.H \
  	dump=dump.H final_RTM_image=img.s.H 

#Calculate the gradient of the local objective function

dumps = dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H

seis.lin.s.H: vel/simpmod-flat-smooth.H $B/AdjRTM2d.x 
	$B/AdjRTM2d.x par=par/simp.flat.P vel=$< \
	src_snap=src.snap.s.H wavelet=ricker.H \
	recv_snap=recv.snap.s.H refl_model=refl.s.f.H \
	seis_record=$@ dump=dump.H 
seis.lin.s.2.H: vel/simpmod-flat-smooth.H $B/AdjRTM2d.x 
	$B/AdjRTM2d.x par=par/simp.flat.P vel=$< \
	src_snap=src.snap.s.2.H wavelet=ricker.H \
	recv_snap=recv.snap.s.2.H refl_model=refl.s.f.H \
	seis_record=$@ dump=dump.H boundary_type=4

seis.lin.s-g.H: vel/simpmod-flat-smooth.H 
	$B/CudaModeling2d.gpux par=par/simp.flat.P vel=$< \
	src_snap=src.snap.s-g.H wavelet=ricker.H \
	recv_snap=recv.snap.s-g.H refl_model=refl.s.f.H \
	seis_record=$@ is_linear_modeling=1 \
	take_2nd_deriv_adjRTM=1 \
	${dumps}

seis.nl.s-g.H: vel/simpmod-flat.H ricker.H 
	$B/CudaModeling2d.gpux par=par/simp.flat.P vel=$< \
	src_snap=snap.s.f.H wavelet=ricker.H \
	recv_snap=dummy.H \
	seis_record=$@ dump=dump.H dump2=dump2.H
seis.lin.s.2nd.H: seis.lin.s.H
	$B/TimeDeriv.x <$< order=2 > $@

#Make simple uniform vel
dx=20

vel/vel1.H:
	Vel n1=140 n2=300 n3=1 o3=0 d3=1 o1=0 o2=0 vc=3000 d1=${dx} d2=${dx} > $@

## %1.0 vel perturbation
vel/vel2.H: vel/vel1.H
	Vel n1=140 n2=300 n3=1 o3=0 d3=1 o1=0 o2=0 vc=3030 d1=${dx} d2=${dx} > $@
## %-1.0 vel perturbation
vel/vel3.H: vel/vel1.H
	Vel n1=140 n2=300 n3=1 o3=0 d3=1 o1=0 o2=0 vc=2970 d1=${dx} d2=${dx} > $@

# %10 perturb
vel/vel4.H:vel/vel1.H
	Vel n1=140 n2=300 n3=1 o3=0 d3=1 o1=0 o2=0 vc=3300 d1=${dx} d2=${dx} > $@
# -%10 perturb
vel/vel5.H:vel/vel1.H
	Vel n1=140 n2=500 n3=1 o3=0 d3=1 o1=0 o2=0 vc=2700 d1=${dx} d2=${dx} > $@


#True vel, with the anomaly
vel/vela1.H:
	Vel n1=140 n2=500 n3=1 o3=0 d3=1 o1=0 o2=0 vc=3000 xan=2500 zan=500 \
	ranx=400 ranz=200 rankx=200 rankz=100 exank=1.0 exan=1.0 dvank=-300 \
	dvan=-300 d1=10 d2=10 > $@
#without the anamoly
vel/vela2.H: 
	Vel n1=140 n2=500 n3=1 o3=0 d3=1 o1=0 o2=0 vc=3000 d1=10 d2=10 > $@
vel/vela3.H:
	Vel n1=140 n2=500 n3=1 o3=0 d3=1 o1=0 o2=0 vc=3000 d1=10 d2=10 > $@

img1.H imgs1.H : seis-simp.H vel/vel1.H ricker.H $B/RTM2d.x par/beamSimp.P
	$B/RTM2d.x par=par/beamSimp.P vel=vel/vel1.H \
	src_snap=snap1.H wavelet=ricker.H \
	seis_record=$< dump=wav2.H final_RTM-image=$@ RTM_image=imgs1.H


#Single Shots
seis1.H: vel/vel1.H ricker.H $B/Modeling2d.x par/beamSimp.P
	$B/Modeling2d.x par=par/beamSimp.P vel=$< \
	src_snap=snap1.H wavelet=ricker.H \
	seis_record=$@ dump=wav2.H

grad1.H: vel/vel4.H ricker.H $B/beamWET_grad.x par/beamSimp.P seis1.H
	$B/beamWET_grad.x par=par/beamSimp.P vel=$< data_obs=seis1.H \
	wavelet=ricker.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=3 moveOutType=2 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H	

#Yi Luo's Method
grad1-2.H: vel/vel3.H ricker.H par/beamSimp.P seis1.H
	$B/beamWET_grad.x par=par/beamSimp.P vel=$< data_obs=seis1.H \
	wavelet=ricker.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=2 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H	

#classical FWI gradient
grad1-3.H: vel/vel3.H ricker.H $B/beamWET.x par/beamSimp.P seis1.H
	$B/beamWET_grad.x par=par/beamSimp.P vel=$< data_obs=seis1.H \
	wavelet=ricker.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=1 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H	

#slow vel
seis3.H: vel/vel1.H ricker.H $B/Modeling2d.x par/beamSimp.P
	$B/Modeling2d.x par=par/beamSimp.P vel=$< \
	src_snap=snap1.H wavelet=ricker.H \
	seis_record=$@ dump=wav2.H

grad3.H: vel/vel4.H ricker.H $B/beamWET.x par/beamSimp.P seis1.H
	$B/beamWET_grad.x par=par/beamSimp.P vel=$< data_obs=seis1.H \
	wavelet=ricker.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=3 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H	

#impluse response Yi Lu
seisImp.H: vel/vel1.H ricker.H $B/Modeling2d.x 
	$B/Modeling2d.x par=par/beamSimpImpulse.P vel=$< \
	src_snap=snap1.H wavelet=ricker.H \
	seis_record=$@ dump=wav2.H

gradImp.H: vel/vel3.H ricker.H $B/beamWET.x 
	$B/beamWET_grad.x par=par/beamSimpImpulse.P vel=$< data_obs=seisImp.H \
	wavelet=ricker.H recv_snap=snap1B.H \
	grad=$@ objFunc=2 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H	
#Multiple Shots
seis2.H: vel/vel1.H ricker.H $B/Modeling2d.x par/beamSimp2.P
	$B/Modeling2d.x par=par/beamSimp2.P vel=$< \
	src_snap=snap1.H wavelet=ricker.H \
	seis_record=$@ dump=wav2.H

#Yi Luo's gradient
grad2-2.H: vel/vel3.H ricker.H $B/beamWET.x par/beamSimp2.P seis2.H
	$B/beamWET_grad.x par=par/beamSimp2.P vel=$< data_obs=seis2.H \
	wavelet=ricker.H recv_snap=snap1B.H \
	grad=$@ objFunc=2 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H	

#Do the FWI for multiple iterations
#YL's method multiple shots
velMod2-2.H: vel/vel3.H ricker.H $B/beamWET.x par/beamSimp2.P seis2.H
	$B/beamWET.x par=par/beamSimp2.P vel=$< data_obs=seis2.H \
	wavelet=ricker.H recv_snap=snap1B.H grad=grad2-2.H \
	objFunc=2 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 finalData=dat_$@ >$@

#single shot
velMod1-2.H: vel/vel3.H ricker.H $B/beamWET.x par/beamSimp.P seis1.H
	$B/beamWET.x par=par/beamSimp.P vel=$< data_obs=seis1.H \
	wavelet=ricker.H recv_snap=snap1B.H grad=grad2-2.H \
	objFunc=2 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 \
	>$@ grad=grads1-2.H finalData=dat_$@

#regular 1D DaveHale's boundary
#There might be coding bugs. It is less effect than boundary_type=1
seis-simp-flat-regular.H: vel/simpmod-flat.H ricker.H $B/Modeling2d.x par/modelingSimpVel.P
	$B/Modeling2d.x par=par/modelingSimpVel.P vel=$< \
	src_snap=snap-simp-flat-regular.H wavelet=ricker.H \
	wavelet_interped=wavInterped.H recv_snap=dummy.H \
	seis_record=$@ boundary_type=3

cropped-image-simp-flat.H: image-simp-flat.H
	< $< Window3d min1=300.0 > $@

imag-stack.H: cropped-image-simp-flat.H
	< $< Stack axis=3 >$@

###  compare Ali

aliN1=120
aliN2=241

vel/v1.H:
	Vel n1=${aliN1} n2=${aliN2} n3=1 o3=0 d3=1 o1=0 o2=0 vc=3000 d1=25 d2=25 > $@
vel/v2.H: vel/v1.H
	Vel n1=${aliN1} n2=${aliN2} n3=1 o3=0 d3=1 o1=0 o2=0 vc=2900 d1=25 d2=25 > $@
vel/v3.H: vel/v1.H
	Vel n1=${aliN1} n2=${aliN2} n3=1 o3=0 d3=1 o1=0 o2=0 vc=3300 d1=25 d2=25 > $@
vel/v4.H: vel/v1.H
	Vel n1=${aliN1} n2=${aliN2} n3=1 o3=0 d3=1 o1=0 o2=0 vc=3500 d1=25 d2=25 > $@
vel/v5.H: vel/v1.H
	Vel n1=${aliN1} n2=${aliN2} n3=1 o3=0 d3=1 o1=0 o2=0 vc=3030 d1=25 d2=25 > $@

wav15.H:
	Wavelet n1=550 d1=0.0038 wavelet=ricker2 phase=0 fund=15 tdelay=0.075 > $@
	echo 'peakFreq=15' >> $@

seisa1.H: vel/v1.H ricker.H $B/Modeling2d.x 
	$B/Modeling2d.x par=par/aliComp.P vel=$< \
	src_snap=snap1.H wavelet=wav15.H \
	seis_record=$@ dump=wav2.H 
seisa1g.H: vel/v1.H $B/CudaModeling2d.gpux 
	$B/CudaModeling2d.gpux par=par/aliComp.P vel=$< \
	src_snap=snap1.H wavelet=wav15.H \
	seis_record=$@ dump1=dump1.H dump2=dump2.H dump3=dump3.H



grada3-2.H: vel/v3.H ricker.H 
	$B/beamWET_grad.x par=par/aliComp.P vel=$< data_obs=seisa1.H \
	wavelet=wav.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=2 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H backProjData=bda3-2.H

grada3-5.H: vel/v3.H ricker.H 
	$B/beamWET_grad.x par=par/aliComp.P vel=$< data_obs=seisa1.H \
	wavelet=wav15.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=5 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H backProjData=bda3-5.H \
	vTrue=3000 v0=3300 corr=corra3.H 
grada3-6.H: vel/v3.H ricker.H 
	$B/beamWET_grad.x par=par/aliComp.P vel=$< data_obs=seisa1.H \
	wavelet=wav15.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=6 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H backProjData=bda3-6.H \
	vTrue=3000 v0=3300 corr=corra3.H


#Do the FWI for multiple iterations !!!!   !!!!!! 

$R/%.v: %.H
	< $< Grey newclip=1 color=j wantscalebar=1 out=$@ title=" " >/dev/null

$R/%.pdf: $R/%.v
	pstexpen $< $R/$*.ps color=y fat=3 fatmult=1.5 invras=n
	epstopdf $R/$*.ps


grada3g.H: vel/v3.H ricker.H $B/CudaBeamWET_grad.gpux 
	$B/CudaBeamWET_grad.gpux par=par/aliComp.P vel=$< data_obs=seisa1.H \
	wavelet=wav15.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=2 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H	
	
grada2.H: vel/v2.H 
	$B/beamWET_grad.x par=par/aliComp.P vel=$< data_obs=seisa1.H \
	wavelet=wav15.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=2 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H	

%_w1.H: %.H
	< $< Window n1=80 f1=20 n2=201 f2=20 > $@
#biondo's method
grada2_o3.H: vel/v4.H ricker.H $B/beamWET.x 
	$B/beamWET_grad.x par=par/aliComp.P vel=$< data_obs=seisa1.H \
	wavelet=wav.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=3 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H	

seisa1t.H: vel/v1.H ricker.H $B/Modeling2d.x 
	$B/Modeling2d.x par=par/aliComp1t.P vel=$< \
	src_snap=snap1.H wavelet=wav.H \
	seis_record=$@ dump=wav2.H

grada3t.H: vel/v2.H ricker.H $B/beamWET.x 
	$B/beamWET_grad.x par=par/aliComp1t.P vel=$< data_obs=seisa1t.H \
	wavelet=wav.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=2 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H 

	
grada3t_o1.H: vel/v4.H ricker.H $B/beamWET.x 
	$B/beamWET_grad.x par=par/aliComp1t.P vel=$< data_obs=seisa1t.H \
	wavelet=wav.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=1 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H	


grada2t.H: vel/v2.H ricker.H $B/beamWET.x 
	$B/beamWET_grad.x par=par/aliComp1t.P vel=$< data_obs=seisa1t.H \
	wavelet=wav.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=2 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H	


###Biondo's  paper
vel/bvel1.H:
	Vel n1=130 n2=552 n3=1 o3=0 d3=1 o1=-150 o2=-200 vc=1000 d1=10 d2=10 > $@

## %0.1 vel perturbation
vel/bvel2.H: vel/vel1.H
	Vel n1=130 n2=552 n3=1 o3=0 d3=1 o1=-150 o2=-200 vc=1001 d1=10 d2=10 > $@
## % %8.0 vel perturbation
vel/bvel3.H: vel/vel1.H
	Vel n1=130 n2=552 n3=1 o3=0 d3=1 o1=-150 o2=-200 vc=1200 d1=10 d2=10 > $@
vel/bvel4.H: vel/vel1.H
	Vel n1=130 n2=552 n3=1 o3=0 d3=1 o1=-150 o2=-200 vc=800 d1=10 d2=10 > $@

rickerb.H:
	Wavelet n1=500 d1=0.004 wavelet=ricker2 phase=0 fund=10 tdelay=0.15 > $@
	echo 'peakFreq=15' >> $@

seisb1.H: vel/bvel1.H rickerb.H $B/Modeling2d.x 
	$B/Modeling2d.x par=par/bioComp.P vel=$< \
	src_snap=snap1.H wavelet=rickerb.H \
	seis_record=$@ dump=wav2.H boundary_type=1
seisb1_bndPar.H: vel/bvel1.H rickerb.H $B/Modeling2d.x 
	$B/Modeling2d.x par=par/bioComp.P vel=$< \
	src_snap=snap1_bndPar.H wavelet=rickerb.H \
	seis_record=$@ dump=wav2.H boundary_type=11

gradb1-3-line.H: vel/bvel1.H $B/beamWET.x 
	$B/beamWET_grad.x par=par/bioComp.P vel=$< data_obs=seisb1.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=3 corr=corrb1-3.H moveOutType=1 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H	
gradb1-3-quad.H: vel/bvel1.H $B/beamWET.x 
	$B/beamWET_grad.x par=par/bioComp.P vel=$< data_obs=seisb1.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=3 corr=corrb1-3.H moveOutType=2 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H	

gradb2-2.H: vel/bvel2.H $B/beamWET.x 
	$B/beamWET_grad.x par=par/bioComp.P vel=$< data_obs=seisb1.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=2 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H	
gradb2-1.H: vel/bvel2.H $B/beamWET.x 
	$B/beamWET_grad.x par=par/bioComp.P vel=$< data_obs=seisb1.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=1 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H	
gradb2-3-line.H: vel/bvel2.H $B/beamWET.x 
	$B/beamWET_grad.x par=par/bioComp.P vel=$< data_obs=seisb1.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=3 corr=corrb2-3.H moveOutType=1 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H	
gradb2-3-quad.H: vel/bvel2.H $B/beamWET.x 
	$B/beamWET_grad.x par=par/bioComp.P vel=$< data_obs=seisb1.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=3 corr=corrb2-3.H moveOutType=2 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 \
	DJDU=DJDU2-3-quad.H backProjData=backData2-3-quad.H

gradb2-4.H: vel/bvel2.H $B/beamWET_grad.x 
	$B/beamWET_grad.x par=par/bioComp.P vel=$< data_obs=seisb1.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=4 corr=corrb2-4.H dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H A2=A2-2-4.H B2=B2-2-4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H	

gradb3-1.H: vel/bvel3.H $B/beamWET.x 
	$B/beamWET_grad.x par=par/bioComp.P vel=$< data_obs=seisb1.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=1 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H	
gradb3-2.H: vel/bvel3.H $B/beamWET.x 
	$B/beamWET_grad.x par=par/bioComp.P vel=$< data_obs=seisb1.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=2 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H backProjData=bdb3-2.H
gradb3-2-g.H: vel/bvel3.H 
	$B/CudaBeamWET_grad.gpux par=par/bioComp.P vel=$< data_obs=seisb1.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=2 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H backProjData=bdb3-2-g.H

gradb3-3-line.H: vel/bvel3.H $B/beamWET.x 
	$B/beamWET_grad.x par=par/bioComp.P vel=$< data_obs=seisb1.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=3 corr=corrb3-3.H moveOutType=1 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H \
	DJDU=DJDU3-3-line.H  H=H3-3-line.H E=E3-3-line.H \
	backProjData=bdb3-3-line.H bd-l=bdb3-3-l.H bd-m=bdb3-3-m.H
gradb3-3-quad.H: vel/bvel3.H $B/beamWET.x 
	$B/beamWET_grad.x par=par/bioComp.P vel=$< data_obs=seisb1.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=3 corr=corrb3-3.H moveOutType=2 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H	\
	DJDU=DJDU3-3-quad.H backProjData=backData3-3-quad.H

gradb3-4.H: vel/bvel3.H $B/beamWET_grad.x 
	$B/beamWET_grad.x par=par/bioComp.P vel=$< data_obs=seisb1.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=4 corr=corrb3-4.H dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H fh=fhb3-4.H gh=ghb3-4.H \
	Emat=Emat3-4.H Fmat=Fmat3-4.H A2=A2-3-4.H B2=B2-3-4.H backProjData=bdb3-4.H \
	trcA=trcA3-4.H trcA1=trcA1-3-4.H trcB1=trcB1-3-4.H Hfg=Hfg3-4.H  \
	Csub=Csub3-4.H C1sub=C1sub3-4.H C2sub=C2sub3-4.H H12=H12-3-4.H 
gradb3-6.H: vel/bvel3.H $B/beamWET_grad.x 
	$B/beamWET_grad.x par=par/bioComp.P vel=$< data_obs=seisb1.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=6 corr=corrb3.H dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H fh=fhb3-4.H gh=ghb3-4.H \
	aCube=aCube3-6.H bCube=bCube3-6.H vTrue=1000 v0=1200 backProjData=bdb3-6.H

#Use DSO
gradb3-7.H: vel/bvel3.H 
	$B/beamWET_grad.x par=par/bioComp.P vel=$< data_obs=seisb1.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H	\
	backProjData=bdb3-7.H objFunc=7 corr=corrb3.H \
	dCdx=dCdxb3.H dCdx2=dCdx2b3.H

gradb3-7-g.H: vel/bvel3.H 
	$B/CudaBeamWET_grad.gpux par=par/bioComp.P vel=$< data_obs=seisb1.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H	\
	backProjData=bdb3-7-g.H objFunc=7 corr=corrb3-g.H \
	dCdx=dCdxb3-g.H dCdx2=dCdx2b3-g.H

## accurat a,b for the end points. a=2.0287e-15, b=3.477e-4
scanParamb3-4.H: vel/bvel3.H $B/beamWET_scanMoveoutParams.x 
	$B/beamWET_scanMoveoutParams.x par=par/bioComp.P vel=$< data_obs=seisb1.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=4 corr=corrb3-4.H dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H \
	o_a=-3.70e-7 d_a=3.70e-8 n_a=21 o_b=-3.14e-4 d_b=3.14e-5 n_b=21 \
	scanCube=$@ aCube=aCubeScan.H bCube=bCubeScan.H	

gradb3-4g.H: vel/bvel3.H 
	$B/CudaBeamWET_grad.gpux par=par/bioComp.P vel=$< data_obs=seisb1.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=4 corr=corrb3-4g.H dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H	\
	A2=A2-3-4g.H B2=B2-3-4g.H

gradb4-3-line.H: vel/bvel4.H $B/beamWET.x 
	$B/beamWET_grad.x par=par/bioComp.P vel=$< data_obs=seisb1.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=3 corr=corrb4-3.H moveOutType=1 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H	
gradb4-3-quad.H: vel/bvel4.H $B/beamWET.x 
	$B/beamWET_grad.x par=par/bioComp.P vel=$< data_obs=seisb1.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=3 corr=corrb4-3.H moveOutType=2 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H	

G = Graph title=" "
gcn = Grey newclip=1 color=j wantscalebar=1 gainpanel=e title=" "

$R/fhb3-4.v $R/ghb3-4.v: fhb3-4.H ghb3-4.H 
	< fhb3-4.H $G out=$R/fhb3-4.v ${dn}
	< ghb3-4.H $G out=$R/ghb3-4.v ${dn}

$R/Emat-3-4-11.v: Emat3-4.H
	< $< Window f3=0 n3=1 | $G out=$@ ${dn}
$R/Emat-3-4-12.v: Emat3-4.H
	< $< Window f3=1 n3=1 | $G  out=$@ ${dn}
$R/Emat-3-4-22.v: Emat3-4.H
	< $< Window f3=2 n3=1 | $G out=$@ ${dn}
$R/Fmat-3-4-11.v: Fmat3-4.H
	< $< Window f3=0 n3=1 | $G out=$@ ${dn}
$R/Fmat-3-4-12.v: Fmat3-4.H
	< $< Window f3=1 n3=1 | $G out=$@ ${dn}
$R/Fmat-3-4-22.v: Fmat3-4.H
	< $< Window f3=2 n3=1 | $G out=$@ ${dn}

$R/trcA3-4-l.v $R/trcA3-4-m.v: trcA3-4.H
	< $< Window n2=1 f2=0 | $G out=$R/trcA3-4-l.v ${dn}
	< $< Window n2=1 f2=1 | $G out=$R/trcA3-4-m.v ${dn}
$R/trcA1-3-4-l.v $R/trcA1-3-4-m.v: trcA1-3-4.H
	< $< Window n2=1 f2=0 | $G out=$R/trcA1-3-4-l.v ${dn}
	< $< Window n2=1 f2=1 | $G out=$R/trcA1-3-4-m.v ${dn}
$R/trcB1-3-4-l.v $R/trcB1-3-4-m.v: trcB1-3-4.H
	< $< Window n2=1 f2=0 | $G out=$R/trcB1-3-4-l.v ${dn}
	< $< Window n2=1 f2=1 | $G out=$R/trcB1-3-4-m.v ${dn}

$R/A2-3-4.v: A2-3-4.H
	< $< $G out=$@ ${dn}
$R/B2-3-4.v: B2-3-4.H
	< $< $G out=$@ ${dn}
$R/H12-3-4-1.v $R/H12-3-4-2.v:: H12-3-4.H 
	< $< Window n3=1 f3=0 | $G out=$R/H12-3-4-1.v ${dn}
	< $< Window n3=1 f3=1 | $G out=$R/H12-3-4-2.v ${dn}

$R/Hfg3-4-l.v $R/Hfg3-4-m.v: Hfg3-4.H 
	< $< Window n2=1 f2=0 | $G out=$R/Hfg3-4-l.v ${dn}
	< $< Window n2=1 f2=1 | $G out=$R/Hfg3-4-m.v ${dn}

$R/Csub3-4-l.v $R/Csub3-4-m.v: Csub3-4.H
	< $< Window n3=1 f3=0 | ${gcn} ${dn} out=$R/Csub3-4-l.v
	< $< Window n3=1 f3=1 | ${gcn} ${dn} out=$R/Csub3-4-m.v
$R/C1sub3-4-l.v $R/C1sub3-4-m.v: C1sub3-4.H
	< $< Window n3=1 f3=0 | ${gcn} ${dn} out=$R/C1sub3-4-l.v
	< $< Window n3=1 f3=1 | ${gcn} ${dn} out=$R/C1sub3-4-m.v
$R/C2sub3-4-l.v $R/C2sub3-4-m.v: C2sub3-4.H
	< $< Window n3=1 f3=0 | ${gcn} ${dn} out=$R/C2sub3-4-l.v
	< $< Window n3=1 f3=1 | ${gcn} ${dn} out=$R/C2sub3-4-m.v

$R/bdb3-3-l.v $R/bdb3-3-m.v: bdb3-3-l.H bdb3-3-m.H
	< bdb3-3-l.H Grey color=j newclip=1 bpclip=0 epclip=100 wantscalebar=1 ${dn} out=$R/bdb3-3-l.v
	< bdb3-3-m.H Grey color=j newclip=1 bpclip=0 epclip=100 wantscalebar=1 ${dn} out=$R/bdb3-3-m.v

figs:$R/fhb3-4.pdf $R/ghb3-4.pdf $R/Emat-3-4-11.pdf $R/Emat-3-4-12.pdf $R/Emat-3-4-22.pdf \
	$R/Fmat-3-4-11.pdf $R/Fmat-3-4-12.pdf $R/Fmat-3-4-22.pdf \
	$R/Csub3-4-l.pdf $R/Csub3-4-m.pdf $R/C1sub3-4-l.pdf $R/C1sub3-4-m.pdf \
	$R/C2sub3-4-l.pdf $R/C2sub3-4-m.pdf $R/trcA3-4-l.pdf $R/trcA3-4-m.pdf \
	$R/trcA1-3-4-l.pdf $R/trcA1-3-4-m.pdf $R/trcB1-3-4-l.pdf $R/trcB1-3-4-m.pdf \
	$R/A2-3-4.pdf $R/B2-3-4.pdf $R/H12-3-4-1.pdf $R/H12-3-4-2.pdf \
	$R/Hfg3-4-l.pdf $R/Hfg3-4-m.pdf $R/gradb3-2.pdf $R/gradb3-3-line.pdf \
	$R/bdb3-3-l.pdf $R/bdb3-3-m.pdf
	echo "make figs"

seisbm.H: vel/bvel1.H rickerb.H $B/Modeling2d.x 
	$B/Modeling2d.x par=par/bioComp2.P vel=$< \
	src_snap=snap1.H wavelet=rickerb.H \
	seis_record=$@ dump=wav2.H
seisb3m-g.H: vel/bvel1.H rickerb.H 
	$B/CudaModeling2d.gpux par=par/bioComp3.P vel=$< \
	src_snap=snap1.H wavelet=rickerb.H \
	seis_record=$@ dump=wav2.H
	

gradbm3-3-line.H: vel/bvel3.H 
	$B/beamWET_grad.x par=par/bioComp2.P vel=$< data_obs=seisbm.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=3 corr=corrbm3.H moveOutType=1 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H	
gradbm3-3-quad.H: vel/bvel3.H
	$B/beamWET_grad.x par=par/bioComp2.P vel=$< data_obs=seisbm.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=3 corr=corrbm3.H moveOutType=2 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H	

gradbm3-3-line-g.H: vel/bvel3.H 
	$B/CudaBeamWET_grad.gpux par=par/bioComp2.P vel=$< data_obs=seisbm.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=3 corr=corrbm3.H moveOutType=1 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H \
	grads=gradsbm3-3-line-g.H

gradbm3-7-g.H: vel/bvel3.H 
	$B/CudaBeamWET_grad.gpux par=par/bioComp2.P vel=$< data_obs=seisbm.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ corr=corrbm-3.H dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H \
	objFunc=7 dCdx=dCdxbm3-7-g.H dCdx2=dCdx2bm3-7-g.H backProjData=bdbm3-7-g.H
gradbm3-8-g.H: vel/bvel3.H 
	$B/CudaBeamWET_grad.gpux par=par/bioComp2.P vel=$< data_obs=seisbm.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ corr=corrbm-3.H dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H \
	objFunc=8 backProjData=bdbm3-8-g.H grads=gradsbm3-8-g.H
gradb3m3-7-g.H: vel/bvel3.H 
	$B/CudaBeamWET_grad.gpux par=par/bioComp3.P vel=$< data_obs=seisb3m-g.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ corr=corrb3m-3.H dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H \
	objFunc=7 dCdx=dCdxb3m3-7-g.H dCdx2=dCdx2b3m3-7-g.H

#Dispersion compare
seisbc1.H: vel/bvel1.H rickerb.H $B/Modeling2d.x 
	$B/Modeling2d.x par=par/bioComp3.P vel=$< \
	src_snap=snap1.H wavelet=rickerb.H \
	seis_record=$@ dump=wav2.H \
	use_4th_order_in_time=0
seisbc2.H: vel/bvel1.H rickerb.H $B/Modeling2d.x 
	$B/Modeling2d.x par=par/bioComp3.P vel=$< \
	src_snap=snap1.H wavelet=rickerb.H \
	seis_record=$@ dump=wav2.H \
	use_4th_order_in_time=1 

%_wb.H: %.H
	< $< Window n1=90 f1=20 n2=512 f2=20 > $@
%-sum.H: %-quad_wb.H %-line_wb.H
	Math file1=$*-quad_wb.H file2=$*-line_wb.H exp='file1+file2' > $@

#### Vel Anomaly, 1000 + 200 m/s
vel/velAb1.H vel/velAb1_0.H: $B/lens_twin.x
	$< nx=552 nz=130 dx=10 dz=10 x11=2760 x13=600 x21=2760 x23=900 \
		width1=500 height1=500 \
	alpha1=-0.2 alpha2=0.0 zref=1000 v0=1000	\
	refl=refl.H bvel=vel/velAb1_0.H > vel/velAb1.H

vel/velAb1w.H: $B/lens_twin.x
	$< nx=552 nz=130 dx=10 dz=10 x11=2760 x13=600 x21=2760 x23=900 \
		width1=3000 height1=500 \
	alpha1=-0.2 alpha2=0.0 zref=1000 v0=1000	\
	refl=refl.H bvel=vel/velAb1_0.H > vel/velAb1w.H
vel/velAb1n.H: $B/lens_twin.x
	$< nx=552 nz=130 dx=10 dz=10 x11=2760 x13=600 x21=2760 x23=900 \
		width1=150 height1=500 \
	alpha1=-0.2 alpha2=0.0 zref=1000 v0=1000	\
	refl=refl.H bvel=vel/velAb1_0.H > vel/velAb1n.H

vel/velAb2.H vel/velAb2_0.H: $B/lens_twin.x
	$< nx=552 nz=130 dx=10 dz=10 x11=2760 x13=600 x21=2760 x23=900 \
		width1=500 height1=500 \
		alpha1=0.2 alpha2=0.0 zref=1000 v0=1000	\
	refl=refl.H bvel=vel/velAb2_0.H > vel/velAb2.H

vel/velAbc.H vel/velAbc0.H: $B/lens_twin.x
	$< nx=552 nz=170 dx=10 dz=10 x11=2760 x13=845 x21=4000 x23=845 \
		width1=400 height1=400 \
		alpha1=-0.1 alpha2=0 zref=1200 v0=1000	\
	refl=refl.H bvel=vel/velAbc0.H > vel/velAbc.H


seisAb1.H: vel/velAb1.H rickerb.H $B/Modeling2d.x 
	$B/Modeling2d.x par=par/bioComp.P vel=$< \
	src_snap=snapAb1.H wavelet=rickerb.H \
	seis_record=$@ dump=wav2.H boundary_type=1
seisAb1w.H: vel/velAb1w.H rickerb.H $B/Modeling2d.x 
	$B/Modeling2d.x par=par/bioComp.P vel=$< \
	src_snap=snap1w.H wavelet=rickerb.H \
	seis_record=$@ dump=wav2.H boundary_type=1
seisAb1w-g.H: vel/velAb1w.H rickerb.H $B/CudaModeling2d.gpux 
	$B/CudaModeling2d.gpux par=par/bioComp.P vel=$< \
	src_snap=snap1w-g.H wavelet=rickerb.H \
	seis_record=$@ dump=wav2.H boundary_type=1
#### multiple shots
seisAb1wm.H: vel/velAb1w.H rickerb.H $B/Modeling2d.x 
	$B/Modeling2d.x par=par/bioComp2.P vel=$< \
	src_snap=snap1wm.H wavelet=rickerb.H \
	seis_record=$@ dump=wav2.H boundary_type=1

#narrow anomaly
seisAb1n.H: vel/velAb1n.H rickerb.H $B/Modeling2d.x 
	$B/Modeling2d.x par=par/bioComp.P vel=$< \
	src_snap=snap1n.H wavelet=rickerb.H \
	seis_record=$@ dump=wav2.H boundary_type=1
seisAb1n-g.H: vel/velAb1n.H rickerb.H $B/CudaModeling2d.gpux 
	$B/CudaModeling2d.gpux par=par/bioComp.P vel=$< \
	src_snap=snap1n-g.H wavelet=rickerb.H \
	seis_record=$@ dump=wav2.H boundary_type=1

gradAb1-3-line-g.H: vel/velAb1_0.H $B/CudaBeamWET_grad.gpux 
	$B/CudaBeamWET_grad.gpux par=par/bioComp.P vel=$< data_obs=seisAb1.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=3 corr=corrAb1-g.H dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H \
	DJDU=B2-Ab1-3-line-g.H  H=H-Ab1-3-line-g.H E=E-Ab1-3-line-g.H \
	backProjData=bdAb1-3-line-g.H moveOutType=1

#Use small window
gradAb1-3-line-beam5-g.H: vel/velAb1_0.H $B/CudaBeamWET_grad.gpux 
	$B/CudaBeamWET_grad.gpux par=par/bioComp.P vel=$< data_obs=seisAb1.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=3 corr=corrAb1-g.H dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H \
	DJDU=B2-Ab1-3-line-beam5-g.H  H=H-Ab1-3-line-beam5-g.H E=E-Ab1-3-line-beam5-g.H \
	backProjData=bdAb1-3-line-beam5-g.H moveOutType=1 beam_size=5
#Wide anamoly
gradAb1w-3-line-g.H: vel/velAb1_0.H $B/CudaBeamWET_grad.gpux 
	$B/CudaBeamWET_grad.gpux par=par/bioComp.P vel=$< data_obs=seisAb1w-g.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=3 corr=corrAb1w-g.H dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H \
	DJDU=B2-Ab1w-3-line-g.H  H=H-Ab1w-3-line-g.H E=E-Ab1w-3-line-g.H \
	backProjData=bdAb1w-3-line-g.H moveOutType=1

gradAb1w-3-line.H: vel/velAb1_0.H 
	$B/beamWET_grad.x par=par/bioComp.P vel=$< data_obs=seisAb1w.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=3 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H \
	DJDU=B2-Ab1w-3-line.H  H=H-Ab1w-3-line.H E=E-Ab1w-3-line.H \
	backProjData=bdAb1w-3-line.H moveOutType=1 corr=corrAb1w.H 

gradAb1wm-3-line.H: vel/velAb1_0.H 
	$B/beamWET_grad.x par=par/bioComp2.P vel=$< data_obs=seisAb1wm.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=3 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H \
	DJDU=B2-Ab1wm-3-line.H  H=H-Ab1wm-3-line.H E=E-Ab1wm-3-line.H \
	backProjData=bdAb1wm-3-line.H moveOutType=1 corr=corrAb1wm.H 
gradAb1wm-3-line-g.H: vel/velAb1_0.H 
	$B/CudaBeamWET_grad.gpux par=par/bioComp2.P vel=$< data_obs=seisAb1wm.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=3 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H \
	DJDU=B2-Ab1wm-3-line-g.H  H=H-Ab1wm-3-line-g.H E=E-Ab1wm-3-line-g.H \
	backProjData=bdAb1wm-3-line-g.H moveOutType=1 corr=corrAb1wm-g.H 

### shot on the edge
seisAb1w-shotL-g.H: vel/velAb1w.H rickerb.H $B/CudaModeling2d.gpux 
	$B/CudaModeling2d.gpux par=par/bioComp.P vel=$< \
	src_snap=snap1w-shotL-g.H wavelet=rickerb.H \
	seis_record=$@ dump=wav2.H boundary_type=1 shot1_ix=50 T=5.0
gradAb1w-3-line-shotL-g.H: vel/velAb1_0.H $B/CudaBeamWET_grad.gpux 
	$B/CudaBeamWET_grad.gpux par=par/bioComp.P vel=$< data_obs=seisAb1w-shotL-g.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=3 corr=corrAb1w-shotL-g.H dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H \
	DJDU=B2-Ab1w-3-line-shotL-g.H  H=H-Ab1w-3-line-shotL-g.H E=E-Ab1w-3-line-shotL-g.H \
	backProjData=bdAb1w-3-line-shotL-g.H moveOutType=1 shot1_ix=50 T=5.0

#narrow anomaly
gradAb1n-3-line-g.H: vel/velAb1_0.H $B/CudaBeamWET_grad.gpux 
	$B/CudaBeamWET_grad.gpux par=par/bioComp.P vel=$< data_obs=seisAb1n-g.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=3  dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H \
	DJDU=B2-Ab1n-3-line-g.H  H=H-Ab1n-3-line-g.H E=E-Ab1n-3-line-g.H \
	backProjData=bdAb1n-3-line-g.H moveOutType=1 corr=corrAb1n-g.H
gradAb1n-3-line.H: vel/velAb1_0.H 
	$B/beamWET_grad.x par=par/bioComp.P vel=$< data_obs=seisAb1n.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=3 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H \
	DJDU=B2-Ab1n-3-line.H  H=H-Ab1n-3-line.H E=E-Ab1n-3-line.H \
	backProjData=bdAb1n-3-line.H moveOutType=1 corr=corrAb1n.H 

seisAb2.H: vel/velAb2.H rickerb.H $B/Modeling2d.x 
	$B/Modeling2d.x par=par/bioComp.P vel=$< \
	src_snap=snap1.H wavelet=rickerb.H \
	seis_record=$@ dump=wav2.H boundary_type=1


gradAb1-2.H: vel/velAb1_0.H 
	$B/beamWET_grad.x par=par/bioComp.P vel=$< data_obs=seisAb1.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=2 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H backProjData=bdAb1-2.H
gradAb1-6.H: vel/velAb1_0.H $B/beamWET_grad.x 
	$B/beamWET_grad.x par=par/bioComp.P vel=$< data_obs=seisAb1.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=6 corr=corrAb1.H dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H fh=fhb3-4.H gh=ghb3-4.H \
	o_a=-3.70e-7 d_a=3.70e-8 n_a=21 o_b=-3.14e-4 d_b=3.14e-5 n_b=21 \
	aCube=aCubeAb1-6.H bCube=bCubeAb1-6.H vTrue=1000 v0=1200 backProjData=bdAb1-6.H
gradAb1-5.H: vel/velAb1_0.H $B/beamWET_grad.x 
	$B/beamWET_grad.x par=par/bioComp.P vel=$< data_obs=seisAb1.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=5 corr=corrAb1.H dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H fh=fhb3-4.H gh=ghb3-4.H \
	o_a=0.0 d_a=1.0 n_a=1 o_b=-2e-4 d_b=0.5e-5 n_b=81 \
	aCube=aCubeAb1-5.H bCube=bCubeAb1-5.H vTrue=1000 v0=1200 backProjData=bdAb1-5.H
gradAb1-4.H: vel/velAb1_0.H $B/beamWET_grad.x 
	$B/beamWET_grad.x par=par/bioComp.P vel=$< data_obs=seisAb1.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=4 corr=corrAb1.H dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H fh=fhb3-4.H gh=ghb3-4.H \
	backProjData=bdAb1-4.H
gradAb1-3-line.H: vel/velAb1_0.H $B/beamWET_grad.x 
	$B/beamWET_grad.x par=par/bioComp.P vel=$< data_obs=seisAb1.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=3 corr=corrAb1.H dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H \
	DJDU=B2-Ab1-3-line.H  H=H-Ab1-3-line.H E=E-Ab1-3-line.H \
	backProjData=bdAb1-3-line.H moveOutType=1

gradAb1-7.H: vel/velAb1_0.H 
	$B/beamWET_grad.x par=par/bioComp.P vel=$< data_obs=seisAb1.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ corr=corrAb1.H dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H \
	backProjData=bdAb1-7.H objFunc=7 dCdx=dCdxAb1.H dCdx2=dCdx2Ab1.H
gradAb1-7-g.H: vel/velAb1_0.H 
	$B/CudaBeamWET_grad.gpux par=par/bioComp.P vel=$< data_obs=seisAb1.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ corr=corrAb1-g.H dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H \
	backProjData=bdAb1-7-g.H objFunc=7 dCdx=dCdxAb1-g.H dCdx2=dCdx2Ab1-g.H

gradAb2-5.H: vel/velAb2_0.H $B/beamWET_grad.x 
	$B/beamWET_grad.x par=par/bioComp.P vel=$< data_obs=seisAb2.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=5 corr=corrAb2.H dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H fh=fhb3-4.H gh=ghb3-4.H \
	o_a=0.0 d_a=1.0 n_a=1 o_b=-2e-4 d_b=0.5e-5 n_b=81 \
	aCube=aCubeAb2-5.H bCube=bCubeAb2-5.H vTrue=1000 v0=1200 backProjData=bdAb2-5.H

seisAbm1.H: vel/velAb1.H rickerb.H $B/Modeling2d.x 
	$B/Modeling2d.x par=par/bioComp2.P vel=$< \
	src_snap=snap1.H wavelet=rickerb.H \
	seis_record=$@ dump=wav2.H
seisAb3m1.H: vel/velAb1.H rickerb.H 
	$B/CudaModeling2d.gpux par=par/bioComp3.P vel=$< \
	wavelet=rickerb.H \
	seis_record=$@ dump=wav2.H

gradAbm1-5.H: vel/velAb1_0.H $B/beamWET_grad.x 
	$B/beamWET_grad.x par=par/bioComp2.P vel=$< data_obs=seisAbm1.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=5 corr=corrAbm1.H dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H fh=fhb3-4.H gh=ghb3-4.H \
	o_a=0.0 d_a=1.0 n_a=1 o_b=-2e-4 d_b=0.5e-5 n_b=81 \
	aCube=aCubeAbm1-5.H bCube=bCubeAbm1-5.H vTrue=1000 v0=1200 backProjData=bdAbm1-5.H

gradAbm1-4.H: vel/velAb1_0.H $B/beamWET_grad.x 
	$B/beamWET_grad.x par=par/bioComp2.P vel=$< data_obs=seisAbm1.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=4 corr=corrAbm1.H dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H fh=fhb3-4.H gh=ghb3-4.H \
	backProjData=bdAbm1-4.H
gradAbm1-3-line.H: vel/velAb1_0.H $B/beamWET_grad.x 
	$B/beamWET_grad.x par=par/bioComp2.P vel=$< data_obs=seisAbm1.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=3 corr=corrAbm1.H moveOutType=1 \
	dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H \
	backProjData=bdAbm1-3-line.H
gradAbm1-3-line-g.H: vel/velAb1_0.H 
	$B/CudaBeamWET_grad.gpux par=par/bioComp2.P vel=$< data_obs=seisAbm1.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=3 corr=corrAbm1-g.H moveOutType=1 \
	dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H \
	backProjData=bdAbm1-3-line-g.H grads=gradsAbm1-3-line-g.H

gradAbm1-7-g.H: vel/velAb1_0.H 
	$B/CudaBeamWET_grad.gpux par=par/bioComp2.P vel=$< data_obs=seisAbm1.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ corr=corrAbm1-g.H dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H \
	backProjData=bdAbm1-7-g.H objFunc=7 dCdx=dCdxAbm1-g.H dCdx2=dCdx2Abm1-g.H grads=gradsAbm1-7-g.H
gradAbm1-8-g.H: vel/velAb1_0.H 
	$B/CudaBeamWET_grad.gpux par=par/bioComp2.P vel=$< data_obs=seisAbm1.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ corr=corrAbm1-g.H dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H \
	backProjData=bdAbm1-8-g.H objFunc=8 grads=gradsAbm1-8-g.H
gradAb3m1-7-g.H: vel/velAb1_0.H 
	$B/CudaBeamWET_grad.gpux par=par/bioComp3.P vel=$< data_obs=seisAb3m1.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ corr=corrAbm1-g.H dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H \
	backProjData=bdAb3m1-7-g.H objFunc=7 dCdx=dCdxAb3m1-g.H dCdx2=dCdx2Ab3m1-g.H \
	grads=gradsAb3m1-7-g.H

gradAbm1-2.H: vel/velAb1_0.H $B/beamWET_grad.x 
	$B/beamWET_grad.x par=par/bioComp2.P vel=$< data_obs=seisAbm1.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=2 corr=corrAbm1.H dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H fh=fhb3-4.H gh=ghb3-4.H \
	backProjData=bdAbm1-2.H
gradAbm1-2-g.H: vel/velAb1_0.H 
	$B/CudaBeamWET_grad.gpux par=par/bioComp2.P vel=$< data_obs=seisAbm1.H \
	wavelet=rickerb.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=2 corr=corrAbm1-g.H dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H \
	backProjData=bdAbm1-2-g.H

seisAbmc-g.H: vel/velAbc.H rickerb.H 
	$B/CudaModeling2d.gpux par=par/bioCompc.P vel=$< \
	wavelet=rickerb.H \
	seis_record=$@ dump=wav2.H
seisAbmc-alias-g.H: vel/velAbc.H rickerb.H 
	$B/CudaModeling2d.gpux par=par/bioCompc-alias.P vel=$< \
	wavelet=rickerb.H \
	seis_record=$@ dump=wav2.H

gradAbmc-7-g.H: vel/velAbc0.H 
	$B/CudaBeamWET_grad.gpux par=par/bioCompc.P vel=$< data_obs=seisAbmc-g.H \
	wavelet=rickerb.H \
	grad=$@ corr=corrAbmc-g.H dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H \
	backProjData=bdAbmc-7-g.H objFunc=7 dCdx=dCdxAbmc-g.H dCdx2=dCdx2Abmc-g.H \
	grads=gradsAbmc-7-g.H

gradAbmc-2-g.H: vel/velAbc0.H 
	$B/CudaBeamWET_grad.gpux par=par/bioCompc.P vel=$< data_obs=seisAbmc-g.H \
	wavelet=rickerb.H src_snap=src_snapAbmc.H recv_snap=recv_snapAbmc.H \
	grad=$@ corr=corrAbmc-g.H dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H \
	backProjData=bdAbmc-2-g.H objFunc=2 \
	grads=gradsAbmc-2-g.H


#### Jan.snap:
	< snap-simp-flat.H Window n3=1 f3=10 > temp.H;
	${gcn} < temp.H | Tube


include Tomo.Makefile

%.ps: %.v
	pstexpen $< $@ color=y fat=${FAT} fatmult=1.5 invras=n
%.pdf: %.ps
	epstopdf $<
%.pdf: %.eps
	epstopdf $<

clean: jclean
	-rm ${BINDIR}/*.* -f
	-rm ${OBJDIR}/*.* -f
	-rm *.mod *.o *.H -f
	-rm $R/*.* -f
#	@(cd $(srcdir);${MAKE} clean)

ss:
	rm ${BINDIR}/* -f

.PHONY: clean ss default1

#Type 'make' to build all the figures
#And then type 'scons' to build the report. 

include $(SEPINC)/SEP.bottom


